# Spinodal decomposition and domain coarsening in a multi-layer Cahn-Hilliard
model for lithium intercalation in graphite

Antoine Cordoba  Univ. Grenoble Alpes, CEA, LITEN, DEHT, 38054 Grenoble,
France  Laboratoire de Physique de la Matière Condensée, CNRS, Ecole
Polytechnique, Institut Polytechnique de Paris, 91120 Palaiseau, France
Marion Chandesris  Univ. Grenoble Alpes, CEA, LITEN, DEHT, 38054 Grenoble,
France  Mathis Plapp<EMAIL_ADDRESS>Laboratoire de Physique de
la Matière Condensée, CNRS, Ecole Polytechnique, Institut Polytechnique de
Paris, 91120 Palaiseau, France

###### Abstract

During the intercalation of lithium in layered host materials such as
graphite, lithium atoms can move within the plane between two neighboring
graphene sheets, but cannot cross the sheets. Repulsive interactions between
atoms in different layers lead to the existence of ordered phases called
“stages”, with stage $n$ consisting of one filled layer out of $n$, the others
being empty. Such systems can be conveniently described by a multi-layer Cahn-
Hilliard model, which can be seen as a mean-field approximation of a lattice-
gas model with intra- and interlayer interactions between lithium atoms. In
this paper, the dynamics of stage formation after a rapid quench to lower
temperature is analyzed, both by a linear stability analysis and by numerical
simulation of the full equations. In particular, the competition between
stages 2 and 3 is studied in detail. The linear stability analysis predicts
that stage 2 always grows the fastest, even in the composition range where
stage 3 is the stable equilibrium state. This is borne out by the numerical
simulations, which show that stage 3 emerges only during the non-linear
coarsening of stage 2. Some consequences of this finding for the charge-
discharge dynamics of electrodes in batteries are briefly discussed.

##  I Introduction

The advent of portable consumer electronics in the 1970s has boosted the
research on lithium batteries, leading in 1991 to the first commercialization
by Sony and Asahi Kasei of a rechargeable lithium-ion battery (LIB), with
graphite at the negative electrode and LiCoO2 at the positive electrode [1,
2]. Beyond the portable electronic devices market, LIBs are now successfully
used for electromobility and are considered as a technology of interest for
stationary energy storage [3]. Current research investigates a large variety
of LIB active materials and postlithium ion concepts to meet the demand for
better performing and safer storage technologies [4]. Nevertheless, graphite
still remains the dominant negative electrode material in commercial cells
[5]. It has the benefit of being abundant and cheap, to have an excellent
cycling capacity, a low thermodynamic potential and a decent specific capacity
of 372 mAh/g. However, it also has some well-known limitations, such as its
slow kinetics, solid-electrolyte interface (SEI) formation upon cycling, and
lithium plating upon fast-charge [6]. A better understanding of its intrinsic
properties is necessary to overcome these problems.

Graphite is a lamellar material made of two-dimensional graphene layers
stacked upon each other and linked by relatively weak van der Waals forces.
Lithium ions can be reversibly intercalated in between adjacent graphene
layers [7, 8, 9], but the diffusion of lithium across the graphene layer is
highly unlikely [10, 11, 12, 13]. For high enough lithium concentration, a
staging phenomenon occurs: lithium ions form ordered structures with filled
layers regularly separated by empty ones [14]. These structures are referred
to by a stage number, stage $n$ meaning that one out of every $n$ galleries is
filled, with the others being empty. The different stages correspond to
ordered phases of the layered structure, and coexistence between different
stages is possible. During charge and discharge, a sequence of stages occurs,
and the voltage-charge curve exhibits plateaux during the replacement of one
stage by another one.

There is a large consensus with regards to the structures of lithiated
graphite at high lithium concentrations [15, 16]. In contrast, for low lithium
concentrations, i.e, stages with a number higher than 2, there is an ongoing
discussion about the precise compositions and structure of the stages [7, 8,
9]. Moreover, a hysteresis between charge and discharge is also reported both
in the entropy measurement [17] and in the structural evolution between the
stages [18, 19, 20, 21]. Indeed, during lithiation the sequence is pure
graphite, dilute stage 1, stage 3L, stage 2 and stage 1 whereas during
delithiation a stage 2L can be observed between stage 2 and 3L. Here, the
denomination L stands for “liquid-like”, in reference to the absence of
lithium ordering in the plane between the graphene layers. For a good
description of these transitions, a model is needed that captures all the
different stages with their respective symmetries, in particular stage 3.

The representation in which a layer is completely full or empty, known as the
Rüddorf-Hofmann model [14], cannot explain the kinetics of staging transitions
between stages with even and odd numers, for example between stages 2 and 3.
This is because, as lithium migration from one layer to another can occur only
at the edges of the graphite particles or at grain boundaries, the Rüddorf-
Hofmann representation would require an entire gallery to be emptied or filled
for this transition to occur. Five decades ago, Daumas and Herold [22]
proposed the domain model, where the lithium ions occupy the galleries as
islands, locally maintaining the staging. This structure in islands is
schematically illustrated in Fig. 1 with a coexistence of two regions, one in
stage 2 on the left and one in stage 3 on the right. The fact that graphite
tends to phase-separate into these Li-rich and Li-poor regions is a key
element to understand the kinetics of stage transitions and will affect both
the lithium transport inside the graphite particles and the lithium insertion
at its surface [23, 24]. A model able to capture the formation of these Li-
rich and Li-poor domains is therefore necessary to study the dynamics of the
transition between stage 2 and 3.

Figure 1: Illustration of graphite staging structure based on the Daumas-
Hérold representation. Black lines represent graphene planes and grey domains
the lithium islands. Stage 2 is represented on the left and stage 3 on the
right together with the interactions between lithium atoms.

Diffuse-interface models based on the physics of first-order phase transitions
[25, 26, 27] have been widely used to study phase-separating systems [28, 29].
They naturally capture domain formation and interface motion, and work on
diffusive time scales. In the field of LIBs, renewed attention has been given
to these approaches for their ability to predict and simulate phase separation
in two-phase materials such as LFP [30, 31, 32, 33, 34, 35]. To capture
staging in intercalation compounds, a multi-layer diffuse-interface model has
been developed in the 1980’s [36, 37]. In this approach, the concentration of
lithium in each gallery is treated by a separate evolution equation. The model
considers attractive in-plane interactions and repulsive interplanar
interactions between intercalants. A mean-field approximation is then used to
obtain a free-energy model. It contains a square-gradient energy term for each
layer that arises from the attractive in-plane interactions.

This multi-layer diffuse-interface framework is an attractive tool to analyze
multi-phase lamellar LIB active materials such as graphite. Hawrylak and
Subbaswamy [38] simulated a six-layer system and demonstrated the intrinsic
tendency of the system to form Daumas-Herold islands and its ability to
transition between stages 2 and 3. More recently, using a bi-layer model,
Bazant et al. [39, 40, 24] have studied the dynamics of intercalation between
three stable phases, stage 1’ (dilute solid solution), stage 2 and stage 1 and
compared it to the direct optical imaging from Guo et al. [40]. They highlight
the occurrence of ”checkerboard” domains which cannot be captured by solid-
solution models based on Fickian diffusion. Such models can be refined in
various ways to yield a more quantitative description of specific systems. For
instance, the free-energy model of Hawrylak and Subbaswamy [38] leads to the
formation of a stage $3/2$ which has two full layers and one empty and is thus
symmetric to the stage 3 (through an exchange of full and empty layers). Such
a stage is, however, not observed in the experiments on graphite. Chandesris
et al. [41], showed that the asymmetry of the graphite phase diagram can be
captured by introducing a continuous screening effect for the repulsive
interaction between second neighbours. The resulting model predicts the
formation of the staged phases 2 and 3, without the occurrence of the 3/2
stage. Agrawal and Bai [42] used the two-layer model of Bazant et al. [39, 24]
in a multi-scale simulation framework to describe charge and discharge of a
composite electrode containing many graphite particles. Given the
computational cost associated to the multi-layered phase-field approach, some
authors [24, 43, 44, 45, 46] also developed reduced phase-field models to
capture some of the features of the phase transforming behavior of lithiated
graphite while studying the coupling with the other transport phenomena at the
electrode scale and the consequences during fast charge.

The latter references go beyond a simple use of multilayer phase-field models.
Indeed, the electrode of a real battery has a complex architecture, and
different transport processes occur in series during the charge-discharge
process: lithium atoms have to be transported through the electrolyte, then
enter the active material, which can have itself a complex and polycrystalline
microstructure, through an intercalation reaction at the surface, and finally
migrate inside the particles of the active material. Each of these transport
processes occurs on a different scale and exhibits different complexities, so
that multiscale and multiphysics modelling are required for a holistic
modeling of battery operation [43, 47, 48, 49, 44, 45, 46]. However, a more
detailed knowledge on each of the sub-processes is also needed for a precise
description.

Whereas the transition between stages 2 and 1 at high lithium concentration is
certainly more relevant for the optimization of battery electrode operation in
the fast charging regime [46], the competition between stages 2 and 3 at low
lithium concentration remains a poorly understood problem that is interesting
from a fundamental point of view and might impact other aging mechanisms
induced by internal mechanical stress and concentration dependent side
reactions. To shed new light on this question, we deliberately choose to focus
on a particularly simple situation: the spinodal decomposition of a
homogeneous unstable initial state, without external driving, in a one-
dimensional stack of graphite without defects or grain boundaries. The great
advantage of this setting is that it can by analyzed analytically by a linear
stability analysis. Hawrylak and Subbaswamy [38] mention such an analysis, but
do not provide any details. Agrawal and Bai [42] perform a linear stability
analysis in the context of open driven systems, taking into account local
interfacial currents. But since the underlying model has only two-layers, the
growth rate of the stage 3 inside the particle cannot be captured. In
conclusion, no complete treatment of this question is as of yet available.

Here, we lay out a general formulation for the stability analysis of multi-
layer staging models without external driving, and apply it to the model of
Chandesris et al. [41]. The analysis is carried out in terms of normal modes,
so that any stage that is compatible with the symmetries of the system can
appear. Therefore, competition between stages 2 and 3 can be studied in
systems with a number of layers that is a common multiplier of 2 and 3; we
choose the smallest such number, which is six. We confirm a statement already
made by Hawrylak and Subbaswamy [38], namely, that stage 2 grows faster than
stage 3, even for initial compositions for which stage 3 is the
thermodynamically stable state. We conduct a systematic parametric study to
ensure that this result is generic to all models that exhibit the
experimentally observed staging sequence. Numerical simulations are then
carried out to verify this prediction, and to study how the system approaches
the correct equilibrium by nonlinear domain-coarsening processes.

In the following, we will first introduce the model in Sec. II, and then
present the linear stability analysis in Sec. III. In Sec. IV, numerical
simulations are performed to test the prediction on the LSA and to study the
transition from stage 2 to stage 3 during domain coarsening. Some possible
consequences of these findings for behavior during charge and discharge are
discusses in Sec. V, followed by the general conclusion.

##  II Model

We consider a system of $N_{z}$ graphene layers stacked along the $z$-axis
(See Fig. 1). Since lithium atoms cannot jump between layers, a separate
concentration field is defined in each layer. These fields are functions of
the continuous coordinates in the plane $x$ and $y$ [38, 24].

###  II.1 Free-energy model

Let $c_{j}(x,y)$ be the local lithium concentration (in units of mole per unit
volume) in layer number $j$. We define the scaled concentration
$\tilde{c}_{j}(x,y)=c_{j}(x,y)/c_{\rm max}$, where $c_{\rm max}$ is the
maximum possible concentration in the layer for the studied material. The
latter can be related to the number of possible intercalation sites per unit
volume, $N_{V}$, by $N_{V}=\mathcal{N}_{A}c_{\rm max}$, where
$\mathcal{N}_{A}$ is the Avogadro number. The total Gibbs free energy of the
system, $G$, may be written in terms of local volumetric free energy density
$\mathcal{G}$, as:

$\displaystyle G$ $\displaystyle=\int_{V}\mathcal{G}~{}dV$
$\displaystyle=\sum_{j=1}^{N_{z}}L_{j}\int_{S}\bigg{[}\dfrac{1}{2}\kappa(\nabla\tilde{c}_{j})^{2}+g(\tilde{c}_{j})$
$\displaystyle+g^{(1)}(\tilde{c}_{j},\tilde{c}_{j\pm
1})+g^{(2)}(\tilde{c}_{j},\tilde{c}_{j\pm 1},\tilde{c}_{j\pm
2})\bigg{]}dx\,dy$ (1)

Here, $V$ is the volume of the material, and $S$ and $L_{j}$ are the surface
area and the thickness of one graphite layer. In this work, we will neglect
the volume dilatation due to intercalation and take $L_{j}=L_{z}$ as a
constant independent of $j$. In the first term of Eq. (II.1), the gradient
operator is applied in the plane (along directions $x$ and $y$), and $\kappa$
is a constant coefficient of dimension energy per unit length. This is the
gradient-square term introduced by Cahn and Hilliard [25]. It describes the
energy penalty linked to any inhomogeneity in concentration and arises from
attractive in-plane interactions between the lithium atoms.

The second term in Eq. (II.1) describes the thermodynamics of a single layer.
The simplest model that leads to phase separation within a layer is the
regular solution model that is a mean-field approximation of the lattice gas
model [29]:

$\displaystyle g(\tilde{c}_{j})=$ $\displaystyle
N_{V}\bigg{\\{}k_{B}T\big{[}\tilde{c}_{j}\ln(\tilde{c}_{j})+(1-\tilde{c}_{j})\ln(1-\tilde{c}_{j})\big{]}$
$\displaystyle+\Omega_{a}\tilde{c}_{j}(1-\tilde{c}_{j})+\tilde{c}_{j}\mu_{\rm
ref}\bigg{\\}}$ (2)

where $k_{B}$ is the Boltzmann constant, and $T$ is the temperature. The first
term in Eq. (II.1) is due to the configurational entropy of full and empty
intercalation sites, while the second term results from the attractive
interaction between lithium atoms within the same plane, with $\Omega_{a}$
being the corresponding interaction parameter. Finally, $\mu_{\rm ref}$ is a
reference potential.

The last two terms of Eq. (II.1) arise from the interactions between lithium
atoms in different planes. The first-neighbor interaction energy $g^{(1)}$ may
be written as:

$\displaystyle g^{(1)}(\tilde{c}_{j},\tilde{c}_{j\pm
1})=N_{V}\dfrac{\Omega_{b}}{2}(\tilde{c}_{j}\tilde{c}_{j+1}+\tilde{c}_{j}\tilde{c}_{j-1})$
(3)

where $\Omega_{b}$ describes the interaction strength. Since this parameter is
positive, it will cost energy to have high lithium concentration in adjacent
layers for a given position $(x,y)$, which will favor the appearance of stage
2.

Interlayer interactions of longer range are necessary to obtain staging
phenomena beyond stage 2. A simple repulsive second-neighbor interaction of
the form $g^{(2)}=\Omega_{c}\tilde{c}_{j}\tilde{c}_{j+2}$, as in [36, 38]
would favor stage 3. However, it was found that this expression leads to a
phase diagram that is symmetric with respect to half-filling
($\tilde{c}_{j}=1/2$): a new stage 3/2 appears, in which two full layers are
separated by a single empty one [37, 41]. Since this stage is not observed in
experiments, we therefore use a more complex expression for the second-
neighbor interactions which includes a continuous screening effect introduced
in a previous work [41]:

$\displaystyle g^{(2)}(\tilde{c}_{j},\tilde{c}_{j\pm 1},\tilde{c}_{j\pm 2})=$
$\displaystyle
N_{V}\dfrac{\Omega_{c}}{2}\bigg{[}\tilde{c}_{j}(1-\tilde{c}_{j+1})\tilde{c}_{j+2}$
$\displaystyle+\tilde{c}_{j}(1-\tilde{c}_{j-1})\tilde{c}_{j-2}\bigg{]}$ (4)

With this expression, high concentrations of intercalants in the next-nearest
neighbor layers are energetically penalized, except if intercalants are also
present in between. The effect of this screening can be assessed considering
two limits:

  * •

$\Omega_{c}(1-\tilde{c}_{j+1})\approx\Omega_{c}$ for low value of
$\tilde{c}_{j+1}$, favoring the formation of stage 3 at low lithium
concentration;

  * •

$\Omega_{c}(1-\tilde{c}_{j+1})\approx 0$ for high value of $\tilde{c}_{j+1}$,
favoring the formation of stage 2 (instead of stage 3/2) at high lithium
concentration.

Interactions of longer range are mentioned in several papers [36, 37, 38]. As
these interactions rapidly decrease with the distance between layers and are
not crucial to analyze transitions between stage 2 and stage 3, they are
neglected in the present study.

###  II.2 Kinetics

The kinetic equations for the concentration fields are formulated using the
thermodynamic theory of irreversible processes [25, 50, 51, 52]. Since the
atoms cannot jump from one layer to another, each concentration field
satisfies a conservation law,

$\displaystyle\dfrac{\partial\tilde{c}_{j}}{\partial t}=-\nabla\cdot J_{j},$
(5)

where $J_{j}$ is the flux of lithium atoms in layer $j$. This flux can be
expressed as

$\displaystyle J_{j}=-\mathcal{M}\tilde{c}_{j}\nabla\mu_{j},$ (6)

where $\mu_{j}$ is the chemical potential in layer $j$. The latter is defined
as the variational derivative of the free energy with respect to the
concentration within this layer:

$\displaystyle\mu_{j}=\dfrac{1}{c_{max}}\dfrac{\delta G}{\delta\tilde{c}_{j}}$
(7)

In Eq. (6), $\mathcal{M}$ is the mobility. Whereas a constant mobility is used
in some papers [38, 41], we choose here an expression that is consistent with
diffusion of atoms on a lattice [24],

$\displaystyle\mathcal{M}(\tilde{c}_{j})=\dfrac{D}{\mathcal{N}_{A}k_{B}T}(1-\tilde{c}_{j}),$
(8)

with $D$ being the tracer diffusivity in the dilute limit. Gathering Eq. (5)
to (8), we have:

$\displaystyle\dfrac{\partial\tilde{c}_{j}}{\partial
t}=\nabla\cdot\bigg{(}M(\tilde{c}_{j})\nabla\dfrac{\delta
G}{\delta\tilde{c}_{j}}\bigg{)}$ (9)

where we have defined

$\displaystyle
M(\tilde{c}_{j})=\dfrac{D}{N_{V}k_{B}T}\tilde{c}_{j}(1-\tilde{c}_{j}).$ (10)

The factor $\tilde{c}_{j}(1-\tilde{c}_{j})$ ensures that there is no mass
transport when there is no lithium, and when the layer is completely full.

##  III Theoretical analysis

###  III.1 Phase diagram

For given interaction parameters, the phase diagram can be determined by
minimization of the free-energy functional under the constraint of mass
conservation. The details of this procedure are described in Ref. [41];
therefore, we will summarize here only the most important points.

For a single layer, the intra-layer free energy $g(\tilde{c}_{j})$ is a
double-well function due to the competition between configurational entropy
and the attractive in-plane interactions of strength $\Omega_{a}$. This leads
to the coexistence of concentrated and dilute domains in the same layer, with
equilibrium compositions given by the common tangent to $g(\tilde{c}_{j})$
[24, 41]. For multiple layers, the formation of the ordered stages is driven
by the repulsive interlayer interactions (illustrated by $\Omega_{b}$ and
$\Omega_{c}$ in Fig. 1). A given stage at equilibrium is characterized by the
co-existence of lithium-poor and lithium-rich domains. Since all the layers
are in contact with the same external medium (the electrolyte), the chemical
potential must be the same in all layers. Therefore, the equilibrium
compositions of a given stage can again be determined by a double-tangent
construction. This also yields the equilibrium potential, which can be
different for each stage. The equilibrium compositions also differ from the
ones in a single layer.

The equilibrium between two different stages occurs when their chemical
potential and their grand potential, defined as $\mathcal{G}(c)-\mu^{eq}c$,
are the same. The coexisting stages have different average compositions
(filling fractions). For filling fractions that fall in between these two
values, the two stages coexist, with volume fractions that are determined by
the lever rule.

In the following, we will use a parameter set that was developed in Ref. [41],
see Table 1. The corresponding phase diagram is presented in Fig. 2. With the
screened second neighbor interaction of Eq. (II.1), the phase diagram is not
symmetric. For mean filling fractions below 0.5, there are two regions of
phase coexistence: one for dilute stage 1/stage 3 and one for stage 3/stage 2.
In contrast, for mean filling fractions above 0.5, the only phase coexistence
is between stage 2 and stage 1. At room temperature, $298~{}K$, the sequence
of stages that occurs for increasing filling fraction is: dilute stage
$1^{\prime}$, coexistence of stages $1^{\prime}$ and $3$, pure stage $3$,
coexistence of stages $3$ and $2$, pure stage $2$, coexistence of stages $2$
and $1$, and finally pure stage $1$. When only one stage is present, the
equilibrium chemical potential is an increasing function of the filling
fraction; in contrast, during coexistence of two stages, the potential remains
fixed to its equilibrium value at phase coexistence, and only the phase
fractions change. Therefore, the curve of potential versus filling fraction
exhibits a plateau whenever two stages coexist. Note that the quantity that is
usually measured is the electric potential, which is proportional to the
negative of the chemical potential. In the model, the value of $\Omega_{a}$
determines the concentration difference between dense and dilute domains and
the width of the voltage plateaus, whereas the interlayer interaction energies
$\Omega_{b}$ and $\Omega_{c}$ will influence the potential differences between
the different plateaus and the limits of phase coexistence. The orders of
magnitude of all the interaction parameters are similar to previously reported
values [30, 39, 24].

The maximum concentration of lithium $c_{\rm max}$ can be evaluated
theoretically from the crystal structure or estimated experimentally from the
reversible capacity. The value used in this work, $c_{\rm max}=30\,000$
mol/m3, falls within the range of averaged values reported in the liiondb
database for graphite [53]. The value of the reference potential $\mu_{\rm
ref}$ was set to zero as it has no impact on the phase diagram. The values of
the gradient penalty coefficient and of the diffusivity will be discussed
later.

$c_{\rm max}$ | $30~{}000\quad\rm{mol}\;\rm{m}^{-3}$
---|---
$\Omega_{a}$ | $64.3\quad\rm{meV}$
$\Omega_{b}$ | $23.1\quad\rm{meV}$
$\Omega_{c}$ | $4.1\quad\rm{meV}$
$\mu_{ref}$ | $0\quad\rm{J}$
$\kappa$ | $3\times 10^{-6}\quad\rm{J}\;\rm{m}^{-1}$
$D$ | $1.25\times 10^{-12}\quad\rm{m}^{2}\;\rm{s}^{-1}$

Table 1: List of model parameters. Figure 2: Phase diagram of lithium
intercalated graphite with the thermodynamic parameters of TABLE 1.

Compared to experimental phase diagrams close to room-temperature [54, 55, 7],
the predicted diagram does not capture all the details at low average filling
fraction and overestimates the mean filling fraction of stage 3. This might be
improved by including longer-range lithium interactions [37] and/or more
complex interactions with the graphite structure [56]. These more complex
approaches, however, still need to demonstrate their ability to fully describe
the transitions occurring at low filling fraction.

###  III.2 Spinodal decomposition

Spinodal decomposition is a spontaneous transition from a single homogeneous
but unstable phase into two phases with lower free energy. It can be induced
by a quench, that is, the system temperature is rapidly changed. Even though
such experiments have, to our knowledge, not been reported for intercalation
materials, it is still instructive to analyze this situation. Indeed, the
initial stages of spinodal decomposition can be understood analytically with
the help of a linear stability analysis.

To conduct such an analysis for our model, we consider a system with periodic
boundary conditions in the $z$-direction. Furthermore, in the following we
simplify the calculations by supposing that the concentration is invariant
along the $y$ direction; the layers are thus one-dimensional. This represents
a cut through a graphite particle in a direction that is normal to the
particle-electrolyte interface and to the graphene planes. The results of the
stability analysis would be unchanged for a full two-dimensional treatment
because of the rotational invariance within the planes. A particle of linear
size $L$ is considered, with no-flux boundary conditions at both sides.

To make notations more compact, the expressions of the interaction energies
are rewritten as follows:

$\displaystyle g^{(1)}(\tilde{c}_{j},\tilde{c}_{j\pm
1})=\dfrac{1}{2}\bigg{[}f^{(1)}(\tilde{c}_{j},\tilde{c}_{j+1})+f^{(1)}(\tilde{c}_{j},\tilde{c}_{j-1})\bigg{]}$
(11)

with
$f^{(1)}(\tilde{c}_{j},\tilde{c}_{j+1})=N_{V}\Omega_{b}\tilde{c}_{j}\tilde{c}_{j+1}$.
Similarly, for the next-nearest neighbor interaction term, we write:

$\displaystyle g^{(2)}(\tilde{c}_{j},\tilde{c}_{j\pm 1},\tilde{c}_{j\pm
2})=\dfrac{1}{2}$
$\displaystyle\bigg{[}f^{(2)}(\tilde{c}_{j},\tilde{c}_{j+1},\tilde{c}_{j+2})$
$\displaystyle+f^{(2)}(\tilde{c}_{j},\tilde{c}_{j-1},\tilde{c}_{j-2})\bigg{]}$
(12)

with
$f^{(2)}(\tilde{c}_{j},\tilde{c}_{j+1},\tilde{c}_{j+2})=N_{V}\Omega_{c}\tilde{c}_{j}(1-\tilde{c}_{j+1})\tilde{c}_{j+1}$.
Furthermore, in the following, $\partial_{\ell mn}f$ corresponds to the
derivation of $f$ $\ell$ times with respect to its first argument, $m$ times
with respect to its second argument and $n$ times with respect to the third
one, and $\partial_{t}$ is the regular time derivative.

With these notations, Eq. (9) combined with the free energy model of Eqs.
(II.1) to (II.1) reads:

$\displaystyle\partial_{t}$ $\displaystyle\tilde{c}_{j}=\nabla\cdot
M\nabla\bigg{[}-\kappa\nabla^{2}\tilde{c}_{j}+g^{\prime}(\tilde{c}_{j})+\partial_{10}f^{(1)}(\tilde{c}_{j},\tilde{c}_{j+1})$
$\displaystyle+\partial_{10}f^{(1)}(\tilde{c}_{j},\tilde{c}_{j-1})+\partial_{100}f^{(2)}(\tilde{c}_{j},\tilde{c}_{j+1},\tilde{c}_{j+2})$
$\displaystyle+\partial_{100}f^{(2)}(\tilde{c}_{j},\tilde{c}_{j-1},\tilde{c}_{j-2})+\partial_{010}f^{(2)}(\tilde{c}_{j-1},\tilde{c}_{j},\tilde{c}_{j+1})\bigg{]}$
(13)

To linearize this expression, we consider that all filling fractions consist
of fluctuations around an average initial composition $\bar{c}$ that is the
same for all layers:

$\displaystyle\tilde{c}_{j}(x,t)=\overline{c}+\delta\tilde{c}_{j}(x,t)$ (14)

Using Taylor series and neglecting all term that are non-linear in the
perturbations ($\delta\tilde{c}_{j}\delta\tilde{c}_{i}$ etc.), the following
system is obtained:

$\displaystyle\partial_{t}$ $\displaystyle\delta\tilde{c}_{j}=\nabla\cdot
M\nabla\bigg{[}-\kappa\nabla^{2}\delta\tilde{c}_{j}+g^{\prime\prime}\delta\tilde{c}_{j}$
$\displaystyle+\partial_{11}f^{(1)}\delta\tilde{c}_{j+1}+2\partial_{20}f^{(1)}\delta\tilde{c}_{j}+\partial_{11}f^{(1)}\delta\tilde{c}_{j-1}$
$\displaystyle+\partial_{101}f^{(2)}\delta\tilde{c}_{j+2}+2\partial_{110}f^{(2)}\delta\tilde{c}_{j+1}+2\partial_{200}f^{(2)}\delta\tilde{c}_{j}$
$\displaystyle+\partial_{101}f^{(2)}\delta\tilde{c}_{j-2}+2\partial_{110}f^{(2)}\delta\tilde{c}_{j-1}+2\partial_{020}f^{(2)}\delta\tilde{c}_{j}\bigg{]}$
(15)

where only the zeroth order term is kept in the mobility, and all the
derivatives with respect to $\tilde{c}$ are evaluated at the unperturbed
concentration $\bar{c}$. This equation can be written in matrix form:

$\partial_{t}\delta\tilde{c}=H\delta\tilde{c}$ (16)

where $\delta\tilde{c}=\begin{pmatrix}\delta\tilde{c}_{1}\\\ \vdots\\\
\delta\tilde{c}_{N_{z}}\end{pmatrix}$ and $H$ is the symmetric circulant
matrix:

$H=\begin{pmatrix}A&B&C&0&\cdots&0&C&B\\\ B&&\ddots&\ddots&\ddots&&0&C\\\
C&\ddots&\ddots&&C&0&&0\\\ 0&\ddots&&A&B&&\ddots&\vdots\\\
\vdots&\ddots&C&B&&&\ddots&0\\\ 0&&0&&&\ddots&\ddots&C\\\
C&0&&\ddots&\ddots&\ddots&&B\\\ B&C&0&\cdots&0&C&B&A\end{pmatrix}$

where

$\displaystyle A=$ $\displaystyle\nabla\cdot
M\nabla\left(-\kappa\nabla^{2}+g^{\prime\prime}+2\partial_{20}f^{(1)}\right.$
$\displaystyle\left.+2\partial_{200}f^{(2)}+\partial_{020}f^{(2)}\right)$ (17)
$\displaystyle B=\nabla\cdot
M\nabla\left(\partial_{11}f^{(1)}+2\partial_{110}f^{(2)}\right)$ (18)
$\displaystyle C=\nabla\cdot M\nabla\partial_{101}f^{(2)}.$ (19)

The circulant structure is due to the invariance of the system with respect to
discrete translation along $z$ by one layer.

The system (16) is a partial differential equation with a time derivative on
the left-hand side and space derivatives on the right-hand side. Separation of
variables then yields that the solution is exponential in time. Furthermore,
the translation invariance in space implies that the spatial solution can be
expanded in terms of plane waves. Therefore, the general solution is of the
form

$\delta\tilde{c}(x,t)=\sum_{m=0}^{N_{z}-1}\sum_{n=-\infty}^{+\infty}\frac{1}{2}\left(a_{mn}e^{ik_{n}x}e^{\omega_{m}(k_{n})t}v_{m}\;+c.c.\right)$
(20)

where $v_{m}$ are the eigenvectors of the matrix $H$, $k_{n}=2\pi n/L$ is the
wavenumber along $x$ (the discrete values are due to the finite system size),
$\omega_{m}(k_{n})$ is the growth rate associated to the eigenvector $v_{m}$
and wavenumber $k_{n}$, and $a_{mn}$ are the amplitudes of each mode. In the
following, we will assume that the system size $L$ is large enough for the
$k_{n}$ to form a quasi-continuum, and drop the index $n$. With these
notations, the growth rates are directly related to the eigenvalues
$\lambda_{m}$ of the matrix $H$:

$\displaystyle\omega_{m}(k)=\lambda_{m}(k)$ (21)

The number of different eigenvectors equals the number of layers. The
eigenvectors of a circulant matrix are the discrete plane waves along $z$:

$\displaystyle v_{m}=\dfrac{1}{\sqrt{N_{z}}}\begin{pmatrix}1\\\
\exp\left(\dfrac{2i\pi m}{N_{z}}\right)\\\ \vdots\\\ \exp\left(\dfrac{2i\pi
m(N_{z}-1)}{N_{z}}\right)\end{pmatrix}.$ (22)

Introducing the solution (20) in the system (16), we can express the
coefficients $A$, $B$ and $C$ of the matrix as function of the wavenumber $k$:

$\displaystyle A(k)$ $\displaystyle=-Mk^{2}\left(\kappa
k^{2}+g^{\prime\prime}+2\partial_{20}f^{(1)}\right.$
$\displaystyle\left.+2\partial_{200}f^{(2)}+\partial_{020}f^{(2)}\right)$ (23)
$\displaystyle
B(k)=-Mk^{2}\left(\partial_{11}f^{(1)}+2\partial_{110}f^{(2)}\right)$ (24)
$\displaystyle C(k)=-Mk^{2}\partial_{101}f^{(2)}$ (25)

As the matrix $H$ is a symmetric circulant matrix, the expression for the
eigenvalues is easy to find and reads

$\displaystyle\lambda_{m}(k)=A(k)+2B(k)\cos{\dfrac{2\pi
m}{N_{z}}}+2C(k)\cos{\dfrac{4\pi m}{N_{z}}}.$ (26)

Although the number of layers in a graphite particle may be very large, we can
consider a much smaller number of layers in the analysis depending on the
number of stages we want to study. In the following, to detail the expressions
of the different growth rates, we consider $N_{z}=6$ layers which is
compatible with the appearance of both stages 2 and 3. It is useful here to
detail the relations between eigenvectors and stages, and to label the growth
rates by the stage number rather than the index of the eigenvalue.

The eigenvector associated to the mode $m=0$ is uniform (all elements are
$1$). The corresponding growth rate therefore describes the evolution of the
local average concentration. Since this has the same symmetry as stage 1 (all
layers behave the same way) we label the corresponding growth rate with index
$1$.

$\displaystyle\omega_{1}(k)$ $\displaystyle=A+2B+2C$
$\displaystyle=-Mk^{2}\bigg{\\{}\kappa
k^{2}+N_{V}\bigg{[}\dfrac{k_{B}T}{\bar{c}(1-\bar{c})}$
$\displaystyle-2\Omega_{a}+2\Omega_{b}-2\Omega_{c}(3\bar{c}-1)\bigg{]}\bigg{\\}}.$
(27)

The modes $m=1$ and $m=5$ are degenerate due to the up-down symmetry. The
corresponding eigenvectors correspond to discrete plane waves, in the $z$
direction, of wavelength $6L_{z}$, propagating in the positive and negative
$z$ direction, respectively. Therefore, these modes would correspond to a
stage 6, which does not appear in the phase diagram. The corresponding growth
rate is

$\displaystyle\omega_{6}(k)$ $\displaystyle=A+B-C$
$\displaystyle=-Mk^{2}\bigg{\\{}\kappa
k^{2}+N_{V}\bigg{[}\dfrac{k_{B}T}{\bar{c}(1-\bar{c})}$
$\displaystyle-2\Omega_{a}+\Omega_{b}-\Omega_{c}(\bar{c}+1)\bigg{]}\bigg{\\}}$
(28)

Similarly, the modes $m=2$ and $m=4$ are also degenerate, and the
corresponding eigenvectors are discrete plane waves of wavelength 3; these
modes are therefore associated with stage 3, and their growth rate is

$\displaystyle\omega_{3}(k)$ $\displaystyle=A-B-C$
$\displaystyle=-Mk^{2}\bigg{\\{}\kappa
k^{2}+N_{V}\bigg{[}\dfrac{k_{B}T}{\bar{c}(1-\bar{c})}$
$\displaystyle-2\Omega_{a}-\Omega_{b}+\Omega_{c}(3\bar{c}-1)\bigg{]}\bigg{\\}}$
(29)

Finally, the last mode, $m=3$, has an eigenvector proportional to $(-1)^{j}$
and corresponds to stage 2, with a growth rate of

$\displaystyle\omega_{2}(k)$ $\displaystyle=A-2B+2C$
$\displaystyle=-Mk^{2}\bigg{\\{}\kappa
k^{2}+N_{V}\bigg{[}\dfrac{k_{B}T}{\bar{c}(1-\bar{c})}$
$\displaystyle-2\Omega_{a}-2\Omega_{b}+2\Omega_{c}(\bar{c}+1)\bigg{]}\bigg{\\}}$
(30)

An example for a stability spectrum (the growth rates of the different modes
as a function of the wavenumber $k$) is displayed in Fig. 3 for the parameters
from Table 1, and for an average composition $\bar{c}=0.3$, at $T=298~{}K$. A
mode is unstable (and a phase decomposition is observed) when its growth rate
is positive. For example, for stage $2$, the system is unstable for
wavenumbers between $-6.34\times 10^{6}~{}m^{-1}$ and $6.34\times
10^{6}~{}m^{-1}$. On the contrary, we can see that stage $6$ is stable for all
wavenumbers for this temperature and average composition. For this mean
filling fraction of $\bar{c}=0.3$, the growth rate of stage $2$ is higher than
the one of stage $3$. This is remarkable because, according to the phase
diagram, stage $3$ should be present at equilibrium. The theoretical analysis
thus predicts that stage $2$ should emerge during the first few instants of a
spinodal decomposition right after a quench, even though the final state must
be stage $3$.

Figure 3: Growth rates for the different stages for an average composition of
$\bar{c}=0.3$, a temperature of $T=298~{}K$ and the thermodynamic parameters
of Table 1.

###  III.3 Characteristic wave numbers and growth rates

For each mode (stage), there are two noteworthy wavenumbers: the one where the
growth rate reaches its maximum value, $k_{\rm max}$; and the marginal one for
which the growth rate is zero, $k_{0}$. To compute these numbers, the growth
rates are rewritten in the general form

$\displaystyle\omega_{\alpha}(k)=-Mk^{2}[\kappa
k^{2}-\Gamma_{\alpha}(\bar{c})]$ (31)

for $\alpha=1$, $2$, $3$ or $6$, where $\Gamma_{\alpha}$

$\displaystyle\Gamma_{1}$ $\displaystyle=$ $\displaystyle
N_{V}\bigg{[}-\dfrac{k_{B}T}{\bar{c}(1-\bar{c})}+2\Omega_{a}-2\Omega_{b}+2\Omega_{c}(3\bar{c}-1)\bigg{]}$
$\displaystyle\Gamma_{2}$ $\displaystyle=$ $\displaystyle
N_{V}\bigg{[}-\dfrac{k_{B}T}{\bar{c}(1-\bar{c})}+2\Omega_{a}+2\Omega_{b}-2\Omega_{c}(\bar{c}+1)\bigg{]}$
$\displaystyle\Gamma_{3}$ $\displaystyle=$ $\displaystyle
N_{V}\bigg{[}-\dfrac{k_{B}T}{\bar{c}(1-\bar{c})}+2\Omega_{a}+\Omega_{b}-\Omega_{c}(3\bar{c}-1)\bigg{]}$
$\displaystyle\Gamma_{6}$ $\displaystyle=$ $\displaystyle
N_{V}\bigg{[}-\dfrac{k_{B}T}{\bar{c}(1-\bar{c})}+2\Omega_{a}-\Omega_{b}+\Omega_{c}(\bar{c}+1)\bigg{]}$

is the part that originates from interactions and entropy in Eqs. (27) to
(29). Computing the $k_{0}$ and $k_{\rm max}$ values only makes sense for
unstable states, when $\Gamma_{\alpha}(\bar{c})>0$. For $k_{0}$, one needs to
find where $\omega_{\alpha}(k)=0$ with $k\neq 0$,

$k_{0}^{(\alpha)}=\sqrt{\dfrac{\Gamma_{\alpha}}{\kappa}}$ (32)

Figure 3 shows that $k_{\rm max}$ is found when the first derivative of
$\omega_{\alpha}(k)=0$ and $k\neq 0$,

$k_{\max}^{(\alpha)}=\sqrt{\dfrac{\Gamma_{\alpha}}{2\kappa}}.$ (33)

The corresponding maximal growth rate is

$\omega_{\max}^{(\alpha)}=M\frac{\Gamma_{\alpha}^{2}}{4\kappa}.$ (34)

The mode corresponding to this maximal growth rate will dominate the emerging
pattern; therefore, the inverse of its growth rate gives a characteristic time
scale for phase separation, which we will denote by $\tau_{\rm dec}$, and
which depends on the initial composition,

$\tau_{\rm dec}=\frac{4\kappa}{M\Gamma_{\alpha}(\bar{c})^{2}}$ (35)

The growth rates of the different stages depend on the interaction energies
($\Omega_{a}$, $\Omega_{b}$ and $\Omega_{c}$) through the coefficients
$\Gamma_{\alpha}$. The condition $\Gamma_{\alpha}(\bar{c})>0$ determines the
composition range in which mode $\alpha$ is unstable. Since the growth rate
increases with $\Gamma_{\alpha}$, the mode with the highest $\Gamma$ grows the
fastest and has the widest unstable composition range. Of particular interest
is the competition between stages $2$ and $3$. The condition
$\Gamma_{2}>\Gamma_{3}$ is equivalent to

$\dfrac{\Omega_{c}}{\Omega_{b}}<\frac{1}{3-\bar{c}}.$ (36)

For the parameters of Table 1, this is the case for all values of $\bar{c}$,
such that stage $2$ always grows the fastest at short times. It can also be
shown that for these parameters, the onset of spinodal decomposition, obtained
by solving $\Gamma_{\alpha}(\bar{c})=0$, occurs for lower composition for
stage 2 than for stage 3. The competition of growth rates between stages 1 and
2 deserves also some attention. The condition $\Gamma_{2}>\Gamma_{1}$ is
equivalent to

$\dfrac{\Omega_{c}}{\Omega_{b}}<\frac{1}{2\bar{c}}.$ (37)

If we consider all the possible filling fraction, this condition reduces to
$\Omega_{c}<2\Omega_{b}$, which is fulfilled for the parameters of Table 1.
This means that stage 2 will always appear before any potential decomposition
between stage 1 and stage 1’.

As will be shown below (Section V), these conclusions remain valid for other
ratios of the inter-layer free-energies $\Omega_{c}$ and $\Omega_{b}$, as long
as the ratio leads to the right staging sequences.

##  IV Numerical simulations

###  IV.1 Numerical method

To study the kinetics of spinodal decomposition beyond the early times during
which the linearization is valid, numerical simulations are needed. The free-
energy model of section II.1 is introduced in the multi-layer formulation of
Eq. (9). As in the theoretical analysis, a system of 6 layers makes it
possible to capture both stage $2$ and stage $3$. Periodic boundary conditions
are used in the $z$ direction. The system is supposed to be invariant along
the $y$ direction. The particle size is fixed to $25\;\rm{\mu m}$, which
corresponds to a typical size of graphite particles in standard electrodes.
The value of $D$ is taken from Ref. [24], in coherence with the dependence of
the mobility on the local lithium filling fraction (see Eq. 10).

An important issue is the choice of the gradient penalty term $\kappa$. In a
mean-field approximation, its value is related to the strength and spatial
range of attrative in-plane interaction between lithium atoms [29]; an order-
of magnitude estimation of $10^{-10}~{}\rm{J}\;\rm{m}^{-1}$ has been given in
Ref. [30]. Since the gradient enegy penalizes spatial variations of the
concentration, it directly influences the thickness of the diffuse interfaces
between concentrated and dilute domains in the same layer. This thickness,
which can be estimated as $\lambda\approx\sqrt{\kappa/(\Omega_{a}N_{V})}$, is
of the order of a few nanometers for the quoted value of $\kappa$. The
thickness of the interfaces yields a practical constraint, since the spacing
of the numerical grid must be smaller than the interface thickness in order to
properly resolve the interfaces. For this reason, a particle of size $25\;\mu
m$ cannot be simulated with this value of $\kappa$ within reasonable
simulation times. Therefore, in the present work, a larger value is used,
$\kappa=3\times 10^{-6}~{}\rm{J}\;\rm{m}^{-1}$. Some consequences of this
choice will be discussed below.

A finite-volume method is used to discretize the equations in the
$x$-direction. The grid spacing is chosen to ensure that it is small enough to
capture the interfaces between lithium-poor and lithium-rich domains. Given
the stiffness of the problem and to keep the computational time reasonable the
time integration is carried out using an implicit method. The nonlinear system
of equations arising from the implicitation is solved using a Newton method
combined with a GMRES solver for the resolution of the linear problem [57,
58]. The algorithm adapted from [57, 58] has been implemented in an in-house
code, originally developed for the study of diffuse interface models in fluid
mechanics [59].

Spinodal decomposition occurs when a homogeneous phase is linearly unstable
and spontaneously separates into two phases. The simulation of this phenomenon
is possible by imposing a zero lithium flux at both $x=0$ and $x=L$ and by
initializing the simulation with an homogeneous filling fraction $\bar{c}$ in
all the layers. Since the decomposition cannot happen if the composition is
perfectly homogeneous, a small perturbation $\delta\tilde{c}_{j}(x)$ is added
to the average initial composition. In the generic case a randomly uniformly
distributed perturbation of up to $\pm 5\%$ was added to each grid point. For
the comparison between simulation and theory the initial conditions are given
below.

Figure 4: Numerical (symbols) and analytical (lines) growth rates for stage 2
(diamonds and dashed red line) and stage 3 (circles and dash-dotted green
line) for an initial composition $\bar{c}=0.3$, a temperature $T=298K$ and the
parameters of TABLE 1.

###  IV.2 Validation: linear regime

A comparison between the theoretical growth rates predicted by the linear
stability analysis and the simulated ones has been performed. Since the linear
stability analysis is valid only for small $\delta\tilde{c}_{j}$, the
comparison is performed at short times after the quench. Two successive
Fourier transforms are applied to the composition fields obtained from the
simulations. The first is a discrete Fourier transform in the $z$ direction,

$\hat{c}_{m}(x,t)=\dfrac{1}{N_{z}}\sum_{j=0}^{N_{z}-1}\tilde{c}_{j}(x,t)e^{iq_{m}j}$
(38)

where $q_{m}=2\pi m/N_{z}$ is the discrete wave number in the $z$-direction.
This post-treatment allows us to distinguish the stages at all times and
depths in the graphite particle. Furthermore, we perform a continuous Fourier
transform along the $x$ direction,

$\bar{c}_{mn}(t)=\dfrac{1}{L}\int_{0}^{L}\hat{c}_{m}(x,t)e^{ik_{n}x}dx$ (39)

where $k_{n}=2\pi n/L$ is the wavenumber along $x$, which takes discrete
values due to the finite system size. The combination of the two transforms
yields the Fourier coefficient of mode number $m$ with wavenumber $k_{n}$
along $x$,

$\bar{c}_{mn}(t)=\dfrac{1}{N_{z}L}\int_{0}^{L}\left(\sum_{j=0}^{N_{z}-1}\tilde{c}_{j}(x,t)e^{iq_{m}j}\right)e^{ik_{n}x}dx$
(40)

In the numerical implementation, Simpson’s rule has been used to compute the
integral. In the regime immediately after the quench, the logarithm of the
amplitude $|\bar{c}_{mn}(t)|$ varies linearly with time, and a fit of the
slope provides the growth rate.

Using the linear superposition principle, it should be theoretically possible
to extract the growth rates of all the modes from a single simulation with a
random initial fluctuation. However, the precision is deteriorated by roundoff
errors, and therefore simulations were performed for each individual mode and
each wavelength. All the simulations have an average initial composition of
$\bar{c}=0.3$, as in the results of Fig. 3. The perturbations corresponds to a
prescribed mode and a prescribed plane wave. For stage two, the initial
compositions are:

$\displaystyle\tilde{c}_{j}(x,0)=\bar{c}+\dfrac{\cos(k_{n}x+j\pi)}{1000}\quad\text{with}~{}n=[0,\cdots,27],$
(41)

while for stage three:

$\displaystyle\tilde{c}_{j}(x,0)=\bar{c}+\dfrac{\cos(k_{n}x+2j\pi/3)}{1000}$
$\displaystyle\text{with}~{}n=[-23,\cdots,23].$ (42)

These expressions correspond, respectively, to the eigenvectors $m=3$ and
$m=2$ of Eqs. (20) and (22), with amplitudes $a_{3n}=a_{2n}=1/1000$. The
number of simulated modes was determined as a function of system size and the
theoretically calculated marginal wavenumber such as to cover all the relevant
range of the stability spectrum. The amplitude of the fluctuations ($1/1000$)
was small enough to ensure a sufficiently long duration of the linear regime
for obtaining good fits of the growth rates. A total number of 75 simulations
have been performed to analyse the growth rate of each individual mode and
each wavelength. Figure 4 shows the comparison of the growth rates obtained
from the numerical simulations and from the linear stability analysis. The
excellent agreement between theory and simulation confirms the validity of the
linear stability analysis as well as the implementation of the numerical
model.

Figure 5: Top: Snapshots of spinodal decompositions for two different initial
mean filling fraction : $\bar{c}=0.5$ (left) and $\bar{c}=0.2$ (right). The
snapshots are taken slighty after the end of the linear regime to achieve a
sufficient contrast between dilute and dense domains. Bottom: Corresponding
linear stability analysis growth rate curves with an emphasis on the wave
number for which the stage 2 is at its maximum, $k_{max}^{(2)}$ (black dashed
lines).

All the plane waves with $k<k_{0}$ will grow exponentially in time, the
fastest one being the closest to $k_{\rm max}$. Therefore, at the end of the
linear stage the size of the lithium islands should be of the order of
$2\pi/k_{\rm max}$. This prediction can be verified by performing simulations
of spinodal decompositions with random initial perturbation for two different
initial mean filling fractions, $\bar{c}=0.5$ and $\bar{c}=0.2$. The stability
spectra are presented in Fig.5 together with snapshots of the simulations at
early times. For $\bar{c}=0.5$, $k_{\rm max}$ is expected to be larger than
for $\bar{c}=0.2$, leading to smaller islands. This is consistent with the
snapshot pictures, where the wavelength of the fastest growing mode,
$2\pi/k_{\rm max}$, is displayed for comparison. The agreement is very
satisfactory.

The linear stability analysis also predicts the growth rate of the
perturbations. At $\bar{c}=0.5$, the growth rate of the leading mode
$\omega_{2}(k_{\rm max})$ is of the order $10.5~{}s^{-1}$, corresponding to a
characteristic decomposition time $\tau_{\rm dec}\approx 0.1s$, in coherence
with the formation of well-developed domains in less than half a second, as
can be seen Fig. 5 on the left. At lower filling fraction, $\bar{c}=0.2$, the
growth rate of the leading mode is much lower, $\omega_{2}(k_{\rm max})\approx
0.034~{}s^{-1}$ corresponding to a characteristic decomposition time
$\tau_{\rm dec}\approx 30s$. Even after 162 seconds, one can barely identify
the domains in the snapshot picture of Fig. 5 on the right. While for
$\bar{c}=0.3$, stage 6 is stable for any wave number (Fig. 3), we can notice
that it is slightly unstable for $\bar{c}=0.5$ as can be noticed on Fig 5
(left). However, since it is not the dominant unstable mode, it will not
easily be observed.

###  IV.3 Coarsening: the non-linear regime

Figure 6: Snapshots of a spinodal decomposition at different times (see also
movie in the Supplemental Material [60]). The initial average composition is
$\bar{c}=0.3$ and the system has $N_{z}=6$ layers. For each subfigure, top:
lithium content in each layer, bottom: coefficients of the Fourier transform
$|\hat{c}_{m}(x,t)|$ (see Eq. 38).

Figure 6 presents snapshots at different times of the spinodal decomposition
starting from an initial filling fraction of $\bar{c}=0.3$ with a random
initial fluctuation. At the top part of each subfigure, the lithium content in
all of the 6 layers and for the whole particle depth can be seen in the form
of a grey-scale heat-map. The bottom part displays the result of the discrete
Fourier transform (Eq. (38)): the amplitudes of the Fourier coefficients can
be used to identify quantitatively the stages present at each position within
the particle and thus to follow the evolution of the system beyond the linear
regime.

The particle starts from a homogeneous state with small fluctuations (a). Very
quickly, all the layers undergo spinodal decomposition. As predicted by the
linear stability analysis, stage $2$ grows the fastest, and a corresponding
checkerboard pattern starts to appear in the heat map, whereas the average
composition still remains homogeneous (b). The amplitude of this pattern
increases, until domains of high (close to $\tilde{c}=1$) and low (close to
$\tilde{c}=0$) composition clearly emerge and form vertical columns (c). The
average composition is now higher in the center of the columns than in the
frontier between two columns. This is a clear indication of a non-linear
interaction between different modes, and thus the linear stability analysis
and Eq. (III.2) do not apply anymore.

Since the simulation was started from a random perturbation, the size of the
islands (columns) is not everywhere the same, and the pattern is slightly
disordered. This subsequently leads to coarsening events: some domains expand
in size, while others disappear (d). During this coarsening, stage 3
progressively appears, even though the global average filling fraction stays
constant (e). This is visible on the greyscale map, but also on the plot of
the local Fourier coefficients. One particular mechanism by which stage 3
appears in regions previously occupied by stage 2 is the following: the
disappearance of one concentrated domain in a single layer is followed by a
sequence of elimination events in the other layers which makes an entire
“diagonal row” of concentrated domains disappear. Subsequently, the domains
located on both sides of the eliminated row grow closer to each other, until
they locally reach a stage 3 pattern. Several such events can be seen in the
movie included in the Supplemental Material [60], for example in the center of
the system at $t\approx 100$s. At the end of the simulation, the entire system
is in stage 3 (f).

The column structure is nicely represented by the Fourier coefficients
$\hat{c}_{m}(x,t)$ defined in Eq. (38). They are complex numbers, with an
amplitude and a phase; only the amplitude is plotted in Fig. 6. In stage 2,
two neighboring columns have the same amplitude, but a phase difference of
$\pi$. In stage 3, two neighboring columns can have a phase difference of $\pm
2\pi/3$, which correspond to the “upward and downward staircases” visible in
Fig. 6; these two different states reflect the two degenerate linear modes
that contribute to stage 3. In the Fourier representation, the boundary
between two columns is thus a “phase jump”. Since, in each layer, the
composition variation in space is controlled by the gradient energy
coefficient, the thickness of such a domain wall should be close to the one of
an interface between dilute and concentrated domains in a single layer.

To follow quantitatively the global presence of the different stages during
the coarseing, the amplitude of the Fourier coefficients $\hat{c}_{m}(x,t)$
can be averaged over the length of the particle from $x=0$ to $x=L$. The time
evolution of these averages, noted
$\langle|\hat{c}_{m}(x,t)|\rangle_{x\in[0,L]}$, are presented in Fig. 7. This
figure clearly shows that stage $2$ grows the fastest at early times in
accordance with the linear stability analysis. After a few seconds, stage $3$
appears and increases to finally dominate after 45s, as expected from the
result of the phase diagram for an initial mean filling fraction of
$\bar{c}=0.3$.

Simulations with the multi-layer Cahn-Hilliard framework are therefore a
powerful tool to analyze the dynamics of stage formation and evolution, and to
bridge the gap between the early-time dynamics predicted by the linear
stability analysis and the final equilibrium state predicted by the phase
diagram.

Figure 7: Time evolution of the global stage fractions
$\langle|\hat{c}_{m}(x,t)|\rangle_{x\in[0,L]}$ during the spinodal
decomposition at $T=298K$, with an initial composition of $\bar{c}=0.3$.

##  V Discussion

The most noteworthy result of the preceding analytical and numerical studies
is that, for the interaction energy parameters listed in Table 1, stage 2
always grows faster than stage 3, even in the composition range where stage 3
is stable at thermodynamic equilibrium. As already stated in Eq. (36), stage 3
would be favored if the ratio $\Omega_{c}/\Omega_{b}$ was large enough.
However, in Ref. [41] it was found that in order to obtain the right sequence
of stages for increasing filling fraction (and, in particular, to avoid the
occurrence of the non-physical stage 3/2), $\Omega_{c}/\Omega_{b}$ has to lie
within a certain range. We have not found a closed-form analytic expression
for this latter condition, but we have numerically investigated the phase
diagram for a large range of $\Omega_{b}$ and $\Omega_{c}$. As an example, we
show in Fig 8 the free-energy curves for all the stages together with their
convex envelopes and the corresponding stage sequence obtained for a fixed
ratio of $\Omega_{a}/\Omega_{b}=2.78$, a fixed temperature of 298K, and
various $\Omega_{c}/\Omega_{b}$ ratio. Stage 3/2 appears before the condition
$\omega_{3}>\omega_{2}$ is satisfied. We have not found any set of values for
which both the equilibrium sequence of stages was correct and stage 3 grew the
fastest after a quench. Therefore, we must conclude that, in the framework of
a free-energy model of the form of Eqs. (3) and (II.1), the fast growth of
stage 2 in the spinodal decomposition of a homogeneously filled sample is a
general phenomenon. This methodology could easily be extended to other intra-
and interlayer interactions energies [24, 36, 38] to determine if this finding
holds for an even wider class of free-energy models.

Figure 8: Sequence of stages occurring with increasing filling fraction, for
various ratio $\Omega_{c}/\Omega_{b}$. The ratio $\Omega_{c}/\Omega_{b}$
increases from left to right. The values of $\Omega_{b}$ and $\Omega_{c}$ are
non-dimensionalized using the prefactor $k_{B}T$ at $298K$. Above the arrow:
free-energy curves of all the stages (black: stage 1 or dilute, red: stage 2,
green: stage 3 and 3/2) together with the convex envelope (cyan). Below the
arrow: resulting scheme of the staging sequence. Ratio $\Omega_{c}/\Omega_{b}$
below 0.06: no stage 3; between 0.06 and 0.3: right sequence, between 0.3 and
0.47: occurrence of stage 3/2; between 0.47 and 0.55: no stage 2 and
occurrence of stage 3/2; above 0.55 no stage 2.

Let us now discuss the consequences of our findings for the situation, more
relevant for comparisons to experiments, in which an initially empty electrode
is progressively charged. The active material initially is in the dilute stage
1. Lithium atoms are inserted into the material at the interface with the
electrolyte, and the distribution of atoms inside the material proceeds by
diffusion. There are two distinct characteristic time scales [61, 62]: the
diffusion time $\tau_{D}\sim L^{2}/D$, and the processing time scale (charging
time) $\tau_{I}\sim Q/iA$, where $i$ is the externally imposed charging
current density, $Q$ is the total capacity of the active material, and $A$ is
the surface area over which the current is applied. For slow charging, such
that $\tau_{I}\gg\tau_{D}$, diffusion is fast enough to distribute the lithium
atoms through the entire particle, and the concentration remains almost
homogeneous. The global concentration thus slowly increases until it reaches
the limit where spinodal decomposition can occur. As already discussed above,
the first mode that becomes unstable is stage 2. Therefore, the system should
undergo a sequence of states that is comparable to our simulations of Fig. 6.
It should be noted that the decomposition time given by Eq. (35) depends on
the concentration and is infinite at the spinodal. Therefore, the
decomposition will become visible at a concentration that is larger than the
critical concentration for instability, and that depends on the charging rate.
Indeed, the instability will occur approximately when the concentration is
such that $\tau_{\rm dec}(\bar{c})=\tau_{I}$.

In the opposite limit, where $\tau_{I}<\tau_{D}$, inhomogeneities in
concentration develop inside the active material. Since we start from a dilute
system in which the diffusion is well described by Fick’s law, the vicinity of
the electrolyte will always be the point with the highest concentration, and
will therefore reach the threshold for instability first. Since the
composition threshold for the onset of spinodal decomposition is the lowest
for stage 2, and its growth rate always remains the highest for compositions
beyond the threshold, the present theory predicts that stage 2 forms at the
surface of the electrode. This was indeed observed in the preliminary
simulations of the charging process in Ref. [41]. Whether this surface layer
of stage 2 persists or disappears in the subsequent evolution is an important
question, because the occurrence of staging at the surface will modify the
effective resistance for the entry of new atoms into the material. To address
this issue, interfacial intercalation kinetics with its potential
nonlinearities must also be taken into account [63, 64]. This issue will be
adressed numerically in a forthcoming publication through detailed simulations
of the charge-discharge dynamics.

Furthermore, these results imply that at high charging rates, stage 3 should
disappear completely. This has not been reported experimentally so far, but
could be confirmed using high speed synchrotron X-ray diffraction [19, 65,
21]. It would require to specifically follow the reflections characteristics
of the stage 3 ordering, i.e. the examination of reflections on planes
separated by about 3 times the graphite interlayer distance d. This is
possible at low charging rates [66], but remains challenging at higher C-rate,
given their low intensity.

As already mentioned, we have used for our simulations a value of the gradient
energy coefficient $\kappa$ that is much larger than the one expected from a
mean-field approximation. This is an issue because $\kappa$ appears in the
expression for the growth rates and for the characteristic wavelengths, Eqs.
(33) through (32). The typical initial size of the domains is $2\pi/k_{\rm
max}$; therefore, a measurement of the domain size at early times after a
spinodal decomposition could be an indirect way to determine the order of
magnitude of $\kappa$ from experimental data. To our knowledge, no such
measurement is available as of yet for lithium intercalated in graphite.

Another point linked to the value of $\kappa$ is the domain wall energy, which
scales as the geometric mean of the gradient energy coefficient and the
barrier height of the double-well potential in the free energy, $\Omega_{a}$.
All of the interaction energies have been obtained by comparison of the model
prediction to measurement of the equilibrium potential [41], and are therefore
physically realistic. As a consequence, the use of a larger value of $\kappa$
also implies that the energies of domain walls between domains and columns are
larger than the physically correct value. We have not investigated the impact
of this fact on the coarsening dynamics. Anyway, our present simulations of
coarsening are limited to small systems (six layers), and to two dimensions
(one-dimensional galleries). It would be interesting but computationally
costly to perform simulations in larger system and in three dimensions in
order to assess the coarsening dynamics more precisely.

##  VI Conclusion

In this paper, the dynamics of spinodal decomposition and domain coarsening in
lithium-intercalated graphite after a quench (a sudden change in temperature)
has been investigated, both by an analytical linear stability analysis and by
numerical simulations of a multi-layer Cahn-Hilliard model. The multi-layer
framework is needed to capture the well-documented phenomenon of staging,
which is caused by repulsive inter-layer interactions between lithium atoms.

A normal-mode analysis of the linearized evolution equations has allowed us to
follow separately the growth dynamics of each stage, and to compute
analytically the associated growth rates. A stability spectrum that is typical
for spinodal decomposition has been obtained for each stage, in the form of a
quartic polynomial in the wave number for plane waves along the graphite
planes. Beyond the linear regime, the decomposition in normal modes,
corresponding to a discrete Fourier transform, remains useful to analyze the
further evolution of the stage distribution, which can only be obtained by
numerical simulations of the full evolution equations.

Our analysis predicts that stage 2 is always the fastest to grow after a
quench, even for average compositions for which stage 3 is the
thermodynamically stable state. Stage 3 emerges only at later times, during
the nonlinear coarsening of stage 2, and progressively fills the entire
system. While quench experiments are difficult to perform, our findings also
yield a prediction for the kinetics of stage formation during the charging of
a graphite particle: at the particle surface, stage 2 should always appear
first. In contrast, if the equilibrium sequence of stages is followed, stage 3
(or higher) should appear first. How this competition between stages at the
surface depends on the charge rate remains to be investigated in more detail.
Further work could also focus on the introduction, in the stage stability
analysis, of electrochemical reactions at the particle surface [63, 64, 42].
This is necessary to understand if insertion reactions can impact the
stability of the different stages and therefore modify the appearance of the
different stages during lithiation/delithiation processes.

## Acknowledgements

The CEA Battery FOCUS program on ”Multi-scale simulation of batteries applied
to electrode materials” is acknowledged for the funding of the PhD program of
Antoine Cordoba.

## References

  * Scrosati [2011] B. Scrosati, J. Solid State Electrochem. 15, 1623 (2011).
  * Xie and Lu [2020] J. Xie and Y.-C. Lu, Nat. Commun. 11, 2499 (2020).
  * Dunn  _et al._ [2011] B. Dunn, H. Kamath, and J.-M. Tarascon, Science 334, 928 (2011).
  * Armand  _et al._ [2020] M. Armand, P. Axmann, D. Bresser, M. Copley, K. Edström, C. Ekberg, D. Guyomard, B. Lestriez, P. Novák, M. Petranikova, W. Porcher, S. Trabesinger, M. Wohlfahrt-Mehrens, and H. Zhang, J. Power Sources 479, 228708 (2020).
  * Andersen  _et al._ [2021] H. L. Andersen, L. Djuandhi, U. Mittal, and N. Sharma, Adv. Energy Mater. 11, 2102693 (2021).
  * Asenbauer  _et al._ [2020] J. Asenbauer, T. Eisenmann, M. Kuenzel, A. Kazzazi, Z. Chen, and D. Bresser, Sustain. Energy Fuels 4, 5387 (2020).
  * Dahn [1991] J. R. Dahn, Phys. Rev. B 44, 9170 (1991).
  * Ohzuku  _et al._ [1993] T. Ohzuku, Y. Iwakoshi, and K. Sawai, J. Electrochem. Soc. 140, 2490 (1993).
  * Billaud  _et al._ [1996] D. Billaud, F. Henry, M. Lelaurain, and P. Willmann, J. Phys. Chem. Solids 57, 775 (1996).
  * Meunier  _et al._ [2002] V. Meunier, J. Kephart, C. Roland, and J. Bernholc, Phys. Rev. Lett. 88, 075506 (2002).
  * Yao  _et al._ [2012] F. Yao, F. Güneş, H. Q. Ta, S. M. Lee, S. J. Chae, K. Y. Sheem, C. S. Cojocaru, S. S. Xie, and Y. H. Lee, J. Am. Chem. Soc. 134, 8646 (2012).
  * Krishnan  _et al._ [2013] S. Krishnan, G. Brenet, E. Machado-Charry, D. Caliste, L. Genovese, T. Deutsch, and P. Pochet, Appl. Phys. Lett. 103, 251904 (2013).
  * Li  _et al._ [2021] N. Li, M.-Z. Sun, S. Hwang, S. Li, H.-Y. Zhao, Y.-P. Du, B.-L. Huang, and D. Su, J. Mater. Chem. A 9, 12080 (2021).
  * Rüdorff and Hofmann [1938] W. Rüdorff and U. Hofmann, Z. Für Anorg. Allg. Chem. 238, 1 (1938).
  * Guerard and Herold [1975] D. Guerard and A. Herold, Carbon 13, 337 (1975).
  * Billaud  _et al._ [1981] D. Billaud, E. McRae, J. Mareche, and A. Herold, Synthetic Metals 3, 21 (1981).
  * Allart  _et al._ [2018] D. Allart, M. Montaru, and H. Gualous, J. Electrochem. Soc. 165, A380 (2018).
  * Hess and Novak [2013] M. Hess and P. Novak, Electrochimica Acta 106, 149 (2013).
  * Yao  _et al._ [2019] K. P. C. Yao, J. S. Okasinski, K. Kalaga, I. A. Shkrob, and D. P. Abraham, Energy Environ. Sci. 12, 656 (2019).
  * Didier  _et al._ [2020] C. Didier, W. K. Pang, Z. Guo, S. Schmid, and V. K. Peterson, Chem. Mater. 32, 2518 (2020).
  * Schmitt  _et al._ [2021] C. Schmitt, A. Kube, N. Wagner, and K. A. Friedrich, ChemElectroChem 9 (2021), 10.1002/celc.202101342.
  * Daumas and Herold [1969] N. Daumas and A. Herold, C. R. Acad. Sci. 268 (1969).
  * Bazant [2013] M. Z. Bazant, Acc. Chem. Res. 46, 1144 (2013).
  * Smith  _et al._ [2017] R. B. Smith, E. Khoo, and M. Z. Bazant, J. Phys. Chem. C 121, 12505 (2017).
  * Cahn and Hilliard [1958] J. W. Cahn and J. E. Hilliard, J. Chem. Physics 28, 258 (1958).
  * Cahn and Hilliard [1959a] J. W. Cahn and J. E. Hilliard, J. Chem. Physics 30, 1121 (1959a).
  * Cahn and Hilliard [1959b] J. W. Cahn and J. E. Hilliard, J. Chem. Physics 31, 688 (1959b).
  * Hohenberg and Halperin [1977a] P. C. Hohenberg and B. I. Halperin, Rev. Mod. Phys. 49, 435 (1977a).
  * Langer [1991] J. S. Langer, in  _Solids far from equilibrium_ , Edition Aléa Saclay, edited by C. Godrèche (Cambridge University Press, Cambridge, UK, 1991) pp. 297–363.
  * Han  _et al._ [2004] B. Han, A. V. der Ven, D. Morgan, and G. Ceder, Electrochim. Acta 49, 4691 (2004).
  * Singh  _et al._ [2008] G. K. Singh, G. Ceder, and M. Z. Bazant, Electrochim. Acta 53, 7599 (2008).
  * Cogswell and Bazant [2012] D. A. Cogswell and M. Z. Bazant, ACS Nano 6, 2215 (2012).
  * Li  _et al._ [2018] Y. Li, H. Chen, K. Lim, H. D. Deng, J. Lim, D. Fraggedakis, P. M. Attia, S. C. Lee, N. Jin, J. Moškon, Z. Guan, W. E. Gent, J. Hong, Y.-S. Yu, M. Gabersček, M. S. Islam, M. Z. Bazant, and W. C. Chueh, Nat. Mater. 17, 915 (2018).
  * Zelič and Katrašnik [2019] K. Zelič and T. Katrašnik, Sci. Rep. 9, 2123 (2019).
  * Daubner  _et al._ [2021] S. Daubner, P. Kubendran Amos, E. Schoof, J. Santoki, D. Schneider, and B. Nestler, Phys. Rev. Mat. 5, 035406 (2021).
  * Safran [1980] S. Safran, Phys. Rev. Lett. 44, 937 (1980).
  * Millman and Kirczenow [1983] S. E. Millman and G. Kirczenow, Phys. Rev. B 28, 3482 (1983).
  * Hawrylak and Subbaswamy [1984] P. Hawrylak and K. R. Subbaswamy, Phys. Rev. Lett. 53, 2098 (1984).
  * Ferguson and Bazant [2014] T. R. Ferguson and M. Z. Bazant, Electrochim. Acta 146, 89 (2014).
  * Guo  _et al._ [2016] Y. Guo, R. B. Smith, Z. Yu, D. K. Efetov, J. Wang, P. Kim, M. Z. Bazant, and L. E. Brus, J. Phys. Chem. Lett. 7, 2151 (2016).
  * Chandesris  _et al._ [2019] M. Chandesris, D. Caliste, D. Jamet, and P. Pochet, J. Phys Chem. C 123, 23711 (2019).
  * Agrawal and Bai [2022] S. Agrawal and P. Bai, Cell Reports Physical Science 3, 100854 (2022).
  * Thomas-Alyea  _et al._ [2017] K. E. Thomas-Alyea, C. Jung, R. B. Smith, and M. Z. Bazant, Journal of The Electrochemical Society 164, E3063 (2017).
  * Ahn  _et al._ [2023] S. Ahn, M. Lagnoni, Y. Yuan, A. Ogarev, E. Vavrinyuk, G. Voynov, E. Barrett, A. Pelli, A. Atrashchenko, A. Platonov, S. Gurevich, M. Gorokhov, D. Rupasov, A. W. Robertson, R. A. House, L. R. Johnson, A. Bertei, and D. V. Chernyshov, ACS Applied Energy Materials 6, 8455 (2023).
  * Lian and Bazant [2023] H. Lian and M. Z. Bazant, preprint (2023).
  * Lu  _et al._ [2023] X. Lu, M. Lagnoni, A. Bertei, S. Das, R. E. Owen, Q. Li, K. O’Regan, A. Wade, D. P. Finegan, E. Kendrick, _et al._ , Nature Communications 14, 5127 (2023).
  * Danner  _et al._ [2016] T. Danner, M. Singh, S. Hein, J. Kaiser, H. Hahn, and A. Latz, Journal of Power Sources 334, 191 (2016).
  * Franco  _et al._ [2019] A. A. Franco, A. Rucci, D. Brandell, C. Frayret, M. Gaberscek, P. Jankowski, and P. Johansson, Chemical Reviews 119, 4569 (2019).
  * Taghikhani  _et al._ [2020] K. Taghikhani, P. J. Weddle, J. R. Berger, and R. J. Kee, Journal of The Electrochemical Society 167, 110554 (2020).
  * de Groot and Mazur [1962] S. R. de Groot and P. Mazur,  _Non-equilibrium thermodynamics_ , edited by N.-H. Co (Intersciences publisher, 1962).
  * Kawasaki [1966] K. Kawasaki, Phys. Rev. 145, 224 (1966).
  * Hohenberg and Halperin [1977b] P. C. Hohenberg and B. I. Halperin, Rev. Mod. Phys. 49, 435 (1977b).
  * Wang  _et al._ [2022] A. A. Wang, S. E. J. O’Kane, F. Brosa Planella, J. Le Houx, K. O’Regan, M. Zyskin, J. Edge, C. W. Monroe, S. J. Cooper, D. A. Howey, E. Kendrick, and J. M. Foster, Prog. Energy 4, 032004 (2022).
  * Fischer  _et al._ [1983] J. E. Fischer, C. D. Fuerst, and K. C. Woo, Synth. Met. 7, 1 (1983).
  * Hwang  _et al._ [1990] D. M. Hwang, G. Kirczenow, L. P., A. Magerl, R. Moret, S. Moss, S. Setton, S. Solin, and H. Zabel,  _Graphite intercalation compounds I: Structure and Dynamics_ , edited by H. Zabel and S. A. Solin (Springer Series in Material Science, 1990).
  * Rykner and Chandesris [2022] M. Rykner and M. Chandesris, J. Phys. Chem. C 126, 5457 (2022).
  * Brown and Saad [1990] P. N. Brown and Y. Saad, SIAM Journal on Scientific and Statistical Computing 11, 450 (1990).
  * Choquet and Erhel [1996] R. Choquet and J. Erhel, International Journal for Numerical Methods in Fluids 23, 177 (1996).
  * Jamet [2005] D. Jamet, “Diffuse interface models in fluid mechanics,” Semantic Scholar, Corpus Id: 18437499 (2005).
  * See Supplemental Material at [URL to be inserted by publisher] See Supplemental Material at (URL to be inserted by publisher) for a movie of the simulation displayed in the snapshots of Fig. 6, .
  * Nadkarni  _et al._ [2018] N. Nadkarni, E. Rejovitsky, D. Fraggedakis, C. V. Di Leo, R. B. Smith, P. Bai, and M. Z. Bazant, Phys. Rev. Mater. 2, 085406 (2018).
  * Fraggedakis  _et al._ [2020] D. Fraggedakis, N. Nadkarni, T. Gao, T. Zhou, Y. Zhang, Y. Han, R. M. Stephens, Y. Shao-Horn, and M. Z. Bazant, Energy Environ. Sci. 13, 2142 (2020).
  * Bazant [2017] M. Z. Bazant, Faraday Discuss. 199, 423 (2017).
  * Bazant [2023] M. Z. Bazant, Faraday Discuss. 246, 60 (2023).
  * Finegan  _et al._ [2020] D. P. Finegan, A. H. Quinn, D. Wragg, A. Colclasure, X. Lu, C. Tan, T. Heenan, R. Jervis, D. Brett, S. Das, T. Gao, D. Cogswell, M. Z. Bazant, M. di Michiel, S. Checchia, P. Shearing, and K. Smith, Energy Environ. Sci. , (2020).
  * Tardif  _et al._ [2021] S. Tardif, N. Dufour, J.-F. Colin, G. Gebel, M. Burghammer, A. Johannes, S. Lyonnard, and M. Chandesris, Journal of Materials Chemistry A 9, 4281 (2021).