A multiphysics model that can capture crack patterns in Si thin films based on their microstructure

Fracture in silicon anodes has fascinated the electrochemistry community for two decades, as it can result in a 80% capacity loss over the first few electrochemical cycles and is the limiting factor in commercializing such high capacity anodes. Although numerous experimental data exist illustrating severe fracture patterns and their dependence on the scale of the microstructure, no theoretical model has been able to re-produce and capture such behaviour. In this article, a multi-physics phase-field damage model is presented that can accurately capture the long standing problem of dry bed-lake crack patterns observed for Si thin film anodes. A promising aspect of the model is that, in addition to accounting for Li-ion diffusion, it can explicitly capture the microstructure, and therefore when applied to a Si film with a thickness below 100 nm no fracture was observed, which is consistent with experiments. As fracture in continuous thin films is random, micron-hole patterned Si films were also fabricated and cycled, resulting in ordered crack patterns. The proposed model was able to capture these elaborate, yet ordered, crack patterns, further validating its efficiency in predicting damage during lithiation of Si. This paves the way to using multiscale modeling for predicting the dimensions that limit and control fracture during lithiation, prolonging hence the electrode lifetime.


Introduction
Extensive experimental research has concluded that Si-based anodes are the most promising candidates for next generation Li-ion batteries. Upon the formation of Li-Si alloys the theoretical capacity can be as high as 4200mAh/g, which is eleven times higher than that of graphitic anodes (372mAh/g) that are used commercially. The limiting factor, however, in employing Si is that its high capacity cannot be retained due to the volume expansions (up to 400%) that occur upon lithiation and the subsequent fracture that develops. The most prominent example being dry bed-lake fracture that has been observed for Si and SiSn thin film anodes [1,2]. However, it should be noted that microstructure plays a key role in the extent of fracture. Experimental studies have shown that Si films below a critical thickness of 100 nm did not exhibit fracture. This was initially assumed in Ref. [3] since a high capacity retention was obtained for such film thicknesses, but was later proven by performing scanning electron microscopy (SEM) on amorphous Si film anodes that were 100 nm and 500 nm after 10 electrochemical cycles, showing that the 100 nm film remained undamaged while the 500 nm film experienced severe fracture [4]. Similar observations exist for Sn [5] and Si [6] which show that damage is significantly reduced as the particle size decreases.
Although this unique fracture mechanism of dry bed-lake fracture of metal anodes has captured the attention of the mechanics community https://doi.org/10.1016/j.jpowsour.2018.07.106 Received 15 February 2018; Received in revised form 15 July 2018; Accepted 27 July 2018 since the early 2000s, a model that can simulate it has yet to be proposed, due to the complexity of the problem. Initial studies used pure mechanical models, without explicitly accounting for Li-ion diffusion, to predict the mechanical stability upon maximum lithiation [7][8][9]. This allowed for the development of design criteria that predicted the most promising configurations, particle sizes and matrix materials, that would allow for an increased mechanical stability and hence longer anode cycle life. Although the predictions were consistent with experimental observations, the physical mechanisms that drove crack growth during cycling remained unexplored, and therefore subsequent mechanics models incorporated Li-ion diffusion and considered the fully coupled problem of diffusion and deformation. Despite the rigor of these models (a summary of which can be found in Ref. [10]), which employed, for example, generalized finite element models of diffusion, finite deformation, plasticity and fracture within the framework of cohesive zone modeling [11,12], they were not able to capture the dry bed-lake fracture of Si film anodes. More recently phase field modeling has also been used to model damage evolution during the Li-insertion process ( [13,14]), however, focus was given mostly on cathodes. When applied to Si film anodes [15] it was not possible to re-produce the experimentally obtained crack patterns.
The purpose of the present work is therefore to propose an experimentally validated multiphysics model that can capture fracture initiation and propagation throughout the continuous Li-insertion and deinsertion process. This is done by using a phase field model which accounts for the mechanical stresses resulting from the volume expansion of Si induced by the Li-ion diffusion. Applying the model to a Si thin film anode allowed to simulate the reknowned dry bed-lake fracture and its dependence on film thickness. To further establish the efficiency of the proposed model it is shown that it can also simulate ordered crack networks that form upon electrochemically cycling a new type of Si thin film that is patterned with micron-holes.

Sample preparation
A Ti thin film with a thickness of 500 nm was deposited on a polished quartz substrate by DC magnetron sputtering as the current collector for the electrochemical tests. The pure Si thin films were deposited on the Ti /quartz substrate by co-sputtering with pure silicon target (99.999% purity) and pure titanium target (99.5%). The thickness of the silicon film was 650 nm. In the case that the Si film was patterned with periodic arrays of micro-holes, a 1 mm-thick layer of positive photoresist S1813 (Shipley microposit S1813, Massachusetts, USA) was spin coated onto the Si film and pre-baked on a hot plate at 115°C for 60 s. After 42 s the samples were etched with a 50 nm per min rate by RIE (A PlasmaLab 80 Plus RIE system, Oxford Instruments Company, UK) with methane trifluoride gases. Finally × 4.5 4.5 μm micropore array silicon films where the distance between the micropores was 5.50 μm and the depth was 650 nm were obtained as illustrated in the Hitachi S-4800 SEM image Fig. 1(a).

Electrochemical testing and sample characterization
A Swagelok-type two-electrode cell was assembled in an Ar -filled glove box (MBraun, H O 2 and < O 0.1 2 ppm). The Si film was used as the working electrode, and Li foil was employed as the counter electrode. The electrolyte was 1 mol/L LiPF6 dissolved in ethylene carbonate (EC ≥ 99%): dimethyl carbonate (DMC ≥ 99%) (1:1, v/v). lithium hexafluorophosphate (LiPF6, ≥ 99.99%). The half cell was cycled between 5 mV and 2 V at room temperature using a Land BA2100A battery tester (Wuhan, China). All cells were galvanostatically discharged and charged at a current density of 0.1C. After the 1 st discharge (lithiation), 1 st charge (delithiation), and 10 th charge (delithiation), the cells were disassembled in an Ar -filled glove box and sealed in a special vacuum transfer box in order to transfer them to the SEM chamber for morphological studies.

Crack patterns observed experimentally in Si thin films
The scanning electron microscopy images of Fig. 1 depict the two types of Si film anodes considered: (i) continuous Si thin films of 500 nm thickness, (ii) Si films patterned with an array of micro holes. In both cases Ti was used as the current collector and therefore the Si films were deposited on it. Fig. 1(I) illustrates the change in thickness of the continuous Si film due to Li-ion insertion and Fig. 1(II) the crack network that formed upon the first insertion. Fig. 1(III), (IV) depict the microstructure upon the 1st delithiation, indicating decohesion at the Ti/Si interface. Inducing a hole pattern in the Si film ( Fig. 1(a) and (b)) allowed for some accommodation of the volume expansion during Li-insertion as the hole could relieve the accumulated stress the Li-ion diffusion exerted. Upon the maximum 1st lithiation, Fig. 1(c) and (d) illustrate that the hole diameter decreases from 4.5 μm to 3 μm and no failure occurs at the Ti/ Si interface. Upon the first delithiation ( Fig. 1(e) and (f)), it is seen that an ordered pattern of diagonal cracks connecting the holes is formed. This crack pattern is retained for the first 9 electrochemical cycles as seen in Fig. 1(g) and (h). However, after the tenth delithiation, a new failure network is obtained, which includes a system of disordered secondary cracks as seen in Fig. 1(i), (j), (k), (l). In the sequel a model is developed that can capture failure in both types of these Si anode microstructures.

Multi-physics model
The multi-physics model elaborated in this paper is based on a phase field damage model [16] coupled with the diffusion equation for Li-ion insertion and de-insertion. Our model exploits the capability of the phase field approach to capture failure in brittle materials. Li-ion insertion/de-insertion produces a volume change in the Si, and the resulting strain is added to the unknown elastic strain to form the total strain derived from the material displacement. Elastic strain is produced through Hooke's classical law for elastic solids, while Fick's law drives the ion-diffusion. Damage evolution is formulated using a variational approach of brittle fracture as suggested in Ref. [16]. A finite element resolution of the coupled problem is adopted for its flexibility to handle arbitrary geometry and boundary conditions.

Ion diffusion
It is assumed that the Li-ion diffusion inside the Si film can be modeled using the standard diffusion equation stating that the rate in Li-ion concentration is proportional to its Laplacian. If c denotes the normalized Li-ion concentration, ranging between 0 and 1, corresponding to full de-insertion and full-insertion, respectively, then the diffusion equation is written using Fick's law as: In this equation, ∂ ∂ c t denotes the rate change in concentration and D the diffusivity parameter. In a fully coupled chemo-mechanical problem, D may depend on the local stress state of the material and a source term proportional to the local stress should also be considered. These terms are not considered herein as they involve material constants that are difficult to measure and because good agreement with the experiments can be obtained even without them, as will be seen in the next section. To account for the saturation in Li-ions during insertion when c reaches 1, a non-linear diffusivity law is adopted. To produce sharp concentration profiles that lead to high stress levels [6], the dependency of D with respect to c is assumed to be: where D o is the initial diffusivity, when = c 0. Using this formula, the initial change of D with respect to c tends to zero and when c goes to 1 then D tends to infinity, and this prevents > c 1. For de-insertion, the same expression is adopted but replacing c with − c 1 . The diffusion equation is complemented by appropriate initial and boundary conditions. The initial condition is = c 0 and = c 1, for Liinsertion and de-insertion respectively. As illustrated in Fig. 2, the Liions are allowed to penetrate the Si-film over the entire surface that is in contact with the electrolyte, according to the physical configuration. A prescribed constant flux of Li-ions is considered over the upper electrode surface. In its present form, the resolution of the diffusion of the Li-ions into the Si-film is not coupled with its mechanical state, meaning that stress does not assist ions diffusion. However, during damage evolution electrolyte penetration through the cracks will be accounted for. Therefore, a flux of Li-ions will be considered in the damage zone. This flux is similar to that between the Si surface and the electrolyte.

Mechanical formulation
The volume changes that occur during Li-ion insertion or de-insertion are accommodated by the anode through a chemically induced strain which can take the form of the diagonal strain tensor: where I is the identity matrix and α a dilatation coefficient that sets the strain induced for full insertion ( = c 1). Even though it is well it is seen that no cracks form, (e) top view after the 1st de-insertion; crack pattern formation is seen; (f) magnification of hole after 1st de-insertion, (g) top view upon 10th maximum Li-insertion, it is seen that the crack pattern remains, (h) magnification of cracks at hole, (i) top view after the 10th de-insertion, (i) magnification of crack after 10th de-insertion, (k) interface between Si film and Ti substrate upon 10th insertion, (l) interface between Si and substrate after 10th de-insertion. established that anodes can encounter large strains, infinitesimal transformations are considered in the proposed approach. This choice is motivated by the level of complexity induced by a finite strain formulation and also because the results shown hereafter demonstrate very good agreement with experiments. In this context, the undeformed configuration is taken as a reference and a negative α is adopted for deinsertion.
From the aforementioned considerations, the total strain ε t of the material is decomposed as the sum of the chemically induced strain ε c and the elastic strain ε e : To model material failure, a variational approach based on regularized discontinuities is elaborated. A smoothed description of the crack discontinuity is introduced through the phase field parameter d and the internal length l c . d allows to account for material degradation locally and l c defines the distance from the crack surface over which the discontinuity is spread. l c should be somehow related to the characteristic length of the material microstructure. In the formulation described in Ref. [17], the variation of d is governed by the following equation: in which g c is the fracture energy of the material and a history variable that is the maximum value over time of the strain energy density associated to the extensive part of the elastic deformation + Ψ defined as where + . is the positive part of a real number. It has been assumed in Equation (6) that for an elastic strain state with a negative volume change ( < tr ε ( ) 0 e ) only the shear part of the strain energy promotes damage [16]. In the model, failure is thus driven by the strain energy density. As obtained from Equation (6) this is a quadratic form of the elastic strain ε e . For this strain contribution to increase there has to be a mismatch between the total strain that can be accommodated by the electrode design and the volume change induced by ion insertion/deinsertion. This occurs either when the volume change is not compatible with the mechanical constraint applied to the electrode (for example when it is bonded on a rigid substrate) or when a sharp concentration gradient occurs. The latter is promoted by a high ion flux (high intensity) and the non-linear diffusivity law (Equation (2)) governing ion concentration that is intrinsically related to a change in the atomic structure of the electrode. It has also been assumed in this equation that the elastic behaviour of Si is isotropic and thus characterized by the Lame coefficients λ and μ, the bulk modulus being defined as = + k λ μ 2 3 . While the initial crystallographic structure of the Si anode meads to orthotropic elastic behaviour, it has been shown in Ref. [18] that successive lithiation and delithiation strongly modify that atomic lattice. It has also been shown [19] that SiLi x alloys have an amorphous structure and therefore isotropic assumptions are justified. The Cauchy stress tensor is now related to the elastic strain ε e and to the field d as Using this formulation, instead of solving a non-linear two-field problem for the mechanical equilibrium and damage evolution, the use of a history variable in Equation (5) allows to recast the problem as two linear problems as in Ref. [17]. This is extremely convenient from the computational view point. Further, the history variable used in this formulation enables for strictly satisfying the irreversibility constrain of the damage process. Note also that even if huge volume changes (total strain ε t ) can be experienced by anodes during lithiation, this is mainly induced by the chemically induced strain ε c while the elastic strain ε e (from which stress arises) remains about one order of magnitude lower.

Results and discussion
The most pronounced fracture for Si anodes has been observed for the thin film configuration showin in Fig. 1 (I-IV). Hence, this was the first case which the model simulated. The area for the simulation is taken to be × 10 10 μm 2 , while two thicknesses will be considered in order to capture the size effect. The Li-ion flux entering the Si film is taken to be 1.4 nm/ms leading to a duration for one lithiation/de-lithiation cycle of about 20 h as in the experiments. On the film surface in contact with the electrolyte a stress free condition applies allowing for the volume expansion to occur, while the other surface that is deposited on the Ti, is rigidly fixed and the displacement there is kept zero. The finite element mesh consists from tetrahedral elements with 50 nm edges. The internal length parameter l c is related to the scale of the microstructure and is taken to be 100 nm. The time step size is set to 400 s and the initial diffusivity coefficient D o to 25 nm 2 /s, so that the cycle duration is the same as in the experiments. The elastic modulus and Poisson's ratio for Si are 130 GPa and 0.2 respectively. The dilatation coefficient α is set to 0.3 to match the experimental observations in terms of film thickness variation.

Preliminary results
In order to analyze how the model couples ions diffusion with mechanical stress, a first simulation is performed without computing damage evolution. The concentration and hydrostatic stress profile along a vertical line in the center of the thin film are plotted at different insertion rates in Fig. 3. It is clearly observed in Fig. 3(a) that the nonlinear diffusivity produces, as expected, sharp concentration gradients propagating from the inserted surface to the substrate. In Fig. 3(b), the subsequent stress generated is shown to increase where the concentration gradient is localized until the state of maximum Li-insertion is reached. From this Figure, it is also observed that a higher stress level is obtained close to the inserted surface and it can be concluded that the model will predict crack initiation at the inserted surface during the first insertion.

Electrochemical cycling of a 500 nm Si thin film that experiences fracture
In order to simulate the crack patterns observed in Fig. 1(III) and (IV), the fracture energy g c is set to 2000 J/m 2 . This value has been specifically selected so the simulation results match the experimental observations. Experimental values for the change of the fracture energy of Si during electrochemical cycling do not exist. During the delithiation process which basically consists of breaking atomic bonds, the fracture energy g c was taken to be five times smaller than that during lithiation. Indeed, the de-insertion of Li-ions produces vacancies in the atomic lattice, hence it is assumed to reduce the failure energy. Also, the material at the beginning of de-lithiation is completely different from the almost pure Si film that is lithiated initially. This is another reason why a different g c value can be adopted for delithiation. The value of the Young's modulus also depends on the Li-ion concentration [20], however, for simplicity it is assumed to be the same for both Si and Li-Si in this case.
Experimentally it is observed that for thin films a crack network initiates during the first lithiation and this network further develops during the subsequent delithiation (see Fig. 1(I-IV)). In Fig. 4, the simulation results obtained for the first electrochemical cycle are presented. In Fig. 4(a), the elements where d exceeds 0.95 are removed to visualize the cracks. The color corresponds to the vertical displacement in nm. A maximum displacement of about 500 nm is obtained meaning that the thickness of the film increases from 500 nm to about 1 μm, in good agreement with experimental observations of Fig. 1(I). In Fig. 4(b) which shows the iso-0.95 contour of d, one can observe how the crack network develops in the thickness of the thin film. It is to be emphasized that during lithiation most of the cracks do not propagate through the thin film. However, during de-lithiation, the results presented in Fig. 5 show that cracks reach the bottom surface of the thin film and then break the interface between the Si thin film and the Ti substrate as observed experimentally in Fig. 1(I-IV). Two movies showing the evolution of the vertical displacement and the local concentration on the cracked region during an electrochemical cycle are provided as Supplementary Material. One can thus appreciate how the simulations lead to the data plotted in Fig. 4 for full lithiation and full de-lithiation. It is observed that the crack network obtained after insertion formed almost instantaneously just before full insertion. A detailed inspection of this data shows that cracks initiated at the thin film surface (as expected from the stress profile shown in Fig. 3(b)) and then propagated through the thickness with a small angle with respect to the vertical direction. Conversely, the new cracks created during delithiation propagated vertically. They initiated from the thin film surface after 25% of the full de-lithiation time. These new cracks reached the interface between the thin film and the substrate after 60% of de-lithiation. During the last 40% of de-lithiation, the interface delaminated progressively until full insertion. One can further observe that close to the intersections between the cracks and the Ti/Si interface the vertical displacement was positive. This gave rise to the typical bending shape of the deformed film observed in Fig. 1(III).
Analyzing these results should also give a different way to explain why a lower g c is used for delithiation. In Fig. 3(b), it is shown that during lithiation, the hydrostatic stress is negative meaning that only the contribution of shear is active in the activation of damage as stated by Equation (6). This shear dominated failure mode explains why the cracks initiated during lithiation are not vertical through the thickness, they are oriented along the direction of maximum shear. Conversely, during de-lithiation, the hydrostatic stress is positive and its contribution is thus active in the development of damage. In Fig. 5(c), one can clearly observe the different orientations that cracks follow through the thickness of the film depending if they form during lithiation (shear dominated) or de-lithiation (tension dominated). Many materials are known to have higher failure energy in shear mode than in tensile mode. It would thus be consistent to adopt a g c value that is higher for lithiation (dominated by shear) than for de-lithiation (tension dominated). Another route to account for this effect would be to weight differently the contribution of shear and tension in the definition of the positive strain energy (Equation (6)).

Electrochemical cycling of a 50 nm Si thin film that does not fracture
In continuing, the model is applied to a Si film anode with a thickness of 50 nm in order to test it against experimental observations [4] which illustrated that fracture of amorphous Si was inhibited when the thickness was below 100 nm. The material parameters used (including the internal length value) were the same as for the previous simulation. The boundary conditions were defined to either maintain the same insertion/de-insertion time or the same ion flux thus leading to 10 times shorter insertion/de-insertion times (because the film thickness is 10 times smaller). In both cases, the model predicts that no failure occurs as illustrated in Fig. 6, in agreement with the results of [4]. In the model, capturing the effect of the thin film thickness relies on the internal length l c which governs the non-local nature of damage evolution governed by Equation (5). Because this internal length defines the distance over which damage is spread, its localization is not predicted for these thinner film cases. Quasi homogeneous damage fields (around = d 0.7) are obtained instead.

Lithiation of a unit cell of the patterned Si-film
As a third case study a unit cell of the Si-film patterned with micronholes is considered. The holes were expected to relieve the stresses produced by the colossal volume expansions that occur upon lithiation. As in the case of the thin film a square domain of × 10 10 μm 2 is considered. To account for the actual shape of the hole, a Virtual Image Correlation (VIC) method as first introduced in Ref. [21] and then further developed in Ref. [22], is used to capture the actual hole contour directly from the scanning electron microscope (SEM) image in Fig. 1(a). Then a unit cell with a hole having the detected shape is meshed and the simulation of lithiation is performed. The mechanical and chemical boundary conditions are the same as for the previous example (continuous film) except that now the chemical flux on the top surface is increased to 1.75 nm/ms in order to keep the same lithiation time as for the thin film (the thickness of the Si film being slightly higher than for the patterned film). The hole surface is considered to be stress free and also free from ion flux. However, the mesh is slightly coarser. It is made of four nodes tetrahedral elements with an average edge length of 150 nm. Different values for l c and g c are also adopted to compensate for the increase of the typical mesh size considered in this  simulation hence l c is 300 nm and g c 8000 J/m 2 . In Equation (6), the strain energy density is compared to g l / c c , this set of parameters is thus consistent with the value adopted for the previous case in Section 4.2.
In Fig. 7(a), the deformed shape of the unit cell after the first Liinsertion is shown. As observed experimentally no failure occurs during this first lithiation, and similarly to the experiments a strong reduction of the hole diameter (about 1 μm) is obtained as illustrated in Fig. 7(b). The hole surface is allowed to deform intensively leading to a release of the stresses induced by Li-ion insertion, inhibiting thus failure during the first lithiation. These results are in agreement with the experimental observations in Fig. 1(c) and (d).

Cycling a × 3 3 hole patterned Si film
To further compare with the experimental crack patterns of the patterned Si films, an area containing a × 3 3 hole lattice is simulated during cycling. The imperfections of the etched holes are accounted by rotating the holes by a small random angle (2.5°standard deviation). The boundary conditions and material constants are the same as in the previous example. The mesh size is also identical. After the first lithiation, the results are the same as for one unit cell. The same reduction of the holes' diameter on the top surface and a similar accommodation of the expansion due to the presence of the holes is obtained. During the first de-insertion, a network of diagonal cracks connecting the holes forms, as seen in Fig. 8(a), (b), (c), and is in agreement with the experiments presented in Fig. 1(e) and (f). As for the thin film case, the contraction of the material due to de-lithiation induces failure of the Si-film but also delamination at the Ti/Si interface. After full delithiation, typical bending of the Si film (positive U z displacement around the holes and along the cracks) in Fig. 8(a) is consistent with the observations in Fig. 1(k) and (l).
The simulation is run for ten electrochemical cycles in order to compare with the experimental results, which show the onset of a secondary crack network during the 10th delithiation. To model the progressive degradation of the material properties of the anode after several cycles, the fracture energy g c is decreased by 9% of its initial value after each delithiation. At the tenth cycle, the fracture energy has been decreased by about 80% of its initial value, and new cracks develop. The pattern presented in Fig. 8(d) is in good agreement with Fig. 1(i) and (j) where mostly horizontal and vertical cracks are observed.
In the results presented above, realistic crack patterns are predicted for two different types of anodes. Despites some limitations (e.g. small strain formulation), the model is successful in capturing crack propagation for various Si anode configurations that correspond to experimental data.

Conclusion
In this paper, a multi-physics phase field damage model was developed to simulate fracture of film anodes during the Li-ion insertion/ de-insertion process. The diffusion of Li-ions was accounted for through Fick's law, while the volume expansions induced during lithiation were accounted for through a separate chemical strain that was proportional to the normalized Li-ion concentration. To simulate material degradation, a phase field damage model was adopted, which allowed the fracture energy of Si to vary throughout each electrochemical cycle. This allowed to capture the dry bed-lake fracture that occurs in Si film anodes, and has been an open question for the past 20 years. The model was further validated against a new type of Si film that was patterned with holes to alleviate the volume changes during cycling, and resulted in an ordered crack pattern. Finally, it was shown that maintaining the same material parameters but decreasing the film thickness to 50 nm inhibited the formation of fracture, which was again consistent with experimental data for Si films below 100 nm. This illustrates that the model can capture the scale of the microstructure which is a key parameter in designing efficient anode materials as experiments have shown that the size of the lithiated particles significantly affects their fracture behaviour; with fracture diminishing as the scale reduces [5]. Hence, this hand-in-hand agreement between the proposed model and experiments presents a new way of dealing with fracture in Si anodes, as this multi-physics model can be used to predict the size of the active materials that will limit and control damage, prolonging the electrode life-time.