Population dynamics of driven autocatalytic reactive mixtures

TL;DR

This paper uses the Fokker-Planck equation to explore how autocatalysis and reaction kinetics influence phase separation and population dynamics in driven reactive mixtures, with implications for batteries and materials science.

Contribution

It introduces a novel application of the Fokker-Planck framework to analyze how autocatalysis affects phase behavior beyond thermodynamic predictions.

Findings

Autocatalysis can induce phase separation in single-phase thermodynamic systems.

Autoinhibitory reactions suppress phase separation.

Asymmetric kinetics lead to different population dynamics upon reaction reversal.

Abstract

Motivated by the theory of reaction kinetics based on nonequilibrium thermodynamics and the linear stability of driven reaction-diffusion, we apply the Fokker-Planck equation to describe the population dynamics of an ensemble of reactive particles in contact with a chemical reservoir. We illustrate the effect of autocatalysis on the population dynamics by comparing systems with identical thermodynamics yet different reaction kinetics. The dynamic phase behavior of the system may be entirely different from what its thermodynamics may suggest. By defining phase separation for a particle ensemble to be when the probability distribution is bimodal, we find that thermodynamic phase separation may be suppressed by autoinhibitory reactions, while autocatalysis enhances phase separation and in some cases induce the ensemble that consists of thermodynamically single-phase systems to segregate into two distinct populations, which we term fictitious phase separation. Asymmetric reaction kinetics also results in qualitatively different population dynamics upon reversing the reaction direction. In the limit of negligible fluctuations, we use method of characteristics and linearization to study the evolution of the standard deviation of concentration as well as the condition for phase separation, in good agreement with the full numerical solution. Applications are discussed to Li-ion batteries and {\it in situ} x-ray diffraction.

Figures

Tables

Table 1. Table 1: Definitions of thermodynamic quantities. ⟨ ⋅ ⟩ delimited-⟨⟩ ⋅ \langle{\cdot}\rangle denotes averaging over states in the summation that comes immediately before its appearance. β = 1 / k B ​ T 𝛽 1 subscript 𝑘 B 𝑇 \beta=1/k_{\text{B}}T .

Ensemble type	Free energy	Entropy
Canonical		$S_N=-k_{\text{B}}{\sum }_{{\nu }_N}{p}_{{\nu }_N}\ln p_{{\nu }_N}$
Grand canonical
Ensemble type	Equilibrium distribution
Canonical	$p_{{\nu }_N}^{\text{eq}}=\exp \left[-\beta \left(E_{{\nu }_N}+A_N^{\text{eq}}\right)\right]$
Grand canonical	$p_{\nu }^{\text{eq}}=\exp \left[-\beta \left(E_{\nu }+{\mu }_{\text{res}}{N}_{\nu }+{\Psi }^{\text{eq}}\right)\right],$           $p_N^{\text{eq}}=\exp \left[-\beta \left(A_N^{\text{eq}}+{\mu }_{\text{res}}N+{\Psi }^{\text{eq}}\right)\right]$

Table 2. Table 2: Maximum standard deviation ratio with the starting and ending average fraction for three different exchange current.

Exchange current	$1\left(\Omega >2\right)$	$\sqrt{c\left(1-c\right)}\left(\Omega >1\right)$	$\left(1-c\right)e^{\mu /2}\left(\Omega >1/2\right)$
${\max }_{c_1,c_2}\frac{{\sigma }_{c_1}}{{\sigma }_{c_2}}$	$\frac{1-\sqrt{1-2/\Omega }}{1+\sqrt{1-2/\Omega }}\exp \left(\Omega \sqrt{1-2/\Omega }\right)$	$\frac{1}{\Omega }{e}^{\Omega -1}$	$\frac{1}{2\Omega }{e}^{2\Omega -1}$
$c_1$	$\frac{1\mp \sqrt{1-2/\Omega }}{2}$	$0$ or $1$	$1$
$c_2$	$\frac{1\pm \sqrt{1-2/\Omega }}{2}$	$1-\frac{1}{\Omega }$ or $\frac{1}{\Omega }$	$\frac{1}{2\Omega }$

Equations

$${\text{M}}_{\text{res}}\to \text{M},$$

$$\frac{\mathrm{d}c}{\mathrm{d}t}=R(c,\mu ,{\mu }_{\text{res}}),$$

$$s=\frac{\partial{R}}{\partial{c}}=\mathopen{}\mathclose{{}\left(\frac{\partial{R}}{\partial{c}}}\right)_{\mu}+\mathopen{}\mathclose{{}\left(\frac{\partial{R}}{\partial{\mu}}}\right)_{c}\frac{\mathrm{d}\mu}{\mathrm{d}c}.$$

$$\frac{\partial{f}}{\partial{t}}+\frac{\partial{\mathopen{}\mathclose{{}\left(fR}\right)}}{\partial{c}}=0,$$

$$S=S^{\prime }+⟨S_N⟩,$$

$$\Psi=\sum_{N}{p_{N}\mathopen{}\mathclose{{}\left(A_{N}^{\text{eq}}-\mu_{\text{res}}N+k_{\text{B}}T\ln{p_{N}}}\right)}.$$

$$\mathcal{L}=\Psi -{\Psi }^{\text{eq}}=k_{\text{B}}T\sum _{\nu }{p}_{\nu }\ln \frac{p_{\nu }}{p_{\nu }^{\text{eq}}}>0.$$

$$\frac{\partial p}{\partial t}+\nabla \cdot \mathbf{J}=0,$$

$$\mathbf{J}=-\mathbf{K}p\nabla \frac{\delta \Psi }{\delta p}=-\mathbf{K}{k}_{\text{B}}Tp\nabla \ln \frac{p}{p^{\text{eq}}},$$

$$\frac{\mathrm{d}\Psi}{\mathrm{d}t}=-\int_{\mathbb{R}^{n}}{\frac{\mathbf{J}^{\dagger}\mathbf{K}^{-1}\mathbf{J}}{p}\,\mathrm{d}\mathbf{x}}+\int_{\Omega}{p\frac{\mathrm{d}}{\mathrm{d}t}\mathopen{}\mathclose{{}\left(E-\mu_{\text{res}}N}\right)\,\mathrm{d}\mathbf{x}},$$

$${\dot{W}}_d=-{\dot{\mu }}_{\text{res}}⟨N⟩.$$

$$\frac{\partial p_N}{\partial t}+\frac{\partial J_N}{\partial N}=0,$$

$$J_N$$

$$\frac{\delta \Psi }{\delta p_N}$$

$$\frac{\partial f}{\partial t}+\frac{\partial\mathopen{}\mathclose{{}\left(\mathopen{}\mathclose{{}\left(\mu_{\text{res}}-\mu}\right)kf}\right)}{\partial c}=\frac{k_{\text{B}}T}{N_{t}}\frac{\partial}{\partial c}\mathopen{}\mathclose{{}\left(k\frac{\partial f}{\partial c}}\right),$$

$$\frac{\partial f}{\partial t}+\frac{\partial j}{\partial c}=0,$$

$$j=\mathopen{}\mathclose{{}\left(\mathopen{}\mathclose{{}\left(\mu_{\text{res}}-\mu}\right)k-k\frac{k_{\text{B}}T}{N_{t}}\frac{\partial}{\partial c}}\right)f=kk_{\text{B}}Tf^{\text{eq}}\frac{\partial}{\partial c}\mathopen{}\mathclose{{}\left(\frac{f}{f^{\text{eq}}}}\right),$$

$$\Psi_{m}=\int_{0}^{1}{\mathopen{}\mathclose{{}\left(a-\mu_{\text{res}}c+\frac{k_{\text{B}}T}{N_{t}}\ln{f}}\right)f\,\mathrm{d}c},$$

$$f^{\text{eq}}(c)=e^{-\beta N_{t}\mathopen{}\mathclose{{}\left(a-\mu_{\text{res}}c-\Psi_{m}^{\text{eq}}}\right)},$$

$$\frac{\mathrm{d}c}{\mathrm{d}t}=R+\sqrt{2\frac{k_{\text{B}}Tk}{N_t}}\xi (t),$$

$$\frac{\partial f}{\partial t}+\frac{\partial\mathopen{}\mathclose{{}\left(fR}\right)}{\partial c}=\frac{\partial}{\partial c}\mathopen{}\mathclose{{}\left(D\frac{\partial f}{\partial c}}\right),$$

$$\frac{\partial f}{\partial t}+\frac{\partial j}{\partial c}=0,$$

$$\frac{\mathrm{d}\sigma_{c}^{2}}{\mathrm{d}t}=\frac{\mathrm{d}}{\mathrm{d}t}\mathopen{}\mathclose{{}\left(\int_{0}^{1}{(c-c_{0})^{2}f\,\mathrm{d}c}}\right)\\ =\int_{0}^{1}{(c-c_{0})^{2}\frac{\partial f}{\partial t}\,\mathrm{d}c}=2\int_{0}^{1}{j(c-c_{0})\,\mathrm{d}c}.$$

$$j=f\mathopen{}\mathclose{{}\left(\sum_{m=0}^{\infty}{R^{(m)}(c_{0})\frac{\mathopen{}\mathclose{{}\left(c-c_{0}}\right)^{m}}{m!}}}\right)\\ -\frac{k_{\text{B}}T}{N_{t}}\mathopen{}\mathclose{{}\left(\sum_{m=0}^{\infty}{k^{(m)}(c_{0})\frac{\mathopen{}\mathclose{{}\left(c-c_{0}}\right)^{m}}{m!}}}\right)\frac{\partial f}{\partial c},$$

$$\frac{\mathrm{d}\sigma_{c}^{2}}{\mathrm{d}t}=2\sum_{m=0}^{\infty}{R^{(m)}(c_{0})\frac{\langle\mathopen{}\mathclose{{}\left(c-c_{0}}\right)^{m+1}\rangle}{m!}}\\ +2\frac{k_{\text{B}}T}{N_{t}}\sum_{m=0}^{\infty}{k_{0}^{(m)}\frac{m+1}{m!}\langle\mathopen{}\mathclose{{}\left(c-c_{0}}\right)^{m}\rangle}.$$

$$\frac{\mathrm{d}{\sigma}^{2}_{c}}{\mathrm{d}t}=2\mathopen{}\mathclose{{}\left(\mathopen{}\mathclose{{}\left.\frac{\partial R}{\partial c}}\right|_{c_{0}}\sigma^{2}_{c}+\frac{k_{\text{B}}T}{N_{t}}k_{0}}\right),$$

$$f^{\text{eq}}(c)\approx e^{-\beta N_{t}\mu^{\prime}(c_{0})(c-c_{0})^{2}/2}e^{-\beta N_{t}\mathopen{}\mathclose{{}\left(a(c_{0})-\mu(c_{0})c_{0}-\Psi_{m}^{\text{eq}}}\right)},$$

$$\sigma^{2}_{c}(t)=\sigma^{2}_{c}(t=0)\exp{\mathopen{}\mathclose{{}\left(2\int_{0}^{t}{\frac{\partial R}{\partial c}\,\mathrm{d}t}}\right)}.$$

$$f(c\mathopen{}\mathclose{{}\left(t,c_{0}}\right),t)=f_{0}(c_{0})\exp{\mathopen{}\mathclose{{}\left(\int_{0}^{t}{-\frac{\partial R}{\partial c}\mathopen{}\mathclose{{}\left(c\mathopen{}\mathclose{{}\left(t,c_{0}}\right),t}\right)dt}}\right)},$$

$$f(c\mathopen{}\mathclose{{}\left(t,c_{0}}\right),t)=f_{0}(c_{0})\frac{R(c_{0})}{R(c\mathopen{}\mathclose{{}\left(t,c_{0}}\right),t)}.$$

$$R_{\text{total}}=\frac{\mathrm{d}⟨c⟩}{\mathrm{d}t}={\int }_0^{1}c\frac{\partial f}{\partial t}\mathrm{d}c={\int }_0^{1}j\mathrm{d}c.$$

Full text

I Introduction

Population dynamics has been studied widely in many particulate systems, such as crystallization, aerosol dynamics, emulsion, cell culture, etc Ramkrishna (2000); Ramkrishna and Singh (2014). A common feature of these systems is the existence of intrinsic instability, either thermodynamic or kinetic in origin, that drives the evolution of the probability distribution. For example, in a multi-phase system, the path toward minimizing interfacial energy causes the mean and the variance of the size of the dispersed phase to increase in time, a well-studied phenomenon called Ostwald ripeningLifshitz and Slyozov (1961). Kinetically driven examples of population dynamics include cell growth, where the autocatalytic cell division rate can be connected to the population growth rate Jafarpour et al. (2018), as well as bistability observed in gene regulatory processesShu et al. (2011).

Our work is motivated by such nonlinear dynamics in the context of an ensemble of particles that exhibit bistable reaction kinetics, in particular electrochemical intercalation reactions. Dreyer et al.Dreyer et al. (2010, 2011a, 2011b, 2015); Guhlke et al. (2018) used Fokker-Planck equation to describe the probability distribution of the state of particles in many-particle storage systems such as lithium ion batteries and interconnected rubber balloons, in particular the phase transition due to regions of instability where the chemical potential decreases with lithium concentration and pressure decreases with volume. The model is capable of predicting voltage hysteresis between slow charge and discharge. Other studiesMielke and Truskinovsky (2012); Dreyer et al. (2011b) have used similar equations to model mechanical systems as an ensemble of bistable units. Extensive study has also been done on the bistable kinetics of CO oxidation Bonnefont et al. (2017); Crespo-Yapur et al. (2013); Johanek et al. (2004), which shows S-shaped negative differential resistance in the current-voltage relationship as well as bistability of surface coverage on Pd nanoparticles. Of particular interest is the population dynamics of coupled microelectrodes Crespo-Yapur et al. (2014), where sequential activation of CO oxidation upon current ramping is observed.

Emergent behavior arises in population dynamics with nonlinear dynamics. Herrmann et al. Herrmann et al. (2012, 2014) studied phase transition based on Dreyer’s model in the case of constant total current. The coupling of particles through the constraint of total current gives rise to a variety of regimes ranging from Kramers type transition, oscillatory bimodal distribution, to unimodal distribution with increasing total current. Other examples of emergent behavior have also been observed through nonlocal interaction via moments of distributionFrank (2005) or globally coupling mean field parameters such as those in the synchronization of Kuramoto coupled oscillatorsStrogatz (2000). Concerted behaviors observed in crowd and biological synchrony Mirollo and Strogatz (1990); Strogatz et al. (2005) and recently in Li-ion batteriesLi et al. (2018) motivate further theoretical modeling.

The reaction kinetics for lithium intercalation in models above is linear with respect to the thermodynamic driving force, lending it the same kinetic stability conditions as its thermodynamic stability. However, asymmetry between charge and discharge of lithium ion battery is observed in particle ensemble experiments such as chronoamperometryBai and Tian (2013), as well as single-particle imaging experiments Li et al. (2014); Lim et al. (2016), which have shown that the level of Li concentration heterogeneity is asymmetric between charge and discharge. Such asymmetry cannot be captured by a symmetric linear kinetics.

BazantBazant (2017) extended the stability of open driven systems to reactions with general reaction kinetics under a thermodynamically consistent framework. The kinetic stability landscape of the system can be altered by an explicitly state-dependent autocatalytic or autoinhibitory reaction kinetics, showing the suppression of thermodynamic instability in one direction and enhancement in the other, or the creation of instability in thermodynamically stable systems. For a system driven far from equilibrium, we demonstrate that the interplay between reaction kinetics and thermodynamics results in a rich collection of phenomenon observed not only for a single system, but also for a reaction-controlled particle ensemble. State-dependent kinetics has also been explored in velocity dependent or space dependent friction chemotaxis and Van der Pol oscillator Schweitzer (2007), which leads to undamped limit cycle motion or several stable states in the phase plane. There have also been studies on temporally and spatially dependent friction in activated rate processes using generalized Langevin equation and Fokker-Planck equation Haynes and Voth (1995); Haynes et al. (1993); Voth (1992); Straus et al. (1993); Schell et al. (1981); Pollak and Berezhkovskii (1993).

The state-dependent reaction kinetics has also been implemented in the population dynamics in realistic electrochemical systems such as lithium iron phosphate (LFP) and graphite, pioneered by Multiphase Porous Electrode Theory (MPET) Smith and Bazant (2017); Burch and Bazant (2009); Ferguson and Bazant (2014, 2012). MPET simulates discrete particles of random sizes in a porous electrode and resolves the spatial dependence including the phase separation within the particles through Allen-Cahn/Cahn-Hilliard reaction modelBazant (2013). It also incorporates transport equations in the electrolyte. For phase separating materials, at low current, discrete particles are filled stochastically, a phenomenon known as mosaic instability. The results of MPET model have been validated experimentally in a direct observation of lithium concentration in an ensemble of LFP nanoparticlesLi et al. (2014). A statistical model has also been applied to explain the chronoamperometric response of LFP electrodesBai and Tian (2013). However, only the fractions of untransformed, active and transformed particles are considered, instead of a continuous probability distribution. MPET includes all relevant physics and has shown predictive power in realistic systems but lacks an accurate statistical description as it is limited by the number of particles used in the simulation. Our theory eliminates any spatial dependence and is capable of solving for the probability distribution directly in the reaction-controlled limit.

In section II, we lay out the theory for population dynamics, starting from its mathematical description in terms of the probability distribution in section II.1. Based on the thermodynamics of an ensemble of particles in section II.2, we construct a thermodynamically consistent theory for population dynamics based on generalized Fokker-Planck equation and Langevin equation in section II.3. Section II.5 explores a variety of reaction kinetics and studies how their particular functional form impacts the dynamics analytically. In the end, section III applies the theory to LFP, for which experimental data are available, and to a model system that is at the critical point.

II Theory

II.1 From single particle dynamics to population dynamics

We begin by considering a particle that undergoes reactions when it is out of chemical equilibrium with its surrounding,

[TABLE]

where the single component of interest M reacts with a reservoir with species $\text{M}_{\text{res}}$. The dynamics is described by the change in the concentration of M in the particle,

[TABLE]

where $R$ is the reaction rate as a function of concentration $c$, chemical potential of M in the particle $\mu$ and of $\text{M}_{\text{res}}$ in the reservoir $\mu_{\text{res}}$. The reaction rate is zero when the chemical reaction is at equilibrium $\mu=\mu_{\text{res}}$.

Following Bazant Bazant (2017), we define the total autocatalytic rate of the reaction to be

[TABLE]

which determines whether the concentration is linearly stable ($s<0$) or unstable ($s>0$) to infinitessimal perturbations. (See Ref.  Bazant (2017) for a general stability theorem for multicomponent driven reactive mixtures, in terms of the chemical diffusivity and autocatalytic rate tensors.) In most classical theories, the reaction rate is linear in the thermodynamic driving force, or affinity of the reaction Kondepudi and Prigogine (2014). More generally, if the reaction rate is linearized near equilibrium, $R=R_{0}\Delta\mu$, for small thermodynamic driving force, $\Delta\mu=\mu_{\text{res}}-\mu$ with constant exchange rate prefactor, $R_{0}$, the autocatalytic rate is given by $s=-R_{0}\,\mathrm{d}\mu/\mathrm{d}c$. In this case of near-equilibrium reactions, we arrive at the Duhem-Jougeut Theorem Kondepudi and Prigogine (2014), which states that kinetic stability ($s<0$) is equivalent to thermodynamic stability ($\mathrm{d}\mu/\mathrm{d}c>0$).

Far from equilibrium, the situation is very different, and kinetic stability is altered by the explicit dependence of the reaction rate on concentration or chemical potential Bazant (2017). The concentration profile tends to be destabilized by a positive “solo-autocatalytic rate”, \mathopen{}\mathclose{{}\left(\partial R/\partial c}\right)_{\mu}>0 or, in the case of a thermodynamically stable system, $\mathrm{d}\mu/\mathrm{d}c>0$, by a negative“differential reaction resistance” \mathopen{}\mathclose{{}\left(\partial R/\partial\mu}\right)_{c}<0 (as in the case of Marcus kinetics for electron transfer reactions Bazant (2013)). The stability of the system is hence determined by both thermodynamics and reaction kinetics. The interplay between the two is the source of inspiration for the following discussion.

The analysis of kinetic stability described here is applicable to lithium intercalation reaction in LFP, a phase-separating electrode material. In particular, using a thermodynamically consistent reaction kinetics, Bazant and coworkers Bai et al. (2011); Cogswell and Bazant (2012); Bazant (2013) showed a linearly stable regime in the spinodal region at large enough reaction rate, leading to the suppression of phase separation. The reaction kinetics has been coupled with phase field and mechanical models to predict the spatial pattern of lithium concentration within nanoparticles Bai et al. (2011); Cogswell and Bazant (2012, 2013); Smith and Bazant (2017); Nadkarni et al. (2018). In the Allen-Cahn reaction (ACR) model Bazant (2013), the Cahn-Hilliard term $-\kappa\nabla^{2}c$ is added to the chemical potential $\mu$ to penalize concentration gradient and allow for a diffuse interface. The assumption behind the ACR model that homogeneous reaction takes place in a particle is shown to be valid when surface reaction dominates while the concentration variation in the direction normal to the surface is negligible.

In the limit of vanishing spatial gradient within particles, the models above for lithium intercalation reduce to Eq. 2. In the limit of infinite number of reaction-controlled particles that are small enough to assume intraparticle homogeneity, we switch from a Lagrangian description of the discrete particle dynamics to an Eulerian description in terms of a probability distribution,

[TABLE]

where $f(t,c)$ is the probability distribution function (PDF).

The characteristics of the linear advection equation in the $c-t$ space is given by the reaction kinetics in Eq. 2 Ramkrishna (2000). When the reaction is autocatalytic, $\partial R/\partial c>0$, characteristics diverge and the distribution expands, whereas when $\partial R/\partial c<0$, characteristics converge and the distribution forms a shock. Fig. 1 illustrates the effect of autocatalysis on the evolution of $f(t,c)$ for two cases where a unimodal distribution evolves into a bimodal distribution over time. The evolution of $f(t,c)$ follows the characteristics, or the gray curves, which are solutions to Eq. 2. Due to autocatalysis, particles ahead of others in the reaction direction accelerate and form another peak in the PDF around a concentration removed from the original peak. The formation of peaks in the distribution is analogous to compression shock waves, while broadening distribution is analogous to rarefaction shock waves studied extensively in convection equations. The following sections establish the foundation of this mathematical description in thermodynamics and kinetics.

II.2 Thermodynamics of a particle ensemble

The model above assumes deterministic reaction rate. Now we seek a statistical mechanical description for an ensemble of particles, where thermodynamic fluctuations are important. For simplicity, we consider single-component systems. In the canonical ensemble, the number of M in the system is fixed to be $N$. In the grand canonical ensemble, the system is in a chemical reservoir of $\text{M}_{\text{res}}$ with chemical potential $\mu_{\text{res}}$. In table 1 we define thermodynamic quantities in canonical and grand canonical ensembles as a function of the probability distribution Chandler (1987). $p_{\nu}$ is the probability of a particular microstate $\nu$, $p_{\nu_{N}}=p(\nu|N)$ is the conditional probability of a microstate given $N$ molecules, and $p_{N}$ is the marginal probability of having $N$ molecules. The probability distribution can be generally in or out of equilibrium. For the canonical ensemble, the subscripts $N$ in the free energy $A_{N}$, average energy $\langle{E}\rangle_{N}$ and entropy $S_{N}$ denote the number of molecules specified for the ensemble. We are particularly interested in how the thermodynamic quantities in the grand canonical ensemble are associated with those in the canonical ensemble. For example, using the property of conditional probability $p_{\nu}=p_{\nu_{N}}p_{N}$, the grand canonical entropy $S$, the canonical entropy $S_{N}$ and the entropy due to mole fluctuations $S^{\prime}$ are related by

[TABLE]

where $\langle{S_{N}}\rangle=\sum_{N}{p_{N}S_{N}}=-k_{\text{B}}\sum_{\nu}{p_{\nu}\ln{p_{\nu_{N}}}}$ is the canonical entropy averaged over $N$, and is also known as conditional entropy in information theory.

The free energies of the two ensembles are also related as shown in table 1. In particular, if all degrees of freedom except for $N$ is at equilibrium, the grand canonical ensemble free energy is

[TABLE]

This may describe the free energy of an ensemble of particles that have reached equilibrium within themselves but not necessarily with the reservoir.

For convenience, we use $A=A_{N}^{\text{eq}}$ to denote the canonical ensemble free energy at equilibrium. From thermodynamics, the system is stable when, $\partial^{2}A/\partial N^{2}>0$, if $N$ is treated as a continuous variable. An example of instability, which leads to phase separation, can arise in the regular solution model for lattice fluid discussed in section II.5. Since the chemical potential of $M$ within the particle is defined to be $\mu=\partial A/\partial N$, the stability criterion can also be written as $\partial\mu/\partial N>0$. This is not to be confused with the fact that, at equilibrium, $\partial^{2}\Psi^{\text{eq}}/\partial\mu_{\text{res}}^{2}=-\partial\langle{N}\rangle^{\text{eq}}/\partial\mu_{\text{res}}=-k_{\text{B}}T\sigma_{N}^{2}<0$, where $\sigma_{N}^{2}$ is the variance of $N$. Therefore, at equilibrium, the average number of M in an grand canonical ensemble always increases with reservoir chemical potential regardless of $A$. Note that $\sigma_{N}^{2}$ at equilibrium is constructed naturally as a result of the entropy $S^{\prime}$ associated with mole fluctuations.

At equilibrium, the probability distribution satisfies the Boltzmann distribution and minimizes the free energy. Therefore we may define a Lyapunov function to be free energy $\Psi$ minus its minimum $\Psi^{\text{eq}}$, which is always greater than or equal to zero and, as will be shown later, monotonically decreases for systems with a constant $\mu_{\text{res}}$. The Lyapunov function is called relative entropy, or Kullback-Leibler divergence, which originates in information theory to measure the difference between two probability distributions.

[TABLE]

II.3 Dynamics

In the classical limit, we use a general continuous variable $\bf{x}$ that maps the microstate of the system $\nu$ to $\mathbb{R}^{n}$ and use $p$ to denote the probability density function. In the grand canonical ensemble, the driving force is the gradient of the potential in the phase space. The potential is defined by the variational derivative $\delta\Psi/\delta p=E-\mu_{\text{res}}N+k_{\text{B}}T\ln{p}$. Note that equilibrium corresponds to $\delta\Psi/\delta p$ being a constant that satisfies the normalization condition for $p$. Assuming the flux in the phase space is proportional to the driving force, we generalize the Fokker-Planck equation, which can be regarded as Wasserstein gradient flowRisken (1989); Herrmann et al. (2012); Jordan et al. (1998),

[TABLE]

where

[TABLE]

where $\bf{K}$ is the mobility tensor.

The following second law of thermodynamics is well knownSekimoto (2010); Reguera et al. (2005); Rao and Esposito (2016); Seifert (2005, 2012); Van den Broeck and Esposito (2010); Frank (2005),

[TABLE]

where the first term is the entropy production rate $-T\dot{S}_{i}$ and the second term is the work of external driving force $\dot{W}_{d}$. When the system is driven chemically,

[TABLE]

The second law of thermodynamics requires that $\dot{S}_{i}>0$, which means $\mathbf{K}$ is positive definite. When the system is nondriven, $\Psi-\Psi^{\text{eq}}$ decreases monotonically, hence it is a Lyapunov function.

If the relaxation to equilibrium in terms of its internal degrees of freedom $\nu_{N}$ is much faster than the exchange of matter with the reservoir, which may take place via chemical reactions, we may assume all internal degrees of freedom are at equilibrium, that is, $p_{\nu}=p_{\nu_{N}}^{\text{eq}}p_{N}$. In this case, we can show that the marginal probability in terms of $N$ can be described by the following equation,

[TABLE]

where

[TABLE]

It is convenient to nondimensionalize the number of molecules $N$ in the system by the its maximum $N_{t}$, which is well defined for intercalation systems. Assuming that $N_{t}$ is fixed, the equation can be nondimensionalized via the mole fraction, $c=N/N_{t}$, the PDF in terms of $c$, $f=P_{N}N_{t}$, and prefactor $k={(\mathbf{K})_{N,N}}/N_{t}$. In contrast to linear kinetics, $k$ can be a function of $c$, $\mu$, and $\mu_{\text{res}}$, and we do not make any assumption about the functional form other than $k>0$.

[TABLE]

where $\mu=\partial A_{N}^{\text{eq}}/\partial N$. Alternatively,

[TABLE]

where the flux $j$ is,

[TABLE]

where $f^{\text{eq}}$ is the equilibrium distribution, given by Eq. 18. Defining the maximum number of molecules $N_{t}$ restricts the fraction $c\in[0,1]$. Therefore we impose the no-flux boundary condition $j=0$ at $c=0$ and 1. To be compatible, we require the chemical potential $\mu\to\infty$ as $c\to 1$ and $\mu\to-\infty$ as $c\to 0$ and is undefined outside (0,1).

Following the nondimensionalization, we define the following free energy in the macroscopic sense

[TABLE]

where $a=A/N_{t}$. Following the analysis of the general case (Eq. 8), $\Psi_{m}$ is the Lyapunov function for Eq. 15. Therefore if not chemically driven, or at constant $\mu_{\text{res}}$, the system is Lyapunov stable. In fact, it can also be shown that the eigenvalues of the operator \mathopen{}\mathclose{{}\left(\mathopen{}\mathclose{{}\left(\mu_{\text{res}}-\mu}\right)k-kk_{\text{B}}TN_{t}^{-1}\,\partial/\partial c}\right) are nonnegative, while the eigenfunction corresponding to the zero eigenvalue is the equilibrium distribution Risken (1989). In other words, the equilibrium distribution corresponds to $\partial j/\partial c=0$, or $j=0$, based on the boundary condition,

[TABLE]

where $\Psi^{\text{eq}}_{m}$ is the potential defined in Eq. 17 when at equilibrium and is the constant that normalizes the distribution.

By comparing Eqs. 4 and 15, we find that the reaction rate is R=\mathopen{}\mathclose{{}\left(\mu_{\text{res}}-\mu}\right)k. In the limit of large systems, $N_{t}\to\infty$, Fokker-Planck equation (Eq. 15) reduces to Eq. 4. Now we are able to correlate mole fluctuations with the thermodynamic quantities in the diffusive term in the Fokker-Planck equation and extend the deterministic dynamics, or Eq. 2, to the following Langevin equation, which is equivalent to the Fokker-Planck equation.

[TABLE]

where the Langevin noise $\xi(t)$ satisfies zero average $\langle\xi(t)\rangle=0$ and is uncorrelated in time $\langle\xi(t)\xi(t^{\prime})\rangle=\delta(t-t^{\prime})$. Eq. 15 can also be written as,

[TABLE]

where $D=k_{\text{B}}Tk/N_{t}$ is an expression of the Einstein relation or Fluctuation-Dissipation Theorem for concentration fluctuations across the particle ensemble. The noise intensity decreases as the size of the particle increases.

In section II.1, we mentioned that when $\partial R/\partial c>0$, characteristics diverge and the probability distribution expands in time, and vice versa. This statement can be made more quantitative in the case of non-negligible fluctuations by considering the evolution of the variance of the distribution. Herrmann Herrmann et al. (2012) analyzed the widening of the distribution in the spinodal region for a short time by linearizing the reaction $R$ with respect to $c$, which specifically entails linearizing the chemical potential $\mu$ since the prefactor $k$ is assumed to be constant. Therefore, in this case, the widening of the distribution is determined by $\partial\mu/\partial c$. The analysis requires that the distribution be sharply peaked.

II.4 Evolution of the concentration variance

With the same assumption, we extend the analysis to a general reaction rate expression and approximate the evolution of the concentration variance, $\sigma_{c}^{2}=\int_{0}^{1}{(c-c_{0})^{2}f(c)\,\mathrm{d}c}$. The distribution $f(t,c)$ satisfies the following Fokker-Planck equation

[TABLE]

where $j=fR-k_{\text{B}}TN_{t}^{-1}k\partial f/\partial c$. Denote the mean $c_{0}=\langle c\rangle$, and let us calculate the evolution of the concentration variance

[TABLE]

Expand the current $j$ about $c_{0}$ and use superscript $(m)$ to denote $m^{\text{th}}$ derivative with respect to $c$,

[TABLE]

Therefore

[TABLE]

If $f$ is concentrated around $c_{0}$, higher order moments in the expression above decay very quickly and are negligible. Therefore we retain only moments up to $\langle(c-c_{0})^{2}\rangle=\sigma_{c}^{2}$. In the limit of sharply-peaked distribution, $k_{\text{B}}T/N_{t}$ and $\sigma_{c}^{2}$ are both small, we arrive at the key result,

[TABLE]

where $k_{0}=k(c_{0})$, and $c_{0}=\langle c\rangle$ is the average mole fraction.

For autoinhibitory reactions, $\partial R/\partial c<0$, the variance tends to decrease in time while fluctuations work against it, and the particle ensemble stays homogeneous. For autocatalytic reactions, the distribution widens and the particle ensemble becomes more heterogeneous. If the equilibrium distribution is unimodal with a sharp peak, it can be approximated by a Gaussian function by expanding the free energy $a$ in terms of $c$ to second order,

[TABLE]

where, $\mu_{\text{res}}=\mu(c_{0})$ and \sigma_{c}^{2}=k_{\text{B}}T/\mathopen{}\mathclose{{}\left(N_{t}\mu^{\prime}(c_{0})}\right). The same conclusion is obtained by setting Eq. 25 to zero, since \partial R/\partial c=-k_{0}\partial\mu/\partial c+\mathopen{}\mathclose{{}\left(\mu_{\text{res}}-\mu}\right)\partial k/\partial c. When the system does not phase separate thermodynamically ($\partial\mu/\partial c>0$), the variance can become greater or smaller than the equilibrium variance when there is a net reaction, depending on the sign of $\partial k/\partial c$ and the direction of the net reaction. It is the focus of this paper to explore the effect of $k(c)$ on population dynamics. The next section describes how different $k(c)$ models predict dramatically different variance evolution.

In the limit of negligible fluctuations, $N_{t}\to\infty$, Eq. 25 gives the evolution of variance over time,

[TABLE]

In the same limit, with the initial condition $f(c_{0},t=0)=f_{0}(c_{0})$, using method of characteristics, the solution Ramkrishna (2000) to the Fokker-Planck equation without the diffusive term (Eq. 4) is

[TABLE]

where c\mathopen{}\mathclose{{}\left(t,c_{0}}\right) follows the characteristics and is the solution to Eq. 2 with the initial condition $c(t)=c_{0}$. We see that in the limit of sharply-peaked distribution, by taking the variance on both sides of Eq. 28, Eq. 27 is recovered. If the reaction rate does not explicitly depend on $t$, then

[TABLE]

II.5 Electrochemical reactions

Electrochemical systems offer the most physically realizable benchmark for the population dynamics model, since in experiments the control of current or voltage corresponds to the control of the total reaction rate or the chemical potential of the reservoir $\mu_{\text{res}}$, respectively. In this section, we focus on the effect of reaction kinetics in the case of constant total reaction rate – one of the most commonly used protocol in electrochemistry. The total reaction rate is given by

[TABLE]

In the limit of a sharply-peaked distribution, the center of the peak $c_{0}$ is equal to the average concentration $\langle{c}\rangle$ and shifts linearly in time. The total reaction rate is equal to the reaction rate at $c_{0}$, that is, $\mathrm{d}c_{0}/\mathrm{d}t=R(c_{0})$. Therefore, the evolution of variance in Eq. 25 can also be written as,

[TABLE]

In the limit of negligible fluctuation or infinitely large reaction rate, when the thermodynamic driving force $\mu_{\text{res}}-\mu$ is large, the growth of the standard deviation of the distribution is given by

[TABLE]

In this section, we consider the following simple reaction kinetics,

[TABLE]

where $g(0)=0$, which includes important cases of electrochemical reactions Bazant (2013, 2017). The corresponding autocatalytic rate is

[TABLE]

With this particular form of reaction kinetics, we have

[TABLE]

Here we focus on the explicit dependence of reaction rate on concentration, in particular how the exchange current, $R_{0}(c)$, affects the stability of the particle ensemble by electroautocatalysis Bazant (2017). Hence we discuss the case of positive differential resistance $g^{\prime}>0$, that is, the reaction rate increases with the driving force, which preserves the contribution of the thermodynamic stability to the overall stability. Although not the focus of the paper, Marcus theory for outer-sphere electron transfer, for example, exhibits an inverted region, or negative differential resistanceBazant (2013, 2017), which can destabilize a thermodynamically stable system and make the reaction autocatalytic. Reaction kinetics of this form can be derived from transition state theory. For example, the famous Butler-Volmer equation has the form aboveBazant (2013); Newman and Thomas-Alyea (2012), where $R_{0}$ is called the exchange current and

[TABLE]

where $\alpha$ is the charge transfer coefficient. In the case of symmetric Butler-Volmer kinetics, $\alpha=0.5$, g^{\prime}=\sqrt{\mathopen{}\mathclose{{}\left(g/2}\right)^{2}+1}. The standard deviation grows in regions where $s>0$, or \mathopen{}\mathclose{{}\left(\ln{R_{0}}}\right)^{\prime}>\sqrt{1/4+\mathopen{}\mathclose{{}\left(R_{0}/R}\right)^{2}}\mu^{\prime}\text{sgn}(R), as shown by Bazant in Ref. Bazant (2017). With Butler-Volmer kinetics, when the driving force $\mu_{\text{res}}-\mu$ is small, $g(x)\to x$, and the reaction rate is slow. Therefore the second term in Eq. 32 dominates and whether the standard deviation increases or decreases is mainly determined by the thermodynamic stability,

[TABLE]

At high total reaction rate, both the exchange current and thermodynamics determine the evolution of standard deviation.

[TABLE]

The exchange current of the reaction kinetics is crucial in determining the stability of the systemBazant (2017); Lim et al. (2016). The asymmetry in exchange current leads to asymmetric behavior between forward and backward reactions. In the case of LFP, a skewed exchange current results in more uniform concentration profiles during intercalation while greater inhomogeneity during de-intercalation. To illustrate how the exchange current alters the kinetic stability, we choose regular solution as the thermodynamic model.

[TABLE]

where the first term arises from the entropy of mixing and the second term comes from enthalpic interaction. When $\Omega>2$, there exists a spinodal gap and the system phase separates into two phases with low and high concentrations.

Table 2 lists, in the case of symmetric Butler-Volmer kinetics ($\alpha=0.5$), the maximum standard deviation ratio at infinitely fast reaction rate for three types of exchange current, which is attained in the interval $[c_{1},c_{2}]$ in which $\partial R/\partial c\geq 0$. The starting and ending concentrations $c_{1}$ and $c_{2}$ are also listed in the table. See section III.1 below for a plot and discussion on the region of instability $\partial R/\partial c\geq 0$ at arbitrary reaction rates. The first two exchange currents are symmetric and thus exhibit symmetric behavior for forward and backward reactions. Constant exchange current does not alter the thermodynamic stability – the kinetic instability ($s>0$) grows exactly within the spinodal region, which exists when $\Omega>2$. Fig. 2 shows that evolution of the standard deviation at infinitely large reaction rate. The third row of the figure shows that the evolution for the constant exchange current follows the trend of the chemical potential. For the second exchange current $\sqrt{c(1-c)}$, which is frequently used in electrochemistry, the kinetic instability grows within $[0,1-1/\Omega]$ for forward reaction and $[1/\Omega,1]$ for backward reaction. The instability region exists when $\Omega>1$. Since at low enough reaction rate the system always follows thermodynamic stability, when $1<\Omega<2$, it only becomes linearly unstable at high reaction rate. As shown in the second row of figure 2, the second exchange current results in a larger maximum standard deviation ratio as well as region of instability than the constant exchange current.

As Bazant has shown Bazant (2017), the third exchange current is predominantly autocatalytic in the backward direction due to the asymmetrically increasing exchange current in that direction. When the direction is reversed, the autoinhibition results in greater homogeneity for the forward reaction. In particular, at infinitely large backward reaction rate, the system is kinetically unstable within the concentration range $[1/2\Omega,1]$, which exists when $\Omega>1/2$. At infinitely large forward reaction rate, the system remains stable regardless of the value of $\Omega$, behaving entirely differently from the symmetric exchange currents. Fig. 2 shows that, in the backward direction, the amplification of standard deviation and the length of the instability interval increases in the order of constant, symmetric $\sqrt{c(1-c)}$ and asymmetric $(1-c)e^{\mu/2}$ exchange current, the reason for the latter being 1/2\Omega<\mathopen{}\mathclose{{}\left(1-\sqrt{1-2/\Omega}}\right)/2<1/\Omega. Fig. 2 also shows the agreement between the analytical expression for the evolution of standard deviation from Eq. 38 and numerical simulation for the same exchange currents. The approximation captures the magnitude of maximum standard deviation and where maximum standard deviation is attained.

When the variance $\sigma_{c}^{2}$ gets large, the assumption of sharply-peaked distribution no longer holds. For example, when the system is driven toward high concentration, autocatalytic effect may cause some particles to accelerate to much higher concentration than average. However, the highest concentration a particle can reach is the equilibrium concentration corresponding to $\mu_{\text{res}}$. The stabilizing effect of thermodynamics when the concentration is near equilibrium ($\partial R/\partial c>0$) causes those particles to form another peak in the distribution. We define the resulting bimodal distribution to be phase separation at the ensemble level.

The key result of this work is to show that phase separation can arise as a result of both kinetic and thermodynamic effects, and mathematically, it is a result of diverging and converging characteristics. Based on the analysis, we may approximate the initiation of phase separation as when the estimated variance $\sigma_{c}=\sigma_{c_{0}}e^{K}$ exceeds a certain threshold, where $\sigma_{c_{0}}$ can be taken at the inception of linear instability or the initial condition, when, as explained previously, the magnitude of standard deviation is proportional to the strength of the random fluctuation, that is, $\sigma_{c_{0}}\propto\sqrt{D}$. In the case of linear kinetics with constant exchange current, that is, R=R_{0}\mathopen{}\mathclose{{}\left(\mu_{\text{res}}-\mu}\right), where $R_{0}$ is a constant, we have $K=(\mu(c)-\mu(c_{0}))R_{0}/R$. Thus, we obtain the following critical current $R_{c}$ above which a thermodynamically phase separating system is kinetically stabilized, consistent with Herrmann’s conclusion Herrmann et al. (2012).

[TABLE]

III Applications

III.1 Li-ion Battery Porous electrodes

We first look at the example of lithium iron phosphate (LFP,) the canonical and most studied phase-separating solid-state lithium intercalation material. The lithiation reaction is $\text{Li}^{+}+\text{e}^{-}+\text{FePO}_{4}\to\text{LiFePO}_{4}$ and the reverse reaction is called delithiation. The chemical potential of lithium in the solid can be modeled by the regular solution model given by Eq. 39 as a function of the lithium fraction in the solid, with $c=0$ being $\text{FePO}_{4}$ and $c=1$ being $\text{LiFePO}_{4}$. We model the population dynamics using (a) a classical exchange current, \sqrt{\mathopen{}\mathclose{{}\left(c(1-c)}\right)}, (b) a thermodynamically consistent exchange current derived from transition state theory ($(1-c)e^{\mu/2}$) and (c) an experimentally measured exchange current Lim et al. (2016), all with Butler-Volmer overpotential dependence.

In this section, we focus on the population dynamics of an ensemble of LFP particles under constant current. Fig. 3 shows the dynamic phase diagram of the LFP model as a function of current and average Li fraction and the evolution of the probability distribution function in time at select currents. With a positive current (reaction rate), the average Li fraction increases, which corresponds to going horizontally from left to right in the main figure, and vice versa for a negative current. In the previous section, we define phase separation to be when the PDF is bimodal. This is shown by the gray area in the main plot. The color indicates the prominence of the second largest peak, which qualitatively measures the degree of phase separation. To clearly delineate the phase separation region, the color map is at saturation when the value is 1. The point of exit from phase separation is oscillatory with respect to total reaction rate. This is due to the varying frequency observed in the regime of oscillatory phase transition, which is addressed extensively previously and is beyond the scope of this workHerrmann et al. (2012, 2014); Dreyer et al. (2011b).

At low currents, regardless of the exchange current, thermodynamics drives the particle ensemble to separate into Li poor and Li rich populations, which correspond to growing and shrinking peaks seen at $i/i_{0}=10^{-3}$. As shown in the dynamic phase diagram, at low currents (but not $R\to 0$, as explained below), the phase separation initiates at around the nearest spinodal point.

As mentioned above, Herrmann et al. Herrmann et al. (2012, 2014); Dreyer et al. (2011b) performed extensive analysis on the various regimes of the population dynamics at different currents for the simple kinetics $R=R_{0}\Delta\mu$. Our more general model approaches it asymptotically as $R\to 0$. In the Kramer’s regime, the reaction rate is proportional to \exp{\mathopen{}\mathclose{{}\left(-b/D}\right)}, where $b$ is the free energy barrier, or the difference between the local maximum and minimum of $a-\mu_{\text{res}}c$. Since the reaction rate is equal to the rate of phase transformation in this regime, as $b/D$ increases, the rate of nucleation decreases exponentially. To reach the quasistationary limit, that is, for the phase separating region to span between the two thermodynamic phases (binodal points), the rate must be smaller than the Kramer’s rate at the largest possible free energy barrier ($\mu_{\text{res}}=0$), which is around $10^{-422}$ for the particular set of parameters chosen ($\Omega=3.4$, $D_{0}=2\times 10^{-4}$), therefore the dynamic phase diagram shown is effectively only for current at or above the Kramer’s regime, and the system must at least reach the spinodal point for the nucleation to occur.

The solid curves denote the boundary of linear instability, $s=0$. In between the two curves, $s>0$. At zero current, the region of linear instability is also the spinodal region. The region of phase separation does not follow the shape of the region of linear instability. In fact, nucleation is delayed with increasing current. We seek to quantitatively understand this trend by approximating phase separation to be when the support of the distribution, estimated by the mean plus or minus one standard deviation $\sigma_{c}=\sigma_{c_{0}}e^{K}$, is contained in $[0,1]$, because when phase separation occurs, the new phase is close to $c=1$ for $R>0$ and close to $c=0$ for $R<0$. In other words, we approximate single phase with the following condition,

[TABLE]

For a qualitative understanding, we choose an order-of-magnitude estimate for the initial variance, $\sigma_{c_{0}}=\sqrt{D}$. $c_{0}$ is the first point where $s=0$ or the initial average concentration used in the simulation (0.05 or 0.95), whichever comes later in the direction of reaction. The predicted region of phase separation is bounded by the dashed lines. This simple approximation captures the general feature of the boundary of phase separation. The form of reaction kinetics, in particular the autocatalytic rate $\partial R/\partial c$, is the dominant factor in determining the phase separating region. It also indicates that the delayed nucleation is mainly due to decreasing $s/R$ with increasing reaction rate.

The regimes studied by Herrmann at higher current no longer hold for certain exchange currents. For the simple kinetics, the distribution stays unimodal above a certain current (around 0.1). The qualitative trend remains true for the exchange current (a). And since exchange current (a) is symmetric with respect to $c=0.5$, the phase behavior is also symmetric between lithiation and delithiation. However, the asymmetric exchange currents (b) and (c) induce autocatalytic effect during delithiation and cause a larger degree of phase separation than lithiation while still suppressing the phase separation above a certain lithiation current. In fact, both (b) and (c) predict that the phase separation persists at large negative current. This demonstrates a new emergent phenomenon that autocatalytic effects add to the population dynamics.

Following Lim Lim et al. (2016), we also define uniformity coefficient, $1-\sigma_{c}/\sqrt{c(1-c)}$, which is 1 when the probability of finding particles at the average concentration is 1, $P(c=\langle{c}\rangle)=1$, and 0 when maximum $\sigma_{c}$ is attained, that is, $P(c=0)=1-\langle{c}\rangle$ and $P(c=1)=\langle{c}\rangle$. Fig. 4 shows that the experimental exchange current (c) results in increasing uniformity with increasing total current and that, due to the skewed exchange current and reasons stated above, lithiation results in more uniform distribution than delithiation. The result is in qualitative agreement with the experimental result from Lim Lim et al. (2016). The theory overpredicts the uniformity for lithiation, which could be attributed to additional noise unaccounted for, such as particle size distribution and spatial inhomogeneity of the reaction kinetics. The asymmetry is also observed in another study by Li et al. Li et al. (2014) where the active particle fraction increases as a function of current and it is higher for lithiation than for delithiation. This is consistent with the asymmetric exchange current since a smaller active particle fraction indicates a higher degree of phase separation and a lower uniformity coefficient.

III.2 Driven phase separation at the critical point

In this section, we present an example of a thermodynamically single-phase system that undergoes fictitious phase separation in certain reaction conditions. We choose the regular solution model with $\Omega=2$, so that the system is at the critical point. As mentioned earlier, \partial R/\partial c=-k_{0}\partial\mu/\partial c+\mathopen{}\mathclose{{}\left(\mu_{\text{res}}-\mu}\right)\partial k/\partial c. At the critical point $c_{0}$, $\mu^{\prime}(c_{0})=0$, \partial R/\partial c(c_{0})=\mathopen{}\mathclose{{}\left(\mu_{\text{res}}-\mu(c_{0})}\right)\partial k/\partial c(c_{0}), whether the reaction is autocatalytic or autoinhibitory is solely determined by the explicit dependence of reaction rate on concentration $k(c)$ and the direction of reaction.

For illustration purposes, we choose Butler-Volmer kinetics with the exchange current, $R_{0}=(1-c)^{2}e^{\alpha\mu}$ Bazant (2013). In a lattice model, the exponent 2 indicates that the transition state excludes 2 empty sites Bazant (2013). The quadratic dependence on $(1-c)$ enhances the autocatalysis and accelerates the backward reaction as $c$ decreases. Fig. 1 is a simulation with such thermodynamic and kinetic models, at $\mu_{\text{res}}=-0.2$ and total reaction rate $-0.1$. It illustrates the diverging characteristics during the backward reaction.

Fig. 5 shows the kinetic phase diagram of the system. Instability and bimodal distribution occur with very small negative reaction rate. Unlike LFP, increasing backward reaction rate leads to the expansion of the phase separation region. During the backward reaction, the distribution diverges so much that the population separates into two distinct “fictitious phases” while the forward reaction shows stable single-phase behavior. The region of phase separation is even more asymmetric than LFP with the asymmetric exchange current due to the second order dependence on $1-c$. The fictitious phase separation region is contained within the linearly unstable region, which is below the solid gray curve. We do not give the approximated phase separation region like the one above for LFP, because when the second peak in the PDF appears, the peak position is not always close to 0 or 1.

This leads to important applications in the interpretation of in-situ and operando X-ray diffraction (XRD) experiments of electrodes, which measures the population density as a function of lattice parameters and thus the state of charge of the particlesZhou et al. (2016); Liu et al. (2014); Singer et al. (2014); Hong et al. . One should be cautious in deriving phases of a material from the XRD spectrum of a dynamical experiment such as charge and discharge, as the XRD peaks may arise as a result of kinetic effects.

The phase separation defined for an particle ensemble can also inform us about the initiation/nucleation of phase separation within particles. In particular, the example presented here suggests that pattern formation may also arise when driven away from equilibrium through autocatalytic reactions even though the reactive mixture is thermodynamically stable. This could have implications in understanding and engineering phase separation in a variety of systems and patterns driven by electron transfer reactions. Negative differential resistance ($\partial R/\partial\Delta\mu<0$) is not studied here in the context of population dynamics but can also lead to fictitious phase separation.

IV Summary

We lay out a general theory of describing the equilibrium and dynamics of an ensemble of spatially homogeneous reactive particles that undergoes material exchange with the reservoir. The theory extends beyond the discrete particle simulation or the simplified reaction kinetics that are studied previously. Under the reaction-controlled condition, reaction kinetics determines the characteristics of the governing equation for the probability distribution of the particle ensemble. In regions of autocatalysis or autoinhibition, the probability distribution exhibits shock wave behaviors. We established the framework of introducing mole fluctuation for any generalized reaction kinetics in a thermodynamically consistent way. The corresponding Fokker-Planck equation gives the equilibrium probability distribution that is consistent with statistical physics.

Through analytical approximation and numerical analysis of several model systems, we studied how reaction kinetics of the system affects the population dynamics and in some cases transforms the entire landscape of phase behavior for systems with the same thermodynamic properties. Suppression of phase separation or the induction thereof in single-phased systems have been observed. The analytical approximation for the boundary of phase separation based on the integral of autocatalytic rate shows good agreement with numerical simulation. This reinforces our understanding that autocatalytic effects lie at the core of such phase separation behavior.

Therefore, knowledge of the reaction kinetics is crucial in understanding the dynamics of reaction-controlled systems. One can also envision engineering reactive surfaces to achieve desirable kinetic properties. The emergent phenomena observed here motivate us to further explore kinetically controlled phase separation in experiments using in-situ analysis such as the aforementioned X-ray diffraction.

Although we did not address any possible spatial correlation or dependency in this paper, the theory can be extended to incorporate the external environment that the particles are subjected to by making the probability space a Cartesian product of the particle state space and physical space. For example, we can introduce other fields such as mass, temperature and electric potential to consider transport phenomena, as well as the particle ensemble’s interaction with the external fields. Particles in a porous medium can also be influenced by its highly localized environment. Any intrinsically variable property or such local variability that follow certain probability measures can also be incorporated in the theory. The population dynamics can also shed light on the spatially dependent behavior within the particles themselves. For example, by switching from concentration field in a physical space to probability distribution, population dynamics offers a spatially-agnostic statistical description of the concentration field such as its mean and the variance for as long as it maintains isotropy and before patterns with sharp concentration gradient take over, therefore potentially serving as a precursory analysis for the inception or even termination of pattern formation.

Acknowledgements.

This work was supported by Toyota Research Institute through the D3BATT Center on Data-Driven-Design of Rechargeable Batteries.