Vulnerable window of yield strength for swelling-driven fracture of phase-transforming battery materials

Despite numerous experimental and theoretical investigations of the mechanical behavior of high-capacity Si and Ge Li-ion battery anodes, our basic understanding of swelling-driven fracture in these materials remains limited. Existing theoretical studies have provided insights into elasto-plastic deformations caused by large volume change phase transformations, but have not modeled fracture explicitly beyond Griffith’s criterion. Here, we use a multi-physics phase-field approach to model self-consistently anisotropic phase transformation, elasto-plastic deformation, and crack initiation and propagation during lithiation of Si nanopillars. Our computational results reveal that fracture occurs within a “vulnerable window” inside the two-dimensional parameter space of yield strength and fracture energy and highlight the importance of taking into account the surface localization of plastic deformation to accurately predict the magnitude of tensile stresses at the onset of fracture. They further demonstrate how the increased robustness of hollow nanopillars can be understood as a direct effect of anode geometry on the size of this vulnerable window. Those insights provide an improved theoretical basis for designing next-generation mechanically stable phase-transforming battery materials undergoing large volume changes.


INTRODUCTION
Increasing demand for portable energy storage has motivated a large research activity focused on high-capacity Li-ion battery anodes. Current carbon-based anodes have limited theoretical capacity (372 mA h g −1 for Li 6 C (ref. 1 )). Silicon and germanium have an order of magnitude larger theoretical capacity (3579 mA h g −1 for Li 15 Si 4 , 4200 mA h g −1 for Li 22 Si 5 (ref. 1 ), 1384 mA h g −1 for Li 15 Si 4 (ref. 2 )) but are prone to fracture due to the high, approximately 300%, volume expansion during lithiation 3,4 , which limits their use. Different designs have been explored to overcome this limitation including silicon nanopillars 5,6 , thin films [7][8][9] , open nano-porous crystalline Si structures with ultra-high interfacial area produced by dealloying of Si-based alloys 10,11 , combinations of these 12,13 , as well as composite designs that embed silicon particles inside a more mechanically stable matrix [14][15][16][17][18] .
Basic studies of the lithiation process have shown that crystalline silicon (c-Si) transforms to an amorphous lithiated alloy (a-Li x Si) 2,19,20 . The kinetics of this large volume change phase transformation is understood to be both interface-reaction limited 20 and highly anisotropic, reflecting the two key observations that the velocity of the c-Si/a-Li x Si interface remains approximately constant during lithiation, and that this velocity depends strongly on crystallographic orientation 19 . Numerical studies of elasto-plastic deformations of c-Si particles (nanowires, nano-/micro-pillars, etc.) have demonstrated that the resulting anisotropic swelling can produce both large shape changes of the particles and, as a non-trivial effect of compressive yielding, tensile stresses on their outer surface that can potentially drive fracture 3,[21][22][23] . Those insights have already proven useful to test new anode designs to mitigate fracture 23 . However, our ability to predict when and how fracture occurs in different large volume change materials (e.g., Si versus Ge) and different anode geometries (e.g., solid versus hollow nanopillars 13 ) is still limited. To date, the onset of fracture has been predicted using analytical solutions for stresses obtained in idealized geometries, assuming purely plastic deformation and isotropic swelling, and by applying a Griffith criterion to predict fracture onset for a flaw size comparable to the particle dimension 24,25 . However, crack initiation and propagation in the setting of large elasto-plastic deformation and anisotropic swelling in phase-transforming materials remain largely unexplored.
Here, we use a multi-physics phase-field approach to simulate both anisotropic swelling and fracture of solid and hollow c-Si nanopillars within a unified theoretical framework and derive from our simulations an understanding of when and how fracture occurs as a function of key materials parameters, including yield strength and fracture energy, and geometric parameters such as nanopillar radius and slenderness. Phase transformation is modeled using a phase-field ψ that distinguishes the c-Si and a-Li x Si phases and is evolved dynamically to describe the interfacereaction-limited anisotropic motion of the c-Si/a-Li x Si interface. Fracture, in turn, is modeled using the well-established variational approach that couples elasticity to a phase-field ϕ, which distinguishes pristine and broken regions of the material 26,27 . To realistically model large volume changes, this variational approach is implemented using a large deformation formulation of elastoplasticity combining neo-Hookean nonlinear elasticity and J 2 plasticity to quasi-statically evolve ϕ together with the material displacement field and the plastic deformation gradient tensor. The phase-field approach offers several advantages in the present context. It provides a self-consistent formulation to model simultaneously anisotropic swelling, large elasto-plastic deformation, and fracture. Furthermore, it can describe the evolution of phase boundaries and cracks of arbitrarily complex shapes, as demonstrated in applications to other phase transformations 28 and fracture problems such as thermal shock fracture 29 , mixed mode fracture 30 , ductile fracture 31,32 , and the simpler chemomechanical fracture of single-phase battery cathode particles [33][34][35][36] , also driven by volume expansion due to Li intercalation but only involving small elastic stresses and no phase change. In addition, in contrast to Griffith theory, the phase-field approach is able to describe crack initiation without pre-existing flaws. This property stems from the fact that ϕ varies smoothly in space on a length scale ξ, thereby enabling crack formation on the scale of the "process zone" where elastic energy is transformed into new fracture surfaces. Hence, directly relevant to the present study, the phase-field approach can quantitatively describe crack initiation from surface imperfections such as U-or V-shape notches by treating ξ treated as a material-dependent parameter 37 . V-shape notches, in particular, bear close similarity to surface shape deformations of lithiated Si particles undergoing elasto-plastic deformation during anisotropic swelling 3,21-23 .
To keep computations tractable, we perform 2D plane-strain simulations (∂ z ≡ 0) on a cross-section of an unconstrained nanopillar (τ zz = 0) lithiated from its surface (i.e., outer boundary for a solid nanopillar and both outer and inner boundaries for hollow nanopillars). Furthermore, to dissect the contributions of multiple interacting physical effects (including compressive and tensile yielding, anisotropic swelling, localization of plastic deformation, and crack initiation and propagation), we carry out three different types of computations of increasing complexity. In a first step, we model stress evolution without fracture by assuming that swelling is isotropic and that the stress fields and plastic hardening parameter α only vary radially and are independent of the azimuthal angle θ as depicted schematically in Fig. 1a. This axisymmetric approximation reduces the 2D problem to a 1D radial problem. Stress evolution in a similar idealized geometry has been previously studied analytically by taking into account only plastic deformation 24,25 . By taking into account here both elastic and plastic deformations, we demonstrate that tensile stresses generated on the particle surface by volume expansion reach a maximum value as a function of yield strength σ y . Even though the critical σ y value corresponding to this maximum is outside the experimentally 8,38,39 or theoretically estimated 40 range σ y0 .5-2 GPa for Si, the existence of this critical yield strength provides a valuable theoretical framework to understand fracture behavior inside this lower estimated range of 0.5-2 GPa. For this reason, we investigate stress evolution over a wide range of σ y that encompasses the entire vulnerable window for fracture. In a second step, we carry out a similar computation, still without fracture, but for the full 2D problem without the axisymmetric approximation in which the stress fields and α can vary both radially and azimuthally inside the pillar cross section. This enables us to asses how anisotropic swelling modifies tensile stresses on the pillar surface. We find that tensile stresses become amplified by localization of plastic deformation, but still exhibit a maximum as a function of increasing σ y . Those 1D and 2D computations demonstrate the existence of a vulnerable window of yield strength inside which pillars are prone to fracture. More crucially, for the relevant experimentally reported yield strengths (σ y = 0.5-2 GPa) of Si, the difference between the material and the most critical yield strength controls the magnification of generated tensile stresses. In the last step, we validate the existence of this window by repeating our 2D computations with fracture, showing that pillars fracture only over an intermediate range of yield strength. We compare the results of full 2D simulations with estimates based on our numerically calculated stresses in the previous steps using the Griffith theory framework. We then use experimental estimates of safe pillar radius (i.e., largest pillar radius without fracture) to quantitatively validate our findings by calculating our estimate of yield strength. Plots of maximum hoop stress reached during complete lithiation on the nanopillar outer (r = R, solid lines) and inner (r = R − t, dashed lines) boundaries versus yield strength σ y for different t/R ratios. Plots predict the existence of a "vulnerable window" of σ y inside which the maximum hoop stress can exceed the threshold for fracture. c Radial profiles at different times of the phase transformation ψ field (right vertical axis, gray lines with the c-Si/a-Li x Si interface located at ψ = 0.5), the Kirchhoff hoop stress ffiffi ffi 3 p τ θθ =2 (left vertical axis, red lines with earlier stages shown using lighter red) and von Mises stress τ eq (left vertical axis, blue lines) for σ y = 3 GPa. Plots show compressive yielding followed by subsequent reversal of the sign of the hoop stress and yielding under tension. d Same as c but for higher yield strength σ y = 7 GPa where tensile yielding does not occur.

Model
We model the swelling-driven deformation of the material using the finite J 2 elasto-plasticity framework [41][42][43] and account for the fracture of the material by coupling it to a phase-field fracture model 26,27 . Furthermore, we model the anisotropic motion of the c-Si/a-Li x Si interface during lithiation using a non-conserved phase field ψ where ψ = 0 in the crystalline phase and ψ = 1 in the amorphous phase. The material properties are approximated using a linear rule of mixture between the crystalline and amorphous phases. We define the deformation gradient tensor as F ij = ∂x i /∂X j = 1 + ∂u i /∂X j where X i are the undeformed coordinates and x i = X i + u i are the deformed coordinates of material points and u is the displacement field. We use a multiplicative decomposition of the deformation gradient tensor such that where J ψ = (1 + βψ) 2 is the phase dependent volumetric expansion due to phase change with linear Vegard expansion coefficient β, F p is the plastic deformation, and F e is the elastic deformation. We use the framework of the phase-field method for fracture 26,27 by introducing the fracture phase field ϕ along with the process zone size ξ and combine it with the standard phasefield approach for phase transformations to prescribe the ψ evolution 44 . Accordingly, we write the total free energy functional as the sum FðF e ; ψ; ϕÞ ¼ F el ðF e ; ψ; ϕÞ þ F ϕ ðϕÞ þ F ψ ðψÞ; of contributions corresponding to the elastic energy stored in the material including the amorphous and crystalline phases (F el ), the energy associated with the creation of new fracture surfaces (F ϕ ), and the free-energy associated with the c-Si to a-Li x Si phase transformation (F ψ ), defined as F el ðF e ; ψ; ϕÞ ¼ Z Ω0 ðϕ 2 þ η ξ ÞW þ ðF e ; ψÞ þ W À ðF e ; ψÞ À Á dx; (3) In the sharp-interface limit ξ/R ≪ 1 where the process zone size ξ is much smaller than the nanopillar radius R, F ϕ ðϕÞ reduces to the standard fracture energy G c of convenctional fracture mechanics, defined as the energy required to create new fracture surfaces per unit extension of the crack, where we have defined For the phase transformation contribution F ψ ðψÞ, f dw (ψ) = 4ψ 2 (1 − ψ) 2 is a standard double-well potential with equal minima at ψ = 0 and ψ = 1, corresponding to a c-Si/a-Si interface with thickness w ¼ ffiffiffiffiffiffiffiffiffiffi σ ψ =h p and excess free-energy γ = 2hw/3. In addition, f tilt (ψ) provides the driving force for the solid-state chemical reaction xLi + + xe − + Si → Li x Si by tilting the double well potential (hf dw (ψ) − Γf tilt (ψ)/3) so as to lower the bulk freeenergy density of a-Li x Si (ψ = 1) relative to c-Si (ψ = 1) by an amount proportional to is the chemical potential of Li (Si) in a-Li x Si at the interface,μ Li and μ c Si are the chemical potential of Li in the counter-electrode and Si in c-Si, respectively, and ΔV and F are the electric potential relative to the counter-electrode and Faraday's constant, respectively. The choice f tilt (ψ) = 15(32ψ 3 /3 − 16ψ 4 + 32ψ 5 /5)/16, with limits f tilt (0) = 0 and f tilt (1) = 1 and vanishing first and second derivatives at ψ = 0 and ψ = 1, preserves the freeenergy minima at ψ = 0 and ψ = 1 for arbitrary Γ. Finally, to construct FðF e ; ψ; ϕÞ, we split the elastic energy into positive (W + (F e , ψ)) and negative (W − (F e , ψ)) volumetric parts to distinguish material regions under tension and compression, respectively, and assume a neo-Hookean hyper-elastic material μðψÞ 2b e kk À 2 otherwise 8 > < > : where J e = det(F e ),b e ij ¼ F e ik F e jk =J e is the isochoric left Cauchy-Green deformation tensor, and the shear modulus μ(ψ) = ψμ a + (1 − ψ) μ c + η ξ μ a and bulk modulus κ(ψ) = ψκ a + (1 − ψ)κ c + η ξ μ a are interpolated between their values in the a-Li x Si (μ a , κ a ) and c-Si (μ c , κ c ) phases with the addition of a residual stiffness η ξ that is introduced to ensure non-singularity of the balance equations and can be chosen arbitrarily small η ξ = 2.5 × 10 −4 without affecting the results. This split allows for fracture to take place under shear and under tensile loads but prohibits it under pure volumetric compression. Moreover, the formulation forbids interpenetration of fracture surfaces under compression while allowing them to slide. Following classic J 2 plasticity we assume that the von Mises equivalent stress τ eq ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ij Àb e kk δ ij =2Þ is the deviatoric part of the Kirchhoff stress tensor τ (τ = Jσ where σ is the Cauchy stress tensor), σ y is the yield strength, K is the isotropic hardening constant, and α is the hardening parameter (at the infinitesimal strain limit the isotropic hardening parameter reduces to the equivalent plastic strain). The equations governing elasto-plastic deformation coupled to fracture driven by the large volume change phase transformation are obtained variationally from (2) and consist of two parts. Firstly, for fixed ψ, the displacement field and fracture phase-field are evolved quasi-statically by minimizing (2) with respect to u i and ϕ δF δu i ¼ 0 s:t: respectively, where δF=δ is the Fréchet derivative of the free energy F with respect to field •. Standard J 2 plasticity is incorporated by the constraint in (8) that plastic deformation occurs when the equivalent von Mises stress reaches the yield surface. The coupling of elasto-plasticity and fracture is selfconsistently captured by the fact that (8) and (9) are derived from the same free-energy functional, which couples the evolution of the displacement field u i and fracture phase-field ϕ. Secondly, the ψ field is evolved in time to model the reaction-controlled anisotropic motion of the c-Si/a-Si interface as wheren ¼ ∇ 0 ψ=j∇ 0 ψj is the direction normal to the interface and MðnÞ is the anisotropic interface mobility modeled using the same functional form as An et al. (Equation (5) in ref. 23 ). In the sharp interface limit w/R ≪ 1 45 , (10) reduces to the simple law of curvature-independent interface motion V n ¼ MðnÞΓ with the interface velocity V n controlled by the thermodynamic driving force for the interface reaction Γ assumed to be spatially uniform.
As discussed previously 23 , this assumption holds as long as Li diffusion is much faster than the reaction kinetics so that the chemical potentials of Li and Si remain spatially uniform inside the a-Si phase, which is expected to be the case for the small nanopillar sizes investigated here. In this limit, interface motion is independent of the stress distribution inside the material and purely geometrically determined by the form of the anisotropic mobility MðnÞ. Hence, elastoplastic deformation and fracture are Vulnerable window of yield strength To understand the basic mechanism of stress generation in these components, as the first step, we performed an exhaustive series of axisymmetric simulations without fracture (i.e., simulations in which stress fields and the plastic flow hardening parameter are assumed to only vary radially). Results of the axisymmetric computations are shown in Fig. 1. Figure 1c, d shows the evolution of the hoop stresses and equivalent von Mises stress during lithiation of a solid nanopillar where we can identify three regimes in these figures. As the phase transformation boundary (ψ = 0.5) invades inside the particle from the outer boundary (r = R), it creates compressive stresses due to the large volumetric expansion of the a-Li x Si phase. The resulting compressive stresses generate plastic flow that caps the von Mises stress at σ y . As the crystalline core shrinks further, the compressive stresses on the outer boundary subside and change sign due to the initial compressive yield. Consequently, the hoop stress on the outer boundary changes sign and becomes tensile (Fig. 1c, d), thereby confirming the knock-on effect of compressive yielding on the creation of tensile stresses that has been hypothesized to cause cracking 3,21-23 . Importantly, for σ y smaller than approximately 4 GPa for the present parameters, the tensile hoop stress reaches the yield strength before the c-Si core has vanished, which results in secondary plastic yielding under tension (Fig. 1c). In this range (σ y ≤ 4 GPa), the maximum hoop stress reached during complete lithiation, maxðτ θθ Þ, increases linearly with σ y as shown in Fig. 1b for t/R = 1 corresponding to a solid pillar; since the outer boundary is traction free (τ rr ≡ 0) and the von Mises stress is capped by σ y , maxðτ θθ Þ ¼ 2σ y = ffiffi ffi 3 p on that boundary for planestrain. In contrast, for larger σ y , compressive yielding requires a larger lithiated fraction, which reduces the amount of volumetric expansion available to create tensile stresses during shrinkage of the remaining c-Si core. Therefore, maxðτ θθ Þ remains below the yield strength and decreases with increasing σ y as shown in Fig.  1b. We should also highlight that although all simulations presented in this article were performed using β = 0.7 that corresponds to~280% volume change at full lithiation, our results show that there exists a universal relationship between the dimensionless maximum hoop stress maxðτ θθ Þ=μ a β and the dimensionelss yield strength σ y /μ a β (see Supplementary Fig. 1), such that our findings can be extended to other materials whose phase-transformation result in smaller volume changes. This universality can be readily understood noticing that at smaller expansion coefficients, smaller stresses are generated; therefore, the knock-on effect of the compressive yielding only takes place at smaller yield strength. Crucially, these results show that the vulnerable window of yield strength is shifted to smaller values of yield strength for smaller expansion coefficient. The resulting hat shape of the maxðτ θθ Þ versus σ y plot in Fig. 1b suggests the existence of a vulnerable window for fracture corresponding to the range of σ y where the maximum hoop stress becomes large enough to initiate fracture. Specifically, existence of a maximum generated hoop stress at a critical yield strength increases the available energetic driving force for fracture at lower yield strength as confirmed below by our full 2D simulations (Figs 3-5). Figure 1b also shows plots of maxðτ θθ Þ versus σ y on the inner (r = R − t) and outer (r = R) boundaries of hollow nanopillars. The maximum hoop stress on the outer boundary still exhibits a maximum. However, since the compliance of the annulus is inversely related to its slenderness t/R, the maximum stresses on both boundaries decrease with increasing slenderness.
Next, we performed 2D simulations of anisotropic swelling without fracture for pillars oriented in two crystallographic directions [001] and [112]. In these simulations, the crystalline core is no longer circular and the anisotropic mobility of the amorphization front creates a crystalline silicon core with sharp corners. During lithiation, those crystalline corners concentrate stresses and localize plastic flow in their vicinity. When the stresses change sign and become tensile on the pillar outer boundary, shear localization produces V-shaped notches at orientations corresponding to these corners for lower yield strength, which allows tensile yielding to occur on the periphery subsequent to compressive yielding, but not larger yield strength where tensile yielding does not occur. This difference can be seen in the pillar morphologies in Fig. 3 for σ y = 1 GPa and σ y = 10 GPa (σ y = 7 GPa for [112] oriented pillar) that did not fracture. Those notches further concentrate stresses, thereby augmenting the magnitude of hoop stresses several fold at those orientations. This magnification is shown in Fig. 2 where we compare maxðτ θθ Þ, defined as before as the maximum hoop stress reached in time during complete lithiation, from 1D axisymmetric computations of isotropic swelling and the present 2D computations of anisotropic swelling. For the latter case, we report maxðτ θθ Þ both on the outer surface (blue and yellow diamonds) and at a position inside the particle close to the outer surface (red circles) where maxðτ θθ Þ reaches its maximum value along a vertical axis that contains the corners of the c-Si core. The maximum hoop stress is seen to be magnified both by localization of plastic deformation during  (Fig. 3) magnifies the magnitude of stresses, thereby enlarging the size of the vulnerable window for fracture. The largest stresses are created along the crystalline corners at a small distance from the surface (red circles).
A. Mesgarnejad and A. Karma compressive yielding, which occurs for all σ y reported, and Vshape notches that form for lower σ y due to tensile yielding.
The 2D results in Fig. 2 confirm the existence of a critical yield strength that generates maximal tensile stresses but were obtained from simulations without fracture. To investigate the effect of stress augmentation on fracture, we repeated a series of simulations with fracture for different σ y and the [001] and [112] crystallographic orientations and for fixed process zone size to radius ratio ξ/R = 0.02. The time evolutions of particle morphologies in Fig. 3 show that nanopillar fracture during lithiation over an intermediate range of σ y inside the vulnerable window centered at the critical yield strength. For σ y = 1 GPa, Vshaped notches are created along the crystalline corner directions but the magnitude of surface tensile stresses are insufficient for crack initiation. For σ y = 3 GPa, cracks initiate due to stress concentration at V-shaped notches and propagate unstably towards the crystalline core. For σ y = 5 GPa, tensile yielding and hence V-shaped notches are absent but tensile stresses grow sufficiently large (at a later stage of charging compared to σ y = 3 GPa) to create two pairs of cracks that propagate unstably after initiation. The horizontal pair subsequently arrests and the vertical pair propagates further due to spontaneous symmetry breaking. Finally, for σ y = 10 GPa (σ y = 7 GPa for [112] oriented pillar), reduced compressive plastic flow in the system suppresses the subsequent increase of tensile stresses, thereby preventing crack initiation. We used a 1% hardening coefficient (K/μ a = 0.01) in all the simulations presented in the main text, consistent with experimental measurements of stress evolution in thin films showing negligible hardening 8 . Incorporation of a small amount of hardening avoids numerical convergence issues associated with modeling perfect plasticity (K = 0). We carried out additional simulations to test the effect of varying the hardening coefficient.
The results show that increasing hardening from 1% to 5% (K/μ a from 0.01 to 0.05) yields a slight reduction of localization of plastic deformation ( Supplementary Fig. 2) and a small increase in maximum tensile stress ( Supplementary Fig. 3), thereby leaving the vulnerable window of yield strength largely unchaned. Therefore our results are weakly dependent on the choice of hardening coefficient up to a 5% value that represents a reasonable upper bound based on measurements of a-Li x Si plastic behavior.
Size effects Experimental observations of Si anodic components show a clear size dependency where, for example, nanospheres with radii smaller than~75 nm 5 and nanopillars with radii smaller than 120 nm 6 do not break during lithiation. Such size dependencies, prevalent in brittle and quasi-brittle materials [46][47][48][49] , are typically characterized by a power law relationship between the stress to fracture and component size, τ c $ 1= ffiffiffi R p , for R much larger than the process zone size. Within the theoretical framework of Linear Elastic Fracture Mechanics (LEFM), which treats crack surfaces as sharp boundaries, this power law is readily obtained from the expression for the energy release rate at the tip of a crack of length a under a spatially homogeneous critical stress τ c , which can be written as G ¼ Caσ 2 ð1 À νÞ=μ for plane-strain where C is a dimensionless constant that generally depends on particle geometry and load configuration. Equating the energy release rate with the fracture energy (G = G c ), we obtain the expression and thus the scaling τ c $ 1= ffiffiffi R p by further assuming that the maximum flaw size a increases proportionally to the particle size (a~R and a ≪ R). In the phase-field model used in this article, which describes the state of the material with a spatially varying scalar field ϕ, crack nucleation is an inherent property of the oriented R = 85 nm nanopillars for different σ y and for fixed process zone size to radius ratio ξ/R = 0.02. The color map depicts the hardening parameter α and the thick black line shows the amorphous-crystalline boundary (i.e., ψ = 0.5 contour line). The red dashed boxes show snapshots just before crack initiation and after unstable penetration towards the crystalline core. The results confirm the prediction of a window of σ y for fracture and distinguishes modes of fracture with (second row) and without (third row) localization of plastic deformation creating a surface V-shaped notch prior to fracture (see also Supplementary Videos 1a, b-6a, b).
A. Mesgarnejad and A. Karma model and occurs via an instability that causes ϕ to develop a local dip (ϕ → 0) when the local stress exceeds a critical value τ c $ ffiffiffiffiffiffiffiffiffiffiffiffi G c μ=ξ p 50 (up to a numerical prefactor that also depends on particle geometry and load configuration). Consequently, by comparing the above scaling expression for τ c to Eq. (11), we can physically interpret ξ as playing an analogous role to the dominant flaw size in the LEFM framework. Furthermore, by using the result of a stability analysis of a 1D stretched strip in the phase-field model, which yields the prediction τ c ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 3G c μ=4ξð1 À νÞ p 50 , we obtain the constant C ¼ 4=3 by comparison with Eq. (11) with a = ξ. This value is close to the standard LEFM value C ¼ π=2 for a crack of length 2a in a uniform stress. One main qualitative difference, however, is that crack nucleation in the phase-field model occurs through an instability of the pristine state in which ϕ is spatially uniform and hence does not require the introduction of a flaw in the form of a finite length seed crack as in LEFM.
We take advantage of this property to investigate the particle size dependence of fracture onset by performing simulations at fixed process zone size to particle size ratio ξ ¼ ξ=R. From the dimensional analysis, it is natural to define the dimensionless fracture energy as G c /(μ a R) where μ a is the shear modulus of the amorphous a-Li x Si phase. The dimensionless fracture energy can be readily interpreted as the ratio of the Griffith length scale G c /μ a and the nanopillar radius. Fig. 4a reports the results of an extensive series of 2D phase-field fracture simulations for [001] oriented nanopillars for ξ ¼ 0:02, which identifies regions of the two-dimensional parameters space G c /(μ a R) and σ y where fracture does (red circles) or does not (blue crosses) occur. These results confirm the existence of a vulnerable window of σ y where cracks initiate and propagate in nanopillars during lithiation. This window is centered around the critical value of σ y generating maximal tensile stresses and shrinks in size with increasing G c /(μ a R).
We now assess if the onset of fracture in 2D phase-field simulations can be predicted within Griffith theory by using Eq. (11) together with values of the maximum hoop stress during lithiation obtained in 1D or 2D phase-field simulations without fracture (Fig. 2). For this, we assume stress-free surfaces (τ rr (R) = τ rθ (R) = 0) and small elastic strains, which allows to use a quadratic approximation for the stored elastic energy W þ ' τ 2 θθ ð1 À ν 2 a Þ=2E a where E a = 9κ a μ a /(3κ a + μ a ), ν a = (3κ a − 2μ a )/(6κ a + 2μ a ) are the elastic modulus and Poisson ratio of the a-Si. We found that this small strain quadratic form provides an accurate description of the elastic energy in phase-field simulations calculated using nonlinear neo-Hookean elasticity (Eq. (6)). We can then rewrite Eq. (11) as a function of the maximum hoop stress as where C ¼ 4=3. Substituting in the above expression the values of maxðτ θθ Þ obtained from 1D axisymmetric simulations, we obtain the gray dashed line in Fig. 4a that falls significantly below the boundary, comprised between red circles and blue crosses, corresponding to the onset of fracture in 2D simulations with fracture. Consequently, Griffith theory with axisymmetric tensile stresses underestimates the critical value of G c /(μ a R) for fracture at fixed σ y and, hence, overestimates the safe pillar radius. This discrepancy can be attributed to the fact that the 1D axisymmetric simulations lack the stress amplification due to plastic localization and instability. We therefore conclude that localization of plasticity caused by anisotropic volumetric expansion plays a significant role in fracture of Si nanopillars. This conclusion is further supported by the finding that the prediction of Eq. (12) is significantly improved when we use values of maximum hoop stresses obtained from 2D simulations without fracture, which exhibit stress concentration at V-shaped notches. Unlike in 1D axisymmetric simulations, where the tensile hoop stress is always maximum at the pillar surface, hoop stresses in 2D simulations with localization of plastic deformation reach their maximal values inside the particle at a short distance away from the V-shaped notch (Fig. 2). Therefore, we can reasonably use Eq. (12) together with values of maxðτ θθ Þ both inside the particle and at the tip of the V-shaped notch (corresponding to the red circles and blue squares in Fig. 2a, respectively) to obtain lower and upper bounds for the fracture boundary in the plane of G c /(μ a R) and σ y . The results are depicted by the gray shaded region in Fig. 4a that is comprised between the lower and upper bounds computed in this fashion. The fracture boundary between red circles and blue crosses in 2D phase-field simulations falls for the most part inside this gray shaded region (in particular over the range σ y~0 .5-2 GPa of experimental relevance), thereby confirming that stress concentration near V-shaped notches is an important mechanism promoting fracture.
We can now relate our numerical findings to experimental observations of safe nanopillar sizes. Figure 4a shows that the safe Fig. 4 Vulnerable window of fracture. a Vulnerable window of fracture plotted for yield strength σ y vs. dimensionless fracture energy (particle size) G c /(μ a R) for [001] oriented nanopillars for fixed process zone size to radius ratio ξ/R = 0.02. Circles depict cases with, and crosses show cases without fracture. The results confirm the existence of a vulnerable window of fracture energy (size) and yield strength for fracture. Lines show the comparison with the closed-form approximation (Eq. (12)) based on maxðτ θθ Þ using 1D axisymmetric results shown in Fig. 1b  (dashed gray line), and maxðτ θθ Þ ¼ 2σ y = ffiffi ffi 3 p (dashed orange line). The gray shaded region comprises lower and upper bounds of the fracture boundary computed using Eq. (12) together with values of maxðτ θθ Þ both inside the particle and at the tip of the V-shaped notch obtained from full 2D simulations (corresponding to the red circles and blue squares in Fig. 2a, respectively). The green shaded region shows the approximate value of the dimensionless fracture energy calculated for the safe particle size 120 nm 6 and fracture energies in the range 5-7 J m −2 (ref. 39 ). b Phase diagram of fracture in the plane of dimensionless process zone size 0.01 ≤ ξ/R ≤ 0.2 and yield strength σ y obtained from full 2D simulation using G c /(μ a R) = 0.01 (corresponding to R = 120 nm and fracture energy G c = 6 J m −2 ). The results are consistent with the experimentally estimated range σ y = 0.5-2 Gpa for fracture. nanopillar radius (where pillars with radii less than the safe value do not break) decreases with increasing σ y over the estimated range σ y = 0.5-2 GPa for a-Li x Si. The experimental range of safe nanopillar radius is highlighted by the green shaded region in Fig. 4a. This region was computed using the experimentally observed safe nanopillar radius (120 nm) 6 and the estimated range of fracture energy (5-7 J m −2 ) from experimental measurements 39 . 2D phasefield simulations predict that, inside this green shaded region of Fig. 4a, fracture occurs for σ y between 1.5 GPa (blue crosses) and 2 GPa (red circles), which falls in the upper part of the range σ y = 0.5-2 GPa estimated from experimental measurements 8,38,39 . Phase-field modeling prediction also depend generally on flaw size through the ratio ξ/R. To test this dependence, we repeated a series of simulations by varying ξ/R over the range 0.01 to 0.2, which encompasses the value 0.02 used in all simulations presented so far. These simulations were carried out at fixed dimensionless fracture energy G c /(μ a R) = 0.01 calculated using the average reported fracture energy for a-Li x Si G c = 6 J m −2 39 and the observed safe nanopillar radius R ≃ 120 nm 6 . The results reported in Fig. 4b show that nanopillars become more vulnerable to fracture with increasing ξ/R as theoretically expected. For the largest value ξ/R = 0.2 studied here, fracture occurs for σ y between 0.5 GPa (blue cross) and 1 GPa (red circle), which falls in the lower part of the range σ y = 0.5-2 GPa estimated from experimental measurements 8,38,39 . While the precise value of ξ/R is not known, its value is presumably much smaller than unity given that nanopillars do not typically exhibit large visible flaw sizes prior to lithiation 6,20 .
We can further confirm the above estimate of σ y by redoing the calculations for isotropic lithiation of amorphous Si. Experimental observations demonstrate that a-Si lithiates isotropically [51][52][53] . It is therefore reasonable to use our 1D axisymmetric simulations to estimate the magnitude of hoop stresses generated during lithiation of a-Si. We can simplify Eq. (12) further using the expression maxðτ θθ Þ ' 2σ y = ffiffi ffi 3 p valid for small σ y (see Fig. 1b), which yields the prediction G c =μ a R ¼ ð4C=3Þξð1 À νÞðσ y =μ a Þ 2 for isotropic lithiation of a-Li x Si. Zhao et al. 24 obtained a similar expression of the form of G c =μ a R $ ðσ y =μ a Þ 2 previously by an analysis of lithiation that only considers plastic deformation and computes the dimensionless prefactor numerically assuming an initial flaw size comparable to R. In contrast, here, the prefactor is obtained analytically from the aforementioned 1D stability analysis of crack initiation in the phase-field model 50 . We can use this isotropic estimate along with the safe experimentally observed safe nanopillar radius 1 μm 53 to calculate an upper bound for its yield strength 0.4-1.2 GPa using the process zone size 0.02 ≤ ξ/R ≤ 0.2. We can use a similar analysis for Ge. Since lithiation of Ge is observed to be isotropic, using the estimates of its shear modulus μ a ≃ 19 GPa we can estimate σ y in the range 1.5-4.6 GPa. We calculated the above range similarly using the process zone size 0.02 ≤ ξ/R ≤ 0.2 and based on observed safe nanopillar radius of 250 nm 54 .
Geometrical effects Finally, to highlight the non-trivial role of geometry beyond size effects, we investigate the fracture of hollow nanopillars that have been shown experimentally to be more resistant to fracture 13 . Our axisymmetric computations predict that this geometrical protective effect is present for large enough yield strength due to a decrease of the maximum hoop stress reached during complete lithiation of crystalline Si as a function of σ y (i.e., the decrease of the peak value of plots in Fig. 1b with increasing annulus slenderness). However, for low yield strength, the maximum hoop stress remains bounded by σ y even for large slenderness. This implies that the hollow nanopillar design can mitigate fracture only for materials with moderately high yield strength. To test these predictions, we modeled the lithiation and fracture of hollow nanopillars with a constant cross-sectional area equal to solid nanopillars with radii 85-170 nm and different slenderness 0.08 ≤ t/R ≤ 1 for σ y = 3 GPa. The results illustrated in Fig. 5 show that increased slenderness has a protective effect for this σ y value. For example, nanopillars break for t/R = 0.4 but only exhibit minor cracking near the interior surface of t/R ≃ 0.29. Even more slender nanopillars (t/R ≃ 0.08) do not fracture. This protective geometrical effect, however, does not persist at lower yield strength for ξ/R = 0.02. Our axisymmetric computations predict that the maximum hoop stress generated in even the thinnest annulus is equal to that of a solid nanopillar for σ y = 1 GPa (Fig. 1b). Our 2D simulations confirm this prediction by showing that the hollow nanopillar for t/R ≃ 0.29, which is protected for σ y = 3 GPa, fractures for σ y = 1 GPa.

DISCUSSION
In summary, we have used a multi-physics phase-field approach to model simultaneously anisotropic phase transformation, elastoplastic deformation, and crack initiation and propagation during lithiation of Si nanopillars. Our results identify a vulnerable window of yield strength inside which pillars fracture during lithiation and distinguish two different modes of fracture inside that window with and without surface localization of plastic deformation prior to fracture for lower and higher yield strength, respectively. Those two modes follow from the existence of a critical yield strength that generates maximal tensile stresses during lithiation. Combined with experimental measurements of fracture energy 39 and observations of size dependent fracture 5,6,53 , our results yield an estimate of yield strength within a range σ y ≃ 0.5-2 GPa consistent with experimental 8,38,39 and theoretical 40 estimates. This range is smaller than the critical yield strength consistent with the observation of localized plastic deformation during lithiation of Si nanopillars 3,5,6 . Over this range, plastic deformation mitigates fracture by energy dissipation but, at the same time, promotes it by the creation of stress-concentrating V-shaped notches that precede quasi-brittle fracture.
Our results also suggest that the observed increased robustness of hollow Si nanopillars 13 is due to a reduction of the critical yield strength generating maximal tensile stresses with increasing slenderness. This interpretation, however, warrants further investigation since this protective effect is only significant in simulations with large enough σ y values. The present study highlights the importance of computationally informed geometric design that takes into account the subtle interplay between material properties and geometry to generate reliable predictions of mechanical stability of high-capacity battery materials, paving the way for designs that exploit more complex geometries such as open nanoporous structures with ultra-high interfacial area 11 .
We conclude with an outlook for future work. In addition to tensile fracture originating from the outer surface of the amorphous phase, there have also been experimental observations of interface delamination between the crystalline and amorphous phases during the lithiation of Si pillars 55 . This process may be driven by shear stresses localized at the amorphous-crystalline boundary near the free surface of pillars where the amorphous phase swells both radially and along the pillar axis. Such shear stresses, which are absent in our 2D plane strain simulations, could potentially initiate interfacial cracks that propagate along the pilar axis to cause delamination. An extension of the present study to 3D could shed light on this mechanism. In addition, while we treated the a-Li x Si phase as a perfect elastic-plastic solid using von Mises theory (J 2 formulation) as in previous studies 23 , other formulations of plastic deformation have been proposed [56][57][58][59] . Those formulations have been motivated by results of DFT-based atomistic simulations showing that stresses generated during lithiation of thin films remain below the yield strength of a-Li x Si computed under uniaxial tension 56 . They include a reaction-flow formulation 57 , which extends von Mises theory to include, in addition to the magnitude of deviatoric stresses, a contribution of the chemical driving force for the solid-state reaction to the yield surface, and a different formulation without explicit flow where Li insertion, also driven by a combination of deviatoric stresses and chemical potential, produces an anisotropic compositional expansion 58,59 . Those more elaborate formulations could be implemented in the computational framework developed in the present study to test their effects on fracture behavior. Since they have been mainly validated by comparison with stress measurements in thin films 56,58 , which are already reasonably well-described by von Mises theory with the yield strength as a fit parameter 60 , it is unclear if they will produce dramatically different fracture behaviors. We expect quantitative rather than qualitative differences. Finally, fracture behaviors may also depend on the choice of phase field formulation of fracture. While the formulation used in this work 26,27 does not explicitly consider surface stresses 61 , it has been shown to reproduce non-trivial crack propagation phenomena in varied contexts 29,30 , consistent with Griffith theory, and to quantitatively model crack initiation at V-shape notches 37 without consideration of surface stresses whose effects in swelling-driven fracture of Si pillars remain to be investigated.

METHODS
All the equations are solved using the Galerkin Finite Element Method. Furthermore, to ensure the robustness of solution, only 1/4 of each geometry was simulated and appropriate boundary conditions are applied on the symmetry axes. Since the phase-field fracture method requires that the spatial resolution of discretization resolves the small regularization length ξ and because of memory requirements for storage of history variables (elements of tensorb e in Eq. (6) stored at each quadrature point), the resulting nonlinear problem is often very large (0.75 million degrees of freedom for 2D simulations presented in this article). This necessitates the use of a parallel programming paradigm that utilizes complex numerical tools, thereby limiting the use of commercially available finite element software. We therefore, use our own code, which also allows us precise control of the solution algorithm and has been previously validated experimental observation of brittle fracture 62 . Our implementation is based on PETSc 63 as the linear algebra backbone and libMesh 64 for finite elements bookkeeping. Equation (8) is solved using a Newton method where we calculate the consistent tangent moduli explicitly at each iteration. Furthermore, Eq. (10) is integrated explicitly using 50 substeps during each time step to ensure the accuracy and stability of integration. Table 1 summarizes the values of different parameters used in our simulations.

DATA AVAILABILITY
The data sets generated during the current study are available from the corresponding author on reasonable request.