An electrochemical cell has two electrodes separated by an electrolyte solution, as schematically shown in Figure 1. The electric potential difference between the two electrodes, denoted V_cell, can be modulated with a potentiostat at one’s disposal. By varying V_cell , one gets a handle on controlling the difference in the electrochemical potential of electrons between the two electrodes. Electrons flow from the electrode with the higher electrochemical potential (lower electric potential), via the potentiostat, to the other side. Let us consider for the moment ideally polarizable electrodes. The resultant electron flow cannot cross the electrode-electrolyte interface (EEI) regardless of V_cell. Therefore, excess or deficient electrons, corresponding to a net negative or positive surface charge, are distributed over a skin layer of several angstroms (Å) thick on the electrode surface, which is a consequence of quantum-mechanical behaviors of electrons. Counterions are attracted to and coions are repelled from the electrode surface via electrostatic interactions, as illustrated in Figure 1. In addition, the solvent molecules adjust their orientation according to the strong electric field generated by the net surface charge. These coupled phenomena occur in a non-electroneutral region of a few nanometer (nm) thick, which is termed as the electric double layer (EDL).

In most cases, the electron flow can cross the EEI, resulting in oxidation or reduction of solution species near the electrode. The electron transfer rate depends on the distributions of the electric potential, the concentrations of reactant and product, and the dielectric polarization of the solvation environment. Therefore, the structure and properties of the EDL are important factors influencing interfacial electron transfer reactions.

The distributions of the electric potential, the ion concentrations, and the solvent orientation in the EDL are dictated by the excess surface charge density, denoted σ_M. In an electrochemical cell, the two electrodes have σ_M of the same magnitude but opposite signs, because the full cell must be electroneutral. σ_M can be varied as a function of V_cell, or as long as a single electrode is considered, σ_M is a function of the electric potential, ϕ_M, of the considered electrode. The relation between σ_M and ϕ_M is called surface charging relation. Determining the σ_M - ϕ_M relation is an essential task in modelling the EDL and electrochemical reactions therein.

A schematic illustration of the surface charging relation of the EDL is given in Figure 2. The electric potential in the bulk solution, ϕ_S, is taken as the reference. The electric potential in the electrode bulk, ϕ_M, can be adjusted at will. There is a potential drop, χ, on the electrode surface, which is caused by electron spillover from the solid electrode into the electrolyte solution^[1]. χ is related to the electron density of the solid electrode, therefore, its value depends on how many electrons are included in the calculation. More specifically, χ is much higher if more electrons of the metal are considered explicitly. For a given metal, χ varies slightly with ϕ_M, and we assume that χ is a constant in the subsequent analysis.

Consider first the case of ϕ_M - ϕ_S< χ, as shown in Figure 2(A). The electric potential in the electrolyte solution near the electrode surface is lower than ϕ_S. Therefore, cations are accumulated in the EDL because they have a lower energy there. Contrarily, anions are depleted. Hence, a positive net charge is stored in the EDL, which is accompanied by excess electrons on the metal surface and σ_M < 0. If ϕ_M is increased, the electric potential in the EDL becomes more positive because we assume a constant χ for the electrode. This means that less cations are accumulated and less anions are repelled in the EDL. Therefore, σ_M becomes less negative and shifts towards more positive values.

For the case of ϕ_M - ϕ_S = χ, as illustrated in Figure 2(B), the electric potential in whole electrolyte solution is zero as we have set the potential reference in the bulk solution. There is no excess charge in the EDL or on the electrode surface, i.e., σ_M = 0. This particular value of ϕ_M is the potential of zero charge (pzc)^[2]. For the case of ϕ_M - ϕ_S > χ, shown in Figure 2(C), anions are accumulated and cations are depleted in the EDL, i.e., σ_M > 0. Overall, the σ_M - ϕ_M relation shows a monotonically increasing trend for ordinary EDLs. From this relation, we obtain the differential double-layer capacitance C_dl,

For EDLs with chemisorption, the surface charging relation can be nonmonotonic. In Figure 2, we give such an example, where chemisorption occurs in the high potential range. The chemisorbates are usually partially charged due to the partial charge transfer^[3]. The charged chemisorbates, together with the compensating charge located on metal surface atoms, give rise to a surface dipole moment, which is termed the chemisorption-induced surface dipole moment, μ_chem. μ_chem introduces an additional contribution to the potential drop on the electrode surface, Δχ. Even at a high potential ϕ_M - ϕ_S > χ, the electrode surface with chemisorbates could be negatively charged due to the additional negative potential drop Δχ. Consequently, the surface charging relation could be nonmonotonic with a second pzc in the high potential region for electrocatalytic interfaces^{[4, 5]}. According to its definition, C_dl is negative in the nonmonotonic region of the σ_M - ϕ_M relation. Note that in the presence of chemisorbates, the definition of σ_M should be varied accordingly, which contains not only the excess charge on the electrode surface, but also the net charged carried on the charged adsorbates. In this regard, it is better to be denoted σ_free.

We wish to provide a tutorial for physical modelling of the EDL under both equilibrium and nonequilibrium conditions. Both time-domain and frequency-domain responses are modelled under nonequilibrium conditions. The latter is in fact the electrochemical impedance response. As a tutorial, this paper includes a systematic exposition of relevant models with mathematical details provided. We also provide simulation scripts of these models in the supporting information.

The remaining parts of this paper are organized as follows. We first provide a brief history of the EDL theory. Then, we introduce EDL models under equilibrium conditions, including the Gouy-Chapman-Stern model, then the Bikerman model that takes ion size effects into account, and finally a modified Bikerman model that considers asymmetric ion size effects. Next, we derive nonequilibrium EDL models from a grand potential of the EDL. Afterwards, we present the basics of EIS, and demonstrate how to derive an EIS model from the nonequilibrium EDL models.

A quantitative determination of the σ_M = f(ϕ_M) relation requires a physicochemical model for the EDL. Figure 3 summarizes milestones in the evolution of EDL modelling and simulations. Helmholtz (1879) viewed the EDL as a planar plate capacitor with a constant double-layer capacitance (C_dl) and a linear potential distribution in the space between two plates^[6]. Different from Helmholtz who assumed a rigid lining up of counterions, Gouy and Chapman considered the diffuse nature of counterions in the electrolyte solution in 1910^{[7, 8]}, ten years before Debye and Hückel. In the Gouy-Chapman model, the electrolyte solution is viewed as a cloud of point ions embedded in a dielectric continuum. The distributions of the electric potential and ion concentrations are governed by the Poisson-Boltzmann equation. The Gouy-Chapman model is limited to very dilute solutions at slightly charged interfaces due to its foundation on the point charge assumption. At highly charged interfaces or when ϕ_M is shifted far away from the pzc, the Gouy-Chapman model results in an unphysically high concentration of counterions near the electrode surface. Under such scenarios, the gap between the electrode surface and the midplane of the diffuse layer, denoted d₀, becomes very small, leading to the phenomenon of capacity catastrophe, namely, C_dl grows toward infinity.

he Helmholtz plane (HP)^[9]. This way, regardless of the magnitude of surface charge density on the electrode, d₀ has a lower limit, thus turning the catastrophic growth of C_dl when ϕ_M is shifted away from the pzc into a leveling off region. Bikerman furthered the consideration of ion size effects using a lattice-gas model^[10]. There is a delicate difference regarding the consideration of ion size in Stern’s and Bikerman’s treatments. Although d₀ is constrained in the Stern model, the Poisson-Boltzmann equation is inherited. This means that the ion concentration at the HP and in the diffuse layer can be infinite in the Stern model. In contrast, in the Bikerman model, the ion concentration at the HP and in the diffuse layer has a finite upper limit determined by the lattice size. Consequently, in the Bikerman model, d₀ first narrows down due to counterion crowding and then expands due to counterion overcrowding, when ϕ_M is shifted away from the pzc. Consequently, a camel-shaped double-layer capacitance profile is usually obtained. In highly concentrated solutions, a bell-shaped double-layer capacitance profile is obtained, as d₀ always increases when ϕ_M is shifted away from the pzc.

Grahame extended Stern’s idea in the presence of specific adsorption of ions^[11]. He divided the HP into an inner HP (IHP) where specifically adsorbed ions reside and an outer HP (OHP) where solvated counterions reside. It was implicitly assumed that the specifically adsorbed ions retain the charge they have in the bulk solution, which was later corrected by the concept of partial charge transfer by Lorenz and Salie in 1961^[12]. Moreover, as a first approximation, the potential distribution in the inner layer is considered to be linear. The potential difference across the inner layer is composed of two contributions. One is caused by the net surface charge on the electrode surface, denoted σ_M. The other is caused by the charge carried by the specifically adsorbed ions. Grahame calculated the differential capacitance of the inner layer (C_IHP) as a function of σ_M^[13]. C_IHP is asymmetric and humped with a maximum at positive σ_M. Grahame’s results aroused wide interests among theorists, suggesting two new lines of EDL modelling, namely, description of water dipoles at the IHP and description of metal electrons, which are detailed below.

Grahame spent his sabbatical year in 1958/59 in Britain and had a fruitful collaboration with Roger Parsons. This British visit disseminated his experimental findings among physicists at Cambridge. Watts-Tobin and Mott tried to interpret the rise of C_IHP for anodic polarizations and the hump of C_IHP at a ϕ_M slightly positive to the pzc^[14]. The former phenomenon had been controversial, plausible causes including adsorption of mercury ions, specifically adsorbed hydroxyl ions, potential varying distance between the OHP and the metal surface, among others^[15]. The latter phenomenon is ascribed to the orientational polarization of interfacial water molecules, which was initially suggested by Grahame and latter modelled by Watts-Tobin^[16]. Watts-Tobin assumed that interfacial water molecules may occupy two states (H-down or O-down). The basic idea is that interfacial water molecules are more polarized at more charged surface, resulting in a decreased permittivity and lower C_IHP. Consequently, the hump of C_IHP is located at the pzc where the permittivity of water molecules is maximum. The deviation of the hump from the pzc observed in experiments is caused by the ‘natural field’ on the metal surface, namely, metal electronic effects, which became a hot topic in 1980s^[17]. The Watts-Tobin model had been refined in several rounds by considering more states of water and the hydrogen-bond network, see a review by Guidelli and Schmickler^[18].

Together with Grahame’s data of C_IHP that gave the first hint of the importance of metal electronic effects, Trasatti’s correlation between C_IHP of simple sp metals taken at the pzc and the metal electron density drove theorists to explicitly consider the metal electronic effects^[19]. In 1980s, Schmickler, Badiali, Kornyshev and their associates introduced the jellium model that had been widely used in the theory of metal surfaces to the EDL theory^[20,21,22]. In the jellium model, the metal is treated as an inhomogeneous electron gas situated against a positive background charge corresponding to metal cationic cores. The electron gas was described using local density approximations, such as the Thomas-Fermi-von Weizsäcker theory^[23]. The positive background charge was later replaced with pseudopotentials to consider metal specific behaviors. The jellium model is able to rationalize the metal and surface-charge dependence of C_IHP^[17].

The next milestone in the EDL modelling is the work of Price and Halley in 1995^[24]. They adopted the Car-Pa-rrinello method of combining molecular dynamics and density functional theory (DFT) to simulate the EDL formed at a copper slab and water molecules. This work opened up a new direction in atomistic modelling of the EDL. Since then, the field has been shifted to extensive computer simulations with the Kohn-Sham DFT at the core. Most DFT-based first-principles simulations have been conducted for electroneutral interfaces. In 2006, Otani and Sugino developed a method to simulate charged interfaces^[25]. They employed the Poisson-Boltzmann equation for the electrolyte solution to screen excess charge on the metal slabs. This so-called implicit solvation method was subsequently advanced by several groups of authors, including Arias et al.^[26,27,28], Jinnouchi and Anderson^[29], Hennig et al.^{[30, 31]}, Nishihara and Otani^[32], Marzari et al.^{[33, 34]}, and Melander et al.^[35] Recent ab initio molecular dynamics (AIMD) simulations were able to handle the charged solid electrodes by introducing ions into the solvent layers^{[36,37,38,39]}. This provides a means to properly treat the electrode potential which is calculated from the work function of the system referenced to a (computational) standard hydrogen electrode (SHE).

Recently, Huang and coworkers developed a hybrid density-potential functional theory (DPFT) for the EDL, combining quantum mechanical treatment of many electrons and classical statistical treatment of charged particles in the solution phase^{[1, 40, 41]}. The DPFT avoids the calculation of Kohn-Sham orbitals which is computationally expensive. Instead, the DPFT is based on so-called orbital-free DFT^[42,43,44]. This feature distinguishes the DPFT model from other DFT-based first principles models, such as the joint DFT model for the EDL developed by Arias et al^{[27, 28]}. In Refs.^{[1, 41]}, the orbital-free DFT for metal electrodes consists of the Thomas-Fermi-von Weizsäcker theory for the electronic kinetic energy, a rudimentary kinetic energy density functional (KEDF), and the Dirac-Wigner theory for the exchange-correlation functional, a rudimentary local density approximation. As for the electrolyte solution, Huang developed a statistical field theory considering asymmetric steric effects, solvent polarization, and ion-specific interactions with the metal^[41]. Combined, a hybrid density-potential functional for the grand potential functional of the EDL is obtained. Variational analysis of this functional yields a grand-canonical EDL model described by two Euler-Lagrange equations in terms of the electron density and the electric potential.

The Gouy-Chapman-Stern (GCS) model is a classical toy model of the EDL. The electric potential distributes linearly from the electrode to the HP, and is described by Poisson-Boltzmann (PB) equation in the diffuse layer and the diffusion layer as shown in the third subfigure in Figure 3. The distance from the electrode to the HP is usually taken as the radius of a hydrated ion.

The PB equation describes the distributions of the electric potential and the ion concentrations in the electrolyte solution. Poisson equation reads,

where $\epsilon_{s}$ is the dielectric permittivity of the bulk solution, ϕ the electric potential referenced to the electric potential in the bulk solution, ϕ_S, z_i the charge number of ion i, F the Faraday’s constant, c_i the concentration of ion i. Boltzmann equation further connects c_i and ϕ,

where c^b is the concentration of ion i in the bulk solution, R the gas constant, T the absolute temperature. For a monovalent electrolyte solution in a one-dimensional space, the PB equation is rewritten as,

where c^b the concentration of total anions (cations) in the bulk solution. The dimensionless form of the PB equation is written as,

with the dimensionless quantities, U = Fϕ/RT, X = x/λ_D, and the Debye length λ_D =$\sqrt{RT\epsilon_{s}/2F^{2}c^{b} }$. Note that the dielectric permittivity depends on the local density of solvent molecules and the local electric field. In an aqueous electrolyte solution, the dielectric permittivity is a multiple of $\epsilon_{0}$ inside the HP, and increases to 78.5 $\epsilon_{0}$ in the bulk solution. Therefore, different values of the dielectric permittivity are used on the two sides of the HP.

The boundary conditions to close Eq.(5), a second-order differential equation, are,

where X = 0 represents the left boundary at the HP, and X = ∞ is the right boundary in the bulk solution. ϕ_HP^[4] (dimensional quantity of U_HP) can be calculated from the electrode side by,

where ϕ is the electric potential at the HP, ϕ_M the electrode potential, ϕ_pzc the potential of zero charge (note that this is not the absolute potential drop at the electrode surface, but the one relative to that of a reference electrode which is a constant), $\epsilon_{HP}$ and δ_HP are the dielectric permittivity and the thickness of the space between the electrode and the HP, respectively. The coefficient _s/$\epsilon_{HP}$ is resultant from the following equality in terms of surface charge density on the electrode surface,

Solving Eq. (5) in the following steps,

we obtain the relationship between the surface charge density and the electric potential at the HP,

Bvp4c is a convenient built-in tool in Matlab for solving boundary value problems described as ordinary differential equations. In accord with the syntax of this tool, Eq. (5) is rewritten as,

where $Y=\frac{F \lambda_{\mathrm{D}}}{R T} \frac{\partial \phi}{\partial x}$ is the dimensionless electric field strength.

Figure 4 shows the typical results of the GCS model (the Matlab script is provided in the supporting information of this article), including the spatial distributions of the electric potential, ϕ, and the cation concentration, c₊, at a series of ϕ_M in (A) and (B), and the relationships between the surface charge density, σ_M, and the differential double-layer capacitance, C_dl, with ϕ_M, in (C) and (D). The spatial range for calculation is 150 nm from the HP. When ϕ_M - ϕ_pzc > 0, we find ϕ_HP > 0, σ_M > 0, and cations are repelled. When ϕ_M - ϕ_pzc = 0, we obtain ϕ_HP = 0, σ_M = 0. When ϕ_M - ϕ_pzc < 0, we get ϕ_HP < 0, σ_M < 0, and cations are attracted. As defined in Eq.(1), C_dl has the minimum at ϕ_pzc.

A simple calculation can illustrate the failure of the GCS model in extreme cases. According to Eq.(3), c_i = c_i^b exp(-z_iFϕ/RT), we obtain c_i = 8.18 × 10¹⁶ mol·m^-3, when ϕ = -1 V, z_i = 1, c_i^b = 1 mol·m^-3 and T = 298 K. Consequently, each cation occupies a volume of 2.03 × 10^-35 cm³. However, even for the smallest bare cation, H⁺, the volume is approximately, d³≈(0.56 Å)³ = 1.76 × 10^-25 cm^{3 [45]}. Thus, it is necessary to consider the finite size of ions in the diffuse layer.

In 1942, Bikerman realized the limitations of neglecting ion size in the GCS model and developed a new model, called Bikerman-Poisson-Boltzmann (BPB) model, as shown in Figure 5^{[10, 46]}. In contrast with the GCS model, the BPB model presents a consistent treatment of the finite size of ions both at the HP and in the diffuse layer.

The BPB model treats the electrolyte solution using the lattice-gas approach. Each ion occupies a volume of d_t³, where d_t is the lattice size. The maximum particle number density is n_t = d_t^-3. The electrochemical potential for ion i reads,

where μ_i⁰ is the chemical potential under standard conditions, e₀ the elementary charge, ϕ the electric potential referenced to that in the bulk solution, ϕ_S, k_B the Boltzmann constant, n_i the number density of ion i, (1-d_t³∑n_i)/d_t³ the number density of solvent molecules. For a monovalent electrolyte solution, we have n₊⁰ = n_-⁰ = n^b, with n^b the number density of total anions (cations) in the bulk solution. Under equilibrium conditions, the electrochemical potential for ion i is uniform in the whole EDL,

The number density of ion i is obtained as,

where the bulk volume fraction of solvated ions is v = 2d_t³n^b. The GCS model assumes v = 0.

Combining Eq. (2) and Eq. (18), the BPB model is described as,

In a one dimensional case, the dimensionless form is,

The surface charge density can be calculated as,

The ‘bvp4c’ function in Matlab is employed to solve Eq. (20) closed with the boundary conditions expressed in Eqs. (6) and (7). Figure 6 shows the typical results of the BPB model, including the spatial distributions of ϕ and the anion concentration, c_-, at a series of ϕ_M in (A) and (B), as well as the relationships between σ_M and C_dl with ϕ_M in (C) and (D). For the purpose of comparison, the results of the GCS model at ϕ_M - ϕ_pzc = 0.7 V are shown in the black solid lines. The distributions of ϕ and c_- calculated using the GCS model are steeper than those calculated using the BPB model. In Figure 6(B), a plateau forms when ϕ_M - ϕ_pzc≥ 0.3 V, signifying the natural formation of the Stern layer due to the overcrowding of counterions. Figures 6(C) and 6(D) display how σ_M and C_dl change with ϕ_M at three values of v. A larger v means either larger ions or higher concentrations or both. At larger v, the relationship between σ_M and ϕ_M is less steep, indicating smaller values of C_dl. Interestingly, the shape of C_dl changes from a camel shape with the minimum at ϕ_pzc to a bell shape with the maximum at ϕ_pzc. Kornyshev gives a critical value of v = 1/3 for the camel-to-bell transition^[47].

With the asymmetric size effect, the electrochemical potential, Eq. (16), is rewritten as,

where γ_i is the size coefficient,

with d_i being the length of the cubic cell occupied by particle i, and d_t the reference size, usually taken as that of solvent molecules.

The number density for a monovalent electrolyte solution reads^[48],

To obtain Eq.(24), γ₊ and γ_- in the exponent are approximated as 1. Combining Eq.(2), the dimensionless form of the BPB equation in a one-dimensional case is,

Typical results of the asymmetric BPB model are presented in Figure 7, showing how σ_M and C_dl change with ϕ_M for the cases of different size coefficients. We set v = 0.05, and the size coefficient of cations, γ₊, or that of anions, γ_- , is equal to 1 or 3. At larger γ_i, the relationship between σ_M and ϕ_M is less steep, indicating smaller values of C_dl. γ₊ has bigger impact than γ_- as ϕ_M - ϕ_pzc < 0 because the concentration of cations dominates in this region due to the electrostatic interaction. γ_- plays an important role as ϕ_M - ϕ_pzc > 0.

In this part we consider dynamics of the EDL brought out of equilibrium. We build nonequilibrium models by using a grand potential approach, considering the size asymmetry effects. The solvent polarization which leads to a field-dependent dielectric permittivity is considered in ref.^[49], but is neglected in the following. Being grand-can-onical, the EDL exchanges electrons freely with the electrode and exchanges ions and solvent molecules freely with the bulk solution. Note that the EDL described here is not limited to a multiple of the Debye length, but could be extended to the bulk solution, because the diffuse layer and the diffusion layer are described by the same set of equations.

Under the conditions of constant electrochemical potential, constant T and a fixed volume V, the grand potential Ω of the EDL is written as,

where U is the internal energy, S the entropy, the electrochemical potential of particle i, n_i the number density of particle i, dV the volume unit.

There are multiple ions and solvent molecules in the electrolyte solution in general. Ions, denoted with a subscript α, have a charge number z_α and a number density n_α. There is a population of ions at (near) the transition state of the ion hopping process^[50], denoted with a superscript ≠. Solvent molecules are denoted with a subscript s. These charged particles, namely, ions and solvent molecules at both ground and excited states, interact via coulombic forces among others. According to field theoretic studies of the coulombic fluid^{[51, 52]}, the internal energy U is expressed as,

The first two terms represent the electrostatic interactions, including the self-energy correction of the electric field, $-\frac{1}{2} \epsilon_{s}$(∇ϕ)², and the electrostatic free energies of the ions, e₀ ϕ∑z_αn_α, and that of the transition-state ions, e₀ϕ∑z_α≠n_α≠. The last term accounts for many-body interactions other than the electrostatic interactions, with H_α being the internal energy except the electrostatic contribution, E_a,α≠ the activation energy of ion hopping.

The total entropy S is calculated from the lattice-gas model^[48],

where P is the number of ways arranging all the particles in the volume unit dV,

where N_α = n_α dV, N_s = n_s dV, and N_α≠ = n_α≠ dV are the particle numbers, and N_t = n_t dV is the total number of lattice cells in the volume unit, with n_t the number density.

The lattice cells are fully occupied without any vacancy, thus N_s is given by,

Note that there are several methods to treat size asymmetry in the lattice-gas model, as recently compared by Zhang and Huang^[48]. What we have used in Eqs. (29) and (30) is Huang’s treatment^[53]. The basic idea is to effectively expand the number of total sites. However, the size asymmetry is not considered in the calculation of P expressed in Eq.(29). As shown by Zhang and Huang, this approach captures major phenomena of the asymmetric steric effects and avoids the artificial sequence effects^[48].

Using the Stirling formula and taking the continuous limit (transforming the summation to a volume integration), we reformulate Eq. (28) as,

Combining Eqs. (26), (27), and (31), we rewrite the grand potential as a volume integration of a volumetric grand potential,

with the volumetric grand potential f given by,

Using the Euler-Lagrange equation,

in terms of X = ϕ, we obtain the Poisson equation,

When applying the Euler-Lagrange equation in terms of X = n_α, we must take notice of the relation that n_s = n_t/γ_s -∑n_αγ_α/γ_s -∑n_α≠γ_α≠/γ_s is also a function of n_α. Therefore, adding an ion α will simultaneously reduce γ_α/γ_s solvent molecules. The consideration leads to,

from which we define the electrochemical potential of the ion-solvent pair,

with the electrochemical potential of solvent molecules expressed as,

and the electrochemical potential of ions α as,

When the fictitious lattice cells are occupied exclusively by ions α, namely, n_t = n_α, the chemical potential turns to the standard chemical potential of ions α, denoted as μ_α⁰,

Applying the Euler-Lagrange equation in terms of X = n_α≠, we obtain the standard chemical potential of transition-state ions α^≠,

Although the transition-state ions are explicitly included in the grand-canonical potential, we will make the approximate that n_α≠≪n_α in the following.

According to the Fick’s second law, the continuity equation for particle α in the i^th cubic cell is written as,

where J_αⁱ is the flow flux of particle α cross the i^th cubic cell, the interface between the i^th cubic cell and the (i+1)^th cubic cell,

Eq.(43) means that the ion transport process is pictured as an ion-solvent exchange reaction, which has been proposed earlier to describe ion transport in solid and concentrated electrolytes^{[50, 54]}, where k_α^i→i+1 and k_α^i+1→i represent the forward and backward rates of ion hopping from the i^th to (i+1)^th cubic cell, respectively. According to transition-state theory and using the Brønsted-Evans-Polanyi (BEP) relationship to associate the activation barrier and the Gibbs free energy change, we write k_α^i→i+1 as,

and $k_{\alpha}^{i+1 \rightarrow i}$ as,

where the standard rate constant k_α⁰ is expressed as,

with τ_α,0 being the time constant of the hopping process. The 1/2 here indicates that ions at the transition state are equally likely to go forward to the new state and backward to the original state.

Substituting Eqs. (44), (45) and (46) into (43), we rewrite the flux as,

where D_α is the diffusion coefficient,

In the near-equilibrium regime, sinh(βd_t∇$\bar{\mu}$_α-s/2)≈βd_t∇$\bar{\mu}$_α-s/2, Eq. (47) can be approximated to,

where ∇$\bar{\mu}_{\alpha-s}$ is expanded as,

Combining Eqs. (42), (49) and (50), we obtain a modified Nernst-Planck equation,

The modified Poisson-Nernst-Planck (PNP) equations in Eq. (35) and Eq. (51) constitute the continuum model for multicomponent mass transport in electrolyte solution, which is derived from the grand potential being a functional of the electric potential and the particle number densities.

The boundary conditions and the initial conditions are necessary for solving the PNP equation, a set of partial differential equations. The boundary conditions are commonly divided into three types, Dirichlet, Neumann and Robin. Dirichlet boundary conditions specify the variable value on the boundary, for example, y = y₀ at x = x₀. Neumann boundary conditions specify the derivative of the variable, for example, $\frac{\alpha y}{\alpha x}$= α at x = x₀. Robin boundary conditions are a combination of Dirichlet and Neumann boundary conditions, for example, y- $\frac{\alpha y}{\alpha x}$= 0 at x = x₀.

The reaction plane, designated as the coordinate origin, x = 0, is the left boundary of which conditions are,

where j represents the current density of the overall reaction, m is the number of transferred electrons in the overall reaction, m_α is the stoichiometric number of the particle α in the overall reaction. If particle α does not participate in the reaction, we use m_α = 0. Eq. (53) shows the electric potential at the HP calculated from the electrode side, which is the same as Eq. (8). Eq. (52) is a Neumann boundary condition, while Eq. (53) is a Robin boundary condition.

The bulk solution is the right boundary, x = x_r , with the following natural boundary conditions,

which are the Dirichlet boundary conditions, where n_α^b is the number density of particle α in the bulk solution.

At t = 0, the initial conditions are shown as,

We consider a proton-coupled electron transfer reaction, A + H⁺ + e^- ↔ B, occurring at the HP, with A and B being neutral species. The current density of the reaction, j, is described by the Frumkin-Butler-Volmer (FBV) equation^[55],

where n_M is the areal number density of the electrode. For example, n_M is calculated by ($\sqrt{3}$a_M²)^-1 for M(111) with the lattice constant a_M. The pre-exponential factor k₀₀ is equal to $\frac{k_{B}T}{h}$ exp (-$\frac{\triangle G_{a}^{00}}{k_{B}T}$), with h being the Planck constant, ΔG_a⁰⁰ the activation energy of the reaction at standard equilibrium state. α is the charge transfer coefficient, taken as 0.5. c_B,HP, c_A,HP and c_H+_,HP are the concentrations of B, A and H⁺ at the HP, respectively. c_B⁰, c_A⁰ and c_H+⁰ are the concentrations of B, A and H⁺ under standard conditions, respectively. η is the overpotential, defined as,

where E⁰⁰ is the equilibrium potential of the reaction under standard conditions, calculated by E⁰⁰ = -ΔG⁰/e₀, with ΔG⁰ being the Gibbs free energy under standard conditions.

We numerically solve the model using the built-in partial differential equations solver, pdepe function, in Matlab, with the parameters listed in Table 1. Matlab script of this model is provided in the supporting information.

The typical results of the PNP equation are shown in Figure 8, including the steady current density, j, as a function with ϕ_M, and the distributions of the concentration of A, c_A, at 0.1 s, 1 s and 5 s at 0.4 V_SHE. The spatial range of the EDL is 100 μm, and the time duration is 5 s. The steady current density is taken at 5 s. As ϕ_M > E⁰⁰, the oxidation reaction occurs and B is consumed. The current density increases near exponentially in the low overpotential region and transitions to the diffusion limiting region when ϕ_M > 0.7 V_SHE, caused by the low concentration of B at the HP. When ϕ_M < E⁰⁰, the reduction reaction occurs. As the electric potential decreases, the current density increases and reaches the diffusion limiting current, which is limited by the low concentration of A at the HP. From Figure 8(B), we see as the reduction reaction occurs, the concentration of A at the HP is lower than that in the bulk solution, and decreases as the reaction continues. At 5 s, c_A becomes almost linear and reaches almost zero at the HP, signifying diffusion limiting effects.

Then we apply some approximations to reduce the modified PNP equation back to the classical PNP equation. Firstly, under the assumption n_α≠≪n_α , Eq. (35) is simplified as,

As for the modified Nernst-Planck equation, neglect of the asymmetric steric effects, that is, γ_i = 1, simplifies Eq. (51) to,

Furthermore, if the electrolyte is sufficiently dilute, that is, n_α≪n_t, and n_t≈n_s, the expression is returned back to the classical Nernst-Planck equation,

At equilibrium state, that is J_α = 0, Eq. (62) turns into Eq. (3), the Boltzmann equation. Similarly, when J_α = 0, Eq. (61) turns into Eq. (18), the equation adopted in the symmetric BPB model. Eq. (51) turns into Eq. (24), the equation used in the asymmetric BPB model.

Electrochemical impedance spectroscopy (EIS) is an in-situ, non-invasive characterization tool that can separate multiple physicochemical processes spanning a wide frequency range. In most cases, the EIS of an EDL is analyzed using the electrical circuit model (ECM) as shown in Figure 9. In fact, as to be discussed in the next paragraph, the ECM has a very clear physical meaning. However, it is not rigorous, theoretically. Instead, it is based on several assumptions which may become invalid in some cases. Therefore, it is of general importance to derive the impedance response of the EDL that is described using the PNP theory—the ‘first-principles’ of con-tinuum modelling of the EDL. We recommend the readers to follow the derivation with paper and pencil. This way, you will grasp the process of building a physical impedance model, acquire the basic mathematical tools, and appreciate the beauty of physicochemical modell-ing.

There are usually three physicochemical processes in the EDL: double-layer charging, charge transfer reactions, and diffusion. The double-layer charging involves redistribution of ions in the EDL, namely, change of the net charge stored in the EDL, under the control of the electric potential. As the EDL is usually only a few nanometers thick, ion transport in the EDL is often considered to be completed immediately. Therefore, charging the EDL is equivalent to charging an interfacial capacitance C_dl. As for the charge transfer reactions, it takes less than 1 ps for an electron to transfer between the electrode and the reactant in solution phase. Consequently, we can safely assume that the reaction current flows immediately when a potential difference is imposed. In other words, the current-electric potential relation of the charge transfer reaction is equivalent to that of a resistance R_ct. These two processes are in parallel because they are controlled by the same potential difference, and the total current is the sum of the double-layer charging part and the charge transfer reaction part. That is why the C_dl and R_ct are in parallel in Figure 9.

The W element in Figure 9 represents the diffusion of species involved in the charge transfer reactions in the electrolyte solution. The elements W and R_ct are in series because the transport process precedes/succeeds the charge transfer reactions. Conscious readers may have noticed a logic flaw: did not we consider ion transport in the EDL twice (one time in C_dl, and the other time in W)? There is another puzzle related to it. Given the fact that ion transport in the EDL and that in the diffusion layer are the same physical process, why do we need two elements? Why are C_dl and W located in different branches in the ECM? The way to resolve these puzzles has to be found via rigorous physicochemical modelling.

Before entering into physics-based impedance modelling, the definition of EIS and the fundamental mathematical tool—Fourier transform—will be introduced. Then, we will work on calculating the impedance response of a basic ECM to express the working mechanism of Fourier transform. Readers may find that EIS is more than a specific kind of classical electrochemical techniques. It provides a powerful mathematical physics approach to solve the electrochemical problems and represents a different look at electrochemical problems.

In this section, we introduce concisely the basics of EIS, including the Fourier transform and the example of a simple electrical circuit.

5.1.1 Fourier transform

The Fourier transform of a function f(t) is,

which transforms a time-domain signal f(t) into a frequency-domain signal F(ω). The inverse Fourier transform is,

derived from Eq.(63) as follows,

where we have used the Dirac’s delta function $\delta\left(\omega-\omega^{\prime}\right)=\frac{1}{2 \pi} \int_{-\infty}^{\infty} \exp \left(-j\left(\omega-\omega^{\prime}\right) t\right) \mathrm{d} t.$

Electrochemical processes are usually described by ordinary or partial differential equations. Therefore, the Fourier transform of the n^th derivative of a function is useful,

given that f ^(n-1)(t) = 0 at the initial state.

We prove for the case n = 1,

where we use the natural boundary conditions, f(∞) = f(-∞) = 0. The case of other orders can be proved by re-peating the manipulation of Eq.(67).

Another often-used property of Fourier transform is the convolution theorem,

where the sign “*” denotes the convolution operator defined as $f(t)^{*} g(t)=\int_{-\infty}^{\infty} f(\tau) g(t-\tau) \mathrm{d} \tau$, F(ω) and G(ω) denote the Fourier transform of f(t) and g(t), respectively. Eq. (68) is proved as follows,

5.1.2 Definition of Impedance

For any electrochemical system at stationary states, we apply an arbitrary current or electric potential excitation of small magnitude to ensure the linearity requirement, and obtain corresponding electric potential or current response. The electrochemical impedance is defined as the ratio of the Fourier transform of the potential to that of the current, that is,

A simple ECM is shown in Figure 10(A), which consists of two resistors and one capacitor. The total electric potential of this circuit is,

where

Applying Fourier transform into Eqs. (71) and (72), we obtain,

Combining Eqs. (70) and (74), the electrochemical impedance reads,

The real part, Z′(ω), and the imaginary part, Z″(ω), of Eq. (75) are,

where τ = RC is the time constant of this circuit. The amplitude and phase angle of this impedance are,

5.1.3 Perturbation Analysis

Considering a dilute, symmetrical and covalent electrolyte solution, we apply a potential perturbation to the system,

where U_M is the dimensionless electrode potential, the superscript “0” denotes the stationary electrode potential, the sign “~” denotes the magnitude of the perturbation potential and ω^nd is the dimensionless angular frequency referenced to D/λ_D². When the perturbation is sufficiently weak ($\widetilde {\phi}$_M < 25 mV, or \widetilde{U}_M < 1), the linear response approximation is valid. Then, all system variables are decomposed into stationary parts and perturbation parts of the same frequency, namely,

Substituting Eqs. (81)-(83) into the PNP equation, and omitting the stationary and high-order parts, we obtain,

Usually, C₊⁰, C_-⁰ and U⁰ are X-varying, making it difficult to solve Eqs. (84)-(86) analytically. Nevertheless, at the pzc, we have, C₊⁰ = C_-⁰ = 1 and U⁰ = 0. Therefore, Eqs. (84) and (85) are reduced to,

Substituting Eq. (86) into Eqs (87) and (88) leads,

Two equations above can be rewritten into a matrix form,

The eigenvalues, λ₁ and λ₂, and eigenvectors, V₁ and V₂, of matrix A are,

Introducing two matrixes P and H, and a vector y,

then we substitute the vector x with y and using the equality, P^-1AP = H, transforming Eq. (91) into,

Notice that H is a diagonal matrix, implying that the two elements of vector y, y₁ and y₂, can be solved separately,

where coefficients α₁, α₂, β₁, β₂ are to be determined by the corresponding boundary conditions. Afterwards, the two elements of the vector x are obtained as,

The boundary conditions of Eqs. (106) and (107) are as follows. In the bulk solution, X = X_b, the electric potential is regarded as the reference, namely, it does not change with the excitation. In addition, all ions have their bulk concentrations which also do not change with the potential perturbation,

At the HP, X = 0, the boundary conditions are rephrased as,

where _de is specifically deduced as follows. Firstly, we define,

where Π⁰ = U_M⁰ - U_HP⁰ - U_eq is the dimensionless stationary overpotential and $\widetilde{\Pi}=\widetilde{U}{ }_{\mathrm{M}}-\widetilde{U} \mathrm{HP}_{\mathrm{HP}}$ the dimensionless perturbation of the potential difference between the electrode surface and the HP.

We consider a reaction of metal ions deposition, M⁺ + e^- ↔ M, occurring at the HP. The current density of the reaction, j_de, is described by the FBV equation,

where c_+,HP is the concentration of M⁺ at the HP, c₊⁰ is the concentrations of M⁺ in the bulk solution, and k₀ is the reaction rate constant.

Substituting Eqs. (81) and (112) into the FBV theory, Eq. (113), gives,

Expanding Eq. (114) into a first-order Taylor series leads to,

Removing the stationary parts and second-order parts, we obtain,

Defining two variables related to the reaction rate,

we can reformulate Eq. (116) as,

5.1.4 Impedance Expression

Substituting Eqs. (106)-(107) into the corresponding boundary conditions expressed in Eqs. (108)-(111), we can solve for the four coefficients introduced in Eqs. (104)-(105),

where $r_{\mathrm{c}}=\frac{\epsilon_{\mathrm{s}} \delta_{\mathrm{HP}}}{\epsilon_{\mathrm{HP}} \lambda_{\mathrm{D}}}$ is the ratio between the Gouy-Chapman capacitance $\left(C_{\mathrm{GC}}=\frac{\epsilon_{\mathrm{s}}}{\lambda_{\mathrm{D}}}\right)$ and the Helmholtz capacit-ance $\left(C_{\mathrm{H}}=\frac{\epsilon_{\mathrm{HP}}}{\delta_{\mathrm{HP}}}\right)$.

Then we obtain the explicit expression of _de expressed in Eq. (118) by substituting Eqs. (119)-(122) into Eqs. (104) and (107). Besides the reaction current density j_de, the EDL current density also needs to be calculated, expressed as,

where q_dl is the total charges stored in EDL. The dimensionless form of j_dl is,

Converting Eq. (124) into frequency domain gives,

The total current density is the sum of the double layer current density and the Faradic reaction current density, namely,

So far, we obtain the dimensionless impedance expression,

where $\sum_{\mathrm{dr}}^{1}$ and $\sum_{\mathrm{dr}}^{2}$ represent the terms related to the coupling of charge transfer reaction and ion transport,

Based on previously defined dimensionless variables, we obtain the impedance reference, $Z_{\mathrm{ref}}=\frac{2 \lambda_{\mathrm{D}}^{2}}{D C_{\mathrm{GC}}}$. Notably, if there is no reaction at the HP, ν₁ = ν₂ = 0, we obtain $\sum_{\mathrm{dr}}^{1}=\sum_{\mathrm{dr}}^{2}=0$, and Eq. (127) is reduced to the impedance of an ideal polarizable electrode that has been given in Ref.^[58].

5.1.5 Simplifications

Although we have obtained the analytical solution expressed in Eq. (127), it is too complex. Therefore, we need to simplify this expression under some reasonable assumptions. Firstly, in a real system, we have X_b ≫1≈r_c, then we obtain,

Substituting Eqs. (130) and (131) into Eq. (127) leads,

Usually, $\omega^{\mathrm{nd}}=\omega \frac{\lambda_{\mathrm{D}}{ }^{2}}{D} \approx \omega \times 10^{-9} \ll 1$, then we obtain,

We notice that the final expression of Eq. (133) has a coefficient of $\frac{1}{2}$, which may look a little weird to readers. Therefore, we redefine the impedance reference as $\frac{\lambda_{D}^{2}}{DC_{GC}}$, then Eq. (133) is reformulated as,

We define, R_s^nd = X_b the dimensionless solution resistance, $R_{ct}^{nd}= \frac{2}{v_{2}} \frac{Dc_{0}}{\lambda_{D}}(1+r_{c})-\frac{v_{1}r_{c}}{v_{2}}$ the dimensionless charge transfer resistance and $W^{\mathrm{nd}}=\frac{\tanh \left(\sqrt{j \omega^{\text {nd }}} X_{\mathrm{b}}\right)}{\sqrt{j \omega^{\text {nd }}}}$ the dimensionless Warburg impedance. Eq.(134) embodies the coupling relationships between charge transfer reaction and EDL charging. Specifically, both R_ct^nd and C_dl^nd have the capacitance ratio term, r_c.

Then Eq. (133) is reformulated as,

From this simplified condition, we define the characteristic frequency of charge transfer reaction as, $\Omega _{ct}^{nd}=\frac{1}{R_{ct}^{nd}C_{dl}^{nd}}$, and the characteristic frequency of diffusion as, $\Omega_{d}=\frac{D}{x_{b}^{2}}$, whose dimensionless form is $\Omega_{d}^{nd}=\frac{1}{X_{b}^{2}} $.

When ω^nd < ω_d^nd, Eq. (135) is simplified to,

whose Nyquist plot is a 45^o-line followed by a semi-circle, as shown in Figure 11(D), representing ion diffusion in bulk solution.

When ω_d^nd < ω^nd < ω_ct^nd, Eq. (136) is reduced to,

whose Nyquist plot is an ideal semi-circle, as shown in Figure 11(C), representing the charge transfer reaction.

Finally, when ω^nd > ω_ct^nd, Eq. (137) is reduced to,

whose Nyquist plot is a straight line representing the EDL charging process, as shown in Figure 11(B).

In this section, we introduce the methods of calculating the impedance from time-domain data, which can be obtained from models and experiments. Firstly, the method of an analytical Fourier transform (AFT) is introduced^[59]. Then it is used to calculate the impedance of the deposition reaction of metal ions. Lastly, the fast Fourier transform (FFT), another often used numerical method, is introduced briefly.

5.2.1 Analytical Fourier transform

Applying linear interpolation to time-domain signal, we obtain

where h(t) is the recorded time-domain signal, σ(t) the normalized step function. Then applying Fourier transform to Eq. (139) gives,

where the specific derivations are detailed below. Integrating the last term of (t) gives,

Similarly, we obtain the frequency-domain expression for each term of Eq. (139),

Adding up all terms from (ω) to _n-1(ω), we obtain Eq. (140).

5.2.2 Application of AFT

Figure 12 compares the AFT-calculated and the analytical impedance expressed in Eq. (127). Notice that the results of AFT have some deviation in the entire frequency range. Especially, the deviation in the limiting high-frequency region is more obvious. The main reason is that the AFT method accumulates error in the numerical calculations. Due to the fact that the AFT calculation is a quite time-consuming task in a low frequency range, current results only show a 45^ο-line without a semi-circle that should occur in the low-frequency region.

Except for the AFT method, another often-used Fourier transform method is the FFT, which is widely used in signal processing^[60]. However, FFT is a completely pure numerical method. Compared with AFT, it lacks stability and has higher requirements for the signal-noise ratio of the time-domain signal.

This paper is designed as a tutorial tool on EDL modelling, including equilibrium models (the GCS model and the BPB model), nonequilibrium models (PNP-like models), and models under AC conditions (the EIS models). Exposition of these models begins with physical insights, followed by detailed mathematical derivation, formal analysis, and then practical numerical implementation with the Matlab scripts provided in the supporting information. A viable attempt to craft a physical model for the specific system under one’s own investigation could start with following the model development procedure presented here with pencil and paper.

This work is financially supported by National Natural Science Foundation of China (21802170) and the Alexander von Humboldt Foundation.

The authors declare no competing financial interests.