Electrochemical Modelling Applied to Intercalation Phenomena Using Lattice Kinetic Monte Carlo Simulations: Galvanostatic Simulations

1. Introduction

Introduction:

Lithium-ion batteries (LIBs) are undoubtedly the energy storage systems that have changed the energy paradigm worldwide. However, the global demand for lithium and the development of technologies such as electric vehicles require batteries with higher storage capacity and shorter charging times. The industry is also looking for alternatives to traditional LIBs, such as sodium-ion batteries [1,2,3], which are already a reality, and other alkaline metal batteries. This is why careful optimization and innovative development of electrode materials are essential.

Any battery based on metal-ion technology basically has the same type of operation. A cathode and an anode where the ions are intercalated during the charge and discharge cycles of the electrochemical cell, a separator that prevents electrical contact between the electrodes and an electrolyte, which is the phase (usually liquid) through which the electrons migrate from one electrode to the other (Figure 1).

Most of the key physical phenomena that determine the charge level and the charging rate of the battery occur in the active materials of the electrodes. This is why, in the last years, nanoelectrochemical measurements have become of special relevance. This type of measurement allows the study of "bare" materials at a single particle level, that is, in the absence of conductive additives or binders that are incorporated in the assembly of battery electrodes. This allows a detailed physicochemical study of the factors that determine the properties of materials to be used as electrodes, and a subsequent more efficient design of them. A comprehensive discussion on single particle measurements and the modeling devoted to understanding these systems has been given very recently in reference [4]. This modeling has been usually based on a continuum approach, that is, the solution of partial differential equations describing the diffusion of lithium-ions within the materials, coupled to a charge-transfer equation describing the insertion/disinsertion of these ions into/out of the active material [5,6,7,8,9,10,11]. Though useful, this model has several approximations [10,11,12], among which an important one is the absence of phase transitions. Thus, comparative calculations with more sophisticated models are necessary to check the impact of this approximation on the predicted behavior. One of these models involves the use of a lattice approach to describe the position of the inserted ions, calculating the evolution of the system according to the Kinetic Monte Carlo (kMC) scheme.

kMC simulations have proven to be a valuable tool for studying the materials of intercalation and surface-level batteries [13,14,15]. Numerous applications have been reported, such as the study of the growth and reaction mechanism of the SEI [16], dynamics of solid and liquid electrolytes [17,18], and post-Li-ion batteries [19,20]. The importance of this methodology lies in the fact that it allows for long-term simulations (on the time scale of the intercalation phenomena) with a microscopic detail of the processes. This is something that is not feasible with other simulation techniques (such as Molecular Dynamics) given the high computational costs. For example, in constant current (galvanostatic) experiments, the charging rate of the electrochemical cell is expressed as the C-rate (C). A value of 1C is the charging rate necessary to charge the cell in one hour, 2C for half an hour, and so on. Thus, the charging of materials, depending on the volume of them, can occur on the scale of seconds to hours. In the literature, few works are performing kMC simulations under constant current conditions [21], which is surprising given that these measurements are routinely used in R&D laboratories for the study of materials. In the present work, we offer a new algorithm for this purpose.

kMC has become a popular technique through the years, especially for simulating diffusion processes for both diffusion on surfaces and in the bulk of materials [22,23,24,25,26,27,28,29,30]. However, if one wants to emulate the physics of intercalation, it is essential to incorporate the kinetics of the insertion/exit of ions at the electrode/electrolyte interface. This is key since the performance of materials depends on the coupling between diffusive and interface phenomena, making this a complex problem to study separately. In previous studies, we have developed a model for this purpose, devoted to the study of lithium-ion intercalation in graphite under potentiostatic and voltammetric conditions [31,32,33]. It has been shown that, in this system, both thermodynamics, which dictates the interactions between inserted ions, and kinetics play a very important role. In that work, to bridge the gap between continuum description and microscopic modeling under potentiostatic signals (validation of kMC simulations), we compared kMC simulations with a theoretical work by Montella [34], showing excellent agreement.

In physics, the use of scalable parameters to reduce a problem from many dimensions to only a few is very popular. With this idea in mind, and seeking to reveal the role of finite diffusion and interface kinetics in the loading capacity and rate of materials, we have recently developed a theoretical framework for the lithiation of materials at the single-particle level under galvanostatic conditions. In essence, the concept of this approach is similar to that of Aoki et al. [35] for voltammetry measurements and that of Montella [34] for potentiostatic measurements. We will go briefly into this theoretical framework in the "Scalable Parameters at Galvanostatic Conditions" Section.

The present work is organized as follows: In the “Theory and Methods” section we show in detail the physics of intercalation for continuum-scale models and for kinetic Monte Carlo, as well as the equations linking the macro and microscopic parameters. The “Results and Discussion” section shows the validity of the kMC model when subjecting the system to different conditions to demonstrate the importance of considering both the statistics and the limitations arising from the system's kinetics in the design of simulations. The galvanostatic scaling model previously developed by us is used as a reference in these tests. Finally, the “Conclusions and Perspectives” section summarizes the most important points of this work and provides an overview of future applications of the proposed methodology.

2. Theory and Methods

The redox reaction taking place in intercalation phenomena can be represented with the following equation.

$$M^{+}+\langle *\rangle +e^{-} \begin{array}{c}{k}_f\\ \rightleftharpoons \\ k_b\end{array} \langle M^{+}\rangle$$

(1)

Here, $M^{+}$ is an ion of an alkaline metal, $\langle *\rangle$ is an intercalation site of the electrode material, located at the electrode/electrolyte interface, $e^{-}$ is the elemental charge, $\langle M^{+}\rangle$ is the intercalated metal ion inside the material and $k_f$ and $k_b$ are the rate constants for the forward and the backward redox process, respectively. This interfacial phenomenon is one of the key factors that influence the intercalation of ions. It depends on several factors, like crystalline structure and exposed crystalline face of the material, type of electrolyte, spurious reactions that may form stable fouling films, the metal ion solvation/desolvation properties, metal ion size.

Once the ion is intercalated, it can move randomly inside the solid without causing substantial changes to its structure. This phenomenon, denoted as ion diffusion, strongly depends on the chemical nature of the material, its crystalline structure, and the kind of intercalated ion. Ion diffusion is characterized by the thermodynamic property denominated diffusion coefficient $D$ and also depends on the number of ions that have been intercalated in the material (occupation) [36]. This shows a complex scenario where interfacial phenomena are coupled to diffusive ones and both are affected by the thermodynamics and kinetics of the problem. This requires adequate modeling of the physics of the processes taking place.

In this sense, kinetic Monte Carlo has proven to be a technique that incorporates the complex framework of the phenomena, allowing the study of the intercalation process in different materials. Curiously, to date, few articles have addressed the intercalation phenomenon incorporating all its complexity. Even more, to the best of our knowledge, only one article has studied the response of the potential of the electrode when applying a constant current signal within the kMC approach, being that this type of experimental signal is commonly used in R&D laboratories for the study of intercalation materials and their application in power sources devices. This is why, in the present work, we focus on the kMC approach to intercalation phenomena in potentiometric experiments and its validation with continuum models.

A complete review of the application of kMC to lithium batteries and other systems can be found in reference [13].

2.1. Continuum Single-Particle Model

Since the results of the kMC simulations will be compared with those of a continuum approach, we will revisit the latter within the denominated single-particle model (SPM). There, the diffusion of ions inside active particles (solid phase) is considered using Fick’s laws for planar diffusion.

$${\displaystyle \frac{\partial x}{\partial t}}=D\left({\displaystyle \frac{{\partial }^{2}x}{{\partial r}^2}}\right)$$

(2)

In this equation $x= {\displaystyle \frac{C_{M^{+}}}{C_{s,max}}}$ is the ion fractional concentration in the solid located ($0\le x\le 1$), where $C_{M^{+}}$ is the concentration of ions and $C_{s,max}$ is the maximum concentration; $D$ is the diffusion coefficient inside the host material and $r$ is the distance from the surface ($r=d$) to the center of the slab ($r=0$).

Von Neumann boundary conditions are considered at the frontier points of the space grid:

$${\left({\displaystyle \frac{\partial x}{\partial r}}\right)}_{r=d}=-{\displaystyle \frac{i_c}{FD{C}_{s,max}}}$$

(3)

$${\left({\displaystyle \frac{\partial x}{\partial r}}\right)}_{r=0}=0$$

(4)

Where $i_c$ is the constant current density applied,

The Butler-Volmer approach, Eq. 5, is used to describe charge transfer kinetics at the electrode/electrolyte interface:

$$i_c=i_0\left\{exp\left[{\displaystyle \frac{\left(1-\alpha \right)F\left(\varphi -{\varphi }^0\left(x_s\right)\right)}{RT}}\right]-exp\left[{\displaystyle \frac{-\alpha F\left(\varphi -{\varphi }^0\left(x_s\right)\right)}{RT}}\right]\right\}$$

(5)

with the exchange current density given by:

$$i_0=F{k}^0{C}_{s,max}{x_s}^{\alpha }{\left(1-x_s\right)}^{1-\alpha }$$

(6)

In these equations, $\varphi$ is the electrode potential, ${\varphi }^0$ the equilibrium potential at the fractional occupation $x_s$, $k^0$ the heterogeneous rate constant for charge transfer, $\alpha$ the transfer coefficient (here considered $\alpha =1/2$ ). We use the subscript $s$ to highlight the fact that the Butler-Volmer equation contains the fractional occupation at the surface of the intercalating material. The symbols $F$, $R$ and $T$ have their usual meaning.

Under constant current conditions, the charging rate ($I_c$) satisfies the following relationship

$$I_c={\displaystyle \frac{Q}{t}}$$

(7)

where $t$ is time and $Q$ is a volumetric capacity in Coulomb. The charging rate is generally expressed as a C-rate ($C_r$), instead of $I_c$, in the analysis of battery materials.

$$C_r=\left|{\displaystyle \frac{I_c{t}_h}{Q_{max}}}\right|$$

(8)

Here $C_r$ is the C-rate, an integer or fractional positive number, $Q_{max}$ is the maximum capacity of the particle (maximum amount of charge that a given material can store) and $t_h$ is the time for an hour. If $Q_{max}$ is expressed in Coulomb units, $t_h$ is considered as 3600 seconds.

As shown, in this model the voltage profiles depend on four parameters $D$, $k^0$, $d$ and $C_r$. It is also important to define the State-of-Charge (SoC). In the present framework, the SoC represents the average ion occupancy within the particle and is calculated as a normalized capacity according to:

$$\mathrm{S}\mathrm{o}\mathrm{C}={\displaystyle \frac{Q}{Q_{max}}}$$

(8b)

A SoC ≈ 0 means an empty particle, while a SoC ≈ 1 corresponds to the maximum occupation of ions in the material.

2.2.1. Scalable Parameters at Galvanostatic Conditions

As described in the previous section, the performance of a SP with planar geometry in the model considered is determined by 4 parameters:$D$, $k^0$, $d$ and $C_r$. This means that the galvanostatic response of a material, illustrated in Figure 2a, depends on the combination of these four parameters. We briefly discuss the nature of this curve in the case of a lithiation taking place close to equilibrium conditions ($i_c$ close to 0 in Eqn. (5)). We represent a negative current transient ($i_c<0,$inset) applied to an empty particle, so that at $t$=0 ($x=Q=\mathrm{S}\mathrm{o}\mathrm{C}=0$), we start to fill the material with Li-ions. The transient starts with a low lattice occupation, at high potentials. As lithiation proceeds, the potential decreases, being the lattice half occupied at the potential ${\varphi }^0\left({\displaystyle \frac{1}{2}}\right)=0.$ Then, when lithiation further proceeds, the electrode potential goes towards minus infinity, as demanded by the Nernst equation in the limit where $x_s\to 1$.

The use of scalable parameters defined from $D$, $k^0$, $d$ and $C_r,$ allows us to reduce the dimensionality of the problem. In the work of reference [7], we proposed two scalable parameters: a kinetic parameter $\mathsf{\Xi }$ (Eq. 9) and a finite-diffusion one $l$ (Eq. 10), which reduces the problem to a two-dimensional one instead of a four-dimensional one.

$$\Xi =k^0{\left({\displaystyle \frac{t_h}{C_{r}D}}\right)}^{1/2}$$

(9)

$$l=d^2\left({\displaystyle \frac{C_r}{D{t}_h}}\right)$$

(10)

Using these parameters, it was shown in reference [4] that the potential profiles under constant current conditions (electrode potential vs SoC representation) depend only on $\mathsf{\Xi }$ and $l$, and on the cut-off potential chosen, ${\varphi }_{cut}$.

We have shown in references [4,7], that for a Langmuir intercalation isotherm and an equilibrium potential ${\varphi }^0=0 \mathrm{V}$, the particle reaches almost the full state of charge, SoC = 0.997 at ${\varphi }_{cut}=-0.15\mathrm{V}$. We will use this ${\varphi }_{cut}$ value in the present simulations.

Master curves or maps for the analysis of single-particle materials used as electrodes in LIBs can be constructed with this approach [7,9,12]. The basis of the methodology is simple: The $\mathsf{\Xi }$ and $l$ space is explored for the same ${\varphi }_{cut}$ and the final SoC reached in each point of this space is saved. The final SoC is called here ${SoC}_{fin}$ (see Figure 2a) and can be represented in a 3D or a color map, where ${SoC}_{fin}=f(\mathsf{\Xi }$, $l)$. A representation of ${SoC}_{fin}$ vs log($\mathsf{\Xi }$) and log($l$) is shown in Figure 2b.

The first thing that stands out is that there is a high-capacity (SoC) zone, greater than 80 % (yellow zone) and a poor SoC zone with less than 10 % (blue zone). In recent work, we have shown that one important feature of this map is that it presents two borderlines that define capacity drop edges due to different phenomena. These borders are marked with red lines in Figure 2. Below log(l)=0.3, a borderline is defined by a straight line with slope ½. This line was related to a SoC drop due to charge transfer limitations at the electrode/solution interface. On the other hand, the vertical line at log(l)=0.3, up to log($\mathsf{\Xi }$)=-1, was related to a SoC drop due to diffusion limitations. This is analyzed in detail in reference work [4].

The use of this map simplifies the analysis of the single particle materials, since we only need two dimensionless parameters to know their final state of charge, at a given cut-off potential. Thus, this map serves as a guide to evaluate how different limitations affect kMC simulations.

2.2. Lattice Kinetic Monte Carlo Simulations

2.2.1. Kinetic Monte Carlo Scheme

A kMC simulation performs a chain of Poisson events that are characterized by a hierarchy of rate constants ($\Gamma$). In this hierarchy, events with higher rates are more likely to occur than those with lower rates [37].

The rate constants calculated for all possible events for a given lattice configuration are stored in a list, from which one is chosen based on a random number ${\xi }_1$​ ranging between 0 and 1. After selecting and allowing the event to proceed, time $t$ is updated by an increment $\Delta t$, determined by an exponential distribution given in Eqn (11) and a second random number ${\xi }_2$ in the same range.

$$\Delta t=-{\displaystyle \frac{ln\left({\xi }_2\right)}{{\sum }_i{\Gamma }_i}}$$

(11)

The rejection-free algorithm used in this work ensures that an event takes place at each time step without any rejection. Once the time is updated, all possible events are calculated from the newly generated configuration, and the procedure is repeated. To finish the simulation, a cut-off criterion is required, which will be detailed later, when we discuss how to adapt the algorithm to potentiometric measurements.

As will be seen in more detail in the next sections, ordinary events considered in an intercalation material are intercalation, deintercalation, and diffusion. Transport processes for the motion of ions through the electrolyte, outside the material, may also be introduced, this is not the case in the present work. Other events or restrictions can be easily added, this is conditioned by the model that one wants to build to emulate a particular chemical physics scenario.

2.2.2. Slab Model and kMC Events

The model of the electrochemical hemi-cell within the SPM is shown in Figure 3. Figure 3a shows a general picture, the slab is located between a current collector and the electrolyte solution. The reference electrode is located at some point in the direction indicated with the dashed blue arrow. One of the faces of the slab is in contact with the electrolyte solution, colored with light grey. A section of the slab is marked with black lines. This section is represented in more detail in Figure 3b. The lattice kinetic Monte Carlo model employed for simulating this section of the slab is detailed there. A 3D lattice-gas was constructed with a simple cubic crystal structure where ions can be placed, this is marked with 1 in Figure 3b and detailed in Figure 3c. The number of intercalation sites at each axis are $N_x$, $N_y$ and $N_z$, and so the total number of sites (that defines the capacity of the section) is calculated as $N_T=N_x x N_y x N_z$. In the case of monovalent ions, the maximum capacity of the slab sectional material is then $Q_{max}=e{N}_T$, where $e$ is the elemental charge in Coulombs. The bulk sites have 6 nearest neighbors, while boundary sites located at y = 0 and y = $L_y$ having 5 nearest neighbors.

The simulations consider diffusion, intercalation and deintercalation events. To mimic the real picture of an electrode material in an electrochemical cell, the following conditions were imposed:

• The electrode was set in contact with a reservoir of ions at the x-z plane located at y = 0, so lithium ions can be intercalated (marked with 2 in Figure 3a) or deintercalated (marked with 3) at the boundary sites of this plane. The concentration of ions in the electrolyte was constant.
• Diffusion of ions was allowed to the first neighboring lattice sites (marked with 4).
• Periodic boundary conditions (PBC) were imposed in the x-direction (marked with 4) and the z-axis (marked with 6).
• The diffusional motion of ions was restricted to the length of the simulation box along the y-axis (marked with 7).
• A metal lithium reference electrode was considered, so that the potential of the working electrode is given by $\varphi =-{\mu }_{Li}/e$ [38].

2.2.3. Rate of Events and Potential Calculation in Galvanostatic kMC Simulations

The rate constant for the events is calculated from the following equation, according to reference [39]:

$$\Gamma =v_0\mathrm{e}\mathrm{x}\mathrm{p} \left[-{\displaystyle \frac{{{\Delta }^{*}}_{\sigma }+\alpha (H_F-H_I)}{k_{B}T}}\right]$$

(12)

In this equation $v_0$ is a pre-exponential factor, $\alpha$ is a symmetry parameter, ${{\Delta }^{*}}_{\sigma }$ is the height of the barrier for the events of particle exchange at the interface $\sigma =int$ and diffusion inside the material $\sigma =diff$. $H_F$ and $H_I$ are the Hamiltonian for the final and initial configurations, respectively; $k_B$ is Boltzmann constant and $T$ temperature.

The Hamiltonian $H$ consists of a potential energy term $U$ and a chemical potential term $\mu$ for the sites located at the interface, i.e., for the intercalation and deintercalation events marked with 2 and 3 in Figure 3b:

$$H=U-\mu$$

(13)

For diffusion events, the Hamiltonian only depends on the potential energy term

$$H=U$$

(13b)

Since, for an intercalation event ${\Gamma }_{in}$, $H_I=0$ and for a deintercalation event ${\Gamma }_{out}$, $H_F=0$, this equation yields:

$${\Gamma }_{in}=v_0\mathrm{e}\mathrm{x}\mathrm{p}\left(-{\displaystyle \frac{{{\Delta }^{*}}_{int}}{k_{B}T}}-{\displaystyle \frac{\alpha (U-\mu )}{k_{B}T}}\right)$$

(14)

and

$${\Gamma }_{out}=v_0\mathrm{e}\mathrm{x}\mathrm{p} \left(-{\displaystyle \frac{{{\Delta }^{*}}_{int}}{k_{B}T}}+{\displaystyle \frac{\alpha (U-\mu )}{k_{B}T}}\right)$$

(15)

Based on charge neutrality, for each ion being intercalated, an electron flows into the system. Conversely, each ion being deintercalated results in an electron getting out of the system. Thus, the constant electronic current applied stems from the balance between the deintercalation and the intercalation of ions. This balance leads to the following equation:

$$I_c=e\left({\sum }_i^{N_{fill}}{\Gamma }_{out}-{\sum }_i^{N_{hole}}{\Gamma }_{in}\right)$$

(16)

where $N_{fill}$ is the number of sites occupied by ions and $N_{hole}$ is the number of holes at the interface. Note that intercalation/deintercalation involves a negative/positive current. This is the usual convention in current electrochemistry. Replacing (14) and (15) in Eqn. (16) yields:

$$I_c=e\left\{{\sum }_i^{N_{fill}}{v}_0\mathrm{e}\mathrm{x}\mathrm{p}\left(-{\displaystyle \frac{{{\Delta }^{*}}_{int}}{k_{B}T}}+{\displaystyle \frac{\alpha \left(U_i-\mu \right)}{k_{B}T}}\right)-{\sum }_i^{N_{hole}}{v}_0\mathrm{e}\mathrm{x}\mathrm{p} \left(-{\displaystyle \frac{{{\Delta }^{*}}_{int}}{k_{B}T}}-{\displaystyle \frac{\alpha \left(U_i-\mu \right)}{k_{B}T}}\right)\right\}$$

(17)

Splitting the exponential terms and taking $v_0$ and the exponential term containing ${{\Delta }^{*}}_i$ outside the curly brackets:

$$I_c=e{v}_0\mathrm{e}\mathrm{x}\mathrm{p} \left(-{\displaystyle \frac{{{\Delta }^{*}}_{int}}{k_{B}T}}\right) \left\{{\sum }_i^{N_{fill}}\left[\mathrm{e}\mathrm{x}\mathrm{p}\left({\displaystyle \frac{\alpha U_i}{k_{B}T}}\right)\mathrm{e}\mathrm{x}\mathrm{p}\left(-{\displaystyle \frac{\alpha \mu }{k_{B}T}}\right)\right]-{\sum }_i^{N_{hole}}\left[\mathrm{e}\mathrm{x}\mathrm{p}\left(-{\displaystyle \frac{\alpha U_i}{k_{B}T}}\right)\mathrm{e}\mathrm{x}\mathrm{p}\left({\displaystyle \frac{\alpha \mu }{k_{B}T}}\right)\right]\right\}$$

(18)

In a Langmuirian system, all the sites are energetically equivalent and we may set $U=0$ . Further assuming $\alpha =0.5$, equation S7 reduces to:

$$I_c=e{v}_0\mathrm{e}\mathrm{x}\mathrm{p} \left(-{\displaystyle \frac{{{\Delta }^{*}}_{int}}{k_{B}T}}\right) \left\{N_{fill}\mathrm{e}\mathrm{x}\mathrm{p}\left(-{\displaystyle \frac{\mu }{2k_{B}T}}\right)-N_{hole}\mathrm{e}\mathrm{x}\mathrm{p}\left({\displaystyle \frac{\mu }{{2k}_{B}T}}\right)\right\}$$

(19)

At this point, a numerical method can be used to find roots to get the $\mu$ value, after calculating the summation of the exponential terms for each time step. We applied the Brent’s method for this aim [40]. Then, the potential of the working electrode is obtained from $\varphi =-\mu /e$.

The algorithm constructed for simulating galvanostatic kMC simulations consists of the following steps, detailed in the workflow of Scheme 1.

1)

Set an initial configuration and define a constant current signal (or C-rate).

2)

Estimate the chemical potential with a root-finding method using Eqn. (18) or (19).

3)

Calculate all rate constant events with Eqn. 12.

4)

Generate a random number ${\xi }_1$.

5)

Choose one event using ${\xi }_1$.

6)

Perform the selected event.

7)

Generate a second random number ${\xi }_2$.

8)

Calculate the time increment with Eqn. 11.

9)

Update the time.

10)

Repeat steps from 2 to 9 until reaching the cut-off electrode potential.

The output quantities for each simulation are averaged for the mean average time between two time points $t_j=(t_i+t_{i+1})/2$ provided that $t_i<t\le t_{i+1}$.

2.3. Calculating Continuum Parameters from Microscopic Quantities

2.3.1. Diffusion Coefficient

In a 3D lattice, the diffusion coefficient may be calculated from [41]:

$$D={\displaystyle \frac{1}{6}}{\Gamma }_{tot}{a}^2$$

(20)

where $a$ is a characteristic diffusion length and

$${\Gamma }_{tot}={\sum }_i^n{\Gamma }_i{\Gamma }_i=v_0\mathrm{e}\mathrm{x}\mathrm{p}\left(-{\displaystyle \frac{{{\Delta }^{*}}_{diff}}{k_{B}T}}\right)$$

(21)

In a simple cubic lattice, the $n$ = 6 nearest neighbors present the same rate for diffusion if no interactions are considered and so Eqn. (20) yields:

$$D=v_0{a}^2\mathrm{e}\mathrm{x}\mathrm{p}\left(-{\displaystyle \frac{{{\Delta }^{*}}_{diff}}{k_{B}T}}\right)$$

(22)

This is the same expression derived by [27,29,42] under Langmuir conditions.

For isotropic systems, such as cubic crystals, the displacements in x, y and z are the same, therefore:

$$\langle x^2\rangle =\langle y^2\rangle =\langle z^2\rangle ={\displaystyle \frac{1}{3}}\langle R^2\rangle$$

(23)

Where $\langle R^2\rangle$ is the mean square displacement. Then, the diffusion coefficient in the y-direction is

$$D_y=D=v_0{a}^2\mathrm{e}\mathrm{x}\mathrm{p}\left(-{\displaystyle \frac{{{\Delta }^{*}}_{diff}}{k_{B}T}}\right)$$

(24)

The continuum diffusion length $d$ was calculated with the number of lattice sites in this direction and with the jump distance, i. e., $d=L_y=N_{y}a$. The surface in contact with the electrolyte was calculated as $S=N_x{N}_z{a}^2.$

2.3.2. Charge Transfer Rate Constant

According to the literature [43], the rate constants of the forward and reverse processes of a redox reaction like that given in Eqn. (1) are derived from the following equations:

$$k_f=A_f\mathrm{e}\mathrm{x}\mathrm{p}\left(-{\displaystyle \frac{{\Delta G}_c^{\ne }}{RT}}\right)$$

(25)

$$k_b=A_b\mathrm{e}\mathrm{x}\mathrm{p}\left(-{\displaystyle \frac{{\Delta G}_a^{\ne }}{RT}}\right)$$

(26)

Where ${\Delta G}_c^{\ne }$ and ${\Delta G}_a^{\ne }$ are the free energy barriers for reduction and oxidation, respectively (Figure 4). The changes in these barriers with the displacements from the equilibrium potential ${\varphi }^0$, can be expressed as:

$${\Delta G}_c^{\ne }={\Delta G}_{0,c}^{\ne }-(1-\alpha )(\varphi -{\varphi }^0)$$

(27)

$${\Delta G}_a^{\ne }={\Delta G}_{0,a}^{\ne }+\alpha (\varphi -{\varphi }^0)$$

(28)

Where ${\Delta G}_{0,c}^{\ne }$ y ${\Delta G}_{0,a}^{\ne }$ corresponds to the free energy barriers when the electrode potential $\varphi ={\varphi }^0$ (Figure 4). Replacing in Eqns. (25) and (26), yield:

$$k_f=A_f\mathrm{e}\mathrm{x}\mathrm{p}\left(-{\displaystyle \frac{{\Delta G}_{0,c}^{\ne }}{RT}}\right)\mathrm{e}\mathrm{x}\mathrm{p}\left({\displaystyle \frac{-\alpha F(\varphi -{\varphi }^0)}{RT}}\right)$$

(29)

$$k_b=A_b\mathrm{e}\mathrm{x}\mathrm{p}\left(-{\displaystyle \frac{{\Delta G}_{0,a}^{\ne }}{RT}}\right)\mathrm{e}\mathrm{x}\mathrm{p}\left({\displaystyle \frac{(1-\alpha )F(\varphi -{\varphi }^0)}{RT}}\right)$$

(30)

and when the working electrode potential is different to that of the equilibrium $\varphi$.

From these las equation emerges that when equilibrium is reached ($\varphi ={\varphi }^0$), $k_f=k_b=k^0$, and so Eqns. (29) and (30) yield:

$$k_f=k^0\mathrm{e}\mathrm{x}\mathrm{p}\left({\displaystyle \frac{-\alpha F(\varphi -{\varphi }^0)}{RT}}\right)$$

(31)

$$k_b=k^0\mathrm{e}\mathrm{x}\mathrm{p}\left({\displaystyle \frac{(1-\alpha )F(\varphi -{\varphi }^0)}{RT}}\right)$$

(32)

From last equations arises that

$$k^0=A_f\mathrm{e}\mathrm{x}\mathrm{p}\left(-{\displaystyle \frac{{\Delta G}_{0,c}^{\ne }}{RT}}\right)=A_b\mathrm{e}\mathrm{x}\mathrm{p}\left(-{\displaystyle \frac{{\Delta G}_{0,a}^{\ne }}{RT}}\right)$$

(33)

This gives the Butler Volmer equation at the macroscopic level [44]:

$$I=SF{c}_{max}{k}^0\left[x\mathrm{e}\mathrm{x}\mathrm{p}\left({\displaystyle \frac{(1-\alpha )F(\varphi -{\varphi }^0)}{RT}}\right)-(1-x)\mathrm{e}\mathrm{x}\mathrm{p}\left(-{\displaystyle \frac{\alpha F(\varphi -{\varphi }^0)}{RT}}\right)\right]$$

(34)

An equivalence can be established with the microscopic description given in Section 2.2.3. From the comparison between Eqns. (18) and (34), it emerges that the term

$${\Gamma }^0=v_0\mathrm{e}\mathrm{x}\mathrm{p}\left(-{\displaystyle \frac{{{\Delta }^{*}}_{int}}{k_{B}T}}\right)$$

(35)

is equivalent to the rate constant $k^0$, but in s^-1 units.

In a previous work [32], we derived an expression to estimate the rate constant of the Butler-Volmer equation from the microscopic rate constants. It is just a matter of multiplying the microscopic rate constant by a volume/area parameter, i.e.:

$$k^0={\Gamma }^0{d}^0={d^{0}v}_0\mathrm{e}\mathrm{x}\mathrm{p}\left(-{\displaystyle \frac{{{\Delta }^{*}}_{int}}{k_{B}T}}\right)$$

(36)

where

$$d^0={\displaystyle \frac{volume}{surface}}$$

(37)

To calculate $d^0$ it should be considered the unit cell volume and the unit cell surface that contains one single site of the lattice.

3. Results and Discussion

The code for kMC simulations was written in C++ and was parallelized with OpenMP. Continuum SPM simulations were performed using the open-source software galpynostatic, with the module simulation, available at https://github.com/fernandezfran/galpynostatic [12]. In all cases, Langmuir intercalation isotherms were considered. In kMC simulations, the diffusion jump distance of the simple cubic lattice was kept constant at $a= 1\mathrm{\AA }$, then $d^0$= $1\mathrm{\AA }$ and the pre-exponential factor was as assumed as $v_0=1\mathrm{x}{10}^{13}{\mathrm{s}}^{-1}$, this is a value commonly used for intercalation materials that was calculated for graphite [45].

The first validation test of the galvanostatic kMC code was to check the expected linear relationship between the SoC and the simulated time (see Eqn. (7) and (8b)), typical for galvanostatic experiments. As observed in Figure 5, kMC simulations (represented with symbols) follow the expected theoretical trends (represented with lines), where an increase in current density, from -100 mA/cm² to -1000 mA/cm², increases the slope of the linear behaviour. In these simulations, only the current density was varied keeping constant the rest of the parameters.

Next, the agreement between kMC and continuum potential/SoC profiles obtained after applying a constant current signal was tested. With this purpose, kMC barriers were directly estimated from the macroscopic parameters at the continuum scale. This was done from the relationships provided in Section 2.3. For example, let us consider two intercalation materials A and B defined by the set of continuum parameters ($D$, $k^0$, $d$ and $C_r$) given in Table 1, with cubic simulation boxes of $N_T=20 \mathrm{x}20\mathrm{x}20$ for system A and $30\mathrm{x}30\mathrm{x}30$ for B. To calculate the energy barriers for diffusion and intercalation, Eqs. (24) and (36) were used, respectively. The values obtained for the barriers ${{\Delta }^{*}}_{int}$ and ${{\Delta }^{*}}_{diff}$ are provided in Table 1.

Figure 6 shows the potential/Soc profiles for materials A and B obtained from numerical simulations (blue line) and kMC calculations (yellow dots) for one kMC run. The kMC simulations show fluctuating behaviour around the continuum value. However, despite the restricted statistics of a single kMC run, the two methods show excellent agreement.

To check that the parameters selected correspond to relevant cases of the problem under consideration, the scalable parameters $\mathsf{\Xi }$ and $l$ for these materials were estimated with Equations (9) and (10). The calculated dimensionless values are detailed in Table 1 and at the top of Figure 6. This allows the materials to be located on the galvanostatic map, introduced in Figure 2b. In the two cases presented, the materials reach a similar SoC when the cut-off potential is reached; however, material A allows a charging rate three times greater than material B, which would be beneficial for having faster-charging batteries.

From now on, we will report the values of the scalable parameters $\mathsf{\Xi }$ and $l$ instead of detailing the whole data set, as well as the information necessary for the analysis.

3.1. Number of Samples and Simulation Box Size

In this section, we will analyze how the number of independent samples (number of kMC runs) and the size of the simulation box impact the potential profiles.

Figure 7 shows the comparison between the potential/SoC profiles obtained from the kMC simulations with different numbers of statistical samples and the continuum simulations for the same $\Xi$ and $l$ point (log($\Xi )$=-0.19 and log($l)$=-1.147). As shown before, a single kMC sample already shows a random dispersion close to the Langmuir prediction, Figure 7a. As the number of samples increases, the disperse behaviour becomes smoother and the results closely coincide. This is better shown in Figure 8d when plotting the residuals for different numbers of samples. As observed, the results with 32 and 64 samples are practically the same. Figure 7b, shows the average curve (black line) with 32 samples. The standard deviation of the potential was calculated (blue shadow area). The deviation increases for larger SoC values. This latter fact is understandable, since as the SoC increases, the number of available sites for particle insertion becomes smaller, which makes the outcome for the potential more dispersed.

Figure 8 shows an analysis of the dispersion for the predicted potential as a function of box size using 32 samples, for the same $\mathsf{\Xi }$ and $l$ point as that of the previous case. It can be seen that increasing the number of sites leads to a decrease in the fluctuations of the potential, making the standard deviation much lower. Therefore, increasing the size of the box helps to reduce the number of samples required to obtain a converged isotherm. This effect becomes more remarkable at larger SoCs, where the number of free surface sites where the ions may be intercalated sinks.

A residual analysis was performed for these box sizes varying the number of samples. The residual between the continuum potential/SoC curves and the kMC curves is shown in Figure 8d, e, and f. It is clear that with increasing box size, the residual values decrease for all sample numbers, leading to improved convergence of the voltage profiles. To help understand the sample and lattice size effects, an error estimation using Equation (38) was made. The error of the kMC predictions relative to the continuum simulations was estimated from the differences:

$$error={\displaystyle \frac{\sum _i\left|{\varphi }_i^{kMC}-{\varphi }_i^{cont}\right|}{\sum _i\left|{\varphi }_i^{cont}\right|}}\mathrm{x}100$$

(38)

where the upper indices $kMC$ and $cont$ correspond to the results of the kMC and continuum simulations and the index $i$ runs over simulated kMC points. Table 2 shows the errors as a function of the number of samples and box size.

3.2. Limitation by the Number of Surface Sites

In most of the cases displayed above, it was observed that the kMC potential profiles present a drop that reaches the potential cut-off (-0.15 V) at a SoC value that is lower than the final SoC of the numerical simulations. Although not obvious in the representations shown before, this drop is evident in Figure 9, for three different discretizations of the surface at the same ($\Xi$,$l$) point, log($\Xi )$=-0.19 and log($l)$=-1.147. At this ($\Xi$,$l$) point, labelled as “a” in Figure 11, the rate-controlling process is charge transfer, according to the details given in Figure 2b. As observed, something is preventing the SoC of kMC from reaching its numerical counterpart.

Figure 9 shows results for three different discretizations of the simulation box , keeping constant the number of sites $N_y=20$, and varying the number of surface sites $N_x$x $N_z$ = 25 (5x5), 100 (10x10) and 400 (20x20), for the same $\mathsf{\Xi }$ and $l$ point. As observed, for a single kMC sample, an increase in the number of sites in contact with the electrolyte causes a shift in the potential profile towards higher SoCs, approaching the results of the continuum simulations. This is understood from the previous comments: if the number of lattice sites at the interface is small, the fluctuations in the electrode potential will be larger with an increasing number of occupied surface sites. This results in a potential jump that exceeds the cutoff potential at a SoC value lower than the continuum value.

It is important to remark that changing the number of sites $N_T$, and so $Q_{max}$, leads to different continuum parameters for the same $\Xi$ and $l$ point, according to the equations seen in Section 2.1. These parameters are given in Table 3.

In summary, there is a discretization effect due to the finite number of sites at the interface where ion exchange occurs. Therefore, the interface between the electrolyte and the electrode must contain enough sites to approach the continuum behavior when modeling galvanostatic profiles with kMC.

Figure 9. Numerical and kMC voltage profiles for three different areas at the electrode/electrolyte interface 5x5 (25 sites), 10x10 (100 sites) and 20x20 (400 sites), for the same number of sites at the y-axis, 20 sites. The point considered was log($\Xi )$=-0.19 and log($l)$=-1.147 in all cases. Dashed lines were added with the same colour of the simulation box considered to represent the sudden decrease in potential that crosses the potential cut-off of -0.15 V.

Figure 9. Numerical and kMC voltage profiles for three different areas at the electrode/electrolyte interface 5x5 (25 sites), 10x10 (100 sites) and 20x20 (400 sites), for the same number of sites at the y-axis, 20 sites. The point considered was log($\Xi )$=-0.19 and log($l)$=-1.147 in all cases. Dashed lines were added with the same colour of the simulation box considered to represent the sudden decrease in potential that crosses the potential cut-off of -0.15 V.

3.3. Limitation by Finite-Size Diffusion

Another source of discrepancy between the kMC potential profiles and the continuum curves appears when mass transport becomes rate-controlling. An analysis of the convergence of the kMC and the theoretical isotherms by varying the number of sites in the y direction, say $N_y$, was carried out. For this, a ($\Xi ,l$) point corresponding to [log($\Xi ),\mathrm{l}\mathrm{o}\mathrm{g}\left(l\right)]$=[-0.19, 0.2] was used, corresponding to the case where the SoC drop is associated with diffusional limitations (see lines in Figure 2 and point c in Figure 11 and discussion in reference [4]).

Figure 10 compares kMC and continuum results. As the number of points used to discretize in the y direction increases from $N_y$ = 20 (Figure 10a) to $N_y$ = 40 (Figure 10b) and $N_y$ = 60 (Figure 10c), for the same surface area sites ($N_x$x $N_z$ = 60 x 60), the relative difference between the curves decreases. The errors between kMC and continuum curves were calculated with Eqn. (38) with the maximum number of kMC samples (64) and are indicated in each of the panels of the figure.

We remark again that changing the number of sites of the simulation box would lead to different continuum parameters for the same $\mathsf{\Xi }$ and $l$ point, according to Equations seen in Section 2.1. This is shown in Table 3.

Table 3. Parameters corresponding to data shown in Figure 9 and Figure 10.

$$\begin{gathered} \mathbf{Box} \mathbf{size} \\ \left({\mathit{N}}_{\mathit{x}}\mathbf{x}{\mathit{N}}_{\mathit{y}}\mathbf{x}{\mathit{N}}_{\mathit{z}}\right) \end{gathered}$$	$\mathit{d}$ [Å]	${\mathit{C}}_{\mathit{r}}$	${\mathit{Q}}_{\mathit{m}\mathit{a}\mathit{x}}$ [C]	${\mathit{D}}_{\mathit{y}}$ [cm2 / s]	${\mathit{k}}^0$ [cm / s]	${\mathsf{\Delta }}_{\mathit{d}\mathit{i}\mathit{f}\mathit{f}}^{\mathit{*}}$ [eV]	${\mathsf{\Delta }}_{\mathit{i}\mathit{n}\mathit{t}}^{\mathit{*}}$ [eV]	$\left[\mathbf{log}\right(\mathit{\Xi }),$$\mathbf{log}\left(\mathit{l}\right)]$
5x20x5	20	22454	8.01x10-17	3.5x10-12	3.02x10-6	0.5	0.62	[-0.19, -1.147$]$
10x20x10	20	22454	3.2x10-16	3.5x10-12	3.02x10-6	0.5	0.62	[-0.19, -1.147$]$
20x20x20	20	22454	1.15x10-15	3.5x10-12	3.02x10-6	0.5	0.62	[-0.19, -1.147$]$
60x20x60	20	203381	1.15x10-14	3.5x10-12	2.23x10-5	0.5	0.57	[-0.19, 0.2$]$
60x40x60	40	50845	2.31x10-14	3.5x10-12	1.14x10-5	0.5	0.59	[-0.19, 0.2$]$
60x60x60	60	22597	3.46x10-14	3.5x10-12	7.42x10-6	0.5	0.60	[-0.19, 0.2$]$

Table 3. Parameters corresponding to data shown in Figure 9 and Figure 10.

$$\begin{gathered} \mathbf{Box} \mathbf{size} \\ \left({\mathit{N}}_{\mathit{x}}\mathbf{x}{\mathit{N}}_{\mathit{y}}\mathbf{x}{\mathit{N}}_{\mathit{z}}\right) \end{gathered}$$	$\mathit{d}$ [Å]	${\mathit{C}}_{\mathit{r}}$	${\mathit{Q}}_{\mathit{m}\mathit{a}\mathit{x}}$ [C]	${\mathit{D}}_{\mathit{y}}$ [cm2 / s]	${\mathit{k}}^0$ [cm / s]	${\mathsf{\Delta }}_{\mathit{d}\mathit{i}\mathit{f}\mathit{f}}^{\mathit{*}}$ [eV]	${\mathsf{\Delta }}_{\mathit{i}\mathit{n}\mathit{t}}^{\mathit{*}}$ [eV]	$\left[\mathbf{log}\right(\mathit{\Xi }),$$\mathbf{log}\left(\mathit{l}\right)]$
5x20x5	20	22454	8.01x10-17	3.5x10-12	3.02x10-6	0.5	0.62	[-0.19, -1.147$]$
10x20x10	20	22454	3.2x10-16	3.5x10-12	3.02x10-6	0.5	0.62	[-0.19, -1.147$]$
20x20x20	20	22454	1.15x10-15	3.5x10-12	3.02x10-6	0.5	0.62	[-0.19, -1.147$]$
60x20x60	20	203381	1.15x10-14	3.5x10-12	2.23x10-5	0.5	0.57	[-0.19, 0.2$]$
60x40x60	40	50845	2.31x10-14	3.5x10-12	1.14x10-5	0.5	0.59	[-0.19, 0.2$]$
60x60x60	60	22597	3.46x10-14	3.5x10-12	7.42x10-6	0.5	0.60	[-0.19, 0.2$]$

From the present analysis, it can be deduced that to obtain the correct voltage profiles, provided that there is some degree of diffusional control, the indicated number of sites in the direction of diffusion must be considered.

Figure 10. Impact of the box size on the kMC convergence in the diffusion control zone. As the box depth increases the error between the curves diminishes. .

Figure 10. Impact of the box size on the kMC convergence in the diffusion control zone. As the box depth increases the error between the curves diminishes. .

3.3. Validation of kMC with the Galvanostatic Map

As discussed in Section 2.1.1, different regions in the galvanostatic map are related to different types of control over the kinetics of intercalation. Therefore, this map serves as a guide to check that the kMC simulations reproduce the behavior of continuum models under different types of kinetic limitations. Thus, an analysis of the kMC behaviour at different [log($\Xi ),$ log($l\left)\right]$ couples was made.

Figure 11f shows different points on the map selected to test the kMC simulations, with the corresponding galvanostatic transient given in Figure 11a–e. The kMC models for each of the points were constructed taking into account the limitations previously evaluated in Section 3.1 and Section 3.2. All kMC voltage/SoC profiles show an excellent agreement with the continuum calculations.

Figure 11. Different kMC curves simulated with 64 samples compared to the continuum simulations and their respective position on the SoC map.

Figure 11. Different kMC curves simulated with 64 samples compared to the continuum simulations and their respective position on the SoC map.

The reference point to analyze the simulations was point a in Figure 11f (log($\Xi )$=-0.19, log($l)$=-1.147). At this point, the maximum SoC is practically reached since there is no kinetic limitation, due to a large $k^0/{C_r}^{1/2}$ ratio, and diffusion kinetics is relatively fast, due to a small $d^2/D$ ratio. The kMC-point-a.mp4 video in the Supplementary Information (SI) shows the evolution of this simulation over time. In all the videos the interphase electrode/electrolyte is located on the left-hand side and the diffusion limit is on the right-hand side of the screen. Furthermore, the ions located at the interface are represented with a different colour.

When the ordinate is fixed at log($\Xi )$=-0.19 and log($l)$ is increased from the reference point, the system enters into the diffusion control zone of the map, transients of Figure 11b and c. The video kMC-point-c.mp4 in the SI shows the evolution of the system at point c over time. It is remarkable how the intercalated ion front does not reach the inner part of the electrode.

On the other hand, when the abscissa is fixed at log($l)$=-1.147 and log($\Xi )$ is decreased from point a, the map approaches the zone regulated by charge transfer, with the transients given in Figure 11d,e. In the file kMC-point-e.mp4 in SI, the evolution of the frames at point e over time shows that although the ions cover the entire length of the slab, the slab is not fully occupied with Li-ions. This is so because there is a slow exchange of ions at the interface.

This is a key result for our work since it demonstrates that the kMC model works correctly under different kinetic regimes of the lithiation process, which consolidates the robustness of the kMC model. It also demonstrates the validity of the equations that connect the parameters of the continuum with those of the microscopic level, on the basis of the galvanostatic map developed by us previously [4,7] as a guide to easily deal with a problem with numerous variables to handle.

4. Conclusions and Perspectives

We have shown that kinetic Monte Carlo simulations can replicate the physics of macro-scale intercalation systems under galvanostatic conditions. This bridges the gap between the electrochemistry of battery materials at macro and micro levels, allowing us to perform computational experiments of potential profiles with microscopic details of the events.

An intercalation Langmuir isotherm was used for both numerical and kMC simulations, as well as the kMC algorithm considering a galvanostatic signal. Also, we have shown how the kMC parameters are linked to those of the continuum equations. In particular, how the microscopic parameters can be controlled from that of the continuum level.

We have shown the importance of performing tests to avoid problems due to the discretization of the systems, as well as to evaluate the number of statistical samples necessary to obtain converged potential/state of charge profiles. The validity of the kMC simulations was tested using a theoretical framework developed previously by us, consisting of maps that use two scalable parameters that determine the galvanostatic profiles. The good agreement between both approaches demonstrated the validity of the Monte Carlo modelling.

This work opens a door to explore the kinetics of alkali-metal ion intercalation materials using more realistic modeling to represent the interaction among inserted ions. In particular, this approach may allow us to address problems that arise due to the limitations of continuum models to cover some physical aspects of intercalation materials, for example, phase coexistence, entropic factors arising from the configurations of the intercalated ions, the solid-state mass transport inside the electrodes, the geometry of the solid material, just to name a few factors.