Sharp-Interface Formation during Lithium Intercalation into
Silicon
Esteban Meca∗Andreas M¨unch†Barbara Wagner∗‡
February 25, 2016
Abstract
In this study we present a phase-field model that describes the process of in-
tercalation of Li ions into a layer of an amorphous solid such as a-Si. The gov-
erning equations couple a viscous Cahn-Hilliard-Reaction model with elasticity in
the framework of the Cahn-Larch´e system. We discuss the parameter settings and
flux conditions at the free boundary that lead to the formation of phase boundaries
having a sharp gradient in ion concentration between the initial state of the solid
layer and the intercalated region. We carry out a matched asymptotic analysis to
derive the corresponding sharp-interface model that also takes into account the dy-
namics of triple points where the sharp interface in the bulk of the layer intersects
the free boundary. We numerically compare the interface motion predicted by the
sharp-interface model with the long-time dynamics of the phase-field model.
Key words. Asymptotic Analysis, Phase-Field Model, Interface Dynamics, Numer-
ical Methods
AMS subject classification 74N20, 35Q74, 35B25, 74S10
1 Introduction
Silicon electrodes for lithium-ion batteries are currently the subject of very intense re-
search to make silicon a practical alternative to graphite. Silicon can store a large
amount of lithium, but the large stresses that these electrodes undergo tend to cause
their fracture and pulverization after a few cycles. Various investigations to circumvent
these problems have led to designs of nanostructured electrodes such as arrays of pillars
or nanowires.
However, in order to carry out systematic and knowledge-based optimizations a fun-
damental understanding of the mechanisms involved in the intercalation process of silicon
∗Weierstrass Institute, Mohrenstraße 39, 10117 Berlin, Germany
†Mathematical Institute, University of Oxford, Andrew Wiles Building, Woodstock Road, Oxford,
OX2 6GG
‡Technische Universit¨at Berlin, Institute of Mathematics, Straße des 17. Juni 136, 10623 Berlin,
Germany
1
Sharp-Interface Formation during Lithium Intercalation into Silicon 2
itself is needed. This has lead to a number of experimental and theoretical investiga-
tions, partly with, seemingly, contradictory results. In particular, the experiments of
Sethuraman et al. [28] on the stress evolution of a silicon electrode during intercalation
have had a great impact on the modeling of silicon electrodes, as they have been taken as
a proof of the surprising result that lithiated amorphous silicon behaves plastically, see
for example the analysis in Bower et al. [2]. Their analysis is based on the assumption
that the film of amorphous silicon is thin enough with respect to the substrate to which
it is attached, so that the stress can be assumed uniform in the whole amorphous silicon
film. However, while the thickness of the film can be controlled, the uniformity of the
stress is an unknown In fact, it has been shown recently in Huang et al. [14] and Ngo et
al.[24]that this assumption uniformity of stress can indeed obscure the interpretation of
the results.
In addition, as has been shown more recently in Levitas et al. [18], the theoretical
yield stress is never reached which has lead to some phenomenological modeling to
explain yielding [33], sometimes also discarding plasticity [18]. On the other hand, the
models used to take into account plasticity have problems in the determination of the
parameters, sometimes suggesting a power law with exponent as high as 50 for the
constitutive law [4], which can be taken as a hint that the model may not be complete.
Further analysis for different electrode geometries involving cylindrical and spherical
silicon particles as well as annular structures have been carried out to investigate stress
evolution and deformation of lithiated silicon, see [9, 6, 7].
On the other hand, regarding the stuctural properties of silicon, it is known that after
the first intercalation cycle, the original crystalline silicon becomes amorphous, see the
review by McDowell et al. [21] and references therein. Moreover, as it has been shown
in McDowell et al.[20] and Wang et al.[30], the first intercalation of crystalline as well as
amorphous silicon occurs through a two-phase mechanism. The question whether this
two-phase process also occurs in amorphous silicon in subsequent intercalation cycles
has been raised in [20, 19] and a recent study by Cubuk et al. [8] based on molecular
dynamics simultations relates the two-phase lithiation of amorphous silicon with a sharp
structural transition in amorphous LixSi for x≈2 between a phase in which the Li atoms
are embedded in a covalent silicon matrix to a phase in which Si atoms are packed in
small clusters with few covalent bonds between them, surrounded by a dense, amorphous
structure of Li atoms, suggesting a transition from a mechanically Si-like phase to a softer
Li-like phase. Moreover, experiments show that this two-phase lithiation is self-limiting
in nanowires [19], presumedbly due to the high stresses generated.
In Meca et al. [22] a phase-field model that couples a viscous Cahn-Hilliard-Reaction
model with elasticity has been derived to describe the process of intercalation of Li-ions
into a layer of an amorphous silicon and investigate conditions leading in the long-time
limit to the formation of sharp phase boundaries between the original silicon phase and
the lithiated phase. The analysis of the long time dynamics of the emerging sharp phase
boundaries is the subject of the present study.
After presenting the phase-field formulation in Section 2, we show in the remainder
of this paper how an intercalation process is nucleated, for example due to some defect
Sharp-Interface Formation during Lithium Intercalation into Silicon 3
on the free absorbing boundary of a layer (or nanowire), and upon reaching a critical
lithium concentration, a phase transition sets in, which in turn leads to a sharp phase
boundary that then moves into the bulk of the layer and along its free boundary.
In Section 3 we derive for this regime a sharp-interface model using matched asymp-
totic expansions, whereby the analysis is divided into three regions, the sharp interface
analysis in the bulk of the layer, the analysis close to the free boundary and the analysis
of the triple points, where the sharp interface in the bulk of the layer intersects the free
boundary. All three regions have to be matched and finally yield the sharp-interface
model. Finally, we investigate in Section 4 numerically the dynamics of the long-time
limit of the phase-field model and compare evolution of the emerging phase boundaries
with expression for the velocities found from matched asymptotic analysis.
2 Formulation of the phase-field model
The system we consider consists of a thin layer of amorphous silicon resting on an
undeformable substrate. Lithium enters in the layer as a consequence of the difference of
electrochemical potential with the electrolyte. As the lithium concentration is increased,
the layer experiences a phase transformation from a poorly lithiated phase (a-Si) to a
heavily lithiated one (LixSi with x>2). Lithium insertion causes a stress-free strain.
For a discussion of the role of phase transformation on stress for this process we have
formulated a mathematical model in Meca et al. [22] and we will briefly introduce the
complete phase-field model here.
To facilitate the discussion of the results we confine our description of elasticity to
linear elasticity and follow the standard approach to the coupling of phase transitions
with linear elasticity, see e.g. [10] for details. From Hooke’s law the elastic energy is
given by
W=1
2Cijkl ij −0
ijkl −0
kl(2.1)
where the fourth order tensor Cijkl is defined as the elasticity tensor or as the stiffness
tensor,ij denotes the strain tensor
ij =1
2(∂jui+∂iuj) (2.2)
defined by the displacement field u, as the difference between the actual position of a
material point and the position in the undeformed material x(reference configuration).
The stress-free strain tensor, that is, the strain due e.g. to composition changes in the
absence of stress, is assumed to grow linearly with the concentration:
0
ij =α(c−¯c)δij (2.3)
where αand ¯care constants, possibly depending on the phase. Note that this choice
implies an isotropic strain change with concentration, but this needs not to be the case
and more general relations are possible.
Sharp-Interface Formation during Lithium Intercalation into Silicon 4
Using the symmetries of Cijkl the choice of elastic energy implies for the stress tensor
σij =Cijkl kl −0
kl,(2.4)
and assuming that the timescale of the elastic relaxation is much faster than that of
diffusion or the phase transformation, elastic equilibrium implies that the divergence of
the stress tensor is zero:
∂jσij = 0.(2.5)
Equation (2.5) will have to be fulfilled separately at each phase in the layer and it can
be written explicitly in terms of the displacement field using
σij =E
1 + νM
ij +ν
1−2νM
kkδij,(2.6)
where Eis Young’s modulus,νis Poisson’s ratio and M
ij =ij −0
ij is the mechanical
strain, the difference between the strain tensor and the stress-free strain tensor.
For the present problem the elastic properties of both phases are assumed different.
We limit this difference to Young’s modulus, since according to Shenoy et al. [29] there
is no clear tendency in the variation of νand obtain the equation:
1
2
E0
E∇c∇u+∇uT−E0
E
1 + ν
1−2να(c−¯c)∇c+E0
E
ν
1−2ν(∇ · u)∇c
+1
2∇2u+1
2(1 −2ν)∇(∇ · u)−α1 + ν
1−2ν∇c= 0
(2.7)
where E0is the derivative of Young’s modulus with respect to concentration c.
For the transport of concentration cwe use the viscous Cahn-Hilliard model
∂tc=∇ · (M(c)∇(µ+χ ε ∂tc)) ,(2.8)
where M(c) is the mobility function, which, in the present study, is taken to be a constant.
The last term is the viscous term, see [25] and χcorresponds to a parameter with
dimensions of viscosity. The chemical potential
µ=1
NΩ
δF
δc =−γε ∇2c+γ
ε
1
2c(1 −c)(1 −2c) + ∂cW(ij, c),(2.9)
is the variational derivative of the free energy
F=NΩZΩ
1
2γε |∇c|2+γ
εf(c) + W(ij, c) (2.10)
where f(c) is the usual double-well potential (free energy per particle) and W(ij, c) is
the elastic energy as defined in Eq. (2.1) for an isotropic elasticity tensor. The constant γ
carries the dimensions of energy times length, and NΩis the (global) number of particles
per surface. The parameter εis related with the interface width and the interfacial
energy between the lithiated and the unlithiated phases.
We nondimensionalize eqs. (2.7) and (2.8) via
Sharp-Interface Formation during Lithium Intercalation into Silicon 5
µ=µ∗γH−1
0, x =x∗H0, y =y∗H0, t =t∗H3
0M−1γ−1, ε =ε∗H0,
denoting with ∗the nondimensionalized variables. The characteristic length scale H0
corresponds to the height of the layer and Mis the constant mobility. For the elastic
variables, we apply the scalings
Cijkl =C∗
ijkl
ESi
2(1 + ν), ui=u∗
iαH0, ij =∗
ijα, σij =σ∗
ij
αESi
2(1 + ν).
After dropping the ∗, the nondimensionalized problem can be written as
∂tc=∇2(µ+εβ ∂tc),(2.11a)
µ=−ε∇2c+1
εf0(c) + δ ∂cW(ij, c),(2.11b)
∂jσij = 0,(2.11c)
σij = 2Gij −0
ij+2ν
1−2νGkk −0
kkδij,(2.11d)
where the constitutive laws for the nondimensional shear modulus G=E(c)/ESi and
stress-free strain 0
ij are specified as
G= 1 + g(c)ELixSi
ESi
−1, 0
ij =h(c),
and the derivative of the nondimensional elastic energy takes the form
∂cW(ij, c) = (1 −ν)G0
1−2ν∂1u2
1+∂2u2
2+1
2G0(∂1u2+∂2u1)2
+2νG0
1−2ν∂1u1∂2u2−2(1 + ν)
1−2ν(h(c)G)0∇ · u
+3(1 + ν)
1−2νh(c)2G0.(2.11e)
Here, h(c) and g(c) are interpolating functions such that g(0) = h(0) = 0 and g(1) =
h(1) = 1. Note that we have slightly generalized our previous choice (2.3) for 0
ij, while
letting ¯c= 0. Also, note that we have not defined a scaling for c, but this should define
0 and 1 to be the two equilibrium concentrations. For the boundaries in contact with
the substrate, we will take a no-flux/no-deformation boundary condition:
u= 0,n· ∇c= 0,n· ∇µ= 0,(2.11f)
where nis the normal vector to the surface. In the case of the boundaries in contact
with the electrolyte, we take a no-traction boundary condition and, following [5], assume
a consistent no-flux condition for c(also known as variational boundary condition), to-
gether with a Butler-Volmer type absorption condition for the chemical potential[32]
σ·n= 0,n· ∇c= 0,n· ∇µ=K(µ) = k1−eQ(µ−µext),(2.11g)
Sharp-Interface Formation during Lithium Intercalation into Silicon 6
xx
cσ
Figure 1: Top: Concentration and σxx at t= 0.027. Bottom: Concentration and σxx at
t= 0.037. We use the constant flux boundary condition and the parameters are: β= 0.1
and δ= 0.1. (nondimensional units, the deformation has been scaled to obtain a better
visualization)
where, µext is scaled like µ. Apart from this Butler-Volmer-type absorbtion condition
we will also consider the constant flux boundary conditions n· ∇µ=k.
The problem thus depends on the following non-dimensional groups:
β=χM
H0
, k =krH2
0
Mγ , δ =H0ESiα2
2(1 + ν)γ, Q =γH−1
0
kBT.(2.12)
with the elastic ratio ELixSi/ESi and ν, together with the parameter µext and . Note
that Qand kplay a similar role in the vicinity of equilibrium, and δis the ratio of elastic
to interfacial energies.
For the numerical simulations we have used ELixSi/ESi = 4/9 and ν= 0.25, in
accordance with the calculations form Shenoy et al. [29]. We also use = 0.005 and
k= 4.0 except where indicated.
3 Sharp-interface asymptotics
For a typical scenario, where small defects on the free interface lead to preferred ab-
sorption sites of the incoming ions, the critical concentration that initiates the phase
transition is reached there first and leads to the local formation of sharp phase bound-
aries and the growth of lithiated regions. In Fig. 1 we show snapshots of such an event
(see Section 4 below and Ref. [22] for more information on the numerics). It shows
the concentration and σxx before and after the formation of the sharp interface in the
presence of a non-uniform flux. In this case, we have taken the flux in a small region
near x= 0 to be approximately two times the value outside of it.
On top of Fig. 1 there is a characteristic concentration pattern just before phase
separation. Near x= 0 the concentration is higher than elsewhere, but this is not so
clearly reflected in the stress field, which is more compressive near the top, where the
Sharp-Interface Formation during Lithium Intercalation into Silicon 7
concentration is the highest. After phase separation takes place and the sharp interface
is created we show near the bottom of Fig. 1 the corresponding cross-sections. Now the
deformation is much more visible, as we have two distinct phases with concentrations
near c= 1 and c= 0. We see how σxx is more negative in the transformed phase,
indicating a strong compression and this compression is highest near the triple junctions.
In this example we have prescribed a constant flux as in [22], and taken lateral periodic
boundary conditions for simplicity.
It will be difficult to analyse systematically these nucleation and growth processes,
where diffusion, elasticity and interfacial effects are present and hard to disentangle using
the present model. However, in the past it has been shown to be instructive to derive
reduced sharp-interface models describing the dynamics asymptotically.
To our knowledge, the first development of this kind for a similar problem (an Allen-
Cahn equation coupled with nonlinear elasticity) was made by [11]. They did a com-
prehensive study and developed a Gibbs-Thomson equation that incorporates the right
Eshelby Traction.
The original reference for the analysis of our problem (i.e. the Larch´e-Cahn system,
understood as a Cahn-Hilliard equation coupled with linear elasticity) is [15]. The
authors compare the results with a boundary integral approach and write explicitly the
correspondence between the sharp and the diffuse interface models. According to their
result, the elastic term in the chemical potential has to have a particular scaling with
the interface width in order to recover instantaneous diffusion away from the interface.
Also, they find out that for that purpose, the interpolating functions of cthat define E
and 0
ij must fulfill some conditions.
Further analysis has been performed since then. For instance, [12] proved some exis-
tence and uniqueness results for the diffuse interface model, followed by an asymptotic
analysis [13], finding the corresponding sharp-interface model. More recently, [1] has
proven rigorously these asymptotic results. However, none of these studies consider
realistic boundary conditions.
3.1 Sharp-interface in the bulk
3.1.1 Outer Expansion
c(x, y, t) = c0+εc1+ε2c2+..., (3.1a)
µ(x, y, t) = ε−1µ−1+µ0+εµ1+..., (3.1b)
ui(x, y, t) = ui,0+εui,1+ε2ui,2+..., (3.1c)
σij(x, y, t) = ε−1σij,−1+σij,0+εσij,1+ε2σij,2+..., (3.1d)
For the Cahn-Hilliard equation we obtain, matching powers of ε
∇2µ−1= 0,(3.2)
∂tc0− ∇2µ0= 0.(3.3)
Sharp-Interface Formation during Lithium Intercalation into Silicon 8
Boundary Layer
0(t)
y=0
Sharp Interface
Substrate
Electroyte
y=1
Inner Region
x
Figure 2: Sketch of the layer. We show in pink the lithiated region and we mark the
different regions that require a separate treatment: The inner region, the boundary layer
and the triple junction point at x0(t).
No additional equations are necessary. The values of the chemical potential are:
µ−1=f0(c0),(3.4)
µ0=f00(c0)c1+δ ∂cW(ui,0, c0).(3.5)
Note the dependence of µ0on the elastic energy. Interpolating functions h(c) and g(c)
can be chosen to avoid these terms. This implies
h0(c±
0) = G0(c±
0) = 0,(3.6)
where c±
0are the values of c0on either side of the interface. We will see below that these
are constant, and correspond to minima of F. Hence, together with the requirement that
the value of Fis at a minimum far from the interface, we obtain that c0is a constant
and µ−1= 0
The elastic equilibrium equation has the same form to all orders:
∇ · σm= 0,(3.7)
where the index mrefers to the order of the expansion. A brief inspection of this
equation for m=−1 (and anticipation of continuity of σ−1across the interface) leads
to the conclusion that σ−1= 0, hence we restrict our attention in the outer solution to
m= 0,1, . . ..
3.1.2 Inner Expansion
For the inner variables near the interface we introduce curvilinear coordinates rand s
via the transformations
x=X(s, t) + rY 0(s, t), y =Y(s, t)−rX0(s, t),(3.8a)
Sharp-Interface Formation during Lithium Intercalation into Silicon 9
with X02+Y02= 1, (primes denote differentiation with respect to s) and scale rthe
coordinate normal to the evolving interface as r=ε ρ. In terms of these inner coordinates
(ρ, s), the operators can be written in the following form
∇2→ε−2∂2
ρ+ε−1K∂ρ−ρK2∂ρ+∂2
s+..., (3.8b)
∂t→ −vn
ε∂η−vt∂s+∂t+..., (3.8c)
where Kis the curvature of the interface and vn=Y0˙
X−X0˙
Yand vt=X0˙
X+Y0˙
Yare
the normal and tangential interface velocities We denote the inner variables with tilde
and introduce asymptotic expansions similarly as for the outer problem
˜c(ρ, s, t) = ˜c0+ε˜c1+2˜c2+..., (3.9a)
˜µ(ρ, s, t) = ε−1˜µ−1+ ˜µ0+ε˜µ1+..., (3.9b)
˜ui(ρ, s, t) = ˜ui,0+ε˜ui,1+ε2˜ui,2+..., (3.9c)
˜σij(ρ, s, t) = ε−1˜σij,−1+ ˜σij,0+ε˜σij,1+ε2˜σij,2+..., (3.9d)
We start with the first order of the elastic equilibrium equation, Eq. (2.11c). From
the formulas of the appendix we obtain, to order ε−2:
(∇ · σ)r=∂ρ˜σrr,−1= 0,(3.10a)
(∇ · σ)s=∂ρ˜σrs,−1= 0,(3.10b)
from which follows that
G(˜c0)∂ρ˜ur,0=Ar,(3.11a)
G(˜c0)∂ρ˜us,0=As,(3.11b)
And hence
˜ur,0=ArZρ
0
1
G(˜c0)dρ0+Br,(3.12a)
˜us,0=AsZρ
0
1
G(˜c0)dρ0+Bs.(3.12b)
We will see below that ˜c0will correspond to a kink-like solution with a bounded range,
and G(˜c0) is also bounded by construction, as it is constructed as a polynomial function
of ˜c0. This implies that both integrals will in general diverge as ρ→ ∞, which means
in turn that it will be impossible to match to the outer solution unless Ar=As= 0.
This will also make sure that the displacement field is continuous. The previous result
implies that
∂ρ˜us,0=∂ρ˜ur,0= 0,(3.13)
and therefore
˜σrr,−1= ˜σrs,−1= ˜σsr,−1= 0.
Sharp-Interface Formation during Lithium Intercalation into Silicon 10
The other component of the stress tensor is also zero to this order (˜σss,−1= 0). Similarly,
the strain tensor to this order is also zero.
In the next order, the elastic equations read as follows:
∂ρ˜σrr,0= 0, ∂ρ˜σrs,0= 0,(3.14)
since the lower order term does not contribute. These equations imply that ˜σrr,0and
˜σrs,0do not depend on ρ. Then, the matching condition
lim
ρ→±∞ ˜σij,0−σ±
ij,0= 0,(3.15)
implies that ˜σrr,0and ˜σrs,0are continuous. This means that, to zero order, the condition
of continuous tractions along the discontinuity is fulfilled. This equation together with
the continuity of the displacements ensures a coherent interface.
Note that equations (3.14), while giving continuity of the tractions σ·nthrough the
interface, they also imply a jump condition on the derivatives of the displacement field.
In terms of the displacement field the equations read as follows:
∂ρ[G(˜c0) ((1 −ν)∂ρ˜ur,1+ν(∂s˜us,0+ ˜ur,0K)−(1 + ν)h(˜c0))] = 0,(3.16a)
∂ρ[G(˜c0) (∂s˜ur,0+∂ρ˜us,1− K˜us,0)] = 0.(3.16b)
These equations can be integrated once and then matched term by term to the outer
solution, here given in terms of curvilinear coordinates
G(c±
0)(1 −ν)∂ru±
r,0+ν(∂sus,0+ur,0K)−(1 + ν)h(c±
0)=K1(s),(3.17a)
G(c±
0)∂sur,0+∂ru±
s,0− Kus,0=K2(s).(3.17b)
By subtracting the negative limit from the positive limit we obtain a closed jump con-
dition for the normal derivatives of the displacement field. Once we have found the con-
tinuity of ui,0and the continuity of tractions for σij,0, we move on to the Cahn-Hilliard
equation. We start by writing down explicitly the form of the chemical potential, by
including also the terms coming from Eq. (A.20) of the Appendix:
˜µ−1=−∂2
ρ˜c0+f0(˜c0),(3.18a)
˜µ0=−K∂ρ˜c0+f00(˜c0)˜c1−∂2
ρ˜c1(3.18b)
+δ
2−2h0(˜c0)δij ˜σij,0+G0(˜c0)
G(˜c0)(˜ij,0−h(˜c0)δij) ˜σij,0,
˜µ1=ρK2∂ρ˜c0−∂2
s˜c0− K∂ρ˜c1+1
2f000(˜c0)˜c2
1+f00(˜c0)c2−∂2
ρc2(3.18c)
+δ
2−2h0(˜c0)δijσij,1+h00(˜c0)˜c1δij ˜σij,0+G0(˜c0)
G(˜c0)˜ij,1−h0(˜c0)˜c1δij˜σij,0
+G0(˜c0)
G(˜c0)(˜ij,0−h(˜c0)δij)σij,1+G00(˜c0)G(˜c0)−G0(˜c0)2
G(˜c0)2˜c1(ij,0−h(˜c0)δij) ˜σij,0,
Sharp-Interface Formation during Lithium Intercalation into Silicon 11
note that for ˜µ0and ˜µ1we only write the non-zero elastic terms. Also notice that the
derivative of the elastic energy might bring terms of order −2to the chemical potential,
but they all cancel. This can be easily checked by expanding Eq. (A.20) in the appendix
and using Eq. (3.13) We now proceed order by order with the Cahn-Hilliard equation:
Order ε−3To this order, the Cahn-Hilliard equation reads
∂2
ρ˜µ−1= 0.(3.19)
This means that from Eq. (3.18a) we have:
∂2
ρ−∂2
ρ˜c0+f0(˜c0)= 0.(3.20)
Following the usual assumptions [27] and making ˜µ−1a constant, we take a traveling-
wave solution of the problem ˜c0(ρ) such that
lim
η→±∞ ˜c0(ρ) = c±
0,(3.21)
i.e. it has two constant limits that match the outer solution, implying that the latter is
constant. If we assume that far from the interface the system is in one of the homogeneous
equilibria, then µ−1= 0 will be the outer solution with the right boundary conditions.
This also implies ˜µ−1= 0.
Order ε−2To this order the Cahn-Hilliard equation reads
K∂ρ˜µ−1−βvn∂3
ρ˜c0+∂2
ρ˜µ0= 0.(3.22)
The first term is zero and the last term can be expanded following Eq. (3.18b):
−βvn∂3
ρ˜c0+∂2
ρ−K∂ρ˜c0+f00(˜c0)˜c1−∂2
ρ˜c1+δ ∂˜c˜
W(˜
,˜c)0= 0,(3.23)
with the obvious shorthand
∂˜c˜
W(˜
,˜c)0=1
2−2h0(˜c0)δij ˜σij,0+G0(˜c0)
G(˜c0)(˜ij,0−h(˜c0)δij) ˜σij,0.(3.24)
Eq. (3.23) can be integrated twice to yield:
−βvn∂ρ˜c0− K∂ρ˜c0+f00(˜c0)˜c1−∂2
ρ˜c1+δ ∂˜c˜
W(˜
,˜c)0=Aρ +B, (3.25)
with Aand Bbeing constants. By matching with the outer solution:
lim
ρ→±∞ (˜c1−c1−ρ∂rc0)=0,(3.26)
and since in curvilinear coordinates, ∂rc0= 0, it follows that ˜c1is bounded if c1is
bounded as ρ→ ±∞, and hence A= 0 since all the other elements on the LHS of (3.25)
are bounded as ρ→ ±∞ (we have used implicitly Eq. (3.6) to show that the elastic
Sharp-Interface Formation during Lithium Intercalation into Silicon 12
terms do not contribute). Alternatively, one can deduce the previous result by matching
˜µ0−β∂ρ˜c0to their outer counterparts. Since ∂ρc±
0=∂rµ±
−1= 0 we obtain the desired
result. To obtain the constant B, equal to ˜µ0−βvn∂ρ˜c0, we notice the following. Eq.
(3.25) has the form
f00(˜c0)˜c1−∂2
ρ˜c1=G.(3.27)
It is easy to prove that ∂ρ˜c0is a solution of the homogeneous problem. Since the operator
is self-adjoint, because of the Fredholm alternative we obtain the following solvability
condition: Z∞
−∞
∂ρ˜c0Gdρ = 0 (3.28)
This means that
(β vn+K)Z+∞
−∞
(∂ρ˜c0)2dρ −δZ+∞
−∞
∂ρ˜c0∂˜c˜
W(˜
,˜c)0dρ +B(c+
0−c−
0)=0.(3.29)
Since ∂ρ˜c0goes to zero as ρ→ ±∞,Bmatches µ0. The first integral can be readily
performed, but the second one requires a more careful analysis. First, we observe that
W+−W−=Z+∞
−∞
∂ρ˜c0∂˜c0˜
W(˜
0,˜c0)dρ +Z+∞
−∞
dρ∂ρ˜ij,0∂˜ij,0˜
W(˜
0,˜c0),(3.30)
which can be obtained from the total derivative of ˜
Wwith respect to ρ(note that
∂˜c˜
W(˜
,˜c)0=∂˜c0˜
W(˜
0,˜c0)). We observe that the first integral is the one that we are
interested in, and the second once can be readily computed. First we notice that
Z+∞
−∞
∂ρ˜ij,0∂˜ij,0˜
W(˜
0, c0)dρ =Z+∞
−∞
∂ρ˜ij,0˜σij,0dρ. (3.31)
The sum runs over all indices except ss, since this component does not depend on ρ, see
Appendix B. By integrating by parts we have that
Z+∞
−∞
∂ρ˜ij,0˜σij,0dρ = [ij,0σij,0]±−Z+∞
−∞
˜ij,0∂ρ˜σij,0dρ, (3.32)
but the last integral is zero because of Eqs. (3.14) (note again that the sum excludes
i=sand j=s). By using the continuity of all σij,0except σss,0and because of the
continuity of ss,0we can write more compactly:
Z+∞
−∞
∂ρ˜c0∂˜c0˜
W(˜
0,˜c0)dρ =W+−W−−σ+
ij,0+
ij,0−−
ij,0,(3.33)
where the sum runs over all indices.
Hence, we obtain the following boundary condition:
µ±
0(c+
0−c−
0) = −(β vn+K)I+δ
2hσ+
ij,0+
ij,0−δijh(c+
0)−σ−
ij,0−
ij,0−δijh(c−
0)i
−δσ+
ij,0+
ij,0−−
ij,0,(3.34)
where Iis the integral in (3.29). This same formula without the kinetic term is found
in [15]. The last term corresponds to the elastic energy required to maintain coherence
at the interface [17, 16].
Sharp-Interface Formation during Lithium Intercalation into Silicon 13
Order ε−1To this order the Cahn-Hilliard equation reads as follows:
−vn∂ρ˜c0=K∂ρ˜µ0−β vn∂3
ρ˜c1+∂2
ρ˜µ1−β ∂2
ρ(vt∂s˜c0−∂t˜c0)−β vnK∂2
ρ˜c0,(3.35)
where we have used explicitly ˜µ−1= 0. We can integrate the previous equation once:
−vn˜c0=K˜µ0−β vn∂2
ρ˜c1+∂ρ˜µ1−β ∂ρ(vt∂s˜c0−∂t˜c0)−β vnK∂ρ˜c0+C, (3.36)
where Cis and integration constant, possibly dependent on s. We can take the limit
and match the previous equation term by term to the outer solution in the usual way.
We obtain, after taking the limit:
−vnc±
0=∂rµ±
0+C, (3.37)
where all the terms match to zero or to constants independent of the side of the inter-
face, hence the redefinition of C. From here it is immediate to obtain the conservation
condition:
c+
0−c−
0vn=−∂rµ+
0−∂rµ−
0.(3.38)
3.2 Asymptotic analysis near the absorption boundary
In the following we will assume that the boundary is at y= 1 and we define an inner
coordinate next to the boundary as
η=y−1
ε,(3.39)
and expand
ˆc(x, η, t) = ˆc0+εˆc1+2ˆc2+..., (3.40a)
ˆµ(x, η, t) = ε−1ˆµ−1+ ˆµ0+εˆµ1+..., (3.40b)
ˆui(x, η, t) = ˆui,0+εˆui,1+ε2ˆui,2+..., (3.40c)
ˆσij(x, η, t) = ˆσij,0+εˆσij,1+ε2ˆσij,2+... (3.40d)
The first two orders for (2.11a) in inner coordinates at the boundary are
∂2
ηˆµ−1= 0,(3.41a)
∂2
ηˆµ0= 0,(3.41b)
Remarkably, the third condition (2.11g) brings in terms to all negative orders when ˆµ−1
is not zero at the boundary. These would be impossible to match unless ˜µ−1is zero at
the boundary. Hence we have this boundary condition, together with the matching to
the outer solution:
lim
η→−∞ ˆµ−1=µ−1= 0.
Together with (3.41a), this implies that
ˆµ−1= 0.(3.42)
Sharp-Interface Formation during Lithium Intercalation into Silicon 14
The boundary condition for (3.41b) then is
∂ηˆµ0= 0,(3.43)
as once ˆµ−1= 0 there are no O(ε−1) terms in the right hand side of (2.11g). It follows
that the solution of Eq. (3.41b) is a constant with respect to η, ˆµ0= ˆµ0(x), Thus ˆµ0(x)
is determined by the outer solution via matching, µy=1 = ˆµ0(x), rather than vice-versa.
Notice that rescaling (2.11d) to inner coordinates introduces terms of order −1.
However, if we include terms to this order in the expansion (3.40d) for ˆσij, they lead
to homogeneous problems for the components of σboth for the PDE and the boundary
value problems, leading to a trivial solution. Hence, the terms of order −1in the
expansion of (2.11d) have to vanish as well, which implies that
∂ηˆu1,0= 0, ∂ηˆu2,0= 0.(3.44)
Order ε−1.To this order, (2.11) becomes
∂2
ηˆµ1+∂2
xˆµ−1= 0,(3.45a)
ˆµ−1=−∂2
ηˆc0+f0(ˆc0),(3.45b)
∂ηˆσiy,0= 0, i =x, y. (3.45c)
where ˆµ−1= 0. The first of these implies that ∂ηˆµ1does not depend on η. The boundary
condition for ˆµ1gives
∂ηˆµ1=k1−eQ(µ0−µext).(3.46)
Using the matching condition
lim
η→−∞ (∂ηˆµ1−∂yµ0|y=1) = 0,(3.47)
leads to
∂yµ0|y=1 =k1−eQ(µ0−µext),(3.48)
that is, ˜µ0inherits the boundary condition from the full problem.
Regarding the concentration, (3.45b) has to be fulfilled. Notice that no contribution
of ∂cWappears in this equation, as the possible O(−1) (and, as matter of fact, O(−2))
contributions vanish due to (3.44). For (3.45b), we have on the one hand that
lim
η→−∞ ˆc0=c±,
and on the other ∂ηˆc0|η=0 = 0, which rules out kink-like solutions. This implies that
ˆc0=c±.(3.49)
Integrating (3.45c), and using the O(1) approximation of the boundary condition for ˆσ
from (2.11g) in inner coordinates, we obtain ˆσiy,0= 0 for i=x, y. Matching then gives
the boundary condition
σiy,0|y=1 = 0, i =x, y, (3.50)
for the outer problem for σ0.
Sharp-Interface Formation during Lithium Intercalation into Silicon 15
x(t)
0
y=1 α
ρ
η
ξ
ξξR /2= =
=
−
−R
1
2
1
R /2
Figure 3: Sketch of coordinate system around a triple point
3.3 Conditions at triple points
We now consider a point where the interface between the two phases meets the boundary
of the Silicon domain with the electrolyte following the approach used by [26]. For the
boundary located at y= 1 so that the layer lies in the region y < 1, the triple point is
assumed to have the coordinates x=x0and y= 1. We thus rescale according to inner
scalings in both Cartesian coordinate directions
ξ=x−x0(t)
ε, η =y−1
ε.(3.51)
Rescaling (2.11) gives
ε2∂tˇc−ε˙x0∂ξˇc=∇2
ξ,η ˇµ, (3.52a)
εˇµ=−∇2
ξ,ηˇc+f0(ˇc),(3.52b)
where ∇2
ξ,η =∂2
ξ+∂2
η. Th leading order problem is:
∇2
ξ,η ˇµ0= 0,(3.53a)
−∇2
ξ,ηˇc0+f0(ˇc0)=0.(3.53b)
The rescaled boundary conditions at η= 1 are
∂ηˇc= 0, ∂ηˇµ=εk 1−eQ(ˇµ−µext),(3.54)
thus, to leading order
∂ηˇc0= 0, ∂ηˇµ0= 0,(3.55)
at η= 1. Notice that the flux from the Butler-Vollmer condition has dropped out and
the problem (3.52), (3.55) is the same problem that has been considered for the triple
point at a solid wall elsewhere in the literature, e.g. [3]. We summarise the key arguments
in the following to recover that the interface between the phases meets the boundary
with the electrolyte orthogonally
Sharp-Interface Formation during Lithium Intercalation into Silicon 16
Multiply (3.53b) with ∂ξˇc0and integrate over the box R≡[−R1/2, R1/2] ×[−R2,0],
to get ZZR
f0(ˇc0)∂ξˇc0dηdρ =ZZR
∂ξˇc0∇2
ξ,ηˇc0dηdρ, (3.56)
or
ZZR
f0(ˇc0)∂ξˇc0−∂ξˇc0∂2
ξˇc0+∂ηˇc0∂ξηˇc0dηdρ =ZZR
∂ξˇc0∂2
ηˇc0+∂ηˇc0∂ξηˇc0dηdρ. (3.57)
We know investigate the left hand side (LHS) and right hand side (RHS) of the last
equation separately,
LHS = Z0
−R2f(ˇc0)−1
2(∂ξˇc0)2+1
2(∂ηˇc0)2R1/2
ξ=−R1/2
dη. (3.58)
Taking R1and R2→ ∞ and using that ∂ξˇc0→0 for R1→ ±∞, we obtain
lim
R1,R2→∞ LHS = Z0
−∞
lim
R1→∞ f(ˇc0) + 1
2(∂ηˇc0)2R1/2
ξ=−R1/2
dη. (3.59)
As ξ→ ±∞, the solution ˇc0has to converge to the boundary layer solution at the
boundary with the electrolyte, which are the same for both limits except for a change in
sign. This, however, does not affect the value of the expression in the square brackets,
so in the limit R1→ ∞ the contributions at ξ=R1/2 and ξ=−R1/2 cancel out, and
lim
R1,R2→∞ LHS = 0.(3.60)
For the right hand side, we have
RHS = ZR1/2
−R1/2h∂ξˇc0∂ηˇc0i0
η=−R2
dξ =ZR1/2
−R1/2
∂ξˇc0∂ηˇc0η=−R2
dξ, (3.61)
where we have used (3.55) at η= 1. We now introduce a rotated coordinate system
that is aligned with the interface between the phases, and denote the angle at which
this interface meets the boundary with the electrolyte by α. The aim is of course to
determine αthrough the matching. The new coordinates are
ρ=−ξsin α−ηcos α
ς=ξcos α−ηsin α,
and this gives
RHS = −Z−R1
2sin α+R2cos α
R1
2sin α+R2cos α
Sdρ (3.62)
with
S= cos α(∂ρˇc0)2−(∂ςˇc0)2+sin2α−cos2α
sin α∂ρˇc0∂ςˇc0.(3.63)
Sharp-Interface Formation during Lithium Intercalation into Silicon 17
Now, ∂ςˇc0→0 as ς→ ∞, and thus
lim
R1,R2→∞ RHS = lim
a→∞
∗
lim
R1,R2→∞ Z−R1
2sin α+R2cos α
R1
2sin α+R2cos α
−S dρ, (3.64)
where lim∗is taken under the condition that |R1cos α+R2sin α|< a,
lim
R1,R2→∞ RHS = −cos αZ∞
−∞
∂ρˇc2
0dρ. (3.65)
Equating the limits (3.60) and (3.65) and using that the integral in (3.65) is positive, we
obtain cos α= 0 i. e. α=π/2, as expected.
3.4 Sharp-interface model
In summary the complete sharp-interface model can be written as follows:
∇2µ0= 0,(3.66a)
∇ · σ0= 0,(3.66b)
together with the constitutive relation for stress:
σij,0= 2G±ij,0−0,±
ij +2ν
1−2νG±kk,0−0,±
kk δij,(3.66c)
where G±=G(c±
0) and 0,±
ij =h(c±
0) are constants. The boundary conditions for the
elasticity equation correspond to continuity for the elastic field and for the tractions:
u+
0=u−
0,(3.66d)
n·σ+
0=n·σ−
0.(3.66e)
For the chemical potential equation we have at the interface away from the boundary:
µ±
0(c+
0−c−
0) = −(β vn+K)I+δ
2hσ+
ij,0+
ij,0−δijh(c+
0)−σ−
ij,0−
ij,0−δijh(c−
0)i
−δ σ+
ij,0+
ij,0−−
ij,0,(3.66f)
c+
0−c−
0vn=−∂rµ+
0−∂rµ−
0,(3.66g)
where I=R+∞
−∞ (∂ρ˜c0)2dρ. For the conditions at the absorbtion boundary we have:
∂yµ0|y=1 =k1−eQ(µ0−µext),(3.66h)
σiy,0|y=1 = 0, i =x, y, (3.66i)
and at the triple points the angle is α=π/2.
Note that in the case of constant flux boundary conditions we would obtain the same
results, except for the exponential in Eq. (3.66h), that would not be present.
Sharp-Interface Formation during Lithium Intercalation into Silicon 18
4 Comparisons to phase-field simulations
We now compare the longtime behaviour of the phase-field model (2.11) with the pre-
dictions of the sharp-interface model (3.66). In particular, we are interested in the
convergence of the sharp-interface solution to the profile of the chemical potential as a
function of the interface thickness ε. For that purpose we have computed numerically
the solution of the phase-field model in the one-dimensional case. We have used the con-
stant flux boundary condition, and the simulations have been solved using a nonlinear
adaptive multigrid algorithm [31]. For details on the simulations see [22].
Fig. 4 summarizes the effect of parameters βand δon the convergence. The solutions
of the sharp interface model to leading order in µare straight lines in one dimension,
meeting at the interface. We show in this figure with dashed lines the straight lines
fitted to µfar enough from the inner region, but not too close to the outer boundaries,
by using a heuristic algorithm. On top, we see that the value of εand the value of δ
have a strong effect on the quality of the fit. On the one hand, for the smaller value
of δwe obtain the expected result, i.e. the inner region becomes narrower with smaller
values of εand the convergence improves, this improvement of the convergence is also
reflected in the value of µat y= 1. On the other hand, for δ= 10.0 we see how, while
there appears to be some convergence, the solution curves are far from being straight
lines. This of course indicates that the sharp interface model is far from valid for such
high values of the parameter δ, and a much smaller εis needed to observe convergence.
On Fig. 4 bottom we see the effect of β. For the higher value of βwe have a good
convergence, and the inner region becomes narrower as εdecreases. Nevertheless, the
effect of βis clear, in that there is a depression near the point where both dashed
lines meet. We believe that this is well represented by the expansion in the inner region.
Indeed, in the inner expansion, the Cahn-Hilliard equation contains at order −2a source
term proportional to β(see Eq. 3.22), and hence this result is to be expected. For β=
0.01, the smallest value, we see how this depression in the inner region has disappeared
and the agreement with the linear fits is extremely good.
In Fig. 5 we quantify how good is the agreement between the sharp interface and the
phase field model by using a sharp interface relation, Eq. (3.66g), and computing how
well it is approximated by the quantities extracted from the phase field. We compute
the gradients of the chemical potential from the fitting lines in Fig. 4 and compute the
velocity of the advancing front also numerically. The ratio of the jump in the derivative
and the velocity should be equal to one. We represent the results for β= 0.01 and
β= 0.1 and the agreement is very good. The parabola that fits the β= 0.01 case also
gives a very good approximation to the β= 0.1 case, except for the highest value of ε,
and has a crossing point with the vertical axis at ε= 0 with a value of 1.00028.
This very good agreement and the quadratic convergence seem to indicate, in accor-
dance with similar results in other systems, e.g. Ref. [23], that the smallest order of the
finite εcorrections to Eq. (3.66g) will be Oε2.
Sharp-Interface Formation during Lithium Intercalation into Silicon 19
δ=1.00
0 0.5 1
y
14
15
16
µ
0 0.5 1
y
1
2
3
µ
0 0.5 1
y
14
15
16
µ
0 0.5 1
y
1
2
3
µ
0 0.5 1
y
14
15
16
µ
0 0.5 1
y
1
2
3
µ
0.005 ε
0.010 0.015
δ=10.0
β=0.01
0 0.5 1
y
0
1
2
µ
0 0.5 1
y
0
1
2
µ
0 0.5 1
y
0
1
2
µ
0 0.5 1
y
0
1
2
µ
0 0.5 1
y
0
1
2
µ
0 0.5 1
y
0
1
2
µ
0.005 ε
0.010 0.015
β=0.10
Figure 4: Convergence to the sharp interface limit. Solid lines correspond to the chemical
potential µon the layer, dashed lines correspond to a linear fit near the interface to
compute the derivatives. Top: Results for two values of δ. With fixed flux boundary
conditions, β= 0.01 (nondimensional units). Bottom: Results for two values of the
kinetic parameter β. With fixed flux boundary conditions, δ= 0.1 (nondimensional
units).
Sharp-Interface Formation during Lithium Intercalation into Silicon 20
0 0.005 0.01 0.015
ǫ
1
1.01
1.02
1.03
(∂yµ+−∂yµ−)/vn
β= 0.01
β= 0.10
Figure 5: Convergence to the sharp interface limit, for two values of the kinetic parameter
β. Symbols correspond to the ratio of the expected velocity to the measured front
velocity. The dashed line is a parabolic fit to the β= 0.01 symbols. With fixed flux
boundary conditions, δ= 0.1. (nondimensional units)
5 Conclusion
Understanding the conditions that cause the deformation and eventual destruction of
silicon electrode particles during intercalation is a prerequisite to optimizing design and
loading conditions. In this regard it is interesting to be able to follow the structure
formation and thus the nonuniform stress evolution of a silicon electrode particle in
detail and on long time scales as a function of the control parameters, such as the ion
flux rate or the ratio of elastic to interfacial energies. In particular, the sharp-interface
model we have derived in this study allows for more analytical insight into the complex
intercalation process and except for extreemly low lithiation rates, the time scale on
which the sharp-interface dynamics takes place is the relevant time scale for comparison
with experimental results on the lithiation process.
Results of the comparison of the long-time dynamics of the phase-field model and
the velocity deduced from the sharp interface model show a good quadratic convergence,
in particular for moderate values of the ratio of elastic to interfacial energies. These are
encouraging results, and provide a means to validate the range of parameters in which
the phase-field model gives physically robust results.
Natural extensions of this work will include the stabilty of the moving lithiation front
and the coarsening mechanisms of the evolving structures. Regarding a more realistic
model we note that effects of charge and electrical potential will have to be included in
our model, as at this stage our aim was first to understand the role of non-homogeneous
elasticity in the stress-lithiation curve, see [22]. In addition, we note that one should
be aware of the fact that the underlying phase-field model is only applicable for regimes
that depart little from equilibrium, and we are actually considering non-equilibrium
phase transitions between what can be at best characterized as metastable states.
Sharp-Interface Formation during Lithium Intercalation into Silicon 21
As of now, comparisons to experiments can only be qualitative since, as we have
mentioned before, we have made the simplified assumption of linear elasticity. However,
it will be straight-forward but rather elaborate to include finite strain effects and this is
planned for our upcoming research.
Acknowlegements
This research was carried out in the framework of Matheon supported by Einstein
Foundation Berlin.
A Curvilinear Coordinates
Change of coordinates:
x=X(s, t) + rY 0(s, t),(A.1)
y=Y(s, t)−rX0(s, t),(A.2)
where the prime denotes differentiation with respect to the arclength parameter sand
X0(s, t)2+Y0(s, t)2= 1.
Natural Basis:
g1=Y0(s)
−X0(s),g2=X0(s) + rY 00(s)
Y0(s)−rX00(s).(A.3)
Physical Basis:
ˆ
r≡er=g1,ˆ
s≡es=g2
kg2k=1
hg2,(A.4)
where h= 1 + rK.
Dual Basis:
g1=g1,g2=g2
h2.(A.5)
A.1 Gradient of a Scalar
∇f=∂kfgk=∂rfˆ
r+1
h∂sfˆ
s(A.6)
A.2 Gradient of a Vector
∇v=∂kvi+vjΓi
jkgi⊗gk(A.7)
Sharp-Interface Formation during Lithium Intercalation into Silicon 22
All of the Γi
jk are zero except for:
Γ2
12 = Γ2
21 =K
h,
Γ1
22 =−Kh,
Γ2
22 =rK0
h.
Therefore, it follows that the components of ∇vin the natural basis are:
(∇v)i
k=
∂rv1∂sv1−v2Kh
∂rv2+v2K
h∂sv2+v1K
h+v2rK0
h
(A.8)
In the physical basis, these elements can be readily computed:
∂rvr1
h(∂svr−vsK)
h∂rvs
h+vsK
h
1
h(∂svs+vrK)
.(A.9)
Note that the physical components are related to the natural components as follows:
vr=v1,vs=hv2.
A.3 Divergence of a Second-Order Tensor Field
∇ · S = ∂iSi
j+ Sl
jΓi
il −Si
lΓl
ijgj(A.10)
(∇ · S)1=∂1S1
1+∂2S2
1+ S1
1
K
h+ S2
1
rK0
h−S2
2
K
h,(A.11)
(∇ · S)2=∂1S1
2+∂2S2
2+hS2
1K.(A.12)
In terms of the physical components, we have:
(∇ · S)r=∂rSr
r+1
h∂sSs
r+ Sr
r
K
h−Ss
s
K
h,(A.13)
(∇ · S)s=h∂r(hSr
s) + h∂sSs
s+hSs
rK.(A.14)
A.4 Tensors from the Theory of Elasticity
From Eq. (A.9) we can compute the strain tensor in these coordinates:
=1
2∇u+∇uT=
∂rur1
21
h∂sur+h∂rus
h
1
21
h∂sur+h∂rus
h 1
h(∂sus+urK)
.
(A.15)
Sharp-Interface Formation during Lithium Intercalation into Silicon 23
Since the new coordinates are orthogonal and are related with the old ones through
a rotation (locally) and the elasticity tensor Cijkl is invariant with respect to rotations,
Eq. (2.11d) is still valid. We can write the elements of the stress tensor explicitly:
σrr =2G
1−2ν((1 −ν)rr +νss −(1 + ν)h(c)) (A.16)
=2G
1−2ν(1 −ν)∂rur+ν
h(∂sus+urK)−(1 + ν)h(c)
σrs =σsr = 2Grs =G1
h∂sur+h∂rus
h (A.17)
σss =2G
1−2ν((1 −ν)ss +νrr −(1 + ν)h(c)) (A.18)
=2G
1−2νν∂rur+1−ν
h(∂sus+urK)−(1 + ν)h(c)
In order to write the equations, it is important to obtain the value of the derivative
of the elastic energy (Eq. (2.1)). It can be written as follows:
∂cW(, c) = 1
2−2h0(c)δijσij +G0(c)
G(c)(ij −h(c)δij)σij(A.19)
where the indices run over r,sand z.
We can compute the derivative explicitly in terms of the displacement field:
∂cW(, c) = (1 + ν)3(h(c)2G(c))0
1−2ν
+(1 −ν)G0(c)
1−2ν(∂rur)2+1
h2(∂sus+urK)2
+1
h
2νG0(c)
1−2ν∂rur(∂sus+urK) (A.20)
−2(1 + ν)(h(c)G(c))0
1−2ν∂rur+1
h(∂sus+urK)
+1
2G0(c)1
h∂sur+h∂rus
h2
B Elasticity on the boundary
In this appendix we discuss the continuity of the different elements of the stress and
strain tensor. We begin by writing explicitly their form for the order zero, in the inner
and in the outer case:
0=
∂ηur,1
1
2(∂sur,0+∂ηus,1− Kus,0)
1
2(∂sur,0+∂ηus,1− Kus,0)∂sus,0+ur,0K
(B.1)
Sharp-Interface Formation during Lithium Intercalation into Silicon 24
σ0=2G(c0)
1−2ν(1 −ν)rr,0+νss,0−(1 + ν)h(c0) (1 −2ν)rs,0
(1 −2ν)rs,0(1 −ν)ss,0+νrr0−(1 + ν)h(c0)
(B.2)
We know from Eqs. (3.11) that the displacement fields (and therefore all their s
derivatives) are continuous at the interface.
˜
±
0=
∂r˜u±
r,0
1
2∂s˜ur,0+∂r˜u±
s,0− K˜us,0
1
2∂s˜ur,0+∂r˜u±
s,0− K˜us,0∂s˜us,0+ ˜ur,0K
.(B.3)
Hence it follows that ˜ss,0is continuous at the interface. Since components rr,rs and sr
of the stress tensor are continuous at the interface we obtain two different jump conditions
for the strain tensor from the components of the stress tensor at the interface:
˜
σ±
0=2G(c±
0)
1−2ν(1 −ν)˜±
rr,0+ν˜ss,0−(1 + ν)h(˜c±
0) (1 −2ν)˜±
rs,0
(1 −2ν)˜±
rs,0(1 −ν)˜ss,0+ν˜±
rr,0−(1 + ν)h(˜c±
0)
(B.4)
From the different jump conditions, we extract the following conditions (no summation):
˜σ+
ij,0−˜σ−
ij,0˜+
ij,0−˜−
ij,0= 0 (B.5)
References
[1] H. Abels and S. Schaubeck. Sharp interface limit for the Cahn-Larche system.
Asymptotic Analysis, 91(3-4):283–340, 2015.
[2] A. F. Bower, P. R. Guduru, and V. A. Sethuraman. A finite strain model of
stress, diffusion, plastic flow, and electrochemical reactions in a lithium-ion half-
cell. Journal of the Mechanics and Physics of Solids, 59(4):804–828, 2011.
[3] L. Bronsard, H. Garcke, and B. Stoth. A multi-phase Mullins-Sekerka system:
matched asymptotic expansions and an implicit time discretisation or the geometric
evolution problem. Proc. Roy. Soc. of Edingburgh: Sec. A Mathematics, 128:481–06,
1998.
[4] G. Bucci, S. P. Nadimpalli, V. A. Sethuraman, A. F. Bower, and P. R. Guduru.
Measurement and modeling of the mechanical and electrochemical response of amor-
phous Si thin film electrodes during cyclic lithiation. Journal of the Mechanics and
Physics of Solids, 62:276–294, 2014.
[5] D. Burch and M. Z. Bazant. Size-dependent spinodal and miscibility gaps for inter-
calation in nanoparticles. Nano Letters, 9(11):3795–3800, 2009. PMID: 19824617.
Sharp-Interface Formation during Lithium Intercalation into Silicon 25
[6] J. Chakraborty, C. P. Please, A. Goriely, and S. J. Chapman. Combining mechanical
and chemical effects in the deformation and failure of a cylindrical electrode particle
in a li-ion battery. International Journal of Solids and Structures, 54:66 – 81, 2015.
[7] J. Chakraborty, C. P. Please, A. Goriely, and S. J. Chapman. Influence of constraints
on axial growth reduction of cylindrical li-ion battery electrode particles. Journal
of Power Sources, 279:746 – 758, 2015. 9th International Conference on Lead-Acid
Batteries {LABAT}2014.
[8] E. D. Cubuk and E. Kaxiras. Theory of structural transformation in lithiated
amorphous silicon. Nano Lett., 14(7):4065–4070, Jul 2014.
[9] Z. Cui, F. Gao, and J. Qu. A finite deformation stress-dependent chemical potential
and its applications to lithium ion batteries. Journal of the Mechanics and Physics
of Solids, 60(7):1280–1295, 2012.
[10] P. Fratzl, O. Penrose, and J. L. Lebowitz. Modeling of phase separation in alloys
with coherent elastic misfit. Journal of Statistical Physics, 95(5-6):1429–1503, June
1999.
[11] E. Fried and M. E. Gurtin. Dynamic solid-solid transitions with phase characterized
by an order parameter. Physica D: Nonlinear Phenomena, 72(4):287–308, 1994.
[12] H. Garcke. On Cahn-Hilliard systems with elasticity. Proceedings of the Royal
Society of Edinburgh: Section A Mathematics, 133(02):307–331, 2003.
[13] H. Garcke and D. J. C. Kwak. On asymptotic limits of Cahn-Hilliard systems with
elastic misfit. In A. Mielke, editor, Analysis, modeling and simulation of multiscale
problems, pages 87–111. Springer, 2006.
[14] Y. Huang, D. Ngo, and A. Rosakis. Non-uniform, axisymmetric misfit strain: in
thin films bonded on plate substrates/substrate systems: the relation between non-
uniform film stresses and system curvatures. Acta Mechanica Sinica, 21(4):362–370,
2005.
[15] P. Leo, J. Lowengrub, and H.-J. Jou. A diffuse interface model for microstructural
evolution in elastically stressed solids. Acta materialia, 46(6):2113–2130, 1998.
[16] P. Leo and R. Sekerka. Overview no. 86: The effect of surface stress on crystal-melt
and crystal-crystal equilibrium. Acta Metallurgica, 37(12):3119 – 3138, 1989.
[17] P. H. Leo and R. Sekerka. The effect of elastic fields on the morphological stability
of a precipitate grown from solid solution. Acta metallurgica, 37(12):3139–3149,
1989.
[18] V. I. Levitas and H. Attariani. Anisotropic compositional expansion in elastoplas-
tic materials and corresponding chemical potential: Large-strain formulation and
application to amorphous lithiated silicon. Journal of the Mechanics and Physics
of Solids, 69:84–111, Sep 2014.
Sharp-Interface Formation during Lithium Intercalation into Silicon 26
[19] X. H. Liu, F. Fan, H. Yang, S. Zhang, J. Y. Huang, and T. Zhu. Self-limiting
lithiation in silicon nanowires. ACS Nano, 7(2):1495–1503, Feb 2013.
[20] M. T. McDowell, S. W. Lee, J. T. Harris, B. A. Korgel, C. Wang, W. D. Nix, and
Y. Cui. In situ tem of two-phase lithiation of amorphous silicon nanospheres. Nano
letters, 13(2):758–764, 2013.
[21] M. T. McDowell, S. W. Lee, W. D. Nix, and Y. Cui. 25th anniversary article:
Understanding the lithiation of silicon and other alloying anodes for lithium-ion
batteries. Advanced Materials, 25(36):4966–4985, 2013.
[22] E. Meca, A. M¨unch, and B. Wagner. Thin-film electrodes for high-capacity lithium-
ion batteries: Influence of phase transformations on stress. Preprint-5-2016, Insti-
tute of Mathematics, Technical University Berlin, 2016.
[23] E. Meca, V. B. Shenoy, and J. Lowengrub. Phase-field modeling of two-dimensional
crystal growth with anisotropic diffusion. Phys. Rev. E, 88:052409, Nov 2013.
[24] D. Ngo, Y. Huang, A. Rosakis, and X. Feng. Spatially non-uniform, isotropic misfit
strain in thin films bonded on plate substrates: The relation between non-uniform
film stresses and system curvatures. Thin Solid Films, 515(4):2220–2229, 2006.
[25] A. Novick-Cohen. On the viscous Cahn-Hilliard equation. Material instabilities in
continuum mechanics, Edinburgh, 1985–1986, 329–342. Oxford Sci. Publ., Oxford
Univ. Press, New York, 1988.
[26] N. Owen, J. Rubinstein, and P. Sternberg. Minimizers and gradient flows for singu-
larly perturbed bi-stable potentials with a Dirichlet dondition. Proc. R. Soc. Lond.
A, 429:505–532, 1990.
[27] R. L. Pego. Front migration in the nonlinear Cahn-Hilliard equation. Proceedings of
the Royal Society of London A: Mathematical, Physical and Engineering Sciences,
422(1863):261–278, 1989.
[28] V. A. Sethuraman, M. J. Chon, M. Shimshak, V. Srinivasan, and P. R. Guduru.
In situ measurements of stress evolution in silicon thin films during electrochemical
lithiation and delithiation. Journal of Power Sources, 195(15):5062–5066, 2010.
[29] V. Shenoy, P. Johari, and Y. Qi. Elastic softening of amorphous and crystalline
Li–Si phases with increasing Li concentration: a first-principles study. Journal of
Power Sources, 195(19):6825–6830, 2010.
[30] J. W. Wang, Y. He, F. Fan, X. H. Liu, S. Xia, Y. Liu, C. T. Harris, H. Li, J. Y.
Huang, S. X. Mao, and et al. Two-phase electrochemical lithiation in amorphous
silicon. Nano Lett., 13(2):709–715, Feb 2013.
[31] S. Wise, J. Kim, and J. Lowengrub. Solving the regularized, strongly anisotropic
Cahn-Hilliard equation by an adaptive nonlinear multigrid method. Journal of
Computational Physics, 226(1):414 – 446, 2007.
Sharp-Interface Formation during Lithium Intercalation into Silicon 27
[32] Y. Zeng and M. Z. Bazant. Phase separation dynamics in isotropic ion-intercalation
particles. SIAM Journal on Applied Mathematics, 74(4):980–1004, 2014.
[33] K. Zhao, G. A. Tritsaris, M. Pharr, W. L. Wang, O. Okeke, Z. Suo, J. J. Vlassak,
and E. Kaxiras. Reactive flow in silicon electrodes assisted by the insertion of
lithium. Nano Lett., 12(8):4397–4403, Aug 2012.
