A Comparative Review of Models for All-Solid-State Li-Ion Batteries

1. Introduction

Numerous instances demonstrate that traditional Li-ion batteries, which rely on liquid electrolytes and porous electrodes, are susceptible to chemo-mechanical degradation and face challenges related to environmental and safety concerns due to the flammability of toxic organic liquid electrolytes. Solid electrolytes are used in all-solid-state batteries (ASSB). They improve thermo-electro-chemical stability and remove the possibility of hazardous liquid electrolyte leaks.

Largely funded international projects [5] and strategic action plans ( as for the European battery 2030+ initiative ) assess how the scientific community trusts and supports the development of storage systems to meet the social quest for decarbonization. ASSBs are among the best next-generation candidates, provided that they achieve ultra-high-performances while meeting sustainability and safety. Among the shortcomings of solid electrolytes cells, the technical ones that currently limit the industrialization of ASSBs appear to be: (i) achieving a substantial ionic conductivity at ambient temperature; (ii) reducing the resistance at the interface between the active material and electrolyte [6,7]; (iii) eliminate the mechanical damage [8,9] and the dendritic growth from anodic lithium foils [10,11,12]. Experimental investigations are presently mainly focused on key materials and structures [13,14,15], such as composite electrolyte/electrode morphologies [16] or thin film glass layers [17].

Computational modeling and simulations, which operate at different scales and complement experimental research, allow a profound understanding of the interconnected physics that rule the response of ASSBs. Digital twins enable the comprehension of processes and the identification of limiting factors. Several advanced mathematical models have been documented in the literature. It is always hard, if not even impossible, to categorize the whole bibliography on a subject and undoubtedly identify those cornerstone papers that changed the flow of the subject itself. Perhaps Newman’s theory might be assigned without presumptions for the electrochemistry of batteries. As in several aspects of human beings, history is a succession of little advancements, in a collective way. In this spirit, we detail and compare here four different approaches [18,19,20,21] because of their significant ideas. Unlike previous studies [1,2,3,4], therefore, a broad screening of the existing literature is not the goal of this review.

Fabre et al. [18] developed a one-dimensional model of a Li/LiPON thin-film micro battery, described in Section 2.1. The ionic transfer is modeled via a single-ion conduction within the solid electrolyte. This results in a uniform distribution of lithium ions (${\mathrm{Li}}^{+}$) throughout the solid electrolyte, with concentrations determined by the initial conditions and established by the overall electroneutrality in the electrolyte1. In single-ion conduction, the negative vacant positions in the lattice are tightly bound and do not flow. Since ionic mobility isn’t governed by a concentration gradient, the model’s main goal is to identify the potential drop, which is primarily determined by Ohm’s law.

An advanced model framework for ASSB was further proposed by Landstorfer and coworkers in [19] and will be analyzed in Section 2.2. While using a single-ion conduction ionic transport model as in [18], the authors enriched the formulation of the interface mechanisms. Standing from a rigorous thermodynamic setting, the interface conditions ultimately lead to non-Butler-Volmer type equations. Capacitance within interfaces was captured also in the third paper accounted for in this review, authored by Raijmakers and co-workers in [20]. Their key contribution lies in an innovative two-mechanisms conduction model, wherein both interstitial lithium and negative vacancies move autonomously. This results in establishing a concentration gradient at a steady state, resembling the distribution seen in liquid electrolytes [24,25].

The double mechanism in ionic transport has been further developed in the last contribution discussed here, authored by Cabras et al. [21]. Since the unbalanced negative charges associated with a vacancy in the LIPON lattice cannot flow, the model proposed in [20] appears to be unphysical. Rather Cabras et al. [21] assume that some positive ions hop to occupy neighboring empty vacancies, while the rest of the Li+ ions transition into a meta-stable interstitial phase of the same kind as described in [20]. The complex explanation of interfaces has also been incorporated from [20].

This review paper is organized as follows. A brief theoretical summary of models and their governing equations for ASSBs is proposed in Section 2. An experimental benchmark, taken from [26], is devised in Section 3. The finite element solution scheme provides numerical approximations for the electric potential profiles, interface currents, fluxes, and concentration profiles. These quantities of interest were used to compare the more recent model in [21] against the antecedent published in [20]. In the Section 4, the article provides an overview of the objectives pursued. It also outlines the future plans that are intended to be implemented. Final remarks that summarize the key points discussed conclude the paper.

2. Models and Their Governing Equations

Different mathematical models, some of which reviewed in [3,4], have been proposed and tailored to the battery microstructure, which simplistically can be categorized into two types, i.e., thin films and porous electrodes. A broad set of models that include pseudo-2-D [27,28], multiscale [29,30] and fine-grained models [31] have been designed for electrodes made of porous materials, accounting for the different phases and attempting at capturing realism either in liquid [32] or solid electrolytes [16]. We will restrict our focus on thin films batteries, equipped with Lithium phosphorous oxy-nitride “LiPON” electrolytes.

A planar geometry is generally accepted in thin films batteries, since the ratio between the lateral dimension and the thickness is large enough for the lateral dimension to be considered as infinite, see Figure 1. Since furthermore electrodes and electrolyte can be well approximated by homogeneous planar materials, one-dimensional mathematical models are customary formulations for thin films batteries.

An all-solid-state electrochemical cell incorporates two electrodes and a solid electrolyte, as depicted schematically in Figure 1. Li-ions are extracted from the cathode (positive electrode) during charge and inserted back during discharge. The opposite holds for the anode (negative electrode). Assuming ${\mathrm{LiCoO}}_2$ (shortened in LCO) to be the positive electrode material, the basic electrochemical charge-transfer reaction writes

$$\begin{array}{c}\hfill {\mathrm{LiCoO}}_2\underset{k_{-1}}{\stackrel{k_1}{\rightleftarrows }}{\mathrm{Li}}_{1-x}{\mathrm{CoO}}_2+x{\mathrm{Li}}^{+}+x{\mathrm{e}}^{-}0\le x\le \frac{1}{2}.\end{array}$$

(1)

If lithium foil serves as the negative electrode material, deposition and extraction at the negative surface is described by the reaction

$$\begin{array}{c}\hfill \mathrm{Li}\underset{k_{-2}}{\stackrel{k_2}{\rightleftarrows }}{\mathrm{Li}}^{+}+{\mathrm{e}}^{-}.\end{array}$$

(2)

Figure 2 depicts the solid electrolyte’s structure. In it, $\mathrm{LiO}$ denotes the lithium bound to the non-bridging oxygen atoms, ${\mathrm{Li}}^{+}$ has motion capabilities (either as transferred to the meta-stable interstitial state or hopping), and ${\mathrm{n}}^{-}$ is the uncompensated negative charge associated with a vacancy formed in the LiPON matrix at the place where lithium was originally located. The maximal concentration of host-sites, denoted with $c_0$, is determined by the stoichiometric composition of the electrolyte material. It is reached in the ideal case of absolute zero temperature, when all available host sites are fully filled with lithium ions and the ionic conductivity vanishes because all ions are immobile. In standard conditions, some of the Li ions are thermally excited and the chemical ionization reaction

$$\mathrm{LiO}\underset{k_b^{\mathrm{ion}}}{\stackrel{k_f^{\mathrm{ion}}}{\rightleftarrows }}{\mathrm{Li}}^{+}+{\mathrm{n}}^{-}$$

(3)

occurs, $k_f^{\mathrm{ion}}$ and $k_b^{\mathrm{ion}}$ being the forward and backward rate constants for the ionization and recombination reaction, respectively.

During charge the lithium ions cross the electrolyte and are reduced into metallic $\mathrm{Li}$ at the lithium foil surface; vice-versa during discharge. Mobile species in the electrolyte are therefore ions and predicting the behavior of electrochemical cells requires quantitative modeling of the kinetics of mobile ionic species. ${\mathrm{LiCoO}}_2$ electrode contains lithium oxide and the lithium exists as ions as part of a lithium salt. We can thus state that $\mathrm{Li}$ intercalates in the cathode as ions as well, which will be denoted with ${\mathrm{Li}}^{\oplus }$ to stress its ionic nature “shielded” by its own electron, while distinguishing it from mobile charges ${\mathrm{Li}}^{+}$ in the electrolyte.

In the following, we will overview four electrochemical models proposed in the literature, showcasing numerical results covering electric potential profiles, interfacial currents, fluxes and concentration profiles. In terms of notation, T will stand for the absolute temperature, F for Faraday’s constant, ${\mathrm{D}|}_{\alpha }$ for the diffusivity of species $\alpha$; symbol $\mathbb{1}$ denotes the identity matrix.

2.1. One-Dimensional Single-Ion Conduction Models [18]

A one-dimensional model of a Li/LiPON thin-film micro battery was developed by Fabre et al. [18] and is depicted schematically in Figure 1a. The x-axis, which points out the characteristic lengths of the model is directed from the negative towards the positive electrode; the interface between the negative electrode and the solid electrolyte hosts its origin.

Authors aimed at keeping their model as simple as possible, by introducing proper assumptions in order to come up with a reduced set of input parameters that can be measured from purposely designed experiments. The model is isothermal with no self-heating, neglects volume alterations during charge/discharge, and assumes that redox activities occur solely on the surface area, which remains unchanged throughout cycling. The negative electrode is a metallic film of lithium with negligible Ohmic losses. The ionic transfer in the solid electrolyte is described by a single ion conduction model: as such, the concentration of lithium ions across the solid electrolyte ${\mathrm{Li}}^{+}$ is uniform. This property straightforwardly comes out as long as electro-neutrality approximation holds, a property of the governing equations largely discussed in [22,24,25].

Accordingly, the model traces lithium ${\mathrm{Li}}^{\oplus }$ diffusion and electron $e^{-}$ migration in the positive electrode, single ions ${\mathrm{Li}}^{+}$ migration in the solid electrolyte, and charge-transfer kinetics at the electrode/electrolyte interfaces. Three unknown fields are required, the concentration of lithium in the positive electrode $c_{{\mathrm{Li}}^{\oplus }}(x,t)$, the electric potentials ${\varphi }_e(x,t)$ and ${\varphi }_c(x,t)$ in the solid electrolyte and in cathode, respectively.

The electric potentials in the positive electrode ${\varphi }_c$ and in the solid electrolyte ${\varphi }_e$ are related to the current densities ${\overrightarrow{i}}_c(x,t)$ and ${\overrightarrow{i}}_e(x,t)$ by means of Ohm’s law, through the electrical $k_c$ and ionic $k_e$ conductivities, respectively:

$$\begin{array}{ccc}\hfill & {\overrightarrow{i}}_e(x,t)=-k_e\nabla {\varphi }_e(x,t)\hfill & \hfill 0\le x\le L_e,\end{array}$$

(4a)

$$\begin{array}{ccc}\hfill & {\overrightarrow{i}}_c(x,t)=-k_c\nabla {\varphi }_c(x,t)\hfill & \hfill L_e\le x\le L_e+L_c.\end{array}$$

(4b)

The third governing equation in the positive electrode is a planar solid-state diffusion equation, which describes the (neutral, in the sense detailed beforehand of shielded positive ion) transport of lithium ${\mathrm{Li}}^{\oplus }$ in the electrode:

$$\begin{array}{ccc}\hfill & \frac{\partial c_{{\mathrm{Li}}^{\oplus }}(x,t)}{\partial t}=\mathrm{div}\left[{\mathrm{D}|}_{{\mathrm{Li}}^{\oplus }}\nabla c_{{\mathrm{Li}}^{\oplus }}(x,t)\right]\hfill & \hfill L_e\le x\le L_e+L_c,\end{array}$$

(4c)

with diffusion coefficient ${\mathrm{D}|}_{{\mathrm{Li}}^{\oplus }}$. A concentration-dependent ionic diffusion coefficient was studied in [18], which resulted in more accurate outcomes.

Denoting with $\overrightarrow{n}$ the outward normal to the surface at boundary, the two boundary conditions required for Equation (4c) are the lithium flux density at the reaction surface and the zero-flux condition for lithium ${\mathrm{Li}}^{\oplus }$ at the electrode/collector interface:

$$\begin{array}{ccc}\hfill & {\mathrm{D}|}_{{\mathrm{Li}}^{\oplus }}\nabla c_{{\mathrm{Li}}^{\oplus }}(x,t)·\overrightarrow{n}=\frac{-i_{e/c}\left(t\right)}{F}\hfill & \hfill x=L_e,\end{array}$$

(5a)

$$\begin{array}{ccc}\hfill & {\mathrm{D}|}_{{\mathrm{Li}}^{\oplus }}\nabla c_{{\mathrm{Li}}^{\oplus }}(x,t)·\overrightarrow{n}=0\hfill & \hfill x=L_e+L_c.\end{array}$$

(5b)

Conditions for ${\varphi }_e$ and ${\varphi }_c$, either at the anode-electrolyte ($a/e$) or at the electrolyte-cathode ($e/c$) interfaces

$$\begin{array}{c}\hfill {\overrightarrow{i}}_e(0,t)·\overrightarrow{n}=-i_{a/e}\left(t\right),{\overrightarrow{i}}_e(L_e,t)·\overrightarrow{n}=-i_{e/c}\left(t\right),{\overrightarrow{i}}_c(L_c,t)·\overrightarrow{n}=i_{e/c}\left(t\right)\end{array}$$

(6)

are modeled via Butler-Volmer equations in Fabre et al. [18], in the form

$$i_s\left(t\right)=i_s^0\left(t\right)\left(e^{\frac{\alpha F}{RT}\eta \left(t\right)}-e^{-\frac{(1-\alpha )F}{RT}\eta \left(t\right)}\right)$$

(7)

with s either $a/e$ or $e/c$ and $\alpha$ the so-called anodic (cathodic) charge transfer coefficient, usually both taken to be equal to $0.5$. The overpotential

$$\eta \left(t\right)=〚\varphi 〛-OCP$$

(8)

is the difference between the jump $〚\varphi 〛$ of the electric potential at the electrolyte/electrode interface (always defined as the electrode potential minus the electrolyte potential) and the open circuit potential (OCP) measured experimentally or calculated theoretically as in [33]. The exchange current $i_{a/e}^0$ is given by:

$$i_{a/e}^0\left(t\right)=F{k}_2{k}_{-2}{\left(c_{{\mathrm{Li}}^{+}}^{sat}-c_{{\mathrm{Li}}^{+}}(0,t)\right)}^{\alpha }{c}_{{\mathrm{Li}}^{+}}{(0,t)}^{1-\alpha },$$

(9)

where $k_2$ and $k_{-2}$ are the forward and backward reaction rate constants for the reaction (2), $c_{{\mathrm{Li}}^{+}}^{sat}$ is the saturation concentration of lithium in the electrolyte. The exchange current $i_{e/c}^0$ is given by:

$$i_{e/c}^0\left(t\right)=F{k}_1{k}_{-1}{\left(c_{{\mathrm{Li}}^{+}}^{sat}-c_{{\mathrm{Li}}^{+}}(0,t)\right)}^{\alpha }{c}_{{\mathrm{Li}}^{+}}{(0,t)}^{1-\alpha }{\left(c_{{\mathrm{Li}}^{\oplus }}^{sat}-c_{{\mathrm{Li}}^{\oplus }}(0,t)\right)}^{\alpha }{c}_{{\mathrm{Li}}^{\oplus }}{(0,t)}^{1-\alpha },$$

(10)

where $k_1$ and $k_{-1}$ are the forward and backward reaction rate constants for the reaction (1) and $c_{{\mathrm{Li}}^{\oplus }}^{sat}$ represents the saturation lithium concentration within the cathode. Replacing $c_{{\mathrm{Li}}^{+}}(0,t)$ with the initial, uniform concentration of lithium in the electrolyte, the exchange current densities $i_{a/e}^0\left(t\right)$ and $i_{e/c}^0\left(t\right)$ simplify as

$$\begin{array}{c}\hfill i_{a/e}^0\left(t\right)=F{k}_{-}^{app},i_{e/c}^0\left(t\right)=F{k}_{+}^{app}{\left(c_{{\mathrm{Li}}^{\oplus }}^{sat}-c_{{\mathrm{Li}}^{\oplus }}(0,t)\right)}^{\alpha }{c}_{{\mathrm{Li}}^{\oplus }}{(0,t)}^{1-\alpha }.\end{array}$$

(11)

where $k_{-}^{app}$ and $k_{+}^{app}$ are apparent rate constant for negative electrode and defined as the multiplication of the terms in Eqs. (9) and (10).

Galvanostatic boundary conditions are eventually imposed,

$$\begin{array}{c}\hfill {\overrightarrow{i}}_c(L_e+L_c,t)·\overrightarrow{n}=i_{bat}\left(t\right),\end{array}$$

(12)

where $i_{bat}$ is the given galvanostatic current flowing across the 1D battery. The potential is set arbitrarily to ${\varphi }_e(0,t)=0$.

2.2. An Advanced Framework for Solid Electrolyte Intercalation Batteries [19]

Landstorfer and coworkers studied in [19] a non-porous electrode and a crystalline solid electrolyte. As in [18], they assumed a solid electrolyte with one mobile species (${\mathrm{Li}}^{+}$) and a uniform concentration of vacancies $c_{{\mathrm{n}}^{-}}$ that remains unaltered in time. The model entails a novel view of the electrode/solid electrolyte interface, which consists of an intermediate layer and a space charge region within the electrolyte. A visual representation of the model is given in Figure 1b.

Once ions randomly intercalate in the lattice structure of graphite or ${\mathrm{LiCoO}}_2$, a simple Fickian law diffusion analogous to Eq. (4c) accounts for ionic transport. The electric potential ${\varphi }_c$ is considered to be uniform within the electrodes, neglecting the ohmic loss (depicted in [18] via Eq. (4b)), while ${\varphi }_e$ influences the ionic transport in the electrolyte, ruled by non-equilibrium thermodynamics. The mass balance law relates the concentration of lithium $c_{{\mathrm{Li}}^{+}}(x,t)$ to the actual flux of cations ${\overrightarrow{h}}_{{\mathrm{Li}}^{+}}(x,t)$

$$\begin{array}{ccc}\hfill & \frac{\partial c_{{\mathrm{Li}}^{+}}(x,t)}{\partial t}=-\mathrm{div}\left[{\overrightarrow{h}}_{{\mathrm{Li}}^{+}}(x,t)\right]\hfill & \hfill 0\le x\le L_e.\end{array}$$

(13)

Using the standard linear relationship of Onsager type2

$$\begin{array}{ccc}\hfill & {\overrightarrow{h}}_{{\mathrm{Li}}^{+}}(x,t)=-\mathbf{M}(x,t)\nabla {\overline{\mu }}_{{\mathrm{Li}}^{+}}(x,t)\hfill & \hfill 0\le x\le L_e,\end{array}$$

(15)

between the flux ${\overrightarrow{h}}_{{\mathrm{Li}}^{+}}$ and the gradient of the electrochemical potential ${\overline{\mu }}_{{\mathrm{Li}}^{+}}$, the Clausis-Duhem inequality is satisfied a priori and thermodynamics consistency is granted, as largely discussed in [24,25]. The chemical potential is the functional derivative of the Gibbs free energy (or Helmholtz free energy according to [24,25]) with respect to the ionic concentration $c_{{\mathrm{Li}}^{+}}(x,t)$. The splitting of chemical and electrical potentials used in [19] is classical, as it is the free energy of mobile guest atoms interacting with a host medium, described by a regular solution model [36,37]. The electrochemical potential was eventually derived as the sum of the chemical and electrostatic potential as follows,

$$\begin{array}{ccc}\hfill & {\overline{\mu }}_{{\mathrm{Li}}^{+}}(x,t)={\mu }_{{\mathrm{Li}}^{+}}(x,t)+F{\varphi }_e(x,t)\hfill & \hfill 0\le x\le L_e,\end{array}$$

(16a)

$$\begin{array}{ccc}\hfill & {\mu }_{{\mathrm{Li}}^{+}}(x,t)={\mu }_{{\mathrm{Li}}^{+}}^0+RTln\frac{{\theta }_{{\mathrm{Li}}^{+}}(x,t)}{1-{\theta }_{{\mathrm{Li}}^{+}}(x,t)}+RT\chi |\left[1-2{\theta }_{{\mathrm{Li}}^{+}}(x,t)\right]\hfill & \hfill 0\le x\le L_e,\end{array}$$

(16b)

with3 ${\theta }_{{\mathrm{Li}}^{+}}=c_{{\mathrm{Li}}^{+}}/c_{{\mathrm{Li}}^{+}}^{sat}$ . The mass flux to be inserted into the mass balance law (13) eventually holds:

$$\begin{array}{ccc}\hfill & {\overrightarrow{h}}_{{\mathrm{Li}}^{+}}(x,t)=\left(2\chi |{\theta }_{{\mathrm{Li}}^{+}}-\frac{1}{1-{\theta }_{{\mathrm{Li}}^{+}}}\right){\mathrm{D}|}_{{\mathrm{Li}}^{+}}\nabla c_{{\mathrm{Li}}^{+}}-\frac{{\mathrm{D}|}_{{\mathrm{Li}}^{+}}F}{RT}{c}_{{\mathrm{Li}}^{+}}\nabla {\varphi }_e\hfill & \hfill 0\le x\le L_e.\end{array}$$

(17)

The classical Nernst-Planck flux4

$$\begin{array}{c}\hfill {\overrightarrow{h}}_{\alpha }(x,t)=-{\mathrm{D}|}_{\alpha }\nabla c_{\alpha }(x,t)-\frac{z_{\alpha }F}{RT}{\mathrm{D}|}_{\alpha }{c}_{\alpha }(x,t)\nabla {\varphi }_e(x,t),\end{array}$$

(18)

is attained if the dilute limit is assumed $1-{\theta }_{{\mathrm{Li}}^{+}}\sim 1$ and the energy of interaction vanishes $\chi |=0$.

The whole electrolyte is thought as consisting of a space charge region and a bulk region, as shown in Figure 1b. In the bulk region electroneutrality is imposed a priori, $c_{{\mathrm{Li}}^{+}}$ is defined by the (given and uniform) concentration of vacancies, and Ohm’s law (4a) allows recovering the electric potential. On the contrary, in proximity of the interfaces, where the concentrations of cations and anions differ, the electric potential is obtained by Gauss law, which provides, after constitutive prescriptions, the following Poisson equations:

$$\begin{array}{ccc}\hfill & -{\nabla }^2{\varphi }_e(x,t)=\frac{F}{{\epsilon }_r{\epsilon }_0}(c_{{\mathrm{Li}}^{+}}(x,t)-c_{\mathrm{Li}})\hfill & \hfill 0\le x\le L_{s{r}_a},\end{array}$$

(19a)

$$\begin{array}{ccc}\hfill & -{\nabla }^2{\varphi }_e(x,t)=\frac{F}{{\epsilon }_r{\epsilon }_0}(c_{{\mathrm{Li}}^{+}}(x,t)-c_{{\mathrm{Li}}^{\oplus }})\hfill & \hfill L_e-L_{s{r}_c}\le x\le L_e,\end{array}$$

(19b)

where ${\epsilon }_r$ denotes the relative permittivity of the electrolyte and ${\epsilon }_0$ the vacuum permittivity. Local electroneutrality is not enforced in the space charge region [19], rather a weak (i.e. global) electroneutrality is prescribed in the whole electrolyte

$$\begin{array}{ccc}\hfill & F{\int }_V(c_{{\mathrm{Li}}^{+}}(x,t)-c_{{\mathrm{n}}^{-}})dV=0,\hfill & \hfill F{\int }_V(c_{{\mathrm{Li}}^{+}}(x,t)-c_{{\mathrm{n}}^{-}})dV=0,\end{array}$$

(20)

thus allowing local deviations between cation and anion concentrations while keeping the overall amount of cations and anions equal.

The interface between the electrodes and the solid electrolyte has been modelled as consisting of two intermediate layers, treated as plate capacitors in terms of potential jumps and flux continuity. Authors borrowed from [42] the constitutive equation for potential jumps (defined as the potential at an electrode minus the one at the electrolyte) as

$$\begin{array}{ccc}\hfill & 〚\varphi {〛}_{anode}=\frac{{\epsilon }_a}{C_a}\nabla {\varphi }_e(0,t)·\overrightarrow{n},\hfill & \hfill 〚\varphi {〛}_{cathode}=\frac{{\epsilon }_c}{C_c}\nabla {\varphi }_e(L_e,t)·\overrightarrow{n},\end{array}$$

(21)

with given surface capacitances $C_a$, $C_c$ and permittivities ${\epsilon }_a$, ${\epsilon }_c$.

Boundary conditions on fluxes at electrode - solid electrolyte interfaces resemble Butler-Volmer equations (7) in a form originally presented in [43] and named generalized Frumkin-Butler-Volmer equations. They read:

$$\begin{array}{cc}\hfill & {\overrightarrow{h}}_{{\mathrm{Li}}^{+}}(0,t)·\overrightarrow{n}=-{\tilde{k}}_2{e}^{\frac{▵G_1-\beta F({\varphi }_a-{\varphi }_e(0,t))}{RT}}{c}_{\mathrm{Li}}(0,t)+{\tilde{k}}_{-2}{e}^{\frac{▵G_{-1}-(1-\beta )F({\varphi }_a-{\varphi }_e(0,t)}{RT}}{c}_{{\mathrm{Li}}^{+}}(0,t),\hfill \end{array}$$

(22a)

$$\begin{array}{cc}\hfill & {\overrightarrow{h}}_{{\mathrm{Li}}^{+}}(L_e,t)·\overrightarrow{n}=-{\tilde{k}}_1{e}^{\frac{▵G_2-\beta F({\varphi }_c-{\varphi }_e(L_e,t))}{RT}}{c}_{{\mathrm{Li}}^{\oplus }}(L_e,t)+{\tilde{k}}_{-1}{e}^{\frac{▵G_{-2}-(1-\beta )F({\varphi }_c-{\varphi }_e(L_e,t)}{RT}}{c}_{{\mathrm{Li}}^{+}}(L_e,t),\hfill \end{array}$$

(22b)

where the reaction rate constants defined in reactions (1) and (2) are taken of Arrhenius type, i.e., $k_n={\tilde{k}}_n{e}^{E_n/RT}$, $n=1,2$. The Gibbs energies of activation $▵G_n$ are further parameters of the model and relate to the OCP in Butler-Volmer equations.

2.3. An Advanced All-Solid-State Li-Ion Battery Model [20]

More complex one-dimensional mathematical models have been proposed in a series of publications from Notten’s group [26,44,45] for a micro-battery $Li/LiPON/LiCo{O}_2$. In these studies, ionic transport in the solid electrolyte involves the ionization reaction (3) of immobile, oxygen-bound lithium $\mathrm{LiO}$ to mobile ${\mathrm{Li}}^{+}$ ions and negatively charged vacancies. Charge-transfer kinetics at both electrode/electrolyte interfaces, diffusion and migration of mobile lithium ions in the electrolyte (${\mathrm{Li}}^{+}$) and positive electrode (${\mathrm{Li}}^{\oplus }$) were accounted for. In their most recent work, [20], additional features have been introduced, such as mixed ionic/electronic conductivity in the positive electrode, electrical double layers occurring at both electrode/electrolyte interfaces representing the space-charge separation phenomena that differ from [19], variable ionic and electronic diffusion coefficients that depend on the lithium concentration inside the positive electrode.

Figure 1c displays a discharge process. The anode consists of a lithium $\mathrm{Li}$ foil, the cathode of a ${\mathrm{LiCoO}}_2$ film, while LiPON is used as the electrolyte material. The current collector is tied from the top of the ${\mathrm{LiCoO}}_2$. As for [18] and [19], also this model is isothermal (no self-heating) and redox processes only occur at the interfaces between the electrolyte and the electrode layers; volume changes of the electrolyte during cycling are neglected and the active surface area does not change over cycling.

As a distinctive feature of this class of models, the ionic transfer in the solid electrolyte is not described by a single ion conduction model and the concentration of lithium ions through the solid electrolyte is generally not uniform even though electro-neutrality approximation holds. Denoting with $c_{{\mathrm{Li}}^{+}}$ the concentration of mobile ${\mathrm{Li}}^{+}$ ions, with $c_{\mathrm{LiO}}$ the concentration of immobile lithium, with $c_{{\mathrm{n}}^{-}}$ the concentration of uncompensated negative charges, and with $\delta$ the fraction of Li at equilibrium, the equilibrium concentration of the charge carriers is $c_{{\mathrm{Li}}^{+}}^{eq}=c_{{\mathrm{n}}^{-}}^{eq}=\delta c_0$ and the concentration of the remaining immobile lithium is $c_{\mathrm{LiO}}^{eq}=(1-\delta )c_0$. The overall rate of the charge carrier generation according to reaction (3) is

$$w^{\left(3\right)}=k_f^{\mathrm{ion}}{c}_{\mathrm{LiO}}-k_b^{\mathrm{ion}}{c}_{{\mathrm{Li}}^{+}}{c}_{{\mathrm{n}}^{-}}.$$

(23)

The ratio

$$K_{\mathrm{eq}}^{\mathrm{ion}}=k_f^{\mathrm{ion}}/k_b^{\mathrm{ion}}$$

(24)

is the equilibrium constant of reaction (3) and is related to the the fraction of Li at equilibrium $\delta$, see [20]. In addition to earlier models, two electrical double layer capacitances $C_{a/e}^{dl}$ and $C_{e/c}^{dl}$ and a geometric capacitance $C_{geo}$ were introduced in [20]. As in [19], double layer capacitances attempt at capturing the response of the space-charge very narrow layers as for electric capacitors. Whereas the concept resembles [19], the implementation is different. Capacitors in [19] are described “in series”, whereas in [20] “in parallel” (compare Figures 1b and 1c). In view of this assumption, the current at electrode/electrolyte interfaces splits into two terms, the faradaic contribution that drives the reduction/oxidation charge-transfer reactions (1)-(2) ($i_s^{ct}$) and the non-faradaic contribution that feeds the double layer ($i_s^{dl}$), with s either $a/e$ or $e/c$.

During charging, the ${\mathrm{Li}}^{+}$ ions released from the positive surface must cross the solid electrolyte and are reduced into metallic Li at anode. Electrons, generated by the charge-transfer reactions (1) and (2), flow across the Li foils and the electronic collector, with potential drops that follow Ohm’s law similar to Eq. (4a). Transport of ionic concentrations5 $c_{{\mathrm{Li}}^{+}}$ and $c_{{\mathrm{n}}^{-}}$ in the solid electrolyte is ruled by the mass continuity equation (13), properly extended as

$$\begin{array}{ccc}\hfill & \frac{\partial c_{{\mathrm{Li}}^{+}}(x,t)}{\partial t}=-\mathrm{div}\left[{\overrightarrow{h}}_{{\mathrm{Li}}^{+}}(x,t)\right]+w^{\left(3\right)}(x,t)\hfill & \hfill 0\le x\le L_e,\end{array}$$

(25a)

$$\begin{array}{ccc}\hfill & \frac{\partial c_{{\mathrm{n}}^{-}}(x,t)}{\partial t}=-\mathrm{div}\left[{\overrightarrow{h}}_{{\mathrm{n}}^{-}}(x,t)\right]+w^{\left(3\right)}(x,t)\hfill & \hfill 0\le x\le L_e,\end{array}$$

(25b)

in order to account for the reaction rate $w^{\left(3\right)}$ (see Eq. (23)). The generic mass flux ${\overrightarrow{h}}_{\alpha }$, with $\alpha ={\mathrm{Li}}^{+},{\mathrm{n}}^{-}$, is constitutively described by the Nernst-Planck law (18), which carries the electric potential as a further unknown field. Coupling with an additional relation is thus mandatory, to model the migration process. The most common selection for such an additional relation in battery modeling is the electro-neutrality condition

$$c_{{\mathrm{Li}}^{+}}(x,t)=c_{{\mathrm{n}}^{-}}(x,t).$$

(26)

By substituting eq. (18) into eq. (25) and subtracting eq.(25a) from eq.(25b), two independent partial differential equations eventually arise:

$$\begin{array}{cc}\hfill & \frac{\partial c_{{\mathrm{Li}}^{+}}(x,t)}{\partial t}=\mathrm{div}\left[{\mathrm{D}|}_{{\mathrm{Li}}^{+}}\nabla c_{{\mathrm{Li}}^{+}}(x,t)+\frac{F{\mathrm{D}|}_{{\mathrm{Li}}^{+}}}{RT}{c}_{{\mathrm{Li}}^{+}}(x,t)\nabla {\varphi }_e(x,t)\right]+w^{\left(3\right)}(x,t)\hfill \end{array}$$

(27a)

$$\begin{array}{cc}\hfill & \mathrm{div}\left[({\mathrm{D}|}_{{\mathrm{n}}^{-}}-{\mathrm{D}|}_{{\mathrm{Li}}^{+}})\nabla c_{{\mathrm{Li}}^{+}}(x,t)\right]-\mathrm{div}\left[({\mathrm{D}|}_{{\mathrm{n}}^{-}}+{\mathrm{D}|}_{{\mathrm{Li}}^{+}})\frac{F}{RT}{c}_{{\mathrm{Li}}^{+}}(x,t)\nabla {\varphi }_e(x,t)\right]=0.\hfill \end{array}$$

(27b)

To be solved they require the initial concentrations for $c_{{\mathrm{Li}}^{+}}$

$$\begin{array}{cc}\hfill c_{{\mathrm{Li}}^{+}}(x,0)& =c_{{\mathrm{Li}}^{+}}^{eq}=\delta c_0\hfill \end{array}$$

(28)

and the Neumann conditions on fluxes at the left and right boundaries of the electrolyte

$${\overrightarrow{h}}_{{\mathrm{Li}}^{+}}(0,t)·\overrightarrow{n}=-\frac{i_{a/e}^{ct}\left(t\right)+i_{a/e}^{dl}\left(t\right)}{F},{\overrightarrow{h}}_{{\mathrm{Li}}^{+}}(L_e,t)·\overrightarrow{n}=-\frac{i_{e/c}^{ct}\left(t\right)+i_{e/c}^{dl}\left(t\right)}{F}.$$

(29)

In the positive electrode a mixed ionic/electronic conductivity is considered. The mass balance equations that characterize the transport of ${\mathrm{Li}}^{\oplus }$ ions and electrons ${\mathrm{e}}^{-}$ are similar to Eqs. (25) for the solid electrolyte and write

$$\begin{array}{ccc}\hfill & \frac{\partial c_{{\mathrm{Li}}^{\oplus }}(x,t)}{\partial t}=-\mathrm{div}\left[{\overrightarrow{h}}_{{\mathrm{Li}}^{\oplus }}(x,t)\right]\hfill & \hfill L_e\le x\le L_e+L_c,\end{array}$$

(30a)

$$\begin{array}{ccc}\hfill & \frac{\partial c_{{\mathrm{e}}^{-}}(x,t)}{\partial t}=-\mathrm{div}\left[{\overrightarrow{h}}_{{\mathrm{e}}^{-}}(x,t)\right]\hfill & \hfill L_e\le x\le L_e+L_c.\end{array}$$

(30b)

The generic mass flux ${\overrightarrow{h}}_{\alpha }$, $\alpha =\{{\mathrm{Li}}^{\oplus },{\mathrm{e}}^{-}\}$ is constitutively described by the Nernst-Planck law (18). This choice of independent motion of electrons and ionic intercalated lithium makes the governing equations different from Fabre’s [18] equations ((4c) to be compared with (30a)).

Additionally, the model from [20] makes the positive electrode’s diffusion coefficients dependent on the concentration of ions. Experimental results show that the local electrochemical environment has a major impact on solid-state diffusion. The model further exploits the electro-neutrality approximation inside the cathode, implying $c_{{\mathrm{Li}}^{\oplus }}(x,t)=c_{{\mathrm{e}}^{-}}(x,t)$, and a space-time proportionality of the diffusion coefficients for ${\mathrm{Li}}^{\oplus }$ and ${\mathrm{e}}^{-}$. Classical mathematical passages allow to retrieve Eq. (4c) as the single PDE required to model the mass transport in the electrode, provided that the diffusivity ${\mathrm{D}|}_{{\mathrm{Li}}^{\oplus }}$ is replaced by a suitable combination of electron and ionic diffusivities.

The initial concentration of the charge carrier at $t=0$, when no concentration profile developed yet, is equal to its equilibrium concentration

$$\begin{array}{cc}\hfill c_{{\mathrm{Li}}^{\oplus }}(x,0)& =c_{{\mathrm{Li}}^{\oplus }}^{eq}\hfill \end{array}$$

(31)

Neglecting the impact of geometric capacitance, the Neumann boundary conditions related to fluxes at both the left and right boundaries of the electrolyte result in:

$$\begin{array}{cc}\hfill {\overrightarrow{h}}_{{\mathrm{Li}}^{\oplus }}(L_e,t)·\overrightarrow{n}& =\frac{i_{e/c}^{ct}\left(t\right)+i_{e/c}^{dl}\left(t\right)}{F}\hfill \end{array}$$

(32a)

$$\begin{array}{cc}\hfill {\overrightarrow{h}}_{{\mathrm{Li}}^{\oplus }}(L_e+L_c,t)·\overrightarrow{n}& =0\hfill \end{array}$$

(32b)

$$\begin{array}{cc}\hfill {\overrightarrow{h}}_{{\mathrm{e}}^{-}}(L_e,t)·\overrightarrow{n}& =0\hfill \end{array}$$

(32c)

$$\begin{array}{cc}\hfill {\overrightarrow{h}}_{{\mathrm{e}}^{-}}(L_e+L_c,t)·\overrightarrow{n}& =-\frac{i_{bat}\left(t\right)}{F}\hfill \end{array}$$

(32d)

with $i_{bat}\left(t\right)$ the given galvanostatic current flowing across the 1D battery, defined in Equation (12).

The non-faradaic current (dis)charging the electrical double layers $i_{\alpha }^{dl}\left(t\right)$, can be defined in derivative form as

$$i_{\alpha }^{dl}\left(t\right)=C_{\alpha }^{dl}\frac{\partial 〚\varphi 〛}{\partial t},$$

(33)

where the jump $〚\varphi 〛$ of the electric potential at the electrolyte/electrode interface is always defined as the electrode potential minus the electrolyte potential and $\alpha =a/e,e/c$. Equation (33) shall be compared with (21) in [19].

The faradaic current proposed in [20] emanates from charge transfer kinetics, in a form that extends Butler-Volmer equation (7) to the mass transfer-influenced conditions [46]. The expression of $i_{a/e}^{ct}$ at the metallic lithium electrode interface reads

$$i_{a/e}^{ct}=i_{a/e}^0\left(\frac{c_{\mathrm{Li}}(0,t)}{{\overline{c}}_{\mathrm{Li}}}exp\left[\frac{{\alpha }_{a}F}{RT}{\eta }_a\left(t\right)\right]-\frac{c_{{\mathrm{Li}}^{+}}(0,t)}{{\overline{c}}_{{\mathrm{Li}}^{+}}}exp\left[-\frac{(1-{\alpha }_a)F}{RT}{\eta }_a\left(t\right)\right]\right),$$

(34a)

where ${\overline{c}}_{{\mathrm{Li}}^{+}}$ is the average bulk concentration of species ${\mathrm{Li}}^{+}$, ${\overline{c}}_{Li}$ is the bulk activity of the metallic Li, ${\alpha }_a$ is the charge transfer coefficient for reaction eq. (2), ${\eta }_a$ is the overpotential (8) of the charge transfer reaction at the negative electrode, and the exchange current $i_{a/e}^0$ is given by:

$$i_{a/e}^0=F{k}_2{\left({\overline{c}}_{{\mathrm{Li}}^{+}}\right)}^{{\alpha }_a}{\left({\overline{c}}_{\mathrm{Li}}\right)}^{1-{\alpha }_a},$$

(34b)

with $k_2$ the standard rate constant for reaction eq. (2) . The expression of $i_{e/c}$ at the positive electrode interface is given by the Butler-Volmer equation (7), with the exchange current density

$$i_{e/c}^0=F{k}_1{c}_{{\mathrm{Li}}^{\oplus }}^{sat}{\left(1-\frac{{\overline{c}}_{{\mathrm{Li}}^{\oplus }}}{c_{L{i}^{\oplus }}^{sat}}\right)}^{{\alpha }_c}{\left(\frac{{\overline{c}}_{{\mathrm{Li}}^{\oplus }}}{c_{L{i}^{\oplus }}^{sat}}\right)}^{1-{\alpha }_c}{\left({\overline{c}}_{{\mathrm{Li}}^{+}}\right)}^{{\alpha }_c},$$

(34c)

with $k_1$ the standard rate constant for reaction eq.(1). The reader may refer to [20] for further details on these equations and on the geometric capacitance.

2.4. Two-Mechanisms Model for All-Solid-State Lithium-Ion Batteries

A model that accounts for two mechanisms of ionic conduction was recently published in [21]. Rooted in the thermo-mechanics of continua, the model builds upon the work of Raijmakers et al. [20] to enhance the description of vacancy replenishment in a LIPON solid electrolyte and applies to LLZO as well [47] according to Figure 3. The equations that depict ionic transfer in [21] are multi-scale compatible, too. This feature seems to be particularly pertinent for composite cathodes, as indicated in [16].

As eloquently detailed in [20], some of the Li ions in a solid electrolyte are thermally excited at standard conditions. The chemical ionization reaction (3) occurs, leaving behind uncompensated negative charges associated to a vacancy in the LIPON matrix at the place formerly occupied by lithium. Raijmakers et al. have depicted vacancies motion with the same conceptual formalism used for negative ions in liquid electrolytes, i.e., as able to move in the solid matter driven by an entropic Brownian motion together with migration within an electric field (see eq. (25b)). This picture is thermodynamically encapsulated in the constitutive law (18).

In [21,48] the dynamic filling of vacancies by neighboring positions has been explicitly modeled. To capture this process some ions, denoted henceforth with ${\mathrm{Li}}_{\mathrm{int}}^{+}$, move in a meta-stable interstitial state after the chemical ionization reaction (3) occurs, whereas the remaining ${\mathrm{Li}}^{+}$ ions hop and fill neighboring vacancies. This dynamic behavior is described by a further reaction

$${\mathrm{Li}}^{+}\underset{k_b^{\mathrm{int}}}{\stackrel{k_f^{\mathrm{int}}}{\rightleftarrows }}{\mathrm{Li}}_{\mathrm{int}}^{+},$$

(35)

where $k_f^{\mathrm{int}}$ and $k_b^{\mathrm{int}}$ are the rate constants for reaction (35). When $k_f^{\mathrm{int}}=0$, no interstitial mechanism is accounted for and the model restricts to a classical single ion conducting solid electrolyte [18].

In a nutshell, reaction (3) makes lithium ions capable to leave the host-sites and move within the complex amorphous LIPON structure, either by filling neighboring vacancies or to flow interstitially. The proportion of ions in these two mechanisms is governed by reaction (35). In this formalism, positive ions are the only moving species, whereby the concentration of negatively charged vacancies is the outcome of the motion process and do not possess any intrinsic motility. In this sense, there is no direct flow ${\overrightarrow{h}}_{{\mathrm{n}}^{-}}(x,t)$ of negative charges, contrary to eq. (25b), and the local concentration of vacancies is altered merely by the chemical ionization reaction eq. (35).

This conceptual picture frames into the following set of mass balance equations, which characterize the immobile lithium $\mathrm{LiO}$, the negative charges ${\mathrm{n}}^{-}$, the transport of the lithium ions ${\mathrm{Li}}_{\mathrm{int}}^{+}$ that go interstitial and the remaining ${\mathrm{Li}}^{+}$ that hop:

$$\begin{array}{ccc}\hfill & \frac{\partial c_{\mathrm{LiO}}}{\partial t}=-w^{\left(3\right)}\hfill & \hfill 0\le x\le L_e,\end{array}$$

(36a)

$$\begin{array}{ccc}\hfill & \frac{\partial c_{{\mathrm{n}}^{-}}}{\partial t}=w^{\left(3\right)}\hfill & \hfill 0\le x\le L_e,\end{array}$$

(36b)

$$\begin{array}{ccc}\hfill & \frac{\partial c_{{\mathrm{Li}}^{+}}}{\partial t}+\mathrm{div}\left[{\overrightarrow{h}}_{{\mathrm{Li}}^{+}}\right]=w^{\left(3\right)}-w^{\left(35\right)}\hfill & \hfill 0\le x\le L_e,\end{array}$$

(36c)

$$\begin{array}{ccc}\hfill & \frac{\partial c_{{\mathrm{Li}}_{\mathrm{int}}^{+}}}{\partial t}+\mathrm{div}\left[{\overrightarrow{h}}_{{\mathrm{Li}}_{\mathrm{int}}^{+}}\right]=w^{\left(35\right)}\hfill & \hfill 0\le x\le L_e.\end{array}$$

(36d)

having set

$$w^{\left(35\right)}=k_f^{\mathrm{int}}{c}_{{\mathrm{Li}}^{+}}-k_b^{\mathrm{int}}{c}_{{\mathrm{Li}}_{\mathrm{int}}^{+}}$$

(36e)

according to reaction (35). The ratio

$$K_{\mathrm{eq}}^{\mathrm{int}}=k_f^{\mathrm{int}}/k_b^{\mathrm{int}}$$

is the equilibrium constant of reaction (35). The set of 4 mass balance equations (36a)-(36d) contains 5 unknowns, the concentrations $c_{\mathrm{LiO}}$, $c_{{\mathrm{n}}^{-}}$, $c_{{\mathrm{Li}}^{+}}$, and $c_{{\mathrm{Li}}_{\mathrm{int}}^{+}}$ plus the electric potential $\varphi$, which is constitutively related to the mass fluxes ${\overrightarrow{h}}_{\alpha }$. The additional required equation is Ampère’s law (with Maxwell’s correction) 6. To conclude the set of balance equations, the usual balance of forces in small strains was accounted for in [48]

$$\mathrm{div}\left[\sigma \right]+\overrightarrow{b}=\overrightarrow{0}.$$

(36f)

The chemo-thermo-elastic strain $\epsilon$ is considered to be made up of two separate contributions: an elastic recoverable part after unloading ${\epsilon }^{el}$ and a swelling contribution due to the insertion of species in the host material ${\epsilon }^s$:

$$\epsilon ={\epsilon }^{el}+{\epsilon }^s.$$

(37)

The swelling contribution ( $\alpha ={\mathrm{Li}}^{+},{\mathrm{Li}}_{\mathrm{int}}^{+}$ )

$${\epsilon }^s=\sum _{\alpha }{\omega }_{\alpha }\left(c_{\alpha }-c_{\alpha }^0\right)\mathbb{1}$$

(38)

is assumed to be volumetric and proportional to the deviation $c_{\alpha }-c_{\alpha }^0$ from the reference concentration $c_{\alpha }^0$ by means of the chemical expansion coefficients ${\omega }_{\alpha }$ of species $\alpha$. They equal one third of the partial molar volumes at a given temperature.

To derive governing equations from the balance Eqs. (36), constitutive laws have been derived from a rigorous thermodynamic setting. The hopping mechanism is thermodynamically quite different from the interstitial motion. Therefore, making recourse to the same thermodynamic description for both mechanisms might be questionable. Aware of this limitation and inspired by [49,50,51] restricted to small strains, elaborating the electromagnetic contribution in the Helmholtz free energy $\psi$ from [24,25], the chemical potential (16b) of species $\alpha ={\mathrm{Li}}^{+},{\mathrm{Li}}_{\mathrm{int}}^{+}$ is extended to

$${\mu }_{\alpha }={\mu }_{\alpha }^0+RTln\left[\frac{{\theta }_{\alpha }}{1-{\theta }_{\alpha }}\right]+RT\chi |\left(1-2{\theta }_{\alpha }\right)+\frac{\partial {\psi }_{el}}{\partial c_{\alpha }},$$

(39)

as detailed in [34,35]. A simple choice for the elastic part of the free energy density ${\psi }_{el}$ in the small strain range is the usual quadratic form

$$\begin{array}{c}\hfill {\psi }_{el}(\epsilon ,c_{\alpha })=\frac{1}{2}K\left(c_{\alpha }\right)\mathrm{tr}{\left[\epsilon -{\epsilon }^s\right]}^2+G\left(c_{\alpha }\right)\left|\right|\mathrm{dev}\left[\epsilon -{\epsilon }^s\right]{\left|\right|}^2,\end{array}$$

(40)

where K, G are the bulk and shear modulus respectively and they are made dependent on species concentrations. The stress tensor ${\sigma }^e(\epsilon ,c_{\alpha })$ descends from thermodynamic restrictions (see [34] for details and extension to temperature dependency)

$$\begin{array}{ccc}\hfill {\sigma }^e& =& \frac{\partial {\psi }_{el}}{\partial \epsilon }=2G\mathrm{dev}\left[\epsilon \right]+K\left(\mathrm{tr}\left[\right]-\left|\right|{\epsilon }^s\left|\right|\right)\mathbb{1}.\hfill \end{array}$$

(41)

Note that the derivative $\partial {\psi }_{el}/\partial c_{\alpha }$, in eq. (39) is the sum of two contributions

$$\begin{array}{c}\hfill \frac{\partial {\psi }_{el}}{\partial c_{\alpha }}=-{\omega }_{\alpha }\mathrm{tr}\left[{\sigma }^e\right]+\frac{1}{2}\frac{\partial K}{\partial c_{\alpha }}\mathrm{tr}{\left[\epsilon -{\epsilon }^s\right]}^2+\frac{\partial G}{\partial c_{\alpha }}\left|\right|\mathrm{dev}\left[\epsilon -{\epsilon }^s\right]{\left|\right|}^2.\end{array}$$

(42)

The first emanates from the swelling part of the strain, and is present even if the material properties are independent on concentration of species. Nernst-Planck equation (18) is extended as follows

$$\begin{array}{cc}\hfill {\overrightarrow{h}}_{\alpha }(x,t)=& -{\mathrm{D}|}_{\alpha }\left[1-2\chi |{\theta }_{\alpha }(1-{\theta }_{\alpha })\right]\nabla c_{\alpha }(x,t)\hfill \\ & -3\mathit{M}\left(c_{\alpha }\right)K{\omega }_{\alpha }\left[3{\omega }_{\alpha }\nabla c_{\alpha }-\nabla \mathrm{tr}\left[\epsilon \right]\right]\hfill \\ & -\frac{{\mathrm{D}|}_{\alpha }F}{RT}{c}_{\alpha }(x,t)\nabla {\varphi }_e(x,t),\hfill \end{array}$$

(43)

with $\alpha ={\mathrm{Li}}^{+},{\mathrm{Li}}_{\mathrm{int}}^{+}$. The mobility tensor reads ( ${\theta }_{\alpha }=c_{\alpha }/c_{\alpha }^{sat}$ )

$$\mathit{M}\left(c_{\alpha }\right)=\frac{{\mathrm{D}|}_{\alpha }}{RT}{c}_{\alpha }^{sat}{\theta }_{\alpha }\left(1-{\theta }_{\alpha }\right)\mathbb{1},$$

(44)

It accounts for saturation and in this differs from Eq. (14). The energetic and entropic contributions in the constitutive law (43) have been already described in eq. (17). The mechanical contribution to the mass flux is driven by the chemical expansion coefficient, and derives from thermodynamic consistency.

Mass balance equations, after inserting (43) into eqs. (36c)-(36d), do not form a complete set, because ionic transport entails movement of mass as well as of charge. In order to build a multiscale compatible theory, the generally accepted electroneutrality assumption cannot be taken, since it prevents to impose the conservation of energy across the scales: this concept has been illustrated with great detail in [24,29] and will not be elaborated here further. Multiscale compatibility is granted by using Ampere’s law (with Maxwell’s correction in the realm of small strains)

$$\begin{array}{ccc}\hfill & \mathrm{div}\left[-\epsilon \mathrm{grad}\left[\frac{\partial {\varphi }_e}{\partial t}\right]+F\left({\overrightarrow{h}}_{{\mathrm{Li}}^{+}}+{\overrightarrow{h}}_{{\mathrm{Li}}_{\mathrm{int}}^{+}}\right)\right]=0\hfill & \hfill 0\le x\le L_e.\end{array}$$

(45)

When multiscale is not invoked, the electroneutrality assumption

$$c_{{\mathrm{Li}}^{+}}(x,t)+c_{{\mathrm{Li}}_{\mathrm{int}}^{+}}(x,t)=c_{{\mathrm{n}}^{-}}(x,t).$$

(46)

can be called for.

A widespread choice for the initial conditions for concentrations and electric potential enforces equilibrium conditions. They hold

$$\begin{array}{cc}\hfill & {\varphi }_e(x,0)=0,\hfill \end{array}$$

(47a)

$$\begin{array}{cc}\hfill & c_{\mathrm{LiO}}(x,0)=c_{\mathrm{LiO}}^{eq}=(1-\delta )c_0,\hfill \end{array}$$

(47b)

$$\begin{array}{cc}\hfill & c_{{\mathrm{n}}^{-}}(x,0)=c_{{\mathrm{n}}^{-}}^{eq}=\delta c_0,\hfill \end{array}$$

(47c)

$$\begin{array}{cc}\hfill & c_{{\mathrm{Li}}^{+}}(x,0)=c_{{\mathrm{Li}}^{+}}^{eq}=\frac{\delta c_0}{1+K_{\mathrm{eq}}^{\mathrm{int}}}=K_{\mathrm{eq}}^{\mathrm{ion}}\left(\frac{1}{\delta }-1\right),\hfill \end{array}$$

(47d)

$$\begin{array}{cc}\hfill & c_{{\mathrm{Li}}_{\mathrm{int}}^{+}}(x,0)=c_{{\mathrm{Li}}_{\mathrm{int}}^{+}}^{eq}=K_{\mathrm{eq}}^{\mathrm{int}}\frac{\delta c_0}{1+K_{\mathrm{eq}}^{\mathrm{int}}}=\delta c_0+K_{\mathrm{eq}}^{\mathrm{ion}}\left(1-\frac{1}{\delta }\right).\hfill \end{array}$$

(47e)

Three independet factors influence the equilibrium concentrations within the system: $c_0$ and the equilibrium constants for reactions (3) and (35). While $c_0$ can be determined with precision, estimating $K_{\mathrm{eq}}^{\mathrm{ion}}$ and $K_{\mathrm{eq}}^{\mathrm{int}}$ through experiments is fraught with substantial uncertainties. These three parameters are interrelated and play a role in defining the fraction of lithium existing in a mobile state at equilibrium, which is referred to as $\delta$, as detailed in reference [20]. $\delta$ holds,

$$\delta =\frac{2}{1+\sqrt{1+4\frac{c_0}{K_{\mathrm{eq}}^{\mathrm{ion}}}\frac{1}{1+K_{\mathrm{eq}}^{\mathrm{int}}}-1}}.$$

(48)

Equation (49) can be readily transformed to express $K_{\mathrm{eq}}^{\mathrm{int}}$ in terms of both $\delta$ and $K_{\mathrm{eq}}^{\mathrm{ion}}$ as follows,

$$K_{\mathrm{eq}}^{\mathrm{int}}=\frac{-c_0{\delta }^2+(1-\delta )K_{\mathrm{eq}}^{\mathrm{ion}}}{(\delta -1)K_{\mathrm{eq}}^{\mathrm{ion}}}.$$

(49)

Interface conditions for this advanced model are eqs. (29), accounting for identity (33) for the non-faradaic current contribution. The faradaic current emanates from charge transfer kinetics, as proposed in [20]. In view of the splitting of lithium flux in interstitial and hopping, constitutively specified by eq. (43), faradaic interface conditions split, too:

$$\begin{array}{ccc}\hfill & i_{a/e}^{ct}={i_{a/e}^{ct}}_{{\mathrm{Li}}^{+}}+{i_{a/e}^{ct}}_{{\mathrm{Li}}_{\mathrm{int}}^{+}},\hfill & \hfill i_{e/c}^{ct}={i_{e/c}^{ct}}_{{\mathrm{Li}}^{+}}+{i_{e/c}^{ct}}_{{\mathrm{Li}}_{\mathrm{int}}^{+}}\end{array}$$

(50)

with the interstitial and hopping contributions to charge transfer current clearly identified. They can be inferred from Butler-Volmer eqs. (7) or (34). The exchange current read:

$$\begin{array}{ccc}\hfill & {i_{a/e}^0}_{{\mathrm{Li}}^{+}}=F{k}_2{\left({\overline{c}}_{{\mathrm{Li}}^{+}}\right)}^{\alpha }{\left({\overline{c}}_{\mathrm{Li}}\right)}^{1-\alpha },\hfill & \hfill {i_{e/c}^0}_{{\mathrm{Li}}^{+}}=F{k}_1{c}_{{\mathrm{Li}}^{\oplus }}^{sat}{\left(1-\frac{{\overline{c}}_{{\mathrm{Li}}^{\oplus }}}{c_{{\mathrm{Li}}^{\oplus }}^{sat}}\right)}^{\alpha }{\left(\frac{{\overline{c}}_{{\mathrm{Li}}^{\oplus }}}{c_{{\mathrm{Li}}^{\oplus }}^{sat}}\right)}^{1-\alpha }{\left({\overline{c}}_{{\mathrm{Li}}^{+}}\right)}^{\alpha }\end{array}$$

(51a)

$$\begin{array}{ccc}\hfill & {i_{a/e}^0}_{{\mathrm{Li}}_{\mathrm{int}}^{+}}={i_{a/e}^0}_{{\mathrm{Li}}^{+}}{\left(\frac{{\overline{c}}_{{\mathrm{Li}}_{\mathrm{int}}^{+}}}{{\overline{c}}_{{\mathrm{Li}}^{+}}}\right)}^{\alpha },\hfill & \hfill {i_{e/c}^0}_{{\mathrm{Li}}_{\mathrm{int}}^{+}}={i_{e/c}^0}_{{\mathrm{Li}}^{+}}{\left(\frac{{\overline{c}}_{{\mathrm{Li}}_{\mathrm{int}}^{+}}}{{\overline{c}}_{{\mathrm{Li}}^{+}}}\right)}^{\alpha }.\end{array}$$

(51b)

where ${\overline{c}}_{{\mathrm{Li}}^{+}}$, ${\overline{c}}_{{{\mathrm{Li}}^{+}}_{\mathrm{int}}}$ are average bulk concentrations of species ${\mathrm{Li}}^{+}$ and ${{\mathrm{Li}}^{+}}_{\mathrm{int}}$, respectively. Note that, differently from [20], those averages are not time independent. Lacking more clear understanding, we assume that the non faradaic current $i_{\alpha }^{dl}\left(t\right)$ as in eq. (33) is proportional to the faradaic splitting, i.e.

$$i_s^{dl}={i_s^{dl}}_{{\mathrm{Li}}^{+}}+{i_s^{dl}}_{{\mathrm{Li}}_{\mathrm{int}}^{+}},\mathrm{with}{i_s^{dl}}_{{\mathrm{Li}}^{+}}=\frac{{i_s^{ct}}_{{\mathrm{Li}}^{+}}}{i_s^{ct}}{c}_s^{dl}\frac{\partial 〚\varphi 〛}{\partial t}\mathrm{and}{i_s^{dl}}_{{\mathrm{Li}}_{\mathrm{int}}^{+}}=\frac{{i_s^{ct}}_{{\mathrm{Li}}_{\mathrm{int}}^{+}}}{i_s^{ct}}{c}_s^{dl}\frac{\partial 〚\varphi 〛}{\partial t}.$$

(52)

where the jump $〚\varphi 〛$ of the electric potential at the electrolyte/electrode interface is always defined as the electrode potential minus the electrolyte potential and $s=a/e,e/c$. The Neumann conditions on fluxes at the left and right boundaries of the electrolyte eventually read:

$$\begin{array}{ccc}\hfill & {\overrightarrow{h}}_{{\mathrm{Li}}^{+}}(0,t)·\overrightarrow{n}=-\left({i_{a/e}^{ct}}_{{\mathrm{Li}}^{+}}+{i_{a/e}^{dl}}_{{\mathrm{Li}}^{+}}\right)/F,\hfill & \hfill {\overrightarrow{h}}_{{\mathrm{Li}}^{+}}(L_e,t)·\overrightarrow{n}=-\left({i_{e/c}^{ct}}_{{\mathrm{Li}}^{+}}+{i_{e/c}^{dl}}_{{\mathrm{Li}}^{+}}\right)/F,\end{array}$$

(53a)

$$\begin{array}{ccc}\hfill & {\overrightarrow{h}}_{{\mathrm{Li}}_{\mathrm{int}}^{+}}(0,t)·\overrightarrow{n}=-\left({i_{a/e}^{ct}}_{{\mathrm{Li}}_{\mathrm{int}}^{+}}+{i_{a/e}^{dl}}_{{\mathrm{Li}}_{\mathrm{int}}^{+}}\right)/F,\hfill & \hfill {\overrightarrow{h}}_{{\mathrm{Li}}_{\mathrm{int}}^{+}}(L_e,t)·\overrightarrow{n}=-\left({i_{e/c}^{ct}}_{{\mathrm{Li}}_{\mathrm{int}}^{+}}+{i_{e/c}^{dl}}_{{\mathrm{Li}}_{\mathrm{int}}^{+}}\right)/F.\end{array}$$

(53b)

Continuity of displacements and normal tractions are interface conditions for the mechanical governing equation (36f). The electrodes governing equations and boundary conditions do not differ from Section 2.3.

3. Benchmark Comparison

The two models illustrated in Section 2.3 and Section 2.4 are validated against the experimental outcomes described in [26]. That paper also provides values of material and geometrical parameters, which are collected in Table 1.

The positive electrode is a layer of ${\mathrm{LiCoO}}_2$ with thickness ${\mathrm{L}}_{\mathrm{c}}=0.32\mu \mathrm{m}$, deposited on a platinum substrate. A Lithium metal foil with thickness ${\mathrm{L}}_{\mathrm{a}}=0.50\mu \mathrm{m}$ is used as negative electrode. The solid electrolyte is a one-micron-thick (${\mathrm{L}}_{\mathrm{e}}=1.00\mu \mathrm{m}$) layer of LiPON. The surface area of the deposited electrodes is $A={10}^{-4}{\mathrm{m}}^2$ and the theoretical storage capacity of the battery is ${10}^{-5}\mathrm{Ah}$.

The electrochemical cell is subject to a galvanostatic process of discharge at different $C-rates$, under a temperature-controlled condition of ${25}^{\mathrm{o}}\mathrm{C}$. The current corresponding to $C-rate=j$ is denoted with $I_{jC}$. For $C-rate=1$, $I_{1C}$ amounts at ${10}^{-5}A$. Initial and boundary conditions are made compatible with thermodynamic equilibrium at $t=0$, tuning the density current $i_{bat}\left(t\right)$ in time as:

$$i_{bat}\left(t\right)=(1-e^{-t})\frac{I_{jC}}{A},$$

(54)

with t in seconds. In view of Eq. (54), the concentrations of ions across all battery components are uniform and at equilibrium at $t=0$ since no profiles have established yet. By enforcing the fraction of mobile lithium in the electrolyte $\delta =0.18$ and a maximum concentration of lithium host sites in the electrolyte $c_0=6.01·{10}^4\mathrm{mol}/{\mathrm{m}}^3$, eq. (47) provides

$$\begin{array}{ccccc}\hfill & c_{\mathrm{Li}}(x,0)=7.60·{10}^4\hfill & \hfill \mathrm{mol}/{\mathrm{m}}^3& \hfill & \hfill -L_a\le x\le 0,\end{array}$$

(55a)

$$\begin{array}{ccccc}\hfill & c_{\mathrm{LiO}}(x,0)=4.93·{10}^4\hfill & \hfill \mathrm{mol}/{\mathrm{m}}^3& \hfill & \hfill 0\le x\le L_e,\end{array}$$

(55b)

$$\begin{array}{ccccc}\hfill & c_{{\mathrm{n}}^{-}}(x,0)=1.08·{10}^4\hfill & \hfill \mathrm{mol}/{\mathrm{m}}^3& \hfill & \hfill 0\le x\le L_e,\end{array}$$

(55c)

$$\begin{array}{ccccc}\hfill & c_{{\mathrm{Li}}^{+}}(x,0)=5.68·{10}^3\hfill & \hfill \mathrm{mol}/{\mathrm{m}}^3& \hfill & \hfill 0\le x\le L_e,\end{array}$$

(55d)

$$\begin{array}{ccccc}\hfill & c_{{\mathrm{Li}}_{\mathrm{int}}^{+}}(x,0)=5.12·{10}^3\hfill & \hfill \mathrm{mol}/{\mathrm{m}}^3& \hfill & \hfill 0\le x\le L_e,\end{array}$$

(55e)

$$\begin{array}{ccccc}\hfill & c_{{\mathrm{Li}}^{\oplus }}(x,0)=c_{{\mathrm{Li}}^{\oplus }}^{\mathrm{eq}}=1.20·{10}^4\hfill & \hfill \mathrm{mol}/{\mathrm{m}}^3& \hfill & \hfill L_e\le x\le L_e+L_c.\end{array}$$

(55f)

The electric potential at the interface between the anode and the solid electrolyte is fixed as:

$$\varphi (0,t)=0\left[V\right]\forall t.$$

(56)

All material parameters used in the simulation are listed in Table 1. The equilibrium constants read

$$\begin{array}{ccc}\hfill & K_{\mathrm{eq}}^{\mathrm{ion}}=\frac{1.125·{10}^{-5}}{0.90·{10}^{-8}}=1250,\hfill & \hfill K_{\mathrm{eq}}^{\mathrm{int}}=\frac{8.10·{10}^{-9}}{0.90·{10}^{-8}}=0.9.\end{array}$$

(57)

for this benchmark comparison.

The governing equations are numerically solved with the finite element method, with in house implementation of weak forms in the commercial numerical software Matlab. The geometry and the unknown fields (see Figures 1c, 1d) are discretized with 61 linear elements; 1 element is sufficient for the anode, since the lithium concentration is uniform ad the electric potential is linear; the outcomes refer to a tessellation of 40 finite elements covering the electrolyte and 20 panels that discretize the cathode. In both cases, the mesh is refined near the electrode/electrolyte interface. The time marching is dealt with the backward Euler method, with fixed time increments of $\Delta t=1.0\mathrm{s}$.

The open circuit potential (OCP) used in the simulations has been reconstructed with splines stemming from experimental evidences in [26]. OCP could have been calculated analytically following the approach in [33]. Results, obtained with the analitically OCP, are reported in Appendix A.

The simulations account for a broad range of rates, from $1.0$ to $51.2$. The corresponding experimental discharging curves have been plot with dots in Figure 4b as a function of time, continuous lines correspond to simulations. Measurements and simulations agree well across the wide range of investigated $C-rates$. Obviously, the extracted charge decreases with increasing $C-rate$ due to the higher over-potentials. In fact lower discharge rates implies a lower rate of lithium insertion in the cathode and a more uniform Li distribution in the positive electrode.

Figure 5 depicts the electric potential $\varphi \left(x\right)$ profile in the battery at different times for two different discharge rates, i.e. $C-rate=3.20$ (which theoretically allows to discharge the battery in $1125s$) and $C-rate=51.2$ (for which the battery in principle completes the discharge process in $70s$). Simulations quit when the concentration of lithium ${\mathrm{Li}}^{\oplus }$ inside the cathode reaches the saturation limit $c_{{\mathrm{Li}}^{\oplus }}^{sat}$, after $1085s$ and $50s$ respectively. Cathodic saturation is indeed the limit factor for the battery operation (see also [52] for an extensive discussion on limiting factors in electrochemical cells, induced by materials and architectures).

At the initial time the electric potential discontinuity at the interfaces make Butler-Volmer currents vanishing. Based on the measured battery OCP, at full charge state $\Delta \varphi =4.2V$, as highlighted in Figure 5. During the discharge, the potential drops in the anode and in the electrolyte. However, the major changes in the potential profile occur at interface between electrolyte and cathode: the battery voltage decreases significantly during discharge as depicted in Figure 5a,b and it results almost uniform within the cathode. A similar behavior is observed for the model in [20] ( see Figure 5c,d).

Figure 6 depicts the evolution of lithium concentration $c_{{\mathrm{Li}}^{\oplus }}\left(x\right)$ in the cathode and in the solid electrolyte for the model in [21] ((a) and (b)) as well as for the one proposed in [20] ((c) and (d)). Since two ionic concentrations are concurrently present in the electrolyte, only their sum $(c_{{\mathrm{Li}}^{+}}+c_{{\mathrm{Li}}_{\mathrm{int}}^{+}})$ has been plotted in Figure 6a,b.

The anode, made of a metallic foil of lithium, is unaffected by the lithiation/delithiation processes and considered as an unlimited reservoir of lithium. The lithium ions intercalate inside the cathode and accumulates near at the electrolyte/cathode interface. The discharge process ends when the lithium in the cathode reaches its saturation limit $c_{{\mathrm{Li}}^{\oplus }}^{sat}=23400\mathrm{mol}{\mathrm{m}}^{-3}$. Saturation at the interface electrolyte/cathode is thus the limiting factor for the performance of this battery.

The concentration of $(c_{{\mathrm{Li}}^{+}}+c_{{\mathrm{Li}}_{\mathrm{int}}^{+}})$ in [20], initially uniform, increases near the anode interface, while decreasing at the cathode interface. This “liquid electrolyte” kind of behavior is justified since in [20] negative charges are allowed to move. In the formulation proposed in [21] negative charges are filled by hopping lithium, which in turn is allowed to intercalate. The time evolution of negative charges and total lithium is driven by the ionization reaction rates $w^{\left(3\right)}$ and $w^{\left(35\right)}$, as already discussed, and is very small in Figure 6a,b.

Concentrations ${\mathrm{c}}_{\mathrm{LiO}}$, ${\mathrm{c}}_{{\mathrm{n}}^{-}}$, ${\mathrm{c}}_{{\mathrm{Li}}^{+}}$, and ${\mathrm{c}}_{{\mathrm{Li}}_{\mathrm{int}}^{+}}$ are individually depicted in Figure 7. The concentration of interstitial lithium qualitatively resembles the cationic behavior observed in [20] across varying C-rates. The number of uncompensated negative vacancies, which are offset by interstitial lithium, is higher at the cathode during discharge. Due to electroneutrality and the assumption of no negative charge flow, the total concentration of lithium cannot vanish, unlike in liquid electrolyte systems and in [20] (compare Figure 6b,d). Nevertheless, the contributions of ${\mathrm{c}}_{{\mathrm{Li}}^{+}}$ and ${\mathrm{c}}_{{\mathrm{Li}}_{\mathrm{int}}^{+}}$ can independently approach zero (refer to Figure 7b). The implications of this occurrence are currently unclear.

The concentration of ${\mathrm{Li}}^{\oplus }$ at the cathode increases with discharge time, consistently with the Li-intercalation reaction described in Eq. (1). At higher C-rates, the expected concentration gradient becomes steeper, potentially leading to diffusion issues within the electrode. At $51.2$, this phenomenon is prominently observed in Figure 6b. Similar behavior is noted in [20] - see Figure 6d.

The total lithium concentration is depicted in Figure 8 for the model [21] ((a) and (b)) as well as for the model in [20] ((c) and (d)), at the anode/electrolyte interface with blue line and at the electrolyte/cathode interface with red line. According to [21], the tangent disappears at $t=0$. This is explained by electroneutrality coupled to the equilibrium condition $w=0$ imposed at $t=0$ in eq. (36b). The concentrations at the two electrodes evolve according to mass balance equations (36c) and (36d). Fluxes tend in time to achieve a uniform value across the electrolyte, i.e. the divergence term becomes less and less important and the evolution of concentrations turns out to be driven by the evolution of $w^{\left(3\right)}$ and $w^{\left(35\right)}$. It is thus expected that, at the same equilibrium constant, increasing the reaction constant (i.e. making the reaction faster) would allow reaching the steady state in due time. This is numerically confirmed in Figure 8 (a) and (b), where the higher the $k_f^{\mathrm{ion}}$ and $k_b^{\mathrm{ion}}$, keeping $K_{\mathrm{eq}}^{\mathrm{ion}}=1250$, the sooner the concentration plateau is reached. For the model in [20] the concentration response is similar to liquid electrolytes [24,25], as expected. Note that steady state in that case is obtained only at low C-rates.

4. Conclusion

This note has provided a review of four models of ionic conduction in solid electrolytes, highlighting the assumptions and the most relevant conceptual implications. All the formulations have been comprehensively examined from a theoretical perspective. The model in [21] and its predecessor [20] have also been compared from a computational standpoint via the finite element method. This assessment has been achieved analyzing interface currents, electric potential, flux patterns, and concentration profiles.

An all-solid-state thin-film $\mathrm{Li}/\mathrm{LiPON}/{\mathrm{LiCoO}}_2$ microbattery was modelled in one dimension by Fabre et al. [18], addressing the kinetics involved in charge-transfer at the interfaces between the electrodes and the electrolyte, as well as the diffusion and electron migration of lithium within the cathode, along with the migration of lithium ions throughout the solid electrolyte. A modified Poisson-Nernst-Planck(PNP) system of equations for the solid electrolyte was established in [19]. The model in [20] accounts for diffusion and migration of mobile lithium ions in the solid electrolyte, the kinetics regarding the charge transfer at the interfaces between both electrodes and the electrolyte, and the diffusion and migration of lithium ions and electrons at the cathode are all included in these processes. The LiPON and LCO layers between the metallic Li electrode and Pt current collector act as electronic conductors, resembling a dielectric barrier. To enhance the model, a geometric capacitor is introduced. Additionally, to account for the space-charge effect, double-layer capacitors are positioned at both interfaces of the electrode and electrolyte. Ionic mobility in [21] is a result of hopping and interstitial diffusion, i.e. two-mechanisms for ionic conduction have been considered.

While the four models dealt within the paper are presented at higher degree of complexity, they share some common assumptions. They are all thin-film based: the first model encapsulates this assumption in the constitutive equations, while the others are in principle fully three-dimensional and hence applicable also to composite electrode/electrolyte systems. All presented models are isothermal (no self-heating), although it is quite well known, at least in liquide electrolytes, that temperature affects the response of the battery [53]. All above assumptions shall be removed in order to replicate realistic behavior of batteries, especially with the aim of capturing the aging processes in ASSB that may preclude their own market success.

With battery cycling, the volume of electrodes changes due to repeated insertion and removal of Li atoms, causing disconnection and reduced contacts. Swelling is expected also in the solid electrolyte lattice, although this phenomena has certainly been less studied experimentally. In a recent work of Tian et al. [9], the effect of imperfect contact area was incorporated into a 1-D Newman battery model, by assuming the current and Li concentration will be localized at the contacted area of the interfaces. Among the presented publications, only [21] does not neglect volume changes during charge/discharge. In our simulations, though, the active surface area, where redox processes occur, were unaltered over cycling. Our future research plans are addressed to handle the mechanical problem in a three dimensional setting, with large strain mechanics as in [49,54,55,56].