Fundamental understanding of the fracture toughness of the Li xSi alloys is crucial for designing of Si based high- capacity and failure- resistant electrodes. In this study, molecul ar dynamics simulation informed continuum chemo- mechanical modeling s with conservation integral s were conducted to derive fracture toughness of Li xSi alloys . Our modeling results show reasonabl e agreement with available experimental data, reveal ing that the fracture toughness of Li xSi alloys with low lithium concentration does not vary significantly with lithium concentration. In addition, we demonstrated that, if lithium redistribution caused by the stress gradient around crack tip need s to be cons idered, a n appropriate chemo -mechanical path -independent J -integral should be used as the classic Rices J -integral is path -dependent. The obtained fracture toughness of the Li xSi alloys here 2 provides guidance for the rational design of Si based electrodes , and the presented approach also sheds light for the evaluation of the fracture toughness of other energy materials at different charging/discharging level s. Keywords : Lithium ion battery; Li xSi alloy; molecular dynamics simulation ; chemo -mechanical model; fracture toughness ; lithium redistribution 1. Introduction Lithium (Li) ion batteries (LIBs), due to their high energy density, have been dominantly used as power sources in portable e lectronics . 1-4 However, the relative low capacity of current ly commercialized LIBs has seriously hindered the ir extensive application in some emerging applications , such as electric vehicles and large energy grid.5- 7 It is well- known that s ilicon (Si) stands as one of the most promising anode material s owing to its highest theoretical capacity (4200 mAh/g) .8-10 Nevertheless, Li insertion/extraction during charging/discharging cyclings will inherently induce huge volume tric change and stress inside the Si anode, leading to fracture, pulverization, electrical disconnectivity, and ultimately capacity loss.11-18 Therefore, a fundamental understanding of the degradation mechanisms in the Si anode during lithiation/delithiation cycli ngs, especially the damage initiation and propagation, is crucial for the rational design 3 of next -generation failure -resistant Si based electrodes. Recent experimental studies evidenced that crystalline Si ( c-Si) undergoes ~300% volumetric change during lithiation/delithiation cycling s, leading to the fracture of Si nano - particles even in the first lithiation cycle when the particle size is larger than ~ 150 nm.11, 15-16, 19 For amorphous Si ( a-Si), the volumetric change is ~280%.20-21 However, the critical fracture size of a-Si particles can increase to ~870 nm .20 In order to inve stigate the fracture behavior of lithiated Si, Pharr et al. estimat ed the fracture energy of lithiated Si thin film by measuring the curvature changes of the film at different lithiation and delithiation stages , concluding that the fracture energy of lithi ated Si is almost independent of Li concentration as the obtained fracture toughness is 28.5 4.3 J m± at low Li concentration (~Li 0.7Si) and has the bounds of 25.4 2.2 J m± to 26.9 1.9 J m± at high Li concent ration (~Li 2.8Si).22 However, the nanoindentation based measurements conducted by Wang et al. revealed a highly Li concentration de pendent fracture behavior of lithiated Si, namely , with the increment of Li concentration, the fracture toughness first decreases from 22.85 0.15 J m± for unlithiated a-Si to a minimum value at Li 0.31Si, and then increases substantially, reaching to 28.54 0.72 J m ± for Li 1.09Si, which indicates a brittle -to-ductile transition of fracture with the increase of Li concentration in Si.23 In consideration of the discrepancy between these tw o experimental studies, it is worthwhile to further investigate of the fracture toughness of LixSi as it ultimately dictate s the cyclability of Si based electrodes. Previous works have demonstrated that molecular dynamics (MD) simulation can be adopted as an appropriate method to computationally assess the fracture properties of 4 materials .24-28 For instance, a perfect match was found between the MD simulation results and experimental data for c-Si struc tures , indicating that MD simulation is capable to calculate the fracture toughness.24-26 However, as the lithiated Li xSi can undergo relatively large plastic deformation before fracture failure compared with unlithiated Si, attentions should be paid to the large scale yielding around crack tip when calculating fracture toughness.29-32 Therefore, i n this study, a MD simulation informed continuum chemo- mechanical model is presented for the calculation of fracture toughness of Li xSi alloys. In particular, MD simulations are conducted for the material properties of Li xSi alloys (i.e., Young’ s modulus, Poisson’ s ratio, and yield stress) and the critical fracture strain under which the crack in the pre -cracked specimen starts to propagate. Based on the obtained material properties and critical fracture strain, a continuum chemo -mechanical model is used to calculate the stress distribution in the same cracked specimen under the corresponding critical strain loading. With the stress distribution information, J -integral is directly calculated for the fracture toughness of Li xSi alloys . 2. Methodology 2.1. MD simulation In our MD simulation, the newly developed second nearest -neighbor m odified embedded atom method interatomic potential (2NN MEAM) of Li -Si alloy is adopted to account for the interactions among Li and Si atoms .33-34 This 2NN -MEAM interatomic potential has been adequately demonstrated to provide the primary properties and lattice 5 structures of both crystalline and amorphous Li xSi alloy successfully.33 Based on this potential, MD simulations are conducted for the ma terial properties of Li xSi alloys and the critical fracture strain under which the crack in the pre -cracked specimen starts to propagate. The atomic structures of different Li xSi alloys are created by melting -quenching method.35 First, large enough pure c -Si block is generated by using LAMMPS MD simulator. Second, according to the ratio x between the number of Li and Si atoms of different Li xSi alloys , Li atoms are randomly inserted into the Si block for the corresponding compositions. Then, with periodic boundary conditions (PBC s) prescribed in all three dimensi ons of the simulation cell and the time step set to be 1 fs, the temperature of the Li xSi block is increased from 0 K to 2500 K during the first 100 ps, and maintained there for another 400 ps, and afterward decreased to 10 K with quenching rate ~ 5.0 K/ps during the following 500 ps. Finally, the obtained fully amorphized Li xSi block is relaxed for another 500 ps using NPT ensemble. To understand the mechanical behaviors of Li xSi alloys , MD simulations at 10 K are conducted to obtain the uniaxial tensile stress -strain relationship s. Columns (30 × 30 × 80 Å in x, y, and z directions ) are cut from the fully relaxed LixSi blocks. Tensile simulations with NPT ensemble are carried out by applying a constant strain rate (75 10 1 s× ) along the axial direction (z direction) of the columns. In additi on, PBCs are prescribed in all three directions of the simulation box and the pressure on the other two non-loading direction s is relaxed to zero to ensure the uniaxial tension of the columns . Based on the obtained stress -strain curves, the mechanical properties of Li xSi alloys, such 6 as Young’ s modulus and yield stress, are estimated. For the study of the fracture behaviors of Li xSi alloys, simulat ion cells with the size of 500 × 350 × 15 Å in x, y, and z directions are cut from the fully relaxed Li xSi blocks. PBCs are applied in both y and z directions. An edge crack with 0a 90= Å initial length and 4 Å width in x-y plane is created in the middle of the cell along x direction, as illustrated in Figure 1 (a). After fully relaxed, the cell is stretched in y direction by changing the length of simulation box with a constant strain rate 75 10 1 s× . All the MD simulations are conducted at 10 K in LAMMPS with 1 fs time step . The strain at the critical moment that the crack starts to propagate is recorded for the next step chemo -mechanical calculations. Figure 1. A typical Li xSi alloy cell for MD simulation with a pre -existed edge crack . (a) Initially, the size of the simulation cell is about 500, 350, and 15 Å in x, y, and z dimensions , respectively. Initially, the crack is about 9 0 Å in length and 4 Å in width. (b) Configuration of the simulation cell at the critical moment that the crack starts to propagate. 2.2. Chemo -mechanical simulation 7 The previously developed finite strain chemo- mechanical model is adopted in this study.18, 32, 36 In this model, the concurrent chemical reaction and diffusion processes are simulated in a unified manner by a non- linear diffusion model, despite the difference in interfacial reaction and bulk diffusion. Therefore, the classical diffusion equation is used to describe Li tran sport in the Si: FF,cDtµ∂= ∇ ⋅ = − ∇∂JJ , (1) where, c is the Li concentration, D is the diffusion constant, FJ is the Li flux that related to the chemical potential of Li, µ, in Si defined as follows:37 0 maxlogg m kkcRT Vccµ µ ησ = +−−. ( 2) The 0 maxloggcRTccµ+− and m kkVησ are the stress- independent and stress -dependent parts of the chemical potential, respectively . 0µ is a constant representing the chemical potential at a standard state, gR is the gas constant, maxc is the concentration corresponding to the saturation state, mV is the molar volume, η is the coefficient of compositional expansion, and kkσ is the trace of stress tensor . For the mechanics part of the model, the finite -strain elasto -plastic framework is used, in which the total stretch rate tensor d is decomposed into three additive parts: epc=++dd d d , (3) 8 where, ed , pd , and cd are the elastic, plastic, and chemical stretch rates, respectively. With the same -size edge -cracked Li xSi samples as used in above MD simulation, t he following equilibrium equation : 0 ∇⋅ =σ ( 4) is solved with corresponding boundary conditions. All the material properties obtained from the above MD simulations, such as Young’ s modulus, Poisson’ s ratio, and yield stress, are adopted as inputs for the modeling. Using the implicit coupled temperature - displacement procedure in ABAQUS/Standard with user subroutines, the corresponding Li concentration and stress- strain fields are updated incrementally. Meanwhile, contour integrations are conducted for the classic Rices J -integral :38 21 1()i eiuJ Wn t dxΓΓ∂= −Γ∂ ∫ (5) where, as shown in Figure 2, it is the traction vector on the contour ∂Ω that around the crack tip, in is the unit outward normal on ∂Ω, iu is the displacement vector, dΓ is the length element along the contour, and eW is the elastic energy density. Figure 2. A crack tip contour used in the J -integral calculation. The planar crack is oriented 9 along the 1x direction, the region Ω around the crack tip is surrounded by the closed contour 12 −+ ∂ Ω=Γ +Γ +Γ +Γ , and n the unit outward normal on ∂Ω. 3. Results and discussion Upon successfully setting the simulation cell with corresponding boundary and loading conditions, uniaxial tensile simulations of Li xSi samples, such as crystalline <111> - Si, a-Si, a-Li0.25Si, a-Li0.50Si, a-Li0.75Si, and a- Li1.00Si, are conducted at 10 K under the strain rate of 75 10 1 s× . For the amorphous Li xSi samples, eight different configurations are simulated for each compositional case to average the obtained results. As shown in Figure 3, t he obtained stress -strain curves of a-LixSi exhibit the typical characteristics of metallic materials.35 Material properties of pure Si and Li xSi alloys, such as Young’ s modulus and yield stress, are estimated from the stress -strain curves . As listed in the Table 1, the Young’ s modulus and yield stress of Li xSi alloys are highly dependent on the concentration of Li, while the Poisson’ s ratio shows much less sensitiv ity to the variation of Li concentration. 10 Figure 3. Uniaxial tensile stress -strain curves of crystalline <111> -Si, a-Si, and a- LixSi at 10 K under the strain rate of 75 10 1 s× . MD simulation s of pre-cracked pure Si and Li xSi alloy samples are conducted at 10 K to investigate their fracture behavior. During the simulation s, interaction between the two pre -exist crack surface is shield ed. Under the constant strain rate loading condition, the crack tip open displacement increases gradually. With the increment of Li concentration in the Li xSi alloy samples, much blunter crack tip can be observed. At the critical moment that the crack starts to propagate, as sho wn in Figure 1(b) as an example , the critical strain is recorded and listed in Table 2. It should be noted that, as all the MD simulations are conducted with an extremely high strain rate (75 10 1 s× ), the simulation time is too short for Li migration. Therefore, no Li redistribution around crack tip is observed, even though previous studie s have already demonstrated that Li atoms are highly mobile and can easily 11 be driven by stress gradient.39-45 Table 1. The Young’ s modulus (E), Poisson’ s ratio (ν), and yield stress (Yσ) of pure Si and LixSi alloys obtained from the uniaxial tensile simulations. E (GPa) ν Yσ (GPa) <111> -Si 189.4 0.17 11.7 a-Si 151.7 0.27 6.7 a-Li0.25Si 105.6 0.33 5.4 a-Li0.50Si 82.8 0.35 4.1 a-Li0.75Si 70.8 0.37 2.9 a-Li1.00Si 65.9 0.34 2.7 Using the material properties obtained in the above MD simulations as inputs, continuum level chemo- mechanical calculations are conducted for the fracture toughness of pure Si and Li xSi alloys . The same- size pre- cracked 2D plane strain samples are used, and the corresponding critical strain is applied as the boundary condition. As no Li redistribution is observed in the above MD simulations, Li distribution in the samples is assumed to be uniform throughout the chemo- mechanical simulation. Therefore, for simplicity of the model, no Li transport needs to be considered in the chemo- mechanical modeling. As a demonstration of the method presented in this study, we first calculated the fracture toughness of <111> -Si as 2 C3.60 J m J= based on the contour in tegration of the classic Rices J -integral , leading to the critical stress intensity factor as 12 IC 111 C 20.84 MPa m1EKJν<>= = − , which is in a good agreement with the reported experimental measurements, namely =0.83~0.95MPa mICK ⋅ .46-47 This calculation indicates that the approach presented here is appropriate for the derivation of fracture toughness . With confidence in the meth od, we next calculated the fracture toughness (CJ) of a- Si and Li xSi alloys, as listed in Table 2. The obtained CJ for the low Li concentration cases investigated in this study matches the experimental results published by Pharr et al reasonably well.22 In add ition, with the increment of Li concentration in the Li xSi alloys, the CJ decreases first to a minimum value around Li 0.25Si and then increases. This trend also agrees with the previous experimental measurement conducted by Wang et al .23 It should be noted that, as high strain rate is used in the above study, Li transport during the chemo -mechanical calculation of fracture toughness can be neglected. However, in reality, the loading time scale in quasi -static fracture toughness measurements is comparable to the Li diffusion time scale.43 Therefore , it is important to consider the Li redistribution caused by the stress gradient around the crack tip.45 Table 2. The critical strain and fracture toughness (CJ) of pure Si and Li xSi alloys obtained from MD simulation and chemo -mechanical modeling, respectively. Critical strain 2 CJm)J ( <111> -Si 0.0245 3.60 a-Si 0.0341 5.80 13 a-Li0.25Si 0.0363 4.73 a-Li0.50Si 0.0441 5.39 a-Li0.75Si 0.0449 5.29 a-Li1.00Si 0.0450 4.80 As shown in Figure 4, when a cracked specimen with width 0500L= Å, initial crack length 0a 90= Å, and initial Li concentration 0c is loaded with a far field stress 0σ which is much less than the yield stress Yσ along the vertical direction, the normalized effective stress eff Y/σσ is shown in Figure 4(a). W ithout considering Li transport, the initial uniform Li distribution is unaltered in the sample , therefore the normalized changed amount of Li concentration due to the far field loading is 0/0cc∆= (Figure 4(b)), and the normalized classic Rices J -integral (Figure 4(c) ) obtained from three different integration contours (as shown in Fig. S1 in the Supporting Information) are overlapping with each other, indicating that the J -integral is path -independent. However , if stress gradient driving Li transportation effect is considered and the Li diffusion constant is set to be 2 14mDs−=1×10 ,43 with the same far field stress loading condition, the crack tip s tress gradient can drive Li to aggregate around the crack tip ( Figure 4(e) ), causing more severe stress concentration at the crack tip ( Figure 4(d) ). In addition, Li concentration gradient around the crack tip can lead to the path- dependent behavior of the classic Rices J -integral, as shown in Figure 4(f). 14 Figure 4. Compared with the normalized effective stress (a), normalized changed amount of Li concentration (b), and normalized classic Rice’ s J-integral (c) obtained from the case without Li transport, under the same far field stress loading conditions, the crack tip st ress gradient can drive Li to aggregate around the crack tip (e), causing more severe stress concentration at the crack tip (d), and leading to the path- dependent behavior of the classic Rices J -integral (f) . In order to overcome the limitation of the classic Rices J -integral when there exists 15 a Li concentration gradient, the newly developed chemo- mechanical J- integral need s to be used for the estimation of fracture toughness .37, 48 The chemo- mechanical J -integral is defined by: 210 1( ) ()cm i ieuJ n t d W c TSxψ ψ µµΓΓ∂= − Γ , = − − −∂ ∫ , (6) where, ψ is the grand potential of the system, T is the absolute temperature, and S is the entropy of mixing per volume of mixture giving by: max max max max maxlog( ) (1 )log(1 )gcc c cS Rccc c c= − +− −   . (7) Using this chemo- mechanical J- integral, the crack driving forces for the case with Li transport are calculated, as shown in Figure 5. It should be noted that, as the contribution of Li is considered in the grand potential term in Eq. (6) , the calculated chemo -mechanical J-integral is normalized by max 0 ( T )aYgRc σ+⋅ rather than 0aYσ⋅ that the classic Rices J - integral is normalized by. With the Li flow caused by crack tip stress gradient, the chemo - mechanical J- integrals obtained from different integration contours differ slightly from each other. Once the Li redistribution reaches a steady state, the chemo -mechanical J- integrals will finally converge to a single value , indicating the path- independent property of the chemo -mechanical J- integral. 16 Figure 5. Normalized chemo -mechanical J- integrals obtained from different integration contours . With the Li flow caused by crack tip stress gr adient, the chemo- mechanical J- integrals obtained from different integration contours differ slightly from each other. Once the Li redistribution reaches a steady state, the chemo -mechanical J- integrals will finally converge to a single value, indicating t he path- independent property of the chemo - mechanical J- integral. 4. Conclusions In conclusion, we have presented a molecular dynamics (MD) simulation informed continuum chemo- mechanical model with conservation integrals for the investigation of fracture toug hness of Li xSi alloys. In this model, MD simulations are used to calculate the material properties and the critical fracture strain of Li xSi alloys. Using the obtained MD 17 simulation results as inputs, continuum level chemo- mechanical simulations are conducted for the stress field and thus the fracture toughness of Li xSi samples with small scale yield ing. Using the presented model, we investigated the fracture toughness (CJ) of Li xSi alloys with low Li concentration. O ur modeling results, which are reasonably consistent with the limited publi shed results derived from different experimental measurements, revealed that the fracture toughness of Li xSi alloys with low Li concentration do es not vary significantly with Li concen tration. In addition, we demonstrated that, when the loading time scale of fracture test is comparable to the Li diffusion time scale, Li redistribution caused by the stress gradient around crack tip need s to be considered, as stress gradient can drive Li to migrate from the far field to crack tip , leading to the aggregation of Li there, and this process is both therm odynamically and kinetically favorable. Under such circumstance, the classic Rices J -integral becomes path- dependent. Therefore, an appropria te chemo -mechanical path -independent J -integral need s to be used for the estimation of fracture toughness . The obtained fracture toughness of the Li xSi alloys here provides guidance for the rational design of Si based electrodes, and the approach presented in this study also sheds light for the evaluation of the fracture toughness of other energy materials at different charging/discharging level s. Acknowledgements : H.Y . acknowledges the National 1000 Talents Program of China tenable at the Huazhong University of Science and Technology (HUST), China. H.Y . and J.M.Q. acknowledge 18 National S cience F oundation support under grants CMMI -1624313. ORCID: Hui Yang: 0000- 0002- 2628- 4676 Conflict of Interest: The authors declare no competing financial interest.