Vulnerable Window of Yield Strength for Swelling-Driven Fracture of Phase-Transforming Battery Materials

TL;DR

This paper uses a phase-field model to identify a critical yield strength window where swelling-driven fractures occur in phase-transforming battery materials, revealing different fracture modes and the impact of geometry on stability.

Contribution

It introduces a multi-physics phase-field approach to explicitly model fracture in phase-transforming battery materials, advancing understanding of fracture mechanisms beyond Griffith's criterion.

Findings

Identification of a vulnerable yield strength window for fracture

Different fracture modes depending on yield strength

Hollow nanopillar geometry enhances fracture resistance

Abstract

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 results reveal the existence of a vulnerable window of yield strength inside which pillars fracture during lithiation. They identify 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, and highlight the importance of taking into account this localization to accurately predict the onset of fracture within Griffith theory. The results 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 mechanically stable phase-transforming battery materials undergoing large volume changes.

Tables

Table 1. Table 1: Material properties of a- \ce Si and c- \ce Si used in the simulations.

Material property	Value	Units
${\kappa }_c$ [23]	108	$\mathrm{GPa}$
${\kappa }_a$ [23]	10.8	$\mathrm{GPa}$
${\mu }_c$ [23]	50	$\mathrm{GPa}$
${\mu }_a$ [23]	5	$\mathrm{GPa}$
$K$	50	$\mathrm{MPa}$
$G_c$ [39]	6	$J{m}^{-2}$
$w$	0.85	$\mathrm{nm}$
$\beta$	0.7	-

Equations

$$F_{ij}=\sqrt{J_{\psi }}{F}_{ik}^e{F}_{kj}^p$$

$$\mathcal{F}(F^e,\psi ,\varphi )=$$

$$+$$

$$W^{+}(F^e,\psi )\{\begin{array}{ll}\frac{\mu (\psi )}{2}\left({\hat{b}}_{kk}^e-2\right)+\frac{\kappa (\psi )}{4}\left({J^e}^2-2\ln (J^e)-1\right)& \mathrm{if}{J}^e\ge 1\\ \frac{\mu (\psi )}{2}\left({\hat{b}}_{kk}^e-2\right)& \mathrm{otherwise}\end{array}$$

$$W^{-}(F^e,\psi )=\{\begin{array}{ll}0& \mathrm{if}{J}^e\ge 1\\ \frac{\kappa (\psi )}{4}\left({J^e}^2-2\ln (J^e)-1\right)& \mathrm{otherwise}\end{array}$$

$$\frac{\delta \mathcal{F}}{\delta u_i}$$

$$\frac{\delta \mathcal{F}}{\delta \varphi }$$

$$\frac{\partial \psi }{\partial t}=-M\left(\frac{{\nabla }_0\psi }{|{\nabla }_0\psi |}\right)\left(\frac{1}{e_0}\frac{\delta \mathcal{F}}{\delta \psi }-w^2|{\nabla }_0\psi |{\nabla }_0\cdot \left(\frac{{\nabla }_0\psi }{|{\nabla }_0\psi |}\right)\right)$$

$${\tau }_c=\sqrt{\frac{G_c\mu }{\mathcal{C}a(1-\nu )}}$$

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

Full text

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\text{\,}\mathrm{mA}\text{\,}\mathrm{h}\text{\,}{\mathrm{g}}^{-1}$ for \ceLi6C [1]). Silicon and germanium have an order of magnitude larger theoretical capacity ($3579\text{\,}\mathrm{mA}\text{\,}\mathrm{h}\text{\,}{\mathrm{g}}^{-1}$ for \ceLi15Si4, $4200\text{\,}\mathrm{mA}\text{\,}\mathrm{h}\text{\,}{\mathrm{g}}^{-1}$ for \ceLi22Si5 [1], $1384\text{\,}\mathrm{mA}\text{\,}\mathrm{h}\text{\,}{\mathrm{g}}^{-1}$ for \ceLi15Si4 [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-LixSi) [19, 20, 2]. 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-LixSi 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 $\psi$ 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 $\phi$, 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 elasto-plasticity combining neo-Hookean nonlinear elasticity and $J_{2}$ plasticity to quasi-statically evolve $\phi$ 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 chemo-mechanical 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 $\phi$ varies smoothly in space on a length scale $\xi$, 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 $\xi$ 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, 22, 23].

To keep computations tractable, we perform 2D plane-strain simulations ($\partial_{z}\equiv 0$) on a cross-section of an unconstrained nanopillar ($\tau_{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 $\alpha$ only vary radially and are independent of the azimuthal angle $\theta$ as depicted schematically in Figure 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 $\sigma_{y}$. Even though the critical $\sigma_{y}$ value corresponding to this maximum is outside the experimentally [8, 38, 39] or theoretically estimated [40] range $\sigma_{y}\sim$0.52\text{,}\mathrm{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.52\text{\,}\mathrm{GPa}$. For this reason, we investigate stress evolution over a wide range of $\sigma_{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 $\alpha$ 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 $\sigma_{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 ($\sigma_{y}=$0.52\text{\,}\mathrm{GPa}) of Si, the difference between the material and the most critical yield strength controls the magnification of generated tensile stresses. In a 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.

Results and discussion

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-LixSi interface during lithiation using a non-conserved phase field $\psi$ where $\psi=0$ in the crystalline phase and $\psi=1$ in the amorphous phase. The material properties are approximated using a linear role of mixture between the crystalline and amorphous phases. We define the deformation gradient tensor as $F_{ij}=\partial x_{i}/\partial X_{j}=1+\partial u_{i}/\partial 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

[TABLE]

where $J_{\psi}=(1+\beta\psi)^{2}$ is the phase dependent volumetric expansion due to phase change with linear Vegard expansion coefficient $\beta$, $F^{p}$ is the plastic deformation, and $F^{e}$ is the elastic deformation. We use the framework of the phase-field method for fracture [27, 26] by introducing the fracture phase field $\phi$ along with the process zone size $\xi$. To forbid interpenetration of the fracture faces similar to [32], we write the free energy with using an isochoric-volumetric split in undeformed coordinates as:

[TABLE]

where ${\nabla_{0}\bullet}=(\partial\bullet/\partial X_{i})$, $G_{c}$ is the fracture energy and $e_{0}$ is the energy cost phase change for unit undeformed volume, and $w$ is the characteristic phase change thickness. Furthermore, we assume a neo-Hookean hyper-elastic material

[TABLE]

[TABLE]

where for $J^{e}=det(F^{e})$, we define isochoric left Cauchy-Green deformation tensor $\hat{b}^{e}_{ij}={F^{e}_{ik}F^{e}_{jk}}/{J^{e}}$. In this model, the shear modulus $\mu(\psi)=\psi\mu_{a}+(1-\psi)\mu_{c}$ and bulk modulus $\kappa(\psi)=\psi\kappa_{a}+(1-\psi)\kappa_{c}$ are extrapolated between the shear and bulk moduli of amorphous ($\mu_{a},\kappa_{a}$) and crystalline phase ($\mu_{c},\kappa_{c}$). Following classic $J_{2}$ plasticity we assume that the von Mises equivalent stress $\tau_{eq}=\sqrt{3s_{ij}s_{ij}/2}\leq(\sigma_{y}+K\alpha)$ where $s_{ij}=\mu(\psi)(\hat{b}^{e}_{ij}-\hat{b}^{e}_{kk}\delta_{ij}/N)$ is the deviatoric part of the Kirchhoff stress tensor $\tau$ ($\tau=J\sigma$ where $\sigma$ is the Cauchy stress tensor), $\sigma_{y}$ is the yield strength and $\alpha$ is the isotropic hardening parameter (at the infinitesimal strain limit the isotropic hardening parameter reduces to the equivalent plastic strain). We write the governing equations for the displacement and fracture phase field as the minimizers of (Model):

[TABLE]

where $\delta\mathscr{F}/\delta\bullet$ is the Fréchet derivative of the free energy $\mathscr{F}$ with respect to field $\bullet$. Moreover, we model the curvature-independent interface-reaction-controlled anisotropic motion of the c-Si/a-Si interface as [44]:

[TABLE]

where $M({\nabla_{0}\psi}/{|\nabla_{0}\psi|})$ is the anisotropic mobility of phase change where we use the same form of the anisotropic mobility $M({\nabla_{0}\psi}/{|\nabla_{0}\psi|})$ as An et al. [23].

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 Figure 1. Figure 1(c-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 ($\psi=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 $\sigma_{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 (Figure 1(c-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, 22, 23]. Importantly, for $\sigma_{y}$ smaller than approximately $4\text{\,}\mathrm{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 (Figure 1c). In this range ($\sigma_{y}\leq$4\text{,}\mathrm{GPa}$$), the maximum hoop stress reached during complete lithiation, $\max(\tau_{\theta\theta})$, increases linearly with $\sigma_{y}$ as shown in Figure 1b for $t/R=1$ corresponding to a solid pillar; since the outer boundary is traction free ($\tau_{rr}\equiv 0$) and the von Mises stress is capped by $\sigma_{y}$, $\max(\tau_{\theta\theta})=2\sigma_{y}/\sqrt{3}$ on that boundary for plane-strain. In contrast, for larger $\sigma_{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 $\sigma_{y}$ as shown in Figure 1b. We should also highlight that although all simulations presented in this article were performed using $\beta=0.7$ that corresponds to $\sim 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 $\sigma_{y}/\mu_{a}\beta$ (see Figure \zrefsfig:stt-vs-sy), 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 $\sigma_{y}$ plot in Figure 1b suggests the existence of a vulnerable window for fracture corresponding to the range of $\sigma_{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 (Figures 3–5). Figure 1b also shows plots of $\max(\tau_{\theta\theta})$ versus $\sigma_{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 Figure 3 for $\sigma_{y}=$1\text{,}\mathrm{GPa} and $\sigma_{y}=$10\text{\,}\mathrm{GPa} ($\sigma_{y}=$7\text{,}\mathrm{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 Figure 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 $\sigma_{y}$ reported, and V-shape notches that form for lower $\sigma_{y}$ due to tensile yielding.

The 2D results in Figure 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 $\sigma_{y}$ and the $[001]$ and $[112]$ crystallographic orientations and for fixed process zone size to radius ratio $\xi/R=0.02$. The time evolutions of particle morphologies in Figure 3 show that nanopillar fracture during lithiation over an intermediate range of $\sigma_{y}$ inside the vulnerable window centered at the critical yield strength. For $\sigma_{y}=$1\text{,}\mathrm{GPa}, V-shaped notches are created along the crystalline corner directions but the magnitude of surface tensile stresses are insufficient for crack initiation. For $\sigma_{y}=$3\text{\,}\mathrm{GPa}, cracks initiate due to stress concentration at V-shaped notches and propagate unstably towards the crystalline core. For $\sigma_{y}=$5\text{,}\mathrm{GPa}, tensile yielding and hence V-shaped notches are absent but tensile stresses grow sufficiently large (at a later stage of charging compared to $\sigma_{y}=$3\text{\,}\mathrm{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 $\sigma_{y}=$10\text{,}\mathrm{GPa} ($\sigma_{y}=$7\text{\,}\mathrm{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.

Size effects

Experimental observations of Si anodic components show a clear size dependency where, for example, nanospheres with radii smaller than $\sim$75\text{,}\mathrm{nm}$$ [5] and nanopillars with radii smaller than $120\text{\,}\mathrm{nm}$ [6] do not break during lithiation. Such size dependencies, prevalent in brittle and quasi-brittle materials [45, 46, 47, 48], 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 $\tau_{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=G_{c}$), we obtain the expression

[TABLE]

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\sim R$ and $a\ll R$). In the phase-field model used in this article, which describes the state of the material with a spatially varying scalar field $\phi$, crack nucleation is an inherent property of the model and occurs via an instability that causes $\phi$ to develop a local dip ($\phi\rightarrow 0$) when the local stress exceeds a critical value $\tau_{c}\sim\sqrt{{G_{c}\mu}/\xi}$ [49] (up to a numerical prefactor that also depends on particle geometry and load configuration). Consequently, by comparing the above scaling expression for $\tau_{c}$ to Eq. (8), we can physically interpret $\xi$ 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{{3G_{c}\mu}/{4\xi(1-\nu)}}$ [49], we obtain the constant $\mathcal{C}=4/3$ by comparison with Eq. 8 with $a=\xi$. 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 $\phi$ 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 $G_{c}/(\mu_{a}R)$ where $\mu_{a}$ is the shear modulus of the amorphous a-\ceLi_xSi phase. The dimensionless fracture energy can be readily interpreted as the ratio of the Griffith length scale $G_{c}/\mu_{a}$ and the nanopillar radius. Figure 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 $G_{c}/(\mu_{a}R)$ and $\sigma_{y}$ where fracture does (red circles) or does not (blue crosses) occur. These results confirm the existence of a vulnerable window of $\sigma_{y}$ where cracks initiate and propagate in nanopillars during lithiation. This window is centered around the critical value of $\sigma_{y}$ generating maximal tensile stresses and shrinks in size with increasing $G_{c}/(\mu_{a}R)$.

We now assess if the onset of fracture in 2D phase-field simulations can be predicted within Griffith theory by using Eq. (8) together with values of the maximum hoop stress during lithiation obtained in 1D or 2D phase-field simulations without fracture (Figure 2). For this, we assume stress-free surfaces ($\tau_{rr}(R)=\tau_{r\theta}(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\mathrm{E}_{a}$ where $\mathrm{E}_{a}=9\kappa_{a}\mu_{a}/(3\kappa_{a}+\mu_{a})$, $\nu_{a}=(3\kappa_{a}-2\mu_{a})/(6\kappa_{a}+2\mu_{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. (3)). We can then rewrite Eq. (8) as a function of the maximum hoop stress as

[TABLE]

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 Figure 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}/(\mu_{a}R)$ for fracture at fixed $\sigma_{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. (9) 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 (Figure 2). Therefore, we can reasonably use Eq. (9) 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 Figure 2a, respectively) to obtain lower and upper bounds for the fracture boundary in the plane of $G_{c}/(\mu_{a}R)$ and $\sigma_{y}$. The results are depicted by the gray shaded region in Figure 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 $\sigma_{y}\sim$0.52\text{,}\mathrm{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 $\sigma_{y}$ over the estimated range $\sigma_{y}=$0.52\text{,}\mathrm{GPa} for a-\ceLi_xSi. The experimental range of safe nanopillar radius is highlighted by the green shaded region in Figure [4](#Sx2.F4)a. This region was computed using the experimentally observed safe nanopillar radius ($120\text{\,}\mathrm{nm}$) [[6](#bib.bib6)] and the estimated range of fracture energy ($57\text{\,}\mathrm{J}\text{\,}{\mathrm{m}}^{-2}$) from experimental measurements [[39](#bib.bib39)]. 2D phase-field simulations predict that, inside this green shaded region of Figure [4](#Sx2.F4)a, fracture occurs for $\sigma_{y}$ between $1.5\text{\,}\mathrm{GPa}$ (blue crosses) and $2\text{\,}\mathrm{GPa}$ (red circles), which falls in the upper part of the range $\sigma_{y}=$0.52\text{\,}\mathrm{GPa} estimated from experimental measurements [8, 38, 39]. Phase-field modeling prediction also depend generally on flaw size through the ratio $\xi/R$. To test this dependence, we repeated a series of simulations by varying $\xi/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}/(\mu_{a}R)=0.01$ calculated using the average reported fracture energy for a-\ceLi_xSi $G_{c}=$6\text{,}\mathrm{J}\text{,}{\mathrm{m}}^{-2} [[39](#bib.bib39)] and the observed safe nanopillar radius $R\simeq$120\text{\,}\mathrm{nm} [6]. The results reported in Figure 4b show that nanopillars become more vulnerable to fracture with increasing $\xi/R$ as theoretically expected. For the largest value $\xi/R=0.2$ studied here, fracture occurs for $\sigma_{y}$ between $0.5\text{\,}\mathrm{GPa}$ (blue cross) and $1\text{\,}\mathrm{GPa}$ (red circle), which falls in the lower part of the range $\sigma_{y}=$0.52\text{,}\mathrm{GPa}$$ estimated from experimental measurements [8, 38, 39]. While the precise value of $\xi/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 [20, 6].

We can further confirm the above estimate of $\sigma_{y}$ by redoing the calculations for isotropic lithiation of amorphous Si. Experimental observations demonstrate that a-Si lithiates isotropically [50, 51, 52]. 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. (9) further using the expression $\max(\tau_{\theta\theta})\simeq 2\sigma_{y}/\sqrt{3}$ valid for small $\sigma_{y}$ (see Figure 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-\ceLi_xSi. 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 model [49]. We can use this isotropic estimate along with the safe experimentally observed safe nanopillar radius $1\text{\,}\mathrm{\SIUnitSymbolMicro m}$ [52] to calculate an upper bound for its yield strength $0.41.2\text{\,}\mathrm{GPa}$ using the process zone size $0.02\leq\xi/R\leq 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 $\mu_{a}\simeq$19\text{,}\mathrm{GPa}$$ we can estimate $\sigma_{y}$ in the range $1.54.6\text{\,}\mathrm{GPa}$. We calculated the above range similarly using the process zone size $0.02\leq\xi/R\leq 0.2$ and based on observed safe nanopillar radius of $250\text{\,}\mathrm{nm}$ [53].

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 $\sigma_{y}$ (*i.e., *the decrease of the peak value of plots in Figure 1b with increasing annulus slenderness). However, for low yield strength, the maximum hoop stress remains bounded by $\sigma_{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 $85170\text{\,}\mathrm{nm}$ and different slenderness $0.08\leq t/R\leq 1$ for $\sigma_{y}=$3\text{,}\mathrm{GPa}. The results illustrated in Figure [5](#Sx2.F5) show that increased slenderness has a protective effect for this $\sigma_{y}$ value. For example, nanopillars break for $t/R=0.4$ but only exhibit minor cracking near the interior surface of $t/R\simeq 0.29$. Even more slender nanopillars ($t/R\simeq 0.08$) do not fracture. This protective geometrical effect, however, does not persist at lower yield strength for $\xi/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 $\sigma_{y}=$1\text{\,}\mathrm{GPa} (Figure 1b). Our 2D simulations confirm this prediction by showing that the hollow nanopillar for $t/R\simeq 0.29$, which is protected for $\sigma_{y}=$3\text{,}\mathrm{GPa}, fractures for $\sigma_{y}=$1\text{\,}\mathrm{GPa}.

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 energy [39] and observations of size dependent fracture [6, 5, 52], our results yield an estimate of yield strength within a range $\sigma_{y}\simeq$0.52\text{,}\mathrm{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, 6, 5]. 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 $\sigma_{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].

Methods

All the equations are solved using 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 axises. Our implementation is based on PETSc [54] as the linear algebra backbone and libMesh [55] for finite elements bookkeeping. Eq. (5) is solved using a Newton method where we calculate the consistent tangent moduli explicitly at each iteration. Furthermore, Eq. (7) 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.

Author contributions

A.M. and A.K. conceived the theoretical study and jointly interpreted the numerical results and wrote the paper. A.M. carried out the numerical study.

Acknowledgements

A.M. and A.K. acknowledge the support of Grant No. DE-FG02-07ER46400 from the U.S. Department of Energy, Office of Basic Energy Sciences. The majority of the numerical simulations were performed on the Northeastern University Discovery cluster at the Massachusetts Green High Performance Computing Center (MGHPCC).