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 Li6C (ref. 1)). Silicon and germanium have an order of magnitude larger theoretical capacity (3579 mA h g−1 for Li15Si4, 4200 mA h g−1 for Li22Si5 (ref. 1), 1384 mA h g−1 for Li15Si4 (ref. 2)) but are prone to fracture due to the high, approximately 300%, volume expansion during lithiation3,4, which limits their use. Different designs have been explored to overcome this limitation including silicon nanopillars5,6, thin films7,8,9, open nano-porous crystalline Si structures with ultra-high interfacial area produced by dealloying of Si-based alloys10,11, combinations of these12,13, as well as composite designs that embed silicon particles inside a more mechanically stable matrix14,15,16,17,18.

Basic studies of the lithiation process have shown that crystalline silicon (c-Si) transforms to an amorphous lithiated alloy (a-LixSi)2,19,20. The kinetics of this large volume change phase transformation is understood to be both interface-reaction limited20 and highly anisotropic, reflecting the two key observations that the velocity of the c-Si/a-LixSi interface remains approximately constant during lithiation, and that this velocity depends strongly on crystallographic orientation19. 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 fracture3,21,22,23. Those insights have already proven useful to test new anode designs to mitigate fracture23. 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 nanopillars13) 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 dimension24,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-LixSi phases and is evolved dynamically to describe the interface-reaction-limited anisotropic motion of the c-Si/a-LixSi 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 material26,27. To realistically model large volume changes, this variational approach is implemented using a large deformation formulation of elasto-plasticity combining neo-Hookean nonlinear elasticity and J2 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 transformations28 and fracture problems such as thermal shock fracture29, mixed mode fracture30, ductile fracture31,32, and the simpler chemo-mechanical fracture of single-phase battery cathode particles33,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 parameter37. V-shape notches, in particular, bear close similarity to surface shape deformations of lithiated Si particles undergoing elasto-plastic deformation during anisotropic swelling3,21,22,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 deformation24,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 experimentally8,38,39 or theoretically estimated40 range σy ~ 0.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.

Fig. 1: Axisymmetric results.
figure 1

Results of axisymmetric simulations of lithiation of hollow cylindrical crystalline nanopillars of outer radius R = 85 nm and variable thickness t. a Schematic representation of the hollow pillar geometry. b 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-LixSi interface located at ψ = 0.5), the Kirchhoff hoop stress \(\sqrt{3}{\tau }_{\theta \theta }/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.

Full size image

Results

Model

We model the swelling-driven deformation of the material using the finite J2 elasto-plasticity framework41,42,43 and account for the fracture of the material by coupling it to a phase-field fracture model26,27. Furthermore, we model the anisotropic motion of the c-Si/a-LixSi 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 Fij = ∂xi/∂Xj = 1 + ∂ui/∂Xj where Xi are the undeformed coordinates and xi = Xi + ui 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

$${F}_{ij}=\sqrt{{J}_{\psi }}{F}_{ik}^{e}{F}_{kj}^{p}$$
(1)

where Jψ = (1 + βψ)2 is the phase dependent volumetric expansion due to phase change with linear Vegard expansion coefficient β, Fp is the plastic deformation, and Fe is the elastic deformation. We use the framework of the phase-field method for fracture26,27 by introducing the fracture phase field ϕ along with the process zone size ξ and combine it with the standard phase-field approach for phase transformations to prescribe the ψ evolution44. Accordingly, we write the total free energy functional as the sum

$${\mathscr{F}}({F}^{e},\psi ,\phi )={{\mathscr{F}}}_{el}({F}^{e},\psi ,\phi )+{{\mathscr{F}}}_{\phi }(\phi )+{{\mathscr{F}}}_{\psi }(\psi ),$$
(2)

of contributions corresponding to the elastic energy stored in the material including the amorphous and crystalline phases (\({{\mathscr{F}}}_{el}\)), the energy associated with the creation of new fracture surfaces (\({{\mathscr{F}}}_{\phi }\)), and the free-energy associated with the c-Si to a-LixSi phase transformation (\({{\mathscr{F}}}_{\psi }\)), defined as

$${{\mathscr{F}}}_{el}({F}^{e},\psi ,\phi )={\int}_{\!{\Omega }_{0}}\left(({\phi }^{2}+{\eta }_{\xi }){W}^{+}({F}^{e},\psi )+{W}^{-}({F}^{e},\psi )\right)\ dx,$$
(3)
$${{\mathscr{F}}}_{\phi }(\phi )=\frac{3{G}_{c}}{8\xi }{\int}_{\!{\Omega }_{0}}\left(1-\phi +{\xi }^{2}| {\nabla }_{0}\phi {| }^{2}\right)\ dx,$$
(4)
$${{\mathscr{F}}}_{\psi }(\psi )={\int}_{\!{\Omega }_{0}}\left(h{f}_{dw}(\psi )-\Gamma {f}_{tilt}(\psi )/3+{\sigma }_{\psi }| {\nabla }_{0}\psi {| }^{2}\right)\ dx.$$
(5)

In the sharp-interface limit ξ/R ≪ 1 where the process zone size ξ is much smaller than the nanopillar radius R, \({{\mathscr{F}}}_{\phi }(\phi )\) reduces to the standard fracture energy Gc of convenctional fracture mechanics, defined as the energy required to create new fracture surfaces per unit extension of the crack, where we have defined ∇0• = (∂•/∂Xi). For the phase transformation contribution \({{\mathscr{F}}}_{\psi }(\psi )\), fdw(ψ) = 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=\sqrt{{\sigma }_{\psi }/h}\) and excess free-energy γ = 2hw/3. In addition, ftilt(ψ) provides the driving force for the solid-state chemical reaction xLi+ + xe− + Si → LixSi by tilting the double well potential (hfdw(ψ) − Γftilt(ψ)/3) so as to lower the bulk free-energy density of a-LixSi (ψ = 1) relative to c-Si (ψ = 1) by an amount proportional to \(\Gamma =x({\tilde{\mu }}_{{\rm{Li}}}-{\mu }_{{\rm{Li}}}^{a})+{\mu }_{{\rm{Si}}}^{c}-{\mu }_{{\rm{Si}}}^{a}-xF\Delta V\ge 0\), where \({\mu }_{{\rm{Li}}}^{a}\) (\({\mu }_{{\rm{Si}}}^{a}\)) is the chemical potential of Li (Si) in a-LixSi at the interface, \({\tilde{\mu }}_{{\rm{Li}}}\) and \({\mu }_{{\rm{Si}}}^{c}\) 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 ftilt(ψ) = 15(32ψ3/3 − 16ψ4 + 32ψ5/5)/16, with limits ftilt(0) = 0 and ftilt(1) = 1 and vanishing first and second derivatives at ψ = 0 and ψ = 1, preserves the free-energy minima at ψ = 0 and ψ = 1 for arbitrary Γ. Finally, to construct \({\mathscr{F}}({F}^{e},\psi ,\phi )\), we split the elastic energy into positive (W+(Fe, ψ)) and negative (W−(Fe, ψ)) volumetric parts to distinguish material regions under tension and compression, respectively, and assume a neo-Hookean hyper-elastic material

$${W}^{+}({F}^{e},\psi )=\left\{\begin{array}{*{20}{l}}\frac{\mu (\psi )}{2}\left({\hat{b}}_{kk}^{e}-2\right)+\frac{\kappa (\psi )}{4}\left({{J}^{e}}^{2}-2\,{\mathrm{ln}}\,({J}^{e})-1\right)&{\rm{if}}\ {J}^{e}\ge 1\\ \frac{\mu (\psi )}{2}\left({\hat{b}}_{kk}^{e}-2\right)&{\rm{otherwise}}\end{array}\right.$$
(6)
$${W}^{-}({F}^{e},\psi )=\left\{\begin{array}{*{20}{l}}0&{\rm{if}}\, {J}^{e}\ge 1\\ \frac{\kappa (\psi )}{4}\left({{J}^{e}}^{2}-2\,{\mathrm{ln}}\,({J}^{e})-1\right)&{\rm{otherwise}}\end{array}\right.$$
(7)

where Je = det(Fe), \({\hat{b}}_{ij}^{e}={F}_{ik}^{e}{F}_{jk}^{e}/{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-LixSi (μ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 J2 plasticity we assume that the von Mises equivalent stress \({\tau }_{eq}=\sqrt{3{s}_{ij}{s}_{ij}/2}\le ({\sigma }_{y}+K\alpha )\) where \({s}_{ij}=\mu (\psi )({\hat{b}}_{ij}^{e}-{\hat{b}}_{kk}^{e}{\delta }_{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 ui and ϕ

$$\frac{\delta {\mathscr{F}}}{\delta {u}_{i}}=0\quad s.t.\ \sqrt{3{s}_{ij}{s}_{ij}/2}\le ({\sigma }_{y}+K\alpha ),$$
(8)
$$\,\frac{\delta {\mathscr{F}}}{\delta \phi }=0,$$
(9)

respectively, where \(\delta {\mathscr{F}}/\delta \bullet\) is the Fréchet derivative of the free energy \({\mathscr{F}}\) with respect to field •. Standard J2 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 self-consistently captured by the fact that (8) and (9) are derived from the same free-energy functional, which couples the evolution of the displacement field ui 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

$$\frac{\partial \psi }{\partial t}=-M\left(\hat{n}\right)\left(\frac{\delta {{\mathscr{F}}}_{\psi }}{\delta \psi }+2{\sigma }_{\psi }| {\nabla }_{0}\psi | {\nabla }_{0}\cdot \hat{n}\right)$$
(10)

where \(\hat{n}={\nabla }_{0}\psi /| {\nabla }_{0}\psi |\) is the direction normal to the interface and \(M(\hat{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 ≪ 145, (10) reduces to the simple law of curvature-independent interface motion \({V}_{n}=M(\hat{n})\Gamma\) with the interface velocity Vn controlled by the thermodynamic driving force for the interface reaction Γ assumed to be spatially uniform. As discussed previously23, 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(\hat{n})\). Hence, elastoplastic deformation and fracture are independent of the value of the driving force Γ that only sets the unit of time, i.e., Γ can be scaled out of the problem by measuring time in units of \(R/({M}_{\max }\Gamma )\) where \({M}_{\max }\) is the maximum value of \(M(\hat{n})\).

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-LixSi 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 cracking3,21,22,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 ({\tau }_{\theta \theta })\), 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 ({\tau }_{\theta \theta })=2{\sigma }_{y}/\sqrt{3}\) on that boundary for plane-strain. 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 ({\tau }_{\theta \theta })\) 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 ({\tau }_{\theta \theta })/{\mu }_{a}\beta\) 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 ({\tau }_{\theta \theta })\) 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 ({\tau }_{\theta \theta })\) 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 ({\tau }_{\theta \theta })\), 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 ({\tau }_{\theta \theta })\) 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 ({\tau }_{\theta \theta })\) 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 compressive yielding, which occurs for all σy reported, and V-shape notches that form for lower σy due to tensile yielding.

Fig. 2: Effect of anisotropy on maximum hoop stress.
figure 2

Comparison of 2D and 1D axisymmetric simulations of lithiation of solid nanopillars of radius R = 85 nm without fracture. Plots of maximum hoop stress τθθ vs. yield stress σy for [001] (a) and [112] (b) oriented nanopillars (r is the radial coordinate in the undeformed frame). In 2D simulations, localization of plastic flow at orientations corresponding to sharp corners of the crystalline Si core and concomitant creation of V-shaped notches (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).

Full size image

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, V-shaped 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. In [001] oriented pillars, crack initiate at the four corners corresponding to [100] and [010] directions where stresses are largest. In [112] oriented pillars, crack initiate at the two corners along \([11\underline{1}]\) directions.

Fig. 3: Phase-field simulations of nanopillar fracture.
figure 3

Phase-field simulations of fracture during the lithiation of [001] (a) and [112] (b) 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).

Full size image

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 hardening8. 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-LixSi plastic behavior.

Size effects

Experimental observations of Si anodic components show a clear size dependency where, for example, nanospheres with radii smaller than ~75 nm5 and nanopillars with radii smaller than 120 nm6 do not break during lithiation. Such size dependencies, prevalent in brittle and quasi-brittle materials46,47,48,49, are typically characterized by a power law relationship between the stress to fracture and component size, \({\tau }_{c} \sim 1/\sqrt{R}\), 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={\mathcal{C}}a{\sigma }^{2}(1-\nu )/\mu\) for plane-strain where \({\mathcal{C}}\) is a dimensionless constant that generally depends on particle geometry and load configuration. Equating the energy release rate with the fracture energy (G = Gc), we obtain the expression

$${\tau }_{c}=\sqrt{\frac{{G}_{c}\mu }{{\mathcal{C}}a(1-\nu )}}$$
(11)

and thus the scaling \({\tau }_{c} \sim 1/\sqrt{R}\) 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 model and occurs via an instability that causes ϕ to develop a local dip (ϕ → 0) when the local stress exceeds a critical value \({\tau }_{c} \sim \sqrt{{G}_{c}\mu /\xi }\)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 \({\tau }_{c}=\sqrt{3{G}_{c}\mu /4\xi (1-\nu )}\)50, we obtain the constant \({\mathcal{C}}=4/3\) by comparison with Eq. (11) with a = ξ. This value is close to the standard LEFM value \({\mathcal{C}}=\pi /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 \(\bar{\xi }=\xi /R\). From the dimensional analysis, it is natural to define the dimensionless fracture energy as Gc/(μaR) where μa is the shear modulus of the amorphous a-LixSi phase. The dimensionless fracture energy can be readily interpreted as the ratio of the Griffith length scale Gc/μ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 \(\bar{\xi }=0.02\), which identifies regions of the two-dimensional parameters space Gc/(μaR) 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 Gc/(μaR).

Fig. 4: Vulnerable window of fracture.
figure 4

a Vulnerable window of fracture plotted for yield strength σy vs. dimensionless fracture energy (particle size) Gc/(μaR) 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 ({\tau }_{\theta \theta })\) using 1D axisymmetric results shown in Fig. 1b (dashed gray line), and \(\max ({\tau }_{\theta \theta })=2{\sigma }_{y}/\sqrt{3}\) (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 ({\tau }_{\theta \theta })\) 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 nm6 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 Gc/(μaR) = 0.01 (corresponding to R = 120 nm and fracture energy Gc = 6 J m−2). The results are consistent with the experimentally estimated range σy = 0.5–2 Gpa for fracture.

Full size image

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}^{+}\simeq {\tau }_{\theta \theta }^{2}(1-{\nu }_{a}^{2})/2{{\rm{E}}}_{a}\) where Ea = 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

$$\frac{{G}_{c}}{{\mu }_{a}R}={\mathcal{C}}\bar{\xi }(1-\nu ){\left(\frac{\max ({\tau }_{\theta \theta })}{{\mu }_{a}}\right)}^{2}$$
(12)

where \({\mathcal{C}}=4/3\). Substituting in the above expression the values of \(\max ({\tau }_{\theta \theta })\) 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 Gc/(μaR) 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 ({\tau }_{\theta \theta })\) 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 Gc/(μaR) 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 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-LixSi. 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 measurements39. 2D phase-field 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 measurements8,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 Gc/(μaR) = 0.01 calculated using the average reported fracture energy for a-LixSi Gc = 6 J m−2 39 and the observed safe nanopillar radius R ≃ 120 nm6. 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 measurements8,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 lithiation6,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 isotropically51,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 ({\tau }_{\theta \theta })\simeq 2{\sigma }_{y}/\sqrt{3}\) valid for small σy (see Fig. 1b), which yields the prediction \({G}_{c}/{\mu }_{a}R=(4{\mathcal{C}}/3)\bar{\xi }(1-\nu ){({\sigma }_{y}/{\mu }_{a})}^{2}\) for isotropic lithiation of a-LixSi. Zhao et al.24 obtained a similar expression of the form of \({G}_{c}/{\mu }_{a}R \sim {({\sigma }_{y}/{\mu }_{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 model50. We can use this isotropic estimate along with the safe experimentally observed safe nanopillar radius 1 μm53 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 nm54.

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 fracture13. 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.

Fig. 5: Fracture of hollow nanopillars.
figure 5

Final cross section and fracture pattern as a function of dimensionless fracture energy Gc/(μaR) and slenderness t/R for hollow nanopillars with equivalent cross-section area as solid nanopillars of radius R = 85, 121 and 170 (Gc = 6 J m−2, σy = 3 GPa, and ξ/R = 0.02). Color map depicts the hardening parameter α. The parameter range where the nanopillar fractured is shown with the red background. One can see that at this yield strength the more slender nanopillar mitigates the failure. Also increasing particle size (decreasing the fracture energy) promotes initiation of cracks in-line with the experimental observations (see Supplementary Videos 7a,b and 8a,b for results of simulations for Gc/(μaR) = 0.01 and t/R ≃ 0.29 and t/R = 0.4).

Full size image

Discussion

In summary, we have used a multi-physics phase-field approach to model simultaneously anisotropic phase transformation, elasto-plastic 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 energy39 and observations of size dependent fracture5,6,53, our results yield an estimate of yield strength within a range σy ≃ 0.5–2 GPa consistent with experimental8,38,39 and theoretical40 estimates. This range is smaller than the critical yield strength consistent with the observation of localized plastic deformation during lithiation of Si nanopillars3,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 nanopillars13 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 area11.

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 pillars55. 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-LixSi phase as a perfect elastic-plastic solid using von Mises theory (J2 formulation) as in previous studies23, other formulations of plastic deformation have been proposed56,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-LixSi computed under uniaxial tension56. They include a reaction-flow formulation57, 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 expansion58,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 films56,58, which are already reasonably well-described by von Mises theory with the yield strength as a fit parameter60, 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 work26,27 does not explicitly consider surface stresses61, it has been shown to reproduce non-trivial crack propagation phenomena in varied contexts29,30, consistent with Griffith theory, and to quantitatively model crack initiation at V-shape notches37 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 tensor \({\hat{b}}^{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 fracture62. Our implementation is based on PETSc63 as the linear algebra backbone and libMesh64 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.

Table 1 Material properties of a-Si and c-Si used in the simulations.
Full size table