Forged by charge: polaron-induced matrix formation in silicon nitride conversion-type anodes for lithium-ion batteries

Abstract

The quest for high-capacity anode materials is vital in developing future lithium-ion battery technologies. While silicon-based anodes offer high theoretical capacity, their commercial realization is hindered by instability associated with large volume changes. Amorphous silicon nitride (a-Si₃N₄) has emerged as a promising alternative, acting as a conversion-type anode where lithium incorporation drives the formation of a structurally robust matrix and active phases. Here, we demonstrate that charge trapping, driven by polaron and bipolaron formation, governs the structural transformation of a-Si₃N₄ during the initial lithiation. These charge-induced modifications lead to the formation of a Li–Si–N matrix that stabilizes the anode framework. Matrix generation is accompanied by the development of Si-rich regions, serving as precursors for the active phase. We identify a progression from electronically active polarons to inactive bipolaron states, establishing a direct link between charge localization and matrix formation. These insights recast charge trapping from a passive consequence to a functional design parameter for optimizing conversion-type anodes.

---

1 Introduction

The development of more efficient battery technologies is a critical enabler of the green transition and the shift towards a circular economy. Metal-ion batteries occupy a central role in this transformation, underpinning energy storage and delivery across a wide range of applications – from renewable energy integration and grid-level storage to electric vehicles and personal electronic devices.^1–3 Looking ahead, two key challenges must be addressed by next-generation metal-ion batteries: increased anode capacity and improved cell lifetime.^4–10 These interdependent performance metrics remain central obstacles in the search for advanced anode materials.

These challenges have driven the search for new anode materials with both higher capacity and improved cycle life.^5,10–13 Since the early development of lithium-ion batteries, graphite has dominated as the anode material of choice due to its stability, despite a limited capacity of only 360 mA h g⁻¹.^4–9 Silicon, an alloying-type anode, has long been considered the natural successor to graphite, offering a theoretical capacity of 3579 mA h g⁻¹.^13–18 However, Si anodes suffer from severe structural degradation caused by a volumetric expansion and contraction of approximately 300% during lithiation and delithiation.^18–21 These extreme volume changes lead to mechanical failure of the composite cell, as other components are unable to accommodate the deformation.^13,22 In response, various structuring strategies have been explored to mitigate mechanical stress in silicon-based systems.^16,17,22 In parallel, attention has shifted towards alternative materials that form a Si-rich active phase in situ, particularly sub-stoichiometric oxides and more recently nitrides.^17,23 In these systems, the first lithiation step induces an irreversible conversion that yields both an electrochemically active Si-rich phase and an inert embedding matrix.^17,20 This composite structure has been shown to improve mechanical resilience and provide stable cycling performance after the initial transformation.²⁴ However, the mechanism underlying this initial conversion remains poorly understood.

Silicon nitride has been a technologically important material long before its consideration as a next-generation anode candidate, being employed in a broad range of applications including wear-resistant coatings,²⁵ electronic devices (ReRAM, MOSFETs, and MEMS),^26–28 high-energy optics,^29–31 and integrated photonics.^32,33 Almost universally, Si₃N₄ is deployed as an amorphous thin film (a-Si₃N₄).^{27,28,34–38} In battery applications, silicon nitride has been studied across a range of structures and stoichiometries, and prepared using a variety of growth methods.^{19,23,39–48} These diverse approaches have yielded a wide range of reported lithiation capacities, from 40 mA h g⁻¹ to 2000 mA h g⁻¹. In general, the highest capacities are observed in silicon-rich SiN_x compositions, although this comes at the expense of cycling stability. The optimal reported stoichiometry is approximately SiN_0.9, which delivers a capacity of 1200 mA h g⁻¹ with stable performance over more than 2000 cycles.^45,46 As stoichiometric a-Si₃N₄ is approached, capacity declines sharply. Conversely, further increases in silicon excess result in rapid capacity fade over just tens to hundreds of cycles.⁴⁵ Many explanations have been proposed for this tradeoff, but a definitive link between macroscopic performance and atomic-scale Li incorporation mechanisms remains an open question.

Recent studies of silicon nitride anodes have revealed key insights into the formation and stability of the matrix that emerges during the first lithiation cycle. Early reports proposed the formation of a Li₃N phase, motivated by thermodynamic considerations and observed improvements in ionic conductivity,¹⁹ though no direct structural characterization was available. More recent work by Ulvestad et al.,⁴⁷ employing Pair Distribution Function (PDF) and Energy Dispersive X-ray (EDX) analyses, suggests that the matrix is better described as a mixed Li–Si–N network, with the best-fit composition approximated by Li₂SiN₂ stoichiometry. Concurrent EDX mapping revealed a compositional segregation of the original SiN_x network into Si-rich domains and nitrogen-enriched amorphous lithium–silicon–nitride regions. In contrast, Kilian et al.⁴³ employed ⁷Li NMR to challenge the picture of a static inert matrix, instead suggesting that silicon nitride forms a redox-active Li–Si–N solid solution with evolving local environments throughout cycling. Lovett et al.⁴⁹ further confirmed the domain behavior and demonstrated both the mechanical robustness of the matrix and its stabilizing influence on the active (Si-rich) phase. As exemplified by this ongoing debate, definitive identification remains challenging due to the amorphous character of the matrix, local compositional variations, and potential electrochemical activity. These factors underscore the importance of elucidating the atomic-scale mechanisms that govern matrix formation and evolution.

Despite these important advances, the underlying mechanisms driving matrix formation and stability remain poorly understood. In particular, the atomic-scale processes that govern the initial lithiation of silicon nitride, and the influence of the local atomic environment on lithium incorporation and storage, are not yet resolved. The evolution of the matrix as a function of lithium concentration, and its dependence on the host network structure, are especially difficult to probe via post mortem measurements, since key features of the initial lithiation state are often obscured by subsequent structural relaxation. Furthermore, the relationship between the microstructure of the silicon nitride precursor – including its nanoporosity and local topology – and the morphology of the emerging matrix warrants deeper investigation. Addressing these questions is essential not only for optimizing silicon nitride-based anodes, but also for advancing the broader understanding of conversion type anodes for next-generation lithium-ion battery materials.

In this work, we present a systematic study of lithium incorporation in a-Si₃N₄, focusing on the initial steps of this process at the atomic scale. To ensure that a statistically meaningful range of incorporation sites are considered, we build on the sampling scheme developed in our previous work.^50,51 Using density functional theory (DFT) simulations, the stability, geometry, and interplay between the Li ions and intrinsic charge trapping⁵⁰ is probed. Our study provides a detailed atomic-scale understanding of the irreversible matrix formation and elucidates concomitant reversible charge trapping.

2 Computational methods

All DFT calculations were performed spin-polarized using the CP2K⁵² code with the HSE06 (ref. 53 and 54) exchange-correlation functional in combination with the auxiliary density matrix method (ADMM),⁵⁵ which reduces the computational cost of hybrid functional calculations. The DZVP-SR-MOLOPT⁵⁶ family of basis sets was employed to describe the valence electrons, together with GTH pseudopotentials^57–59 for the core electrons. Energy cutoffs were set after convergence testing to 650 Ry and 70 Ry for the relative cutoff, yielding a precision of 0.1 meV per atom. All geometry optimisations were performed using the BFGS algorithm,^60–63 with convergence criteria of 10⁻⁷ eV for energy differences and 0.001 eV Å⁻¹ for forces.

Using the standard formalism of Zhang and Northrup⁶⁴ at the hybrid functional level, average defect formation energies per lithium atom incorporated into a-Si₃N₄ were calculated according to

Here, E_a-Si3N4 and E_nLi:a-Si3N4 are the DFT total energies of the a-Si₃N₄ simulation cell without and with n inserted Li atoms. The chemical potential of Li (μ_Li) is taken from bulk lithium in its ground state body-centred cubic phase. Other common references that can be found in the literature are EC/DMC solvated Li⁺ ions or Li⁰ atoms in vacuum. Due to reduced Li–Li interaction the corresponding values for μ_Li are ≈1.0 eV and ≈1.8 eV lower, respectively.⁶⁵ E_form as defined in eqn (1) can also be read as average incorporation energy for the n inserted Li atoms relative to bulk Li. An alternative approach would be to calculate the formation energy in a “step-wise” manner where the formation energy for Li_n is referenced to Li_n−1. While this has no impact on the reported trends, the energy range is expanded as there is no averaging over Li sites in this approach (as in each case n = 1), this treatment is shown in Fig. S6 in the SI.

The primary focus of this work is on Li incorporation in amorphous silicon nitride (a-Si₃N₄), where the vast configuration space necessitates the multi-stage sampling workflow shown in Fig. 1. To contextualize the results for the range of sites found in the amorphous system, calculations for Li₁ were performed on crystalline β-Si₃N₄ to serve as an ordered reference model. The following sections detail the generation and sampling of the amorphous structures, while a full description of the crystalline calculations can be found in the SI (see S1 and S2).

Fig. 1 Computational workflow with the key stages highlighted: 1 and 2 are the structure generation via an MD melt quench,⁵⁰ highlighted in red. Stage 3 represents the statistical convergence of the ensemble with respect to a given descriptor space, and stage 4 and 5 (highlighted with blue numbers) represent the site sampling via Voronoi tesselation.⁵¹ Stage 6 and 7 (in black) iteratively repeat an exhaustive PBE screening followed by a GMM sampling and finally the HSE06 production run for a given Li concentration.

The amorphous host network is first generated via a classical melt-quench molecular dynamics trajectory. A 280-atom β-Si₃N₄ supercell was melted and then cooled at a rate of 1 K ps⁻¹ under an isothermal–isobaric (NpT) ensemble using the MG2 interaction potential,⁶⁶ following the protocol in ref. 50 (Fig. 1, stages 1–2 highlighted in red). As structure sampling represents a vital consideration for amorphous systems, an ensemble is generated that is statistically converged with respect to a defined descriptor space. From this ensemble a single, statistically representative 280-atom a-Si₃N₄ cell was selected. This approach disentangles the problem of site sampling from the broad range of trap sites previously described in a-Si₃N₄ (Fig. 1, stages 3–5 highlighted in blue).^50,51

With the host structure defined, the statistical sampling scheme previously developed for intrinsic and H-related defects^50,51 was extended to Li-incorporation. An exhaustive set of initial Li insertion sites was generated via a Voronoi tesselation centered on each atom of the host network (Fig. 1, stages 4 and 5). For each site, the Li atom was placed on the Voronoi polyhedron face with the largest area, while maintaining a minimum distance of 2.0 Å from all neighboring atoms to prevent spurious steric repulsion in the initial configuration. These 5844 distinct starting geometries were then fully relaxed at the PBE level of theory^67,68 (Fig. 1, stage 6). During these relaxations, the Li atom is unconstrained, allowing it to settle into a local minimum within the “cage” imposed by the amorphous network; its final position is therefore not necessarily proximate to the initial reference atom. This exhaustive screening yields a comprehensive dataset of 5844 relaxed, Li-doped structures, providing a robust statistical basis for subsequent analysis. A full breakdown by Li concentration is provided in the SI (Table S2).

A Gaussian Mixture Model (GMM) was leveraged to distill the 5844-structure PBE dataset into a compact yet comprehensive subset for the HSE06 production calculations. This process ensures that the selected structures fully represent the configurational diversity of the original dataset (Fig. 1, stage 6). The GMM is particularly well-suited for this application due to its inherent flexibility in modeling the complex, multimodal distributions of local atomic environments characteristic of amorphous materials.^69–73 For each Li concentration, the optimal number of Gaussian components, was determined by evaluating the Bayesian (BIC) and Akaike (AIC) Information Criteria, which balance model fidelity against complexity.^74–77 From the fitted GMM, a representative ensemble was constructed by sampling configurations from the center of each Gaussian component (central modes) as well as from the tails of the distribution (high-variance outliers). This strategy ensures that the full diversity of local environments is captured. The procedure reduces the full set of 5844 PBE structures to 280 representative configurations, which then serve as the initial geometries for the HSE06 production runs. The process is then repeated across the desired compositional range Li₁–Li₁₀. The simulated concentrations (Li_n, n = 1–10) correspond to gravimetric capacities of 4.8 mA h g⁻¹ to 47.2 mA h g⁻¹. The complete mapping of n, capacity and stoichiometry is given in Table S2 of the SI.

The statistical framework underlying this workflow is deliberately chosen to reflect the non-equilibrium nature of the initial lithiation of a-Si₃N₄. Li incorporation is an irreversible, path-dependent reaction that breaks ergodicity: once inserted, the network becomes confined to a hierarchy of evolving local minima and the a priori probability of accessing any particular minimum is unknown. Standard equilibrium (Boltzmann-weighted) statistics are therefore inapplicable.^78–82 For every physical property we report the full distribution without applying any configurational weights. The GMM retains the character of the original sampling by weighting each identified structural motif in proportion to its population, ensuring that the final ensemble reproduces the diversity of the full energy landscape. The approach is analogous to the “inherent-structure” formalism widely used in glass physics.^83,84

The entire workflow was implemented in Python. The ASE library⁸⁵ was used for structure manipulation and database management, SciPy⁸⁶ for Voronoi analysis, and scikit-learn⁸⁷ for GMM fitting and sampling. Visualizations were generated with matplotlib,⁸⁸ atomic structures with VMD⁸⁹ and Voronoi polyhedra with Ovito.⁹⁰

To characterize local distortions in the a-Si₃N₄ network induced by lithium incorporation, the associated strain field has been analyzed in detail. For each individual configuration containing inserted Li, the displacement vector of every Si and N atom is computed with respect to its position in Li_n−1. Averaging these displacements over all configurations belonging to a given Li_n ensemble yields an average displacement vector for each atom in the structure. To visualize anisotropy in the ensemble-average strain, a continuous scalar strain density field is then constructed by superimposing asymmetric Gaussian functions centered at the undistorted atomic positions. The resulting field is defined as

where R_i is the undistorted position of atom i and Σ_i is the covariance matrix constructed from the direction and magnitude of the associated displacement vectors. To reduce visual noise, displacements with magnitudes less than 0.1 Å are excluded from the strain density evaluation. For each Li_n ensemble the strain-density field is projected onto the crystallographic plane that shows the largest in-plane displacement, as determined from the full three-dimensional strain field. Strain contributions are averaged along its normal direction to yield a scalar field that is rendered as a heat map, with colour intensity proportional to the local strain magnitude. Projections onto other planes—whether at different z values or along alternative orientations—were found to add no further insight and are therefore omitted. The maps are displayed in the original simulation coordinates without re-centering relative to periodic boundary conditions.

3 Results and discussion

3.1 Initial incorporation: Li₁

We first study the incorporation of a single Li in a-Si₃N₄. After relaxation, the incorporation environments can be classified based on the local coordination within the amorphous network and the associated electron trapping sites (Fig. 2). Li incorporation results in a wide range of structural modifications, leading to a broad distribution of formation energies. As shown in Fig. 2a, they span from −2.81 eV to 1.82 eV (referenced to lithium bulk metal–see eqn (1)). For each insertion site, the mean bond length of Li with the neighboring atoms in its first coordination shell is depicted in Fig. 2b. Typically, Li is coordinated by two or three N atoms with a bond length of approximately 2 Å. However, the incorporation of Li often results in one or more bonds being significantly extended by 10% to 20% to accommodate the Li atom. In a small fraction of cases (9.8%), steric crowding forces Li into an extended coordination environment that includes both N and Si atoms, where steric constraints prevent Li from relaxing into a more favorable configuration (Fig. 2c). As observed previously for H defects in a-Si₃N₄,⁵¹ there is a weak correlation between steric repulsion and Li formation energy, with higher formation energies associated with smaller Voronoi volumes (Fig. 2c). It is important to note that the relationship between Li geometry, formation energy, and the morphology of the a-Si₃N₄ network is complex and cannot be fully captured by a simple descriptor.

Fig. 2 (a) Distribution of Li-formation energies in a-Si₃N₄. The total distribution (black outline) is subdivided based on coordination (pink and turquoise) and trapping site (crosshatched). (b) Li-bond length distribution in the first coordination shell, with the same subdivision as in (a). Histogram bars represent the probability–density estimate ; the distribution is normalised by ∫p(x) dx = 1. (c) Relationship between the Voronoi volumes of the incorporated Li atoms and the formation energies. (d) Schematic representation of the main Li-incorporation configurations in the a-Si₃N₄ network.

The incorporation of Li in a-Si₃N₄ is governed by the reaction at the anode: Li⁰ → Li⁺ + e⁻. The interaction of the resulting electron with the amorphous network significantly influences the mode of Li storage. This behavior is in stark contrast to that in crystalline β-Si₃N₄, where the electron does not localize but instead gives rise to a shallow donor state (see SI S1 and S2). This interaction allows the Li sites to be categorized into two main classes. The geometries corresponding to each trapping type are illustrated schematically in Fig. 2d, highlighting both similarities and key differences. In 73% of the sites, the electron is accommodated at an intrinsic trap site within the a-Si₃N₄ network, previously characterized and independent of Li except as the electron source.⁵⁰ These sites are labeled as intrinsic traps in Fig. 2a–c and schematically depicted in Fig. 2 subpanel d : i. The remaining 27% of configurations involve electron traps that are induced by distortions in the network, driven by the presence of Li. In these cases, the electron is trapped within the coordination sphere of Li. These induced traps can be further subdivided based on whether Li induces the distortion without directly interacting with the trap site (Fig. 2dii) or directly interacts with the induced trap site (Fig. 2diii). The formation energy for Li with induced traps is significantly higher than for those with intrinsic traps, dominating configurations with energies above −0.5 eV (Fig. 2a). There is no straightforward geometric predictor to distinguish between sites that result in intrinsic versus induced trapping, as both types are observed across the full range of bond lengths (Fig. 2b), angles, and steric environments (Fig. 2c).

Returning to the geometries of the trapping sites, for the intrinsic traps (Fig. 2di), the electron is decoupled from the Li⁺ site. The Li⁺ is typically coordinated by 2 to 3 N atoms, forming the Li coordination shell. The Li–N bond lengths vary according to the local environment, but generally, 2 to 3 of the bonds are shorter than 2 Å, while 1 to 2 are longer, exceeding 2.5 Å. This variation results in a range of site symmetries, including C₂, C₃, T_d, and their broken symmetry variants. The lowest energy configurations occur when Li⁺ interacts with 2-coordinated N atoms, due to their negative polarization relative to 3-coordinated N atoms. These low-energy sites are rare – both in experiments^34,35 (6%) and our computational setup (1.67%).

For induced traps (Fig. 2dii and iii), the situation is more complex, as the electron-induced distortion and Li incorporation are interdependent. In Fig. 2dii, Li interacts with N atoms connected to, but not directly part of, the trap site. The presence of Li induces a distortion that creates an electron trap at a nearby Si atom. The Li coordination and Li–N separation in this case are similar to those in intrinsic traps, with 2 to 3 N atoms as nearest neighbors. However, due to the link to the trap site, a Si atom is always found as the next neighbor at a distance of approximately 2.5 Å, contrasting with the intrinsic trap scenario. Finally, in the third type of trap (Fig. 2diii), there is a direct interaction between Li and the trapping Si atom, forming a Li–Si interaction. Here, the nearest neighbor is a Si atom at a separation of 2.1 Å, with 2 to 3 N atoms at a similar distance.

The distortion in the amorphous network is primarily driven by the electron trapping process, whether intrinsic or Li-induced. All of the traps shown in Fig. 3 share the same fundamental motif—an excess electron localised on a four-coordinated, strained Si. The distinction lies only in origin: intrinsic traps being an inherent part of the amorphous network (Fig. 3a), whereas induced traps appear only after the network is distorted by a proximate Li ion (Fig. 3b and c). The states occupied and their relationship to Li are shown in Fig. 3d–f. In the case of the intrinsic traps the Li and trap are uncorrelated (Fig. 3d), with the same trap site being occupied regardless of Li position. Whereas for the induced traps the Li is in close proximity, in the next-neighbour or next-nearest-neighbour shell of trap site (Fig. 3e and f), and a broad range of trap sites are occupied. Induced traps are further divided depending of whether the Li drives the distortion only (Fig. 3b, e and h) or whether it directly interacts with the trap site (Fig. 3c, f and i).

Fig. 3 Schematics for the relaxation of each of the of electron traps showing the intrinsic trap (a), induced trap non-interacting Li (b), and the induced trap where Li is interacting (c). The corresponding occupied trap states are visualised in (d) for the intrinsic trap, (e) for the induced trap with non-interacting Li, and (f) induced trap with interacting Li. The density of states are shown in (g) for the intrinsic trap, (h) for the induced trap with non-interacting Li, and (i) induced trap with interacting Li. The filled states in (g)–(i) are indicated by the shaded region below the Fermi energy. For clarity the Li peak has been increased by a factor of 25 as indicated in the legends.

In both intrinsic and non-interacting induced traps, the density of states (DoS) shows an occupied state at (Fig. 3g) or near (Fig. 3h) the valence band maximum (VBM), and a related unoccupied state at (Fig. 3h) or just below (Fig. 3g) the conduction band minimum (CBM), that are predominantly Si-character trap states. The position of the states is dictated by the extent of local relaxation that can be accommodated by a given local environment. However, in the case of induced traps where Li directly interacts with the trap site, both the relaxation and DoS are distinct (Fig. 3c, f and i). This is characterized by the rearrangement of a small number of atoms representing the immediate coordination shell, with the DoS showing trap states of mixed Li and Si character. The different trapping sites exhibit markedly different Mulliken charges associated with the Li sites. As expected, Li remains cationic in all cases. In intrinsic and non-interacting induced traps, the Mulliken charges range from 0.6 e⁻ to 0.9 e⁻. However, when Li directly interacts with the trap site, there is a significant reduction in Mulliken charges, ranging from 0.2 e⁻ to 0.4 e⁻, indicating electron sharing between Si and Li and confirmed by the DoS (Fig. 3i).

3.2 Li₂, Li₃ and Li₄ incorporation

The incorporation of a second Li atom is driven by electrostatic repulsion, maximizing the distance between the two like-charged Li⁺ ions, as might reasonably be expected. The resulting energetic landscape is strongly influenced by charge trapping: the average formation energy per Li atom becomes significantly more favourable due to bi-polaron formation (Fig. 4). In contrast to the Li₁ case, where multiple trap sites are accessible (Fig. 2d), all Li₂ configurations relax to a single dominant trap in which a second electron is trapped at the original intrinsic trap site (Fig. 4a). As a consequence no Li-induced traps are observed. The associated structural relaxation is pronounced: the trap-forming Si atom back-projects to form a Si–Si bond with a neighboring Si atom (Fig. 4c and d). The ensemble-averaged strain map (Fig. 4c) confirms that this distortion is consistent across configurations and decoupled from local Li-induced relaxation. Structural variations specific to individual cells are suppressed in the averaging, leaving only the common bi-polaron signature. It should be noted that the range of E_form displayed in Fig. 4a reflects the ease with which two Li ions can be incorporated into the a-Si₃N₄ network; with the occupied bi-polaron state constant across all configurations, the remaining variation arises solely from the local Li⁺ incorporation environment.

Fig. 4 (a) Distribution of Li₂ formation energies in a-Si₃N₄. As in Fig. 2a, the total distribution (black outline) is subdivided based on coordination and trapping site. Histogram bars represent the probability–density estimate ; the distribution is normalised by ∫p(x) dx = 1. (b) DoS showing the occupied bi-polaron state, the filled states are indicated by the shaded region below the Fermi energy. For clarity the Li peak has been increased by a factor of 25. (c) Average 2D projection of the strain induced as a result of Li₂ incorporation and the associated bi-polaron formation. (d) Doubly occupied bi-polaron state forming a Si–Si bond with a length of 2.32 Å. Blue spheres shows nitrogen, and yellow silicon. Red iso-surface shows the spin up channel, and light blue the spin down channel.

This reorganization modifies the electronic structure (Fig. 4b). In the single polaron case, the trapped electron localizes on a Si dangling bond, yielding occupied and unoccupied defect states in the band gap (Fig. 3d). In contrast, the bi-polaron states are approximately 0.7 eV below the valence band maximum, and no occupied states remain in the band gap. The unoccupied Li-derived states remain in the conduction band, rendering the system electronically inactive. The Si–Si bond persists after removal of the Li atoms (and/or the associated electrons), indicating that the structural rearrangement is irreversible. This mirrors the behavior previously reported for hydrogen incorporation, where a single-electron polaron relaxes irreversibly to form a bi-polaron state on the addition of a second electron.⁵¹

The Li₃ system shares several important similarities with the Li₁ case, with the key difference being that the lowest energy bi-polaron state in the network is occupied and therefore unavailable. This has two significant consequences. First, the distinction between intrinsic and induced traps becomes less clear-cut, as the network distortion is influenced by both the presence of Li and the bi-polaron. To maintain consistency, the following definition is applied: if trapping occurs within the extended coordination sphere (including both the N and Si coordination shells) of a Li atom, it is considered induced. Otherwise, it is classified as intrinsic. Second, the energy difference between induced and intrinsic traps becomes less pronounced. As shown in Fig. 5a and b, both types are now observed across the full range of formation energies rather than being confined to higher energy configurations as seen in Fig. 2a. The lower energy configurations continue to favor Li coordination by 3 or 2 N atoms, with an increased occurrence of Li sites exhibiting coordination numbers greater than 3. The intrinsic and induced traps can still be divided into the broad categories described for Li₁. The main difference is that instead of a single intrinsic trap site, there are now three approximately iso-energetic sites. It is important to note that the formation energy results from both electron trapping and the accommodation of Li⁺ within the lattice.

Fig. 5 (a) Shows a schematic picture of the bi-polaron and the subsequent polaron, and (b) the distribution of Li₃ formation energies. The total distribution (black outline) is divided by coordination (pink and turquoise) and the nature of the trapping site (crosshatched). Histogram bars represent the probability–density estimate ; the distribution is normalised by ∫p(x) dx = 1.

In each case, trapping occurs at an existing wide bond angle within the amorphous network, where the electron localizes. The induced traps exhibit the same behaviors, with Li either inducing a distortion that leads to electron trapping within the Li coordination sphere or directly interacting with the trap site. Both of these configurations have already been illustrated in Fig. 5a and are the same as those shown in Fig. 2d. The presence of the existing electron trap in the amorphous network does not impact Li storage or electron trapping beyond the geometric distortions it induces. Electronically, the Si–Si bi-polaron states remain within the valence band (Fig. 4b) and are thus electrically inactive with respect to both Li⁺ and charge redistribution within the network.

The Li₄ system exhibits several important similarities with Li₂, with the lowest energy bi-polaron already filled, and the next lowest energy configuration also occupied. As before, a single bi-polaron forms independently of the Li geometry and its position relative to the trap site. However, unlike in the Li₂ case, there is no significant decrease in E_form; rather, the energy distribution of Li₄ configurations is tighter, with a greater number of low-energy states (Fig. 6a). The second bi-polaron shows similar relaxation behavior, although it is constrained by a different local environment. Specifically, one Si center back-projects to form an Si–Si interaction, accompanied by an outward relaxation of the neighboring N atoms to accommodate this change. This relaxation is depicted in Fig. 6c, which shows the Si centers moving towards each other and the displacement of the adjacent N atoms. The DoS for Li₄ (Fig. 6b) also reveals a defect free band gap, similar to Li₂, with the bi-polaron states situated below the valence band maximum (VBM). The doubly occupied trap state, shown in Fig. 6d, features the same Si–Si bonding interaction observed in the Li₂ system.

Fig. 6 (a) Distribution of Li₄ formation energies in a-Si₃N₄, the total distribution (black outline) is divided by coordination and trapping site. Histogram bars represent the probability-density estimate ; the distribution is normalised by ∫p(x) dx = 1. (b) DoS showing the occupied bi-polaron state, the filled states are indicated by the shaded region below the Fermi energy. For clarity the Li peak has been increased by ×25. (c) A 2D projection of the strain induced as a result of Li₄ incorporation and the associated bi-polaron formation. (d) The doubly occupied bi-polaron state forming a Si–Si bond with a length of 2.4 Å. Blue spheres shows nitrogen, green lithium, and yellow silicon. Red iso-surface shows the spin up channel, and light blue the spin down channel.

3.3 Increasing Li concentration

As the Li concentration increases, the pattern observed in Li₁ through Li₄ is continued, with a diverse range of single-electron traps giving way to bi-polaron states. The consequences for the formation energies are shown in Fig. 7: initially, there is a strong energetic driver for bi-polaron formation, with a marked drop in formation energy from Li₁ to Li₂, driven by the collapse of the broad variety of electron traps described in Section 3.1 into a single bi-polaron configuration. A similar trend is seen from Li₃ to Li₄, favoring bi-polaron formation, although the magnitude of this stabilization is greatly reduced.

Fig. 7 Progression in formation energy as a function of number of Li in the range Li₁ to L₁₀. For each concentration, the violin plot shows the underlying distribution resulting from the structural ensemble. The red envelope function provides a guide for the eye. The inset zooms into the concentration range from Li₅ to Li₁₀.

Li₄ represents the filling of the last intrinsic bi-polaron site in the host lattice. For Li₅ and above, the trap states are a consequence of Li-incorporation and the resulting structural modifications (resulting from the Li and electron trapping). This shift from filling intrinsic to creating induced trap states is accompanied by an increase in energy from Li₄ to Li₅. The trend then echoes the behavior of the intrinsic trap states (Li₁ to Li₄): a wide distribution of polaron states collapses into a single bi-polaron configuration (at Li₆). In each case the bi-polaron states impart an irreversible change on the lattice, creating Si–Si bonded motifs.

Beyond Li₇, a single dominant bi-polaron configuration no longer exists; instead, there are three to four geometrically distinct configurations. These all exhibit similar relaxation behavior, forming an Si–Si bond with lengths in the range of 2.2 Å to 2.4 Å, and represent the onset of more complex, polyfurcation of bi-polaron formation. Looking at the electronic structure, while all bi-polaron states reside below the VBM, there is a clear shift up towards the VBM as the Li concentration increases.

4 Discussion and summary

The initial incorporation of Li into a-Si₃N₄ induces significant structural modifications through the formation of various one- and two-electron trap states, as governed by the reaction:

nLi₀ → nLi⁺ + ne⁻_trap (3)

Despite the wide range of local environments within the amorphous network, common structural features emerge, as shown in Fig. 8. A broad distribution of Li–N bonds are found ranging from 1.8 Å to 2.5 Å (Fig. 8a), while Li–Si bond lengths, which are spread over a wider range of 2.0 Å to 3.6 Å (Fig. 8b), are more sensitive to the Li-concentration (Fig. 8b). The secondary Li–Si feature arises from Si atoms in the second coordination shell, with bond lengths above 3.0 Å, seen in both induced and intrinsic traps. Finally, the Si–Si bond lengths show a broad distribution centered at 2.45 Å at Li₁ (Fig. 8c), representing the previously described polaron states. These are significantly shorter than the Si–Si first-shell maximum at 3.5 Å (see SI Fig. S3b). At Li₂ and above, the distribution becomes bi-modal, with the bi-polaron states at ≈2.25 Å and a distribution of precursor states below 2.5 Å (Fig. 8c). As Li concentration increases, these features evolve, as shown in Fig. 8. Notably, the Li–N peak broadens significantly after Li₄, coinciding with an upward shift of the formation-energy distribution (Fig. 7), indicating a transition from intrinsic to induced trapping, where Li-induced network distortions facilitate electron trapping. Concurrently, the Li–Si bond length range broadens as induced traps become more prevalent, driving the skew to shorter Li–Si separations at higher Li concentrations. The bi-polaron becomes more pronounced as more Li atoms are incorporated, transitioning from single polaron to bi-polaron states.

Fig. 8 Bond length distributions as a function of Li concentration for (a) Li–N and (b) Li–Si. (c) Distributions for the short Si–Si pair distances that result from bi-polaron formation and the precursor polaron states. The dashed line reflects the full Si–Si pair distribution as shown in the SI S3b. All histograms are empirical, equal-weight distributions of the sampled configurations (no re-weighting by energy).

The progression from intrinsic to induced trapping follows a clear sequence. Initially, those intrinsic traps, which are naturally present in the network, get filled, followed by the formation and occupation of induced traps. The trapping of 2e⁻ drives an irreversible network relaxation, forming Si–Si bonds, which in turn drives the sequential occupation of all intrinsic traps in the system. In the reference cell two intrinsic bi-polarons are able to relax. Once these intrinsic traps are saturated, the induced traps are filled – both polarons and bi-polarons. Although induced traps have higher formation energies than intrinsic traps, they remain significantly favored with respect to the delocalised state, especially at higher Li concentrations. This irreversible structural modification—marked by the creation of Si–Si formation—mirrors earlier findings in a-Si₃N₄ under hydrogen incorporation.⁵¹ This stepwise process provides the atomistic basis for the formation and stability of the Si-rich regions during the initial lithiation. The observed relationship between stoichiometry and Li storage can be understood through these trap states. It is well-established that Si-rich a-Si₃N₄ increases the concentration of electron trap states,^91–93 a trend further enhanced by the presence of Si…Si precursor states capable of trapping electrons.^28,94 This characteristic likely contributes to the capacity fade observed as sub-stoichiometry increases, where the network's ability to accommodate distortions diminishes. However, the models used in this study are all stoichiometric. Therefore, caution should be exercised when extrapolating these findings to significantly sub-stoichiometric compositions unless similar structural features are present.

The structural ensembles generated in this study align well with previously reported experimental data.^43,47,49 Through deconvolution of the PDF, with residuals from reported fits for both lithiated and delithiated samples, showing good agreement with the Li structural features given by the superimposition of Fig. 8a and b. These features include a broad peak between 1.9 Å to 2.2 Å, a sharp peak with negative skew at 2.5 Å, a secondary peak at 2.8 Å, and a broad peak centered around 3.0 Å. The presence of these features in both lithiated and delithiated (along with their absence in the pristine) samples suggests a multi-stage matrix-forming process. Initially, Li incorporation populates intrinsic electron trap states in the network; once these states are filled, a shift to induced trapping occurs. In both cases, the formation of bi-polaron states results in Si–Si rich regions, leading to irreversible changes in the amorphous network. This behavior is in agreement with the experimental EDX which shows a clustering and redistribution of Si post lithiation.⁴⁷ As additional Li is incorporated, the final matrix forms via a number of intermediate compositions, as illustrated schematically in Fig. 9.

Fig. 9 Schematic representation of the initial Li-incorporation, progressing via a number of reversible (polaron) and irreversible (bi-polaron) states.

The electronic structure of a-Si₃N₄ during Li incorporation is characterized by the formation of distinct polaron and bi-polaron states, which play an important role in determining the material's electronic and electrochemical properties. Initially, single polaron states dominate, with occupied states close to the VBM (0 eV to 0.8 eV), and empty states close to the CBM (0 eV to −0.7 eV). These give way to form bi-polaron configurations with occupied states in the valence band. As the Li-concentration further increases, these bi-polaron states increase in energy, becoming near-degenerate with the VBM. This evolution progresses from isolated traps to an interconnected polaron networks underscoring the complex interplay between Li incorporation, local environment, network relaxation, and the resulting electronic properties. It is this coupling that resolves the apparent contradiction of a matrix that is simultaneously redox-active and inert: polaron states enable reversible, electrochemically active behavior, while bi-polaron states drive irreversible structural reorganization and are electronically inactive (Fig. 9).

Finally, the parallels between a-Si₃N₄, a-SiO₂ and other sub-stoichiometric oxides are noteworthy. Both exhibit intrinsic charge trapping, leading to electron polaron and bi-polaron formation,^95,96 and have similar Li incorporation geometries. This suggests a potentially general mechanism for conversion-type anodes based on oxides and nitrides, where intrinsic charge trapping drives bi-polaron formation, causing irreversible modifications to the amorphous network and initiating the matrix formation process. Extending these findings beyond a-Si₃N₄ is currently under investigation, however, the observed parallels provide a strong foundation for further study.

5 Conclusions

This study provides important insights into the structural and electronic behavior of a-Si₃N₄ during the initial Li incorporation, highlighting the formation and evolution of polaron and bi-polaron states. These states play a pivotal role in determining the material's electrochemical properties, particularly in the context of lithium-ion battery anodes. The transition from single polaron to bi-polaron states and the associated changes in structure are vital in generating the initial Si-rich regions of the network. Our findings suggest that managing the distribution and energy of these trap states could be key to optimizing a-Si₃N₄-based anodes for improved capacity and cycling stability.

Conflicts of interest

There are no conflicts to declare.

Data availability

All structures generated in this study are available via https://www.zenodo.org https://doi.org/10.5281/zenodo.15297234, together with input parameters and specifications for the compilation of the CP2K package.

The supporting information contains Li-incorporation in β-Si₃N₄, properties of the pristine cell, and characterisation of Li incorporation sites. See DOI: https://doi.org/10.1039/d5ta04013b.

Acknowledgements

This work was conducted at the Advanced Research Center for Nanolithography, a public–private partnership between the University of Amsterdam (UvA), Vrije Universiteit Amsterdam (VU), Rijksuniversiteit Groningen (RUG), the Netherlands Organization for Scientific Research (NWO), and the semiconductor equipment manufacturer ASML. E. O. is grateful for a WISE Fellowship from the NWO. J. M. and E. O. acknowledge support via Holland High Tech through a public–private partnership in research and development within the Dutch top sector of High-Tech Systems and Materials (HTSM). The use of the national computer facilities in this research was subsidized by NWO Domain Science (grant no. 2024.029) and supported by the SURF Cooperative (grant no. EINF-5451 and EINF-8981). The authors thank SURF (https://www.surf.nl) for the support in using the National Supercomputer Snellius.

References

1. J. B. Goodenough and M. H. Braga, Dalton Trans., 2018, 47, 645–648.
2. H. Au, M. Crespo-Ribadeneyra and M.-M. Titirici, One Earth, 2022, 5, 207–211.
3. J. B. Goodenough, Energy Environ. Sci., 2014, 7, 14–18.
4. Y. Li, Y. Lu, P. Adelhelm, M.-M. Titirici and Y.-S. Hu, Chem. Soc. Rev., 2019, 48, 4655–4687.
5. E. Olsson, J. Yu, H. Zhang, H. Cheng and Q. Cai, Adv. Energy Mater., 2022, 12, 2200662.
6. M. S. Dresselhaus and G. Dresselhaus, Adv. Phys., 2002, 51, 1–186.
7. E. Olsson, J. Cottom, H. Au, M.-M. Titirici and Q. Cai, Carbon, 2021, 177, 226–243.
8. E. Olsson, J. Cottom, H. Au, Z. Guo, A. C. S. Jensen, H. Alptekin, A. J. Drew, M.-M. Titirici and Q. Cai, Adv. Funct. Mater., 2020, 30, 1908209.
9. M. Noel and R. Santhanam, J. Power Sources, 1998, 72, 53–65.
10. X. Yu, S. Chen, B. Tang, X. L. Li, J. Zhou, Y. Ren, J. Wei, C. Yang, Y. Guo, Z. Zhou and S. H. Bo, ACS Energy Lett., 2024, 9, 4441–4449.
11. D. Hua Wu and Z. Zhou, Front. Phys., 2011, 6, 197–203.
12. L. Ma, J. Wu, G. Zhu, Y. Lv, Y. Zhang and H. Pang, J. Mater. Chem. A, 2021, 9, 5232–5257.
13. M. N. Obrovac and V. L. Chevrier, Chem. Rev., 2014, 114, 11444–11502.
14. H. Huo and J. Janek, ACS Energy Lett., 2022, 7, 4005–4016.
15. Y. Zhang, N. Du and D. Yang, Nanoscale, 2019, 11, 19086–19104.
16. D. E. Galvez-Aranda, A. Verma, K. Hankins, J. M. Seminario, P. P. Mukherjee and P. B. Balbuena, J. Power Sources, 2019, 419, 208–218.
17. S. Bai, W. Bao, K. Qian, B. Han, W. Li, B. Sayahpour, B. Sreenarayanan, D. H. S. Tan, S.-Y. Ham, Y. S. Meng, S. Bai, B. Sayahpour, S.-Y. Ham, Y. S. Meng, W. Bao, K. Qian, B. Han, W. Li, B. Sreenarayanan and D. H. Tan, Adv. Energy Mater., 2023, 13, 2301041.
18. S. Y. Ham, E. Sebti, A. Cronk, T. Pennebaker, G. Deysher, Y. T. Chen, J. A. S. Oh, J. B. Lee, M. S. Song, P. Ridley, D. H. Tan, R. J. Clément, J. Jang and Y. S. Meng, Nat. Commun., 2024, 15, 1–9.
19. C.-Y. Wu, C.-C. Chang and J.-G. Duh, J. Power Sources, 2016, 325, 64–70.
20. M. N. Obrovac and L. Christensen, Electrochem. Solid-State Lett., 2004, 7, A93–A96.
21. P. Li, G. Zhao, X. Zheng, X. Xu, C. Yao, W. Sun and S. X. Dou, Energy Storage Mater., 2018, 15, 422–446.
22. H. Huo, M. Jiang, Y. Bai, S. Ahmed, K. Volz, H. Hartmann, A. Henss, C. V. Singh, D. Raabe and J. Janek, Nat. Mater., 2024, 23, 543–551.
23. A. Ulvestad, H. F. Andersen, J. P. Mæhlen, Ø. Prytz and M. Kirkengen, Sci. Rep., 2017, 7, 13315.
24. C. Cao, I. I. Abate, E. Sivonxay, B. Shyam, C. Jia, B. Moritz, T. P. Devereaux, K. A. Persson, H. G. Steinrück and M. F. Toney, Joule, 2019, 3, 762–781.
25. A. Okada, J. Eur. Ceram. Soc., 2008, 28, 1097–1104.
26. V. Y. Doo, D. R. Nichols and G. A. Silvey, J. Electrochem. Soc., 1966, 113, 1279.
27. S.-J. Tsai, C.-L. Wang, H.-C. Lee, C.-Y. Lin, J.-W. Chen, H.-W. Shiu, L.-Y. Chang, H.-T. Hsueh, H.-Y. Chen and J.-Y. Tsai, et al., Sci. Rep., 2016, 6, 28326.
28. V. Gritsenko, V. Kruchinin, I. Prosvirin, Y. N. Novikov, A. Chin and V. Volodin, J. Exp. Theor. Phys., 2019, 129, 924–934.
29. P. T. Törmä, H. J. Sipilä, M. Mattila, P. Kostamo, J. Kostamo, E. Kostamo, H. Lipsanen, N. Nelms, B. Shortt, M. Bavdaz and C. Laubis, IEEE Trans. Nuc. Sci., 2013, 60, 1311–1314.
30. P. T. Törmä, J. Kostamo, H. Sipilä, M. Mattila, P. Kostamo, E. Kostamo, H. Lipsanen, C. Laubis, F. Scholze, N. Nelms, B. Shortt and M. Bavdaz, IEEE Trans. Nuc. Sci., 2014, 61, 695–699.
31. S. Cornaby and D. H. Bilderback, J. Synchrotron Radiat., 2008, 15, 371–373.
32. T. Sharma, J. Wang, B. K. Kaushik, Z. Cheng, R. Kumar, Z. Wei and X. Li, IEEE Access, 2020, 8, 195436–195446.
33. C. Xiang, W. Jin and J. E. Bowers, Photonics Res., 2022, 10, A82–A96.
34. T. Aiyama, T. Fukunaga, K. Niihara, T. Hirai and K. Suzuki, J. Non-Cryst. Solids, 1979, 33, 131–139.
35. M. Misawa, T. Fukunaga, K. Niihara, T. Hirai and K. Suzuki, J. Non-Cryst. Solids, 1979, 34, 313–321.
36. K. W. K. Wakita, H. H. H. Hayashi and Y. N. Y. Nakayama, Jpn. J. Appl. Phys., 1996, 35, 2557.
37. S. V. Deshpande, E. Gulari, S. W. Brown and S. C. Rand, J. Appl. Phys., 1995, 77, 6534–6541.
38. B. S. Sahu, F. Delachat, A. Slaoui, M. Carrada, G. Ferblantier and D. Muller, Nanoscale Res. Lett., 2011, 6, 1–10.
39. R. C. de Guzman, J. Yang, M. Ming-Cheng Cheng, S. O. Salley and K. Y. S. Ng, J. Mater. Chem. A, 2014, 2, 14577–14584.
40. X. Huang, X. Gan, F. Zhang, Q. Huang and J. Yang, Electrochim. Acta, 2018, 268, 241–247.
41. J. Yang, R. C. de Guzman, S. O. Salley, K. S. Ng, B.-H. Chen and M. M.-C. Cheng, J. Power Sources, 2014, 269, 520–525.
42. H. Cheng, D. Li, B. Xu, Y. Wei, H. Wang, B. Jiang, X. Liu, H. Xu and Y. Huang, Energy Storage Mater., 2022, 53, 305–314.
43. S. O. Kilian, B. Wankmiller, A. M. Sybrecht, J. Twellmann, M. R. Hansen and H. Wiggers, Adv. Mater. Interfaces, 2022, 9, 2201389.
44. D. Pandel, J. Neises, S. O. Kilian, H. Wiggers and N. Benson, ACS Appl. Energy Mater., 2022, 5, 14807–14814.
45. A. Ulvestad, J. P. Mæhlen and M. Kirkengen, J. Power Sources, 2018, 399, 414–421.
46. A. Ulvestad, H. F. Andersen, I. J. Jensen, T. T. Mongstad, J. P. Mæhlen, Ø. Prytz and M. Kirkengen, Sci. Rep., 2018, 8, 8634.
47. A. Ulvestad, M. O. Skare, C. E. Foss, H. Krogsæter, J. F. Reichstein, T. J. Preston, J. P. Mæhlen, H. F. Andersen and A. Y. Koposov, ACS Nano, 2021, 15, 16777–16787.
48. A. W. Nemaga, S. Y. Lai, T. Nguyen, A. Ulvestad and C. E. L. Foss, ACS Omega, 2025, 10, 2608–2615.
49. A. J. Lovett, M. Füredi, L. Bird, S. Said, B. Frost, P. R. Shearing, S. Guldin and T. S. Miller, ACS Electrochem., 2025, 1, 962–973.
50. L. Hückmann, J. Cottom and J. Meyer, Adv. Phys. Res., 2024, 3, 2300109.
51. J. Cottom, L. Hückmann, E. Olsson and J. Meyer, J. Phys. Chem. Lett., 2024, 15, 840–848.
52. T. D. Kühne, M. Iannuzzi, M. Del Ben, V. V. Rybkin, P. Seewald, F. Stein, T. Laino, R. Z. Khaliullin, O. Schütt, F. Schiffmann, D. Golze, J. Wilhelm, S. Chulkov, M. H. Bani-Hashemian, V. Weber, U. Borštnik, M. Taillefumier, A. S. Jakobovits, A. Lazzaro, H. Pabst, T. Müller, R. Schade, M. Guidon, S. Andermatt, N. Holmberg, G. K. Schenter, A. Hehn, A. Bussy, F. Belleflamme, G. Tabacchi, A. Glöβ, M. Lass, I. Bethune, C. J. Mundy, C. Plessl, M. Watkins, J. VandeVondele, M. Krack and J. Hutter, J. Chem. Phys., 2020, 152, 194103.
53. J. Heyd, G. E. Scuseria and M. Ernzerhof, J. Chem. Phys., 2003, 118, 8207–8215.
54. J. Heyd, G. E. Scuseria and M. Ernzerhof, J. Chem. Phys., 2006, 124, 219906.
55. M. Guidon, J. Hutter and J. VandeVondele, J. Chem. Theory Comput., 2010, 6, 2348–2364.
56. J. VandeVondele and J. Hutter, J. Chem. Phys., 2007, 127, 114105.
57. S. Goedecker, M. Teter and J. Hutter, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54, 1703–1710.
58. C. Hartwigsen, S. Goedecker and J. Hutter, Phys. Rev. B: Condens. Matter Mater. Phys., 1998, 58, 3641–3662.
59. M. Krack, Theor. Chem. Acc., 2005, 114, 145–152.
60. D. R. J. Boyd, J. Chem. Phys., 1955, 23, 922–926.
61. R. Fletcher, Comput. J., 1970, 13, 317–322.
62. D. Goldfarb, Math. Comput., 1970, 24, 23.
63. D. F. Shanno, Math. Comput., 1970, 24, 647.
64. S. B. Zhang and J. E. Northrup, Phys. Rev. Lett., 1991, 67, 2339–2342.
65. X. Fan, X. Ji, L. Chen, J. Chen, T. Deng, F. Han, J. Yue, N. Piao, R. Wang, X. Zhou, X. Xiao, L. Chen and C. Wang, Nat. Energy, 2019, 4, 882–890.
66. C. M. Marian, M. Gastreich and J. D. Gale, Phys. Rev. B: Condens. Matter Mater. Phys., 2000, 62, 3117–3124.
67. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865–3868.
68. J. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1997, 77, 3865–3868.
69. T. W. Anderson, An Introduction to Multivariate Statistical Analysis, Wiley, New York, 1958, vol. 2.
70. G. J. McLachlan and D. Peel, Finite Mixture Models, Wiley, New York, 2000.
71. G. J. McLachlan, S. X. Lee and S. I. Rathnayake, Annu. Rev. Stat. Appl., 2019, 6, 355–378.
72. A. P. Dempster, N. M. Laird and D. B. Rubin, Stat. Methodol., 1977, 39, 1–22.
73. S. K. Ng, T. Krishnan and G. J. McLachlan, in The EM Algorithm, ed. J. E. Gentle, W. K. Härdle and Y. Mori, Springer Berlin Heidelberg, Berlin, Heidelberg, 2012, pp. 139–172.
74. G. Schwarz, Ann. Stat., 1978, 6, 461–464.
75. A. A. Neath and J. E. Cavanaugh, WIREs Comp. Stat., 2012, 4, 199–203.
76. H. Akaike, IEEE Trans. Autom. Control, 1974, 19, 716–723.
77. H. Akaike, in Information Theory and an Extension of the Maximum Likelihood Principle, ed. E. Parzen, K. Tanabe and G. Kitagawa, Springer New York, New York, NY, 1998, pp. 199–213.
78. P. G. Debenedetti and F. H. Stillinger, Nature, 2001, 410, 259–267.
79. D. J. Wales, Annu. Rev. Phys. Chem., 2018, 69, 401–425.
80. J. C. Mauro and M. M. Smedskjaer, J. Non-Cryst. Solids, 2014, 396–397, 41–53.
81. W. Kauzmann, Chem. Rev., 1948, 43, 219–256.
82. R. Palmer, Adv. Phys., 1982, 31, 669–735.
83. P. K. Roy and A. Heuer, J. Chem. Phys., 2022, 157, 174506.
84. C. Scalliet, L. Berthier and F. Zamponi, Nat. Commun., 2019, 10, 5102.
85. A. H. Larsen, J. J. Mortensen, J. Blomqvist, I. E. Castelli, R. Christensen, M. Dułak, J. Friis, M. N. Groves, B. Hammer, C. Hargus, E. D. Hermes, P. C. Jennings, P. B. Jensen, J. Kermode, J. R. Kitchin, E. L. Kolsbjerg, J. Kubal, K. Kaasbjerg, S. Lysgaard, J. B. Maronsson, T. Maxson, T. Olsen, L. Pastewka, A. Peterson, C. Rostgaard, J. Schiøtz, O. Schütt, M. Strange, K. S. Thygesen, T. Vegge, L. Vilhelmsen, M. Walter, Z. Zeng and K. W. Jacobsen, J. Phys.: Condens. Matter, 2017, 29, 273002.
86. P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, S. J. van der Walt, M. Brett, J. Wilson, K. J. Millman, N. Mayorov, A. R. J. Nelson, E. Jones, R. Kern, E. Larson, C. J. Carey, İ. Polat, Y. Feng, E. W. Moore, J. VanderPlas, D. Laxalde, J. Perktold, R. Cimrman, I. Henriksen, E. A. Quintero, C. R. Harris, A. M. Archibald, A. H. Ribeiro, F. Pedregosa and P. van Mulbregt, Nat. Meth., 2020, 17, 261–272.
87. F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot and E. Duchesnay, J. Mach. Learn. Res., 2011, 12, 2825–2830.
88. J. D. Hunter, Comput. Sci. Eng., 2007, 9, 90–95.
89. W. Humphrey, A. Dalke and K. Schulten, J. Mol. Graphics, 1996, 14, 33–38.
90. A. Stukowski, Modell. Simul. Mater. Sci. Eng., 2009, 18, 015012.
91. J. Robertson, Appl. Phys. Lett., 1991, 59, 3425–3427.
92. D. T. Krick, P. M. Lenahan and J. Kanicki, Phys. Rev. B: Condens. Matter Mater. Phys., 1988, 38, 8226–8229.
93. W. L. Warren, J. Kanicki, J. Robertson, E. H. Poindexter and P. J. McWhorter, J. Appl. Phys., 1993, 74, 4034–4046.
94. Y. N. Novikov and V. Gritsenko, J. Non-Cryst. Solids, 2020, 544, 120186.
95. A.-M. El-Sayed, M. B. Watkins, T. Grasser, V. V. Afanas’ev and A. L. Shluger, Phys. Rev. Lett., 2015, 114, 115503.
96. A.-M. El-Sayed, Y. Wimmer, W. Goes, T. Grasser, V. V. Afanas’ev and A. L. Shluger, Phys. Rev. B: Condens. Matter Mater. Phys., 2015, 92, 014107.

---