Chemo-mechanical failure mechanisms of the silicon anode in solid-state batteries

Silicon is a promising anode material due to its high theoretical specific capacity, low lithiation potential and low lithium dendrite risk. Yet, the electrochemical performance of silicon anodes in solid-state batteries is still poor (for example, low actual specific capacity and fast capacity decay), hindering practical applications. Here the chemo-mechanical failure mechanisms of composite Si/Li6PS5Cl and solid-electrolyte-free silicon anodes are revealed by combining structural and chemical characterizations with theoretical simulations. The growth of the solid electrolyte interphase at the Si|Li6PS5Cl interface causes severe resistance increase in composite anodes, explaining their fast capacity decay. Solid-electrolyte-free silicon anodes show sufficient ionic and electronic conductivities, enabling a high specific capacity. However, microscale void formation during delithiation causes larger mechanical stress at the two-dimensional interfaces of these anodes than in composite anodes. Understanding these chemo-mechanical failure mechanisms of different anode architectures and the role of interphase formation helps to provide guidelines for the design of improved electrode materials.


Supplementary information: Computational details
The atomistic simulations were conducted using the density functional theory (DFT) framework as implemented in the Vienna Ab Initio Simulation Package (VASP).The projector augmented-wave pseudopotentials were used to describe the interaction between ions and electrons, and the exchange-correlation effects were treated using the Perdew-Burke-Ernzerhof (PBE) functional under the generalized gradient approximation (GGA). 1 Herein, the electronic configurations for the PAW potentials were 1s 2 2p 0 for Li, 3s 2 3p 3 for P, 3s 2 3p 4 for S, 3s 2 3p 5 for Cl and 3s 2 3p 2 for Si.The VESTA package was used to visualize the various structures. 2

Stable structure of LPSCl
Li atoms partially occupy the 48h site in LPSCl, and LPSCl exhibits anti-site S/Cl disorder. 3The Supercell code was utilized to tackle the issue of partial occupancy and S/Cl disorder.5][6] An exhaustive search over all possible configurations performed by the Supercell code for medium-size systems makes it possible to describe in-depth the distribution of local composition and resulting local geometrical distortions, with random as well as non-random types of disorder. 7After enumerating the possible structures, Coulomb energy calculations were additionally performed to sort out 10 structures for LPSCl.The structure optimizations based on DFT calculations were then performed, and the structure with the lowest energy was selected for the subsequent calculations.The most stable LPSCl structure was screened in our previous study, 8 which was directly selected to calculate the adsorption energy.

Adsorption of Li atoms on Si and LPSCl surfaces
The most stable (100) and (111) surfaces were selected for LPSCl and Si in the subsequent calculations, respectively.The vacuum space in the z direction was about 15 Å to avoid interaction between neighboring slabs.The energy cut-off is set to 500 eV for the plane-wave expansion of the projector-augmented waves in the self-consistent calculations, whereas the Gamma grid of 1 × 1 × 1 was used.Each system was fully relaxed with residual forces smaller than 0.01 eV Å −1 and the total-energy change was less than 10 −5 eV during the structural optimization.The binding energy of a lithium atom on a surface was determined using the expression,  B =  Li +  surface −  total (1) where  B ,  Li ,  surface and  total represent the binding energy, the energy of a Li atom, of the Si (111) and LPSCl (100) surfaces, and of the adsorbate, respectively.

Melt-and-quench simulations
The "melt-and-quench" simulation scheme was conducted by carrying out a series of ab-initio molecular dynamics (AIMD) simulations. 9The relaxed supercell equilibrated at 300 K for 2 ps, and was then heated to a temperature (2100 K) in 5 ps.The system was allowed to equilibrate at 2100 K for 2000 steps (each step = 1 fs).The process was then reversed to quench the system back to 300 K at rate of 300 K ps −1 .At each step, the structure obtained from the previous AIMD simulation was used as the starting point for the next one.For each AIMD simulation, the system was given 2 ps to reach its thermal equilibrium in order to eliminate its correlation to the previous structure.

AIMD simulations for ionic conductivities of LixSi alloys
The diffusivity of Li + ions in amorphous Li-Si alloys was determined based on AIMD simulations with non-spin polarized calculations.Because AIMD simulations are computationally expensive, a compromise was made between computation efficiency and accuracy, resulting in the use of a kinetic energy cut-off of 300 eV and a gamma-point sampling. 10The initial temperature of the amorphous Li-Si alloys was set to be 100 K, and the velocity of ions was set according to the Boltzmann distribution.Then, the samples were heated to the desired temperature (500 to 1200 K) using a velocity scaling thermostat at a heating rate of 1 K fs −1 .The total simulation time was set to 30 ps.The integration of Newton's equation was treated based on the Verlet algorithm, as implemented in VASP. 1 At the assigned temperature, the AIMD simulation was performed in the NVT ensemble 11 in conjunction with a Nose-Hoover thermostat. 11,12The diffusion coefficient of the Li + ions was calculated by the mean square displacement (MSD) over time as where d is the dimensionality of diffusion, N is the total number of diffusing ions,   () is the displacement of the i-th ion at time t.The activation energy (Ea) is determined from Arrhenius plots of the diffusion coefficient.The Li + ion conductivity was then calculated according to the Nernst-Einstein relationship, where n is the charge carrier concentration of mobile ions as number density, q is the ionic charge (i.e., e0), kB is the Boltzmann constant and T is the temperature. 10

BoltzTrap calculations for electronic conductivities of LixSi alloys
BoltzTraP is a modern implementation of the smooth Fourier interpolation algorithm for electronic bands, which forms the basis of the original and widely used BoltzTraP code. 13The electron conductivity was calculated using BoltzTrap, which was carried out through Fourier expansion of energy bands and special function processing to maintain the symmetry of the space group.The amorphous structures of LiSi3, LiSi, Li12Si7, Li13Si4, Li7Si2, Li13Si4, Li7Si2 and Li15Si4 were all selected for these calculations.Although the relaxation time τ can be defined by estimating the electron-phonon coupling effect, the calculations of the full electron-phonon interactions are complicated.The carrier (electron) relaxation time has been demonstrated to be ~ 10 −14 s in the Materials Project. 14The electron concentration has been also determined from the Bader charge analysis, which has been widely used in calculating the charge transfer in various materials.The electron concentration (n) was calculated by: Where n is the electronic concentration,  is the total transferred charge and  is the volume of the amorphous structure.

Supplementary information: Chemo-mechanical phase-field model
The schematic in Supplementary Fig. 20 shows the model-type electrode for the chemo-mechanical phase-field model, consisting of the silicon anode and the LPSCl electrolyte, with the chemical reaction for the charge/discharge process taking place at the interface between the Si and the SE.The insertion of lithium into the Si causes a volume change, leading to the emergence of stress and elasto-plastic deformation.Hence, the free energy landscape for the chemo-mechanically coupled electrode scenario can be expressed as follows: [15][16][17] where  che represents the chemical free energy that accounts for the energy contribution due to the chemical reaction, i.e., the lithiation/delithiation reaction. ela is the elastic strain energy due to the elastic lattice deformation that originates from the injection and removal of the lithium, and  d denotes the damage density energy for the crack formation and propagation.Given that Li + ions diffuse through the SE phase and are inserted (together with the necessary number of electrons) into the anode, then their contribution to the chemical free energy can be defined as follows: [18][19][20] with R and T being the gas constant and temperature, and  max is the maximum concentration of lithium stored in the Si anode, with the unit of mol/ 3 . Li and  Li + are the normalized concentrations (normalized by C max ) of lithium in the Si anode and Li + ions in the SE, respectively, which are considered to assume an ideal solution state as reflected by Eq. 6 for the configurational entropy.Moreover, it is assumed that the electrons are continuously provided throughout the reaction, and thus their energy contribution is disregarded in this investigation.The potential change resulting from the change in electric neutrality due to the movement of anions and cations is not within the scope of this study and will be addressed in future research.The elastic free energy is given by: where  is the Cauchy stress tensor and  is the identity tensor,  and  e are the total strain and the elastic strain, respectively.  is the concentration of the i-th species, and  ref is its reference concentration, Ω is the partial molar volume.ℂ is the elasticity tensor which depends on both the concentration of the species as well as the damage state variable .The damage density energy is read as follows: 21 where   is the fracture energy, and  0 is the length scale parameter of the crack.In this work,  = 1 represents the fully damaged case and  = 0 denotes the unbroken fully coherent case, respectively.
The equations that govern the diffusion of each species, the linear momentum balance, and the propagation of cracks are expressed as follows: where   is the normalized concentration (divided by   )of i-th species, and its chemical potential is given by   = , which contains the contribution from both the diffusion and mechanical deformation.  is the diffusion coefficient of i-th species.R is the gas constant, and T is the temperature.ṙ is the reaction source term, and   is the stoichiometric coefficient of i-th species.
An isotropic J2 plasticity model is employed for the inelastic plastic response of system, with the yield condition and plastic flow expressed as follows: where   and   represent the yield stress and equivalent plastic strain, respectively. is the unknown plastic multiplier to be solved during the simulation.[] is the deviatoric part of the Cauchy stress.

The crack propagation is governed by
where  is the mobility for the crack propagation, and the variational derivative is defined as: In this work, the kinetics of the lithium insertion and extraction is given by the Butler-Volmer equation as follows: with  0 being the reaction coefficient, and  being the overpotential.The parameters for the phase-field model are listed in Supplementary Table 8.
It is important to note that the presented chemo-mechanical phase-field model is a general model, and can be applied to both 2D and 3D modeling.However, in this study, we have focused on the use of 2D geometry to examine the differences between different types of interfaces, such as the 2D and 3D interfaces.
We like to note that real-world interfaces are rarely true 2D interfaces, and that one might consider them rather as quasi-2D interfaces -in particular when it comes to SEI formation.Distinguishing "2D" from "3D", we like to highlight that a SE-free Si anode shows only one more or less planar SEI layer between Si and SE, while the Si composite anode shows an interconnected 3D network of the SEI.As Supplementary Fig. 1a shows, the SEI grows on the SE-free Si anode as a thin layer into the LPSCl separator, forming a "2D SEI" during cycling.Thus, with "2D" we like to note the planar interphase character of the SEI.In contrast, the "3D SEI" describes the SEI grown in an interconnected 3D network as part of the Si/LPSCl/C composite (Supplementary Fig. 1b).
shows mixed conduction.It does not only conduct electrons, but also conducts ions together with the LPSCl phase.The resulting complex electronic/ionic paths are not well described by the TLM.Therefore, a simplified model in Supplementary Fig. 7b was

Supplementary Note 3
The partial ionic/electronic conductivities obtained by DFT simulations are consistent with the  ̃Li values and the thermodynamic factor obtained by GITT measurements, further confirming the accuracy of the  ̃Li values from the GITT measurements. ̃ can be calculated by the following equations: Li is the lithium (component) diffusion coefficient, i.e., the self-diffusion coefficient of the neutral component Li.  is the thermodynamic factor.  + and   -are the partial ionic and electronic conductivities of the considered phase.R, T, and F are the gas constant, temperature, Faraday's constant, respectively.We can safely assume for the Li-Si alloys that where  Li + ,  e -, and  Li are the concentrations of Li + , e -, and Li (i.e., the neutral component Li 0 ) in the LixSi alloy.
The average value of  ̃Li for a given concentration range results as where  ̅ = 3.0,  ̅ Li + = 1.5×10 -3 S cm -1 , and  ̅ e -= 4.4×10 -4 S cm -1 , respectively (data taken from results in Fig. 3d, 3e, and Supplementary Fig. 19c).We also assume where  m and  Li are molar volume and molar number of Li2Si, respectively.Then, which is in excellent agreement with the  ̅Li = 1.0×10 -8 cm 2 s -1 as obtained from the GITT measurements.We note that this analysis of chemical diffusion of Si, including the comparison of theoretical and experimental data has not been demonstrated before to the best of our knowledge.However, it is well in line with the classical analysis of chemical diffusion in other lithium alloys (e.g.Li3Sb, and Li3Bi), as originally presented by Weppner et al. 23 chemical diffusion and storage) and mechanical processes (elasto-plastic deformation).
To quantitatively understand the different elasto-plastic expansion/shrinkage of Si and the fracture behavior of LPSCl in 3D composites and at 2D LPSCl|Si interfaces, a fully coupled chemo-mechanical phase-field fracture model was developed to simulate the 1 st lithiation and delithiation steps (Supplementary Fig. 20).The model is designed to be versatile and capable of modeling both 2D (circular particles) and 3D (spherical particles) systems.The 2D circular particle model has been used to demonstrate the differences between interfaces in the 3D composite and the case of the planar interface in the SE-free anode.In this model, both the Si and SE particles undergo elasto-plastic deformation, i.e., they exhibit distinct yield stresses and inelastic mechanical response together with linear elastic properties (i.e., Young's modulus).For further details, please refer to our supplementary material.The 3D interface, as shown in Supplementary Fig. 21a and Fig. 21b, exhibits a smooth connection between Si and SE.The smooth interface contact (the smooth circular surface) between the particle and SE is the reason behind the uniform distribution of stress (ranging from -0.2 to 0.2 GPa during lithiation and -0.1 to 0.1 GPa during delithiation) within and around the particle, thereby preventing the formation of noticeable cracks.This confirms that the lithiation process and the associated volume expansion of Si results in the disappearance of interface voids due to a "compression" effect and densification of the LixSi microstructure.Consequently, the lithium distribution within the Si particle is almost uniform throughout the (de)lithiation process (the lithium uniformly increases or decreases within the Si particle), with no significant change in the lithium concentration gradient from its core to the surface being captured, as the SoC increases from 5% to 60% and subsequently decreases from 60% to 10%, as shown in Supplementary Fig. 21a and Fig. 21b.Furthermore, during lithiation, the positive peak value (tensile stress around 0.2 GPa) of the maximum principal stress ( sp 1 ) is primarily observed at the surface of the particle due to the lithium insertion induced volume expansion.In contrast, during delithiation, the top and bottom of the particle are subjected to tensile stress around 0.1 GPa caused by the external pressure (compression from SE).In addition, a section of the SE matrix (located at the top and bottom parts attached to the particle) undergoes plastic deformation, resulting in an equivalent plastic strain of up to 1%.Throughout the entire cycle of lithiation and delithiation, cracks are not prominently visible (the damage percentage d ≪ 100%) as demonstrated in Supplementary Fig. 21a and Fig. 21b, owing to the low levels of stress.This indicates a good interconnection within the LixSi microstructure, as observed in Fig. 4.
The geometry of the 2D interface depicted in Fig. 5 creates stress concentration (due to the geometry singularity) features with stress peaks of up to 300 MPa (0.3 GPa) already upon the first lithiation at the anode/SE interface.The process of lithiation involves the insertion of lithium into the Si particle (the SoC reaches 100% at the end of the lithiation process), causes large volume expansion and exacerbates the stress at the singular interface, as shown in Supplementary Fig. 21c.As lithiation progresses from 5% to 100% (see Supplementary Fig. 21a), the principal stress  sp 1 , which builds up due to the growing misfit at the interface between the adjacent solid phases (i.e.LPSCl SE), reaches a range of approximately (-0.3~0.5 GPa).This value is higher than that observed in the 3D composite and can initiate the development of cracks and plastic flow of the lithiated Si.We expect that a stronger densification effect will occur at the 2D interface.Subsequently, upon delithiation, the removal of lithium (from 100% to 0%) causes the particles to shrink, as shown in Supplementary Fig. 21d.This, in turn, augments the stress at the 2D interface, ultimately leading to the propagation of cracks from the LPSCl into the Si particle, as depicted in Supplementary Fig. 21d.This further confirms the submicron cracks as observed in Fig. 4. The occurrence of fractures leads to a decrease in stress (stress relaxation) in the influenced zone (with a zero stiffness) as the cracks progress and stress levels are lowered to a zero-stress state.This results in a significant difference in stress levels between the damaged region and the undamaged domain.Towards the end of the delithiation process, the plastic deformation is further accumulated, therefore,  sp 1 can climb up to -0.3~0.8GPa, which is higher than that observed during the lithiation process.As a result, several large cracks are observed at the end of the delithiation process in the 2D interface, which is the primary cause of the large voids (around 2μm) depicted in Fig. 4. Furthermore, the surface of the Si particle is subjected to a negative principal stress (compressive stress indicated by the blue region), whereas the delithiation process can cause a positive principal stress (tensile stress) at the surface of the particle.Hence, as a result of the lithiation process, a considerable quantity of plastic strain is accumulated at the interface, which can be further increased by the delithiation process, as depicted in the bottom column of Supplementary Fig. 21c and Fig. 21d, by approximately 10%.This could be a significant contributing factor to the relatively stable porosity of the SE-free Si anode (~26.9%) after the first delithiation, as depicted in Fig. 4.Moreover, the value of the plastic strain is significantly greater than the value observed in the 3D composite (1%).
Based on the findings in this section, it can be inferred that the 2D Si/LPSCl interface geometry shows stress accumulation, leading to the growth of cracks (voids) and undermines the cycling stability.

Supplementary Fig. 8 .
used, where Rion and Rel are combined into Rcomp.Since we focus on the Rint change to study the SEI growth rate, the simplified model works well for both 3D and 2D Si/LPSCl interfaces.Comparison of the binding energy.Schematic structures of (a) bulk Si and (b) bulk LPSCl.(c) Sketch of the model for the calculation of the binding energy between Si (111) and Li.(d) Sketch of the model for the calculation of the binding energy between LPSCl (100) and Li.Supplementary Fig. 9. Specific capacities of different Si/LPSCl anodes.Lithiation curves of different Si/LPSCl anodes at 0.1 C.

Table 1 .
EIS fitting results of the pristine Si/LPSCl composite based on the transmission line model.

Table 2 .
EIS fitting results of the Si/LPSCl composite during resting at the open-circuit voltage after being discharged to E = -0.6V, i.e., to E = 0.02 vs. Li + /Li.

Table 3 .
14ttice parameters of unit cells of crystalline LixSi alloys (materials ID from Materials Project).14

Table 4 .
Total Bader charge of Li atoms in amorphous structures along with the corresponding volume of different LixSi amorphous structures and electron concentration.Electron conductivity was obtained from BoltzTrap calculations.

Table 5 .
EIS fitting results of the SE-free Si anode during the relaxation of GITT.

Table 6 .
EIS fitting results of the SE-free Si anode during resting at the open-circuit voltage after being discharged to E = -0.6V, i.e., to E = 0.02 vs. Li + /Li.

Table 7 .
EIS fitting results of the SE-free Si anode after various cycles.

Table 8
Parameters for the phase-field model