## Abstract

In this study, we revisit experiments by Sethuraman *et al.* (2010 *J. Power Sources*, **195**, 5062–5066. (doi:10.1016/j.jpowsour.2010.02.013)) on the stress evolution during the lithiation/delithiation cycle of a thin film of amorphous silicon. Based on recent work that show a two-phase process of lithiation of amorphous silicon, we formulate a phase-field model coupled to elasticity in the framework of Larché-Cahn. Using an adaptive nonlinear multigrid algorithm for the finite-volume discretization of this model, our two-dimensional numerical simulations show the formation of a sharp phase boundary between the lithiated and the amorphous silicon that continues to move as a front through the thin layer. We show that our model captures the non-monotone stress loading curve and rate dependence, as observed in recent experiments and connects characteristic features of the curve with the structure formation within the layer. We take advantage of the thin film geometry and study the corresponding one-dimensional model to establish the dependence on the material parameters and obtain a comprehensive picture of the behaviour of the system.

## 1. Introduction

In recent years, interest in lithium-ion batteries has surged. Their high energy density and their slow loss of charge make them ideal for applications ranging from portable electronics to electric cars [1,2]. Much research is being devoted to improving its characteristics, e.g. their capacity or their charging time.

A particularly active area of research is the development of new electrodes. Typically, Li-ion batteries consist of an anode of graphite and a cathode of a lithium compound, but other materials are being proposed. In particular, silicon is heralded as a very promising alternative to graphite. When charging, lithium is stored in the anode, and silicon electrodes can store as much as 10 times more lithium than their graphite counterparts [3].

Nevertheless, the use of silicon as an electrode material presents a number of challenges [3]. To begin with, its volume increases by about 300% when fully lithiated [4,5], and hence the material is subject to enormous stresses that cause the mechanical failure and the destruction of the electrode after a small number of charge–discharge cycles. To overcome this difficulty, structures such as nanopillars [6] or nanowires [7] are being investigated, with promising results [8].

In order to optimize the geometry and nanopatterning of the electrodes, it is essential to understand the material properties of silicon under heavy lithiation and this has been the subject of extensive experimental and theoretical research in recent years. It seems to be well established that after the first lithiation–delithiation cycle, crystalline silicon becomes amorphous [3] and hence the relevant material to study is amorphous silicon. In that respect, it has been proved that the first lithiation of crystalline as well as amorphous silicon occurs initially through a two-phase mechanism [9,10] for potentiostatic lithiation, but it is still controversial whether this two-phase mechanism is also present for amorphous silicon in subsequent lithiations [9,11]. A recent first-principles study [12] relates the two-phase lithiation of amorphous silicon with a structural transition that would occur for a lithium molar fraction *x*≈2. Wang *et al.* [10] show that the growing phase has *x*≈2.5 and once phase transformation is completed it can still be further lithiated to *x*≈3.75. Experiments show that this two-phase lithiation is self-limiting in nanowires [11], presumedly due to the high stresses generated.

Moreover, there is some experimental evidence [13–16], as well as first-principles studies [17] that show a plastic behaviour of amorphous silicon when lithiated, which seems surprising for a material which also fractures in a brittle manner. These unusual characteristics have sparked a minor controversy. The maximum stress measured in the experiments does not reach the predicted yield stress for amorphous silicon [18], and this seems to forbid any plastic behaviour, against the experimental evidence. There have been some attempts at explaining this. On the one hand, there is the idea of reactive flow [19], by which lithiated silicon could flow plastically below the yield point. In accordance with this idea, a phenomenological yield stress function that incorporates the effects of lithium concentration has been proposed [19]. While this theory is able to reproduce approximately the hysteresis loop for stress, it still has problems to accommodate the linear relationship observed between the observed stress and the charging rate [16].

An alternative theory [20,21] involves the dependency of the chemical potential on deviatoric stresses. Using a generalized chemical potential, which changes discontinuously with the sign of the rate of change of the lithium molar fraction, good agreement with the experimental curve has been demonstrated, even though this comparison is obtained through a fit.

Additional problems arise when considering the modelling of plasticity. As of now, only simple viscoplastic flow laws have been considered. These are based on power laws of stress and give exponents ranging from 5 to 50 [14,22], suggesting that a piece may be missing from the model. Also, the description of damage in the electrodes should be incorporated in the models, as it has been recently done by Zuo & Zhao [23,24].

In this article, we revisit the experiments determining the plastic behaviour [13,15,16] and study the early stages of lithiation up to the formation of a lithiated phase in amorphous silicon. Results by Wang *et al.* [10] and also by Cubuk & Kaxiras [12] suggest that strong gradients and even phase separation may be present in the experiments. We argue here that a minimal model that couples phase separation dynamics to linear elasticity reproduces many features found in experiments, even though a complete description of the system would require nonlinear elasticity and a model of the plastic behaviour. We thus discuss qualitatively the value of the average stress that would be found as a function of concentration using our phase-field model that exhibit characteristic properties of the loading curve that have been previously identified with yielding. Our approach also brings charging-rate dependence into the problem, as well as hysteresis, as this is implicit in any phase-separation phenomenon. This suggests that the inhomogeneity of the electrode might be behind many of the observed features.

We model the thin electrode from Sethuraman *et al.* [13] as a thin layer attached to a fixed substrate. Our modelling approach starts from the Cahn–Hilliard reaction equation [25] commonly used for lithium intercalation dynamics, and includes linear elasticity, following the Larché-Cahn [26] prescription. We note that while we have made this choice for simplicity, as it is one of the simplest models that couples consistently phase separation dynamics with elasticity, it will allow only for qualitative comparison with experimental results.

In §2, we introduce the model and the assumptions behind it, as well as the relevant non-dimensional parameters. In §3, we give a detailed study of the dependency of the stress loading curve on the different parameters and in §4, we discuss the results and present our conclusions.

## 2. The model

We are interested in the dynamics of the system described in [13]; a schematic is shown in figure 1. Lithium is absorbed from an electrolyte by a thin layer of amorphous silicon (a-Si) sitting on a non-deformable substrate. When the lithium concentration is high enough, the layer separates into a Li-rich phase (a-Li_{x}Si) and the a-Si phase with no or very low lithium concentration. For the purpose of formulating the model equations, we assume that the layer is bounded in the *x* direction, and we work within the plane strain approximation. The equations can, therefore, be formulated effectively in a two-dimensional rectangular domain *Ω*, although later we will also consider the one-dimensional case. Coordinates *x*≡*x*_{1}, *y*≡*x*_{2} and *z*≡*x*_{3} are introduced as indicated in the schematic and *t* represents time. The indices of the tensors run from 1 to 3 and summation is implied over repeated indices.

We assume that lithium ions are slowly absorbed in an electrode held at constant temperature. Lithium ions become neutral on a thin boundary layer, when entering the electrode from the electrolyte. The insertion of lithium causes an isotropic stress-free strain *c*. We assume that the stress-free strain can be written as *α* is the maximum stress-free strain and *c*_{max} is the lithium concentration in the a-Li_{x}Si phase, and *h* interpolates between *h*(0)=0 and *h*(1)=1. We will typically use a linear function for *h*.

The medium is considered to be amorphous, and hence assumed to be isotropic. We use linear elasticity for simplicity, so that the elastic energy density has the form:
*C*_{ijkl} is the fourth-order elasticity tensor and *ϵ*_{ij} is the strain tensor, defined in terms of the deformation

The previous elastic energy implies the following definition of stress:
*ϵ*_{i3}=0.

In order to model the coupled system, we include the elastic energy (2.1) in a free energy functional of the following form:
*f* has a double well form (a specific choice will be made later) and *W*(*ϵ*_{ij},*c*) is the elastic energy density as defined in equation (2.1). The constant *γ* carries the dimensions of energy times length, *N*_{Ω} is the (global) number of particles in *Ω* and the parameter *ϵ* is related to the interface thickness.

We introduce the chemical potential,
*M*(*c*) can in general be a function of the concentration, but here we will only consider the case where it is a constant, *M*(*c*)≡*M*. The last term is the viscous term [27], and *χ* corresponds to a parameter with dimensions of viscosity. This term, which is not commonly used in Cahn–Hilliard-like models, is associated with non-local equilibrium interfacial kinetics. Gurtin [28] obtained such a term by introducing *ϵ* of each term has been chosen to make sure that we obtain the correct asymptotic sharp-interface limit, see Meca *et al.* [30] for details.

Equations (2.5) and (2.6), together with the mechanical equilibrium condition
*F*) boundary conditions at the electrolyte. The latter corresponds to a constant current flowing through the electrode (i.e. the galvanostatic regime). Boundary conditions are formulated explicitly in non-dimensional form at the end of this section.

### (a) Non-dimensionalization

We adapt the non-dimensionalization of Leo *et al.* [31]. For characteristic length scale *L*_{0}, we choose the thickness of the amorphous silicon layer, concentrations are scaled by *c*_{max}, *G*_{0}=*E*_{Si}/[2(1+*ν*)] is the shear modulus of pure amorphous silicon and is used to scale the shear modulus which is also assumed to depend on the concentration. The constants and *E*_{Si} and *ν* are Young’s modulus for pure amorphous silicon and Poisson’s ratio, respectively. The scalings are
*x* represents all spatial dimensions. With these scalings, all variables become non-dimensional. For convenience, we drop the tildes from the variables from now on. We obtain
*a*
*b*
*c*in the bulk. The constitutive law for the stress takes the form,
*d*where
*e*
*f*The constant *E*_{LixSi} is Young’s modulus for fully lithiated amorphous silicon and *g*(*c*) interpolates between *g*(0)=0 and *g*(1)=1. In this paper, we specifically choose *g* to be linear. We also have
*g*at the electrode-substrate boundary, and so-called variational boundary conditions [32] at the top electrode–electrolyte boundary
*h*with a constant, non-dimensional flux *Y* . In the present form, we neglect the dependency of the surface energy on concentration. At the two lateral boundaries of the layer, we impose the same conditions except that we set the flux to zero.

We have six dimensionless parameters, namely the ratio of the Young moduli *E*_{LixSi}/*E*_{Si}, Poisson’s ratio *ν*, the scaled interface thickness *ϵ* and
*X* is the viscous parameter *χ* non-dimensionalized with *M* and *L*_{0}, and hence it can be identified with the kinetic parameter, i.e. a parameter that quantifies the importance of Li attachment–detachment kinetics at the interface between the lithiated and the non-lithiated phases. This can be clearly seen in the sharp interface limit of the equations [29,30]. Parameter *Y* is simply the flux *F* non-dimensionalized again with *M* and *L*_{0}. It is, therefore, directly related with the total current that flows through the electrode, and it is the only parameter that can be experimentally controlled. Parameter *Z* is the ratio of the elastic energy to the interfacial energy, following the notation of Leo *et al.* [31]. Higher values of this parameter imply a more important role of the elasticity. Below this is discussed in more detail.

## 3. Results

We numerically study the model described by equations (2.9) in one and two dimensions via an adaptive nonlinear multigrid algorithm and the solver BSAM [33]. We use a Crank–Nicolson scheme for the time stepping, and discretize the equations with a finite-volume scheme to enforce the conservation of Li.

Regarding the numerical scheme, using finite volumes in a rectangular grid results in a conservative finite-difference scheme. The resulting scheme is well known, involving the standard 5-point stencil for the Laplacian and a more complicated 9-point stencil for the divergence of the stress. Nevertheless, the traction-free boundary conditions are challenging, as they introduce derivatives along the boundary that couple all the ghost points, and hence at every step of the multigrid algorithm and for every single grid a linear system of equations has to be solved in order to find the correct ghost points that give the right boundary conditions for the displacement vector.

Approximate values for some of the parameters can be obtained from the literature. From Shenoy *et al.* [18], we obtain *ν*=0.25 for Poisson’s ratio. We assume that we can neglect the dependency of *ν* on the lithium concentration. The interface thickness is taken to be 1.25 nm, which is in line with the sharp interface observed in experiments [10]. For a layer thickness *L*_{0}=250 nm, we obtain *ϵ*=0.005. By assuming that the high concentration of Li during phase separation corresponds to a-Li_{2.5}Si, as it has been reported by Wang *et al.* [10], we obtain a value of *α*=0.62 following Bower *et al.* [14]. The remaining parameters are characteristic of this transition that remains mostly unknown, and hence we cannot fix the values of *X*, *Y* and *Z*. For this reason, we study the effect of varying these parameters, see the study of the one-dimensional system below.

Finally, we choose *h*(*c*)=*g*(*c*)=*c* and *f*(*c*)=*c*^{2}(1−*c*)^{2}/4 for simplicity. Below we show that the results will necessarily depend quantitatively on the form of *f*, but we expect to capture the qualitative behaviour. We discuss below the bounds of the non-dimensional numbers based on this form of *f*.

To illustrate the behaviour of the system, we solved the equations for a layer with a moderate aspect ratio. The results are displayed in figure 2.

The stress loading curve (referred to as loading curve from now on) shows the average *σ*_{xx} versus the average concentration on *Ω* as the electrode is loaded. At the beginning, the stress grows linearly with the concentration (i.e. it becomes more negative, compressive), but at a certain concentration, marked as point (*a*) the regime changes abruptly, until it reaches a maximum, marked as point (*b*).

To understand this behaviour, we display as inset subfigures in figure 2 the concentration field and the *σ*_{xx} field at the points (*a*) and (*b*). As we see, at point (*a*) phase separation is beginning, starting from the sides (a two-dimensional effect), and negative stress is localized in the highly lithiated phase. As the concentration is increased from point (*a*) to point (*b*), phase separation is completed and we are left with two definite layers with an abrupt change of concentration between them.

Note that the loading curve has a marked minimum, but unlike in the experiments by Sethuraman *et al.* [13], this minimum does not correspond to the yield point of the material, but rather to simple stress localization: as the system phase-separates, lithium and hence stress are concentrated in a thin layer with a smaller value of Young’s modulus than a-Si. Hence, average stress becomes less compressive, but in fact the stresses are much higher in this thin layer than they were before the onset of phase separation.

This behaviour can be understood further by considering figure 3. In this figure, we display the average concentration and *σ*_{xx} on a cross section at *x*=0. It shows that at an early stage, when the average concentration is somewhat smaller than at (*a*), we have an almost flat curve, with a slightly higher value near *y*=1, and the stress is almost constant. The concentration is almost uniform because diffusion is prevalent in this regime, i.e. the flux is small enough for diffusion to dominate. At point (*a*), the concentration increases much more rapidly near *y*=1 than in other parts of the layer, thus signalling the start of a phase-separation process. The stress becomes very negative near *y*=1, as the concentration is also localized there. At this point the lithiated phase is nucleating, and is absorbing lithium not only from the electrolyte but also from the inside of the a-Si layer. This can be seen in the figure, where this Li depletion is marked by a slight depression in the concentration near *y*=0.9. Finally, this nucleation process culminates at point (*b*), and we can observe a very well-defined layer with high concentration and very large negative stress. Now, the a-Si phase is (almost) fully depleted of lithium, i.e. phase separation is complete.

So far, we have presented the results for a single set of parameters. In the following section, we turn to the one-dimensional case to investigate more easily if these features are generic or depend strongly on the parameters.

### (a) Parameter study for the one-dimensional system

In the experiment of Sethuraman *et al.* [13], the layer of a-Si that acts as an electrode has a thickness of 250 nm and the wafer has a diameter of 3 inches (and the curvature is small enough, so we can consider it locally flat). Hence, the problem is one-dimensional to very good approximation. In fact, this is true even for the results in figure 3 despite the moderate aspect ratio used there. Also, the one-dimensional formulation is more fundamental. It is independent of the specific model for the interaction with the electrolyte, as the flux is imposed externally and has no lateral variations.

The model equations (2.9) can be reduced to their one-dimensional form by dropping all dependencies on *x* and assuming that the lateral displacement *u*_{x}=0. This immediately leads to *σ*_{xy}=0. By using equation (2.9c) and boundary conditions (2.9h), we obtain *σ*_{yy}=0 and
*h*(*c*)=*g*(*c*)=*c* we can express the chemical potential as

The qualitative study of (3.3) shows that the effect of the coupling is to leave the energy of pure a-Si unchanged, and raise the energy associated with the fully lithiated state. This suggests that higher values of *Z* will result in higher values of the concentration *c* at which phase separation occurs. It is easy to compute the stability of a uniform profile of concentration in the limiting case of low flux (see appendix A). The main result of this stability analysis confirms indeed that the coupling with elasticity delays phase separation, and there is a linear relation between the concentration at the onset and *Z* for *ϵ*≪1 (equation (A.5)).

Figure 4 displays several loading curves for different values of *Z*. For small values of the concentration, the curve follows very closely the curve that would be expected for a uniform layer, equation (3.2). Nevertheless, at a certain point the system begins to phase-separate and goes into a different regime. In this regime, phase separation is complete (cf. point (*b*) in figure 2). Then, as most of the stress is located in the lithiated layer and it is approximately constant, the average stress will be this value of stress times the ratio of the thickness of the lithiated layer over the total thickness, i.e. the average concentration. The stress on the lithiated layer can be found from equation (3.2) for *c*=1, and the average stress is this value times the average concentration. This relation is depicted in figure 4 as the upper dashed line.

To recapitulate, the lower dashed line and the upper dashed lines in figure 4 correspond to the two possible lithiation regimes. The lower curve represents single-phase lithiation, and its form is given by equation (3.2) with *c* being substituted by the average concentration, which is valid for *c* approximately uniformly distributed. The upper curve corresponds to two-phase lithiation, in which all the lithium and therefore the stress is concentrated in the upper layer, and the average stress is computed by multiplying that value of stress by the relative thickness of the layer, which in turn is given by the average concentration, hence giving the linear law displayed. Our model shows a transition between both regimes.

For *μ* to remain physically meaningful, the value of *Z* cannot be too high, to ensure that the term coming from the elastic energy remains small compared with *f*. Again, the particular value of *Z* at which the results will not be meaningful will depend on the particular form of *f*. We see in figure 4 that phase separation occurs at higher values of the concentration with increasing *Z*, as expected from equation (A.5). We note that the curve for *Z*=5 clearly behaves differently, meaning that for this value the perturbation caused by elasticity to *μ* cannot be considered small anymore. We also remark that the values of the remaining parameters in figure 4 are not important to illustrate this dependence as it is observed for all values.

Figure 5 shows the variation of the loading curve with non-dimensional flux *Y* . Phase separation occurs more abruptly and at a higher concentration as *Y* is decreased. For reference we show in figure 5 the concentration at which a uniform concentration would be unstable and undergo spinodal decomposition for *Z*=1.0 (appendix A). Clearly, the smaller the flux is the more uniform the concentration in the electrode will be, and therefore, the closer the critical concentration will be to the critical concentration for the uniform system. This is a form of rate dependence of the loading curve, i.e. dependence on the rate of charge.

The parameter *Y* is the only one that in principle is controllable in experiments. For a small enough value of the flux, we will observe the behaviour described in figure 5, as there is an absolute limit of the concentration at which the uniform system is unstable, hence the behaviour described is robust for at least an interval of values of *Y* . Also, if the flux is too high, the system jumps immediately into the phase separated regime. Again, the qualitative outcome does not change if we vary the other parameters. Higher *Z* shifts the value of the transition point towards higher concentrations, but the tendency will be the same.

To study the effect of the kinetic parameter *X*, we also varied its value four orders of magnitude (figure 6). The overall effect of increasing *X* is to delay phase separation. When its value is high enough, the extrema of the loading curve disappear (see the line for *X*=50).

The delay of phase separation for increasing kinetic parameter values is to be expected. On the one hand, the growth rate of the instability of the uniform case is reduced for larger values of *X* (appendix A). On the other hand, the delay of phase separation can be expected on more general grounds. In order to grow in the presence of interface friction, which in this model comes from the inclusion of the viscous term, a driving force is needed to counteract this process. This is provided by the inequality of the chemical potentials of the diffusing species, and this implies that the concentration of the growing phase has to be higher than it would be in local equilibrium conditions [34]. This implies that a higher average concentration is needed for phase separation, which will therefore be delayed.

The correct value of the kinetic parameter is unknown, as it is not directly accessible experimentally. Nevertheless, like *Z*, its value should not be too high for the present model to be valid. We see in figure 6 already that for *X*=50 we obtain an extreme change of behaviour, which indicates again that for this value the viscous term is far from being a small perturbation.

## 4. Discussion

We have found that taking into account the possibility of phase-separation dynamics for a lithiated electrode introduces new phenomena that could help explain some of the puzzling features found in experiments. Some of the consequences of two-phase lithiation of a-Si have indeed been found in experiments [9,10].

Our model is a valuable tool to study the effect of phase separation in our system qualitatively. We see that it can lead to a non-monotone behaviour of the loading curve without requiring plasticity, to a dependency on the rate of lithiation and to a hysteresis cycle.

For illustration purposes, we depict in figure 7 a hysteresis loop, which is present, as expected, in the phase-separating system. The electrode is charged with a positive value of the flux (*Y* =4.00) until it phase separates. Then the flux is reversed, and the stress curve goes parallel to the upper dashed line, which corresponds to a fully phase-separated electrode. Before reaching concentration zero the system becomes nearly homogeneous again, thus closing the loop.

Our one-dimensional parameter study shows that phase separation is almost always present, except perhaps for very high values of *Z* or *X* (figures 4 and 6). These high values would require in any case a more accurate model for *f* and the kinetics.

The value or the dependence of the mobility *M* on *c* could also be of importance. In our model it is constant, but there is some theoretical evidence that the diffusivity varies strongly depending on the phase [35]. In our case, while a different diffusivity would clearly change the value of the timescales, we can assume that lithiation is slow enough to allow for instance for diffusion in the a-Si phase to take place. Even if these timescales were not met in practice, the effect that we are describing would still be observable.

In addition, the large volume changes, estimated to be about 280% [5], are beyond the limit of applicability of linear elasticity and finite-strain effects need to be taken into account, that is the main reason for us to limit our results to early stages of lithiation. Yang *et al.* [36] point out that for a cylindrical geometry small-strain elasticity gives a uniform stress on the lithiated shell that is not realistic, which might be relevant to model the damage processes. In our case, for thin-film geometry, we do not expect such effects to challenge the mechanism of stress localization that we have described. Even if the stress was not uniform in the lithiated phase (thus contradicting an assumption made in most experiments), its average would still be comparatively smaller than before phase separation, due to the very different value of Young’s modulus in the lithiated phase. To summarize, the mechanism that we describe should be understood as a minimal model for the phenomenon of stress localization, necessarily qualitative but also independent of the specific model used to describe elasticity.

Moreover, it is important to point out that we cannot discard plasticity as an important effect in this system. On the contrary, stress localization would indeed produce very large compressive stresses that could reach the yield point of the material. The challenge is to explain the fact that the value of the yield stress [13] is significantly lower than the one expected [20]. This could be due to the fact that the experiments [13,15,16] measure the stress indirectly, i.e. through a procedure involving the Stoney equation [37] that considers the stress to be uniform in the thin layer. Therefore, large gradients and two-phase lithiation may have to be taken into account to understand this problem.

The phase-field model in this paper does not attempt to be quantitatively correct and capture all the physical effects. The quantitative results depend, for example, on the precise form of the free energy *f*. We have, however, taken care so that we remain consistent with the correct sharp interface limit [30], and have focused on results and trends that are robust and do not depend on the form of *f* or similar details. In particular, if phase separation is indeed present in the lithiation/delithiation process, stress localization would be a necessary consequence. Regarding phase separation, which is the cause of the surprising uphill diffusion shown on figure 3 and has been found in other phase-changing electrodes such as LiFePO_{4} [38], we stress that we take it as a working hypothesis. While experiments do not show it explicitly, they show indeed two-phase lithiation [9,10] in a potentiostatic regime and theoretical calculations show an abrupt change in material properties as lithium concentration is increased [12]. Hence, even in the absence of phase separation dynamics, two-phase lithiation and the corresponding stress localization should be taken into account in the interpretation of the experiment of Sethuraman *et al.* [13]. On the other hand, this experiment is performed in a galvanostatic regime, and there is no concrete evidence of two-phase lithiation of a-Si in that regime, although this has been observed directly for c-Si [39]. We hope that new experiments will be performed with a-Si that will help clarify this point.

Finally, we consider phase separation between amorphous phases, such as the ones found in amorphous silicon for high pressure [40] or in bulk metallic glasses [41] The results for our model (2.9) predict that the coupling with elasticity hinders phase separation (see appendix A). This observation is a natural result of how a coherency strain (i.e. a strain due to the deformation of the lattice) is developed in a solid solution for a crystalline material [42]. Nevertheless, it is not clear that this effect of the strain on the phase separation is necessarily present in amorphous systems [43], where the picture is more complicated. This topic will be left to future research.

In summary, we make the case that phase-separation dynamics and two-phase lithiation could play an important role in the interpretation of the discussed experiments. and we hope this will spur further experimental investigations. In order to have a more realistic representation of the experiments, we are presently implementing a model that incorporates nonlinear elasticity, plasticity and a more detailed electrochemical modelling.

## Authors' contributions

E.M. in collaboration with B.W. and A.M. developed the model. E.M. performed the numerical simulations and defined the structure of the manuscript. All authors contributed equally to the revision and the writing of the definitive version submitted.

## Competing interests

We have no competing interests.

## Funding

This research was carried out in the framework of Matheon supported by the Einstein Foundation Berlin.

## Appendix A. Linear stability of the one-dimensional model

In order to help the discussion, we develop in this appendix a linear stability analysis of the one-dimensional problem, which is defined by equation (2.9a) with the chemical potential given by equation (3.3). Here we reproduce the classic result by Cahn [42] for spinodal decomposition of a solid solution, slightly augmented with the effect of kinetics. For a general reference on spinodal decomposition and first-order phase transitions, see Gunton *et al.* [44].

We consider a constant solution of equation (3.3), *c*=*c*_{0}. This solution is in principle attainable in the limit of zero flux. With this solution, equation (3.3) can be easily linearized about it:

For an ansatz of the form

We see already that *X* will decrease the growth rate. The critical concentration does not depend on *X* but will depend on the particular form of the potential. Clearly, for *Z*=0 we obtain the standard result that the instability occurs for *f*′′(*c*_{0})<0, beginning with high wavelengths. Setting *λ*=0 in equation (A.3), using *f*(*c*)=*c*^{2}(1−*c*)^{2}/4, and taking *k*=*π* (lowest lying mode for the Neumann problem in [0,1]) we obtain the following value for the critical concentration:
*a*=(*E*_{LixSi}/*E*_{Si}−1) and *b*=(1+*ν*)/(1−*ν*). The previous equation gives two points on the so-called strain or coherent spinodal [42] for this problem. We can prove easily for small values of *ϵ* that the coupling with elasticity hinders phase separation, as expected. If we take the smallest value in equation (A.4) (corresponding to the minus sign) and expand it in powers of *ϵ* we obtain:
*f*′′(*c*)=0. We can write down the first-order term as a function of the original parameters:

- Received February 8, 2016.
- Accepted September 1, 2016.

- © 2016 The Author(s)

Published by the Royal Society. All rights reserved.