Incorporation of charge discreteness and ion correlations into lattice models of ionic liquids

Lattice-based mean-field models of ionic liquids neglect charge discreteness and ion correlations. To address these limitations, we propose separating the short-range and long-range parts of the electrostatic interaction by truncating the Coulomb potential below a fixed distance that is equal to or slightly larger than that between neighboring ions. Interactions and correlations between adjacent ions can then be modeled explicitly, whereas longer-ranged electrostatic interactions are captured on the mean-field level. We implement this approximation into the framework of modeling a compact, solvent-free ionic liquid by, first, considering terms up to the fourth order of the operator that represents the truncated Coulomb potential and, second, by accounting for electrostatic correlations between pairs of neighboring ions on the level of the quasi-chemical approach. A set of boundary conditions for the resulting self-consistent fourth-order differential equation follows from functional minimization of the free energy. The differential capacitance of an ionic liquid in contact with a planar electrode is calculated analytically up to quadratic order in the electrode’s surface charge density by solving the linearized model and applying a perturbation approach valid beyond the linear regime. We demonstrate that charge discreteness enhances the differential capacitance, whereas electrostatic correlations between ion–ion pairs drive the transition from a bell-shaped to a camel-shaped profile of differential capacitance. Our approach offers a systematic way to further improve the treatment of charge discreteness, account for short-range electrostatic and non-electrostatic interactions, and include higher-order ion–ion correlations.

1 Introduction

Ionic liquids consist of mixtures of cations and anions, usually of asymmetric size, that can remain in a liquid phase at temperatures typically below 100°C (Hayes et al., 2015). After initially being described as “water-free salts” (Walden, 1914), ionic liquids continue to gain attention due to their unique properties, including low volatility, high thermal stability, and excellent solvating capabilities (Sowmiah et al., 2009; Eyckens and Henderson, 2019). Given the diversity of molecular structures that form ionic liquids (Hayes et al., 2015) and their potential applications in energy conversion and storage (Armand et al., 2009; Wishart, 2009; Torimoto et al., 2010; Watanabe et al., 2017), it is not surprising that optimizing function–structure relationships is both complex and highly rewarding (Silva et al., 2020; Wang et al., 2020; Khlyupin et al., 2023; Nesterova et al., 2025). The challenges extend even to the level of understanding underlying thermodynamic and dynamic properties of ionic liquids (Weingärtner, 2008; Silva et al., 2020; Nordness and Brennecke, 2020). Compared to the ions in ordinary electrolytes, those of ionic liquids are densely packed and thus subject to ion–ion correlations due to electrostatic and non-electrostatic interactions (Del Pópolo and Voth, 2004).

Classical mean-field modeling (Kornyshev, 2007; May, 2019; Goodwin et al., 2023), which is a powerful conceptual tool to characterize ions in a solvent at sufficient dilution, fails for ionic liquids. Consequently, significant efforts have been undertaken in the past to incorporate ion correlations due to electrostatic and non-electrostatic interactions into the modeling of ionic liquids (Kondrat et al., 2023). Examples include the use of field theory (Démery et al., 2012; Nakamura, 2015; Budkov, 2020; Chao and Wang, 2020), density functional theory (Wu et al., 2011; Henderson et al., 2011; Ma et al., 2014), a microscopic cluster expansion model (Avni et al., 2020), molecular dynamics (Wang et al., 2007; Fedorov and Kornyshev, 2008; Coles et al., 2020; Zeman et al., 2021; Ye and Wang, 2022) and Monte Carlo (Kondrat et al., 2011; Bhuiyan et al., 2012) simulations, as well as phenomenological extensions of mean-field approaches (Frydel, 2016; Goodwin et al., 2017; Downing et al., 2018a). To allow for oscillations in ionic concentration profiles, Bazant, Storey, and Kornyshev (BSK) proposed a phenomenological Landau–Ginzburg-like model that leads to a fourth-order modified Poisson–Boltzmann equation with a correlation length that accounts for short-range ion–ion interactions (Bazant et al., 2011). The BSK approach has been used widely (Yochelis, 2014; Shalabi et al., 2019; Varner and Wang, 2022; Nesterova et al., 2025) and linked to microscopic approaches (Blossey et al., 2017; Avni et al., 2020). Yet, understanding of the physical basis of higher-order Poisson–Boltzmann approaches and their relation to even the most basic molecular properties such as charge discreteness, ion polarizability, non-electrostatic interactions, and electrostatic ion–ion correlations remains incomplete.

The present work aims at incorporating both charge discreteness and electrostatic correlations between ion pairs into a lattice-based model of an ionic liquid that is in contact with a single planar electrode. To account for charge discreteness, we separate the short-range and long-range parts of the electrostatic interaction by truncating the Coulomb potential below a distance $r_0$ to the neighboring ions. The operator that results from this truncation involves an infinite number of higher-order derivatives, yet it becomes equivalent to the BSK model in the limit of small $r_0$. We focus on this limit, which results in an approximate incorporation of electrostatic interactions on the level of continuum electrostatics except for nearest neighbors. In contrast, electrostatic interactions among nearest neighbor ions are modeled explicitly, thereby accounting for their discrete nature. To account for electrostatic correlations among neighboring ions, we adopt the quasi-chemical approximation (QCA) approach (Sher et al., 1987; Davis, 1996) when expressing the configurational entropy of the lattice gas. In principle, QCA can account for a correlation among any number of ions, but in the present work, we only include two-body correlations among nearest neighbor ions. Hence, we propose a method that accounts for both charge discreteness and electrostatic correlations between ion pairs in the modeling of ionic liquids. We recognize that the present work assumes $r_0$ is sufficiently small and neglects ion correlations that involve more than two ions. Yet, we also point out that, in principle, this approach can be extended to account for larger $r_0$ and higher-order correlations. In addition to presenting the method and deriving the self-consistency relationship (a fourth-order differential equation) plus its boundary conditions, we calculate analytic results for the differential capacitance up to quadratic order in the surface charge density of the electrode and numerical results for larger surface charge densities. Our results are compared with limiting cases, such as the absence of electrostatic correlations between ion pairs and the negligence of ion discreteness. We find that charge discreteness enhances the magnitude of the electrostatic potential, and thus also the differential capacitance, whereas electrostatic correlations between ion pairs induce a transition from a bell-shaped to a camel-shaped profile of the differential capacitance.

2 Theory

Consider a compact, solvent-free ionic liquid containing monovalently charged cations and anions of the same size. We represent the ionic liquid by a lattice of coordination number $z$, volume $\nu$ per lattice site, and lattice spacing $b\sim {\nu }^{1/3}$, with each site being occupied by either a cation or an anion. ${\varphi }_1={\varphi }_1(\mathbf{r})$ and ${\varphi }_2={\varphi }_2(\mathbf{r})$ denote the local mole fraction of cations (index “1”) and anions (index “2”), respectively, at position $\mathbf{r}$. The absence of solvent is expressed by the condition ${\varphi }_1+{\varphi }_2=1$ at every position $\mathbf{r}$. In the present work, we aim at separating contributions to the electrostatic energy of the ionic liquid that originate from the short-range and long-range parts of the Coulomb potential. To introduce this method, we write for the electrostatic energy due to the long-range part.

$U=\frac{1}{2{\nu }^2}\int d^3\mathbf{r}\int d^3{\mathbf{r}}^{\prime }\eta \left(\mathbf{r}\right)u\left(|\mathbf{r}-{\mathbf{r}}^{\prime }|\right)\eta \left({\mathbf{r}}^{\prime }\right),(1)$

where we have defined the difference $\eta ={\varphi }_1-{\varphi }_2$ in mole fractions between the cations and anions. Note that $e\eta /\nu$ specifies the local volume charge density, where $e$ denotes the elementary charge, and that $\eta$ is positive in regions with an excess of cations over anions. Equation 1 accounts for all pairwise interactions among the ions within the ionic liquid: $u(r)$ for cation–cation and anion–anion pairs, and $-u(r)$ for cation–anion and anion–cation pairs, given that the two ions are separated by a distance $r=|\mathbf{r}-{\mathbf{r}}^{\prime }|$. We introduce a dimensionless potential $\mathrm{\Psi }=\mathrm{\Psi }(\mathbf{r})$ through:

$\mathrm{\Psi }\left(\mathbf{r}\right)=\frac{1}{\nu }\int d^3{\mathbf{r}}^{\prime }\eta \left({\mathbf{r}}^{\prime }\right)u\left(|\mathbf{r}-{\mathbf{r}}^{\prime }|\right),(2)$

where, here and in the following, we express energies in units of $k_{B}T$, where $k_B$ is the Boltzmann constant and $T$ is the absolute temperature. In order to recast Equation 2 into a local equation for the potential $\mathrm{\Psi }$, we introduce a yet unknown differential operator $\mathcal{A}$ such that $\mathcal{A}u(r)=-\delta (r)$, where $\delta (r)$ is the Dirac delta function. This indeed implies the local equation:

$\mathcal{A}\mathrm{\Psi }=-\frac{\eta }{\nu },(3)$

and corresponding electrostatic energy

$U=-\frac{1}{2}\int d^3\mathbf{r}\mathrm{\Psi }\mathcal{A}\mathrm{\Psi }.(4)$

Only after specifying a concrete interaction potential $u(r)$, we can determine the corresponding differential operator $\mathcal{A}$. This is addressed in the following section.

2.1 Truncated Coulomb potential

To capture the long-range part of the electrostatic interaction, we choose a truncated Coulomb potential

$u\left(r\right)=\frac{l_B}{2r}\left(\frac{r}{r_0}+1-\left|\frac{r}{r_0}-1\right|\right),(5)$

where $l_B$ is the Bjerrum length and $r_0$ is the truncation length. Truncated Coulomb potentials have been employed in the past to cut off the long-range part with the goal to facilitate efficient computer simulations (Linse and Andersen, 1986; Baker et al., 1999). In contrast, we propose to cut off the short-range part of the Coulomb interaction, as shown in Figure 1: the black solid line in both diagrams of Figure 1 shows the truncated potential $u(r)$ according to Equation 5. The operator $\mathcal{A}$ can be identified conveniently in Fourier space, where the Fourier transformation

$\tilde{u}\left(k\right)=\frac{4\pi }{k}\underset{0}{\overset{\infty }{\int }}drr\sin \left(kr\right)u\left(r\right)=4\pi l_B\frac{\sin \left(k{r}_0\right)}{k^3{r}_0},$

together with $\tilde{\mathcal{A}}\tilde{u}(k)=-1$, gives rise to

$\begin{array}{ccc}\hfill \mathcal{A}& =\hfill & \frac{{\nabla }^2}{4\pi l_B}\frac{r_0\nabla }{\sinh \left(r_0\nabla \right)}\hfill \\ \hfill & =\hfill & \frac{{\nabla }^2}{4\pi l_B}{\displaystyle \sum _{i=0}^{\infty }}\frac{\left(2-4^i\right)B_{2i}}{\left(2i\right)!}{\left(r_0\nabla \right)}^{2i}=\frac{{\nabla }^2}{4\pi l_B}\left[1-\frac{1}{6}{\left(r_0\nabla \right)}^2+\frac{7}{360}{\left(r_0\nabla \right)}^4\mp \cdots \right].\hfill \end{array}(6)$

In Equation 6, $\nabla$ and $B_i$ denote the nabla differential operator and the ith Bernoulli number, respectively. Recall that the set of Bernoulli numbers can be defined via a series expansion of the generating function $x/(e^x-1)={\sum }_{i=0}^{\infty }{B}_i{x}^i/i!$. If the truncated Coulomb potential is approximated by the $j$-th order operator,

${\mathcal{A}}_j=\frac{{\nabla }^2}{4\pi l_B}{\displaystyle \sum _{i=0}^j}\frac{\left(2-4^i\right)B_{2i}}{\left(2i\right)!}{\left(r_0\nabla \right)}^{2i},(7)$

then the corresponding potential is given by the following equation:

$u_j\left(r\right)=\frac{2}{\pi }\frac{l_B}{r}\underset{0}{\overset{\infty }{\int }}dk\frac{\sin \left(kr\right)}{k}\frac{1}{{\displaystyle \sum _{i=0}^j}\frac{\left(2-4^i\right)B_{2i}}{\left(2i\right)!}{\left(-k^2{r}_0^2\right)}^i}.(8)$

The function $u_j(r)$ in Equation 8 is the $j$-th order approximation of the truncated Coulomb potential. For $j=0$, Equation 8 gives rise to the Coulomb potential $u_0(r)=l_B/r$. In the limit $j\to \infty$, the truncated Coulomb potential $u_{j\to \infty }(r)=u(r)$, as specified in Equation 5, is recovered. As has been pointed out by Lee et al. (2015), the first-order approximation of the truncated Coulomb potential is as follows:

$u_1\left(r\right)=\frac{l_B}{r}\left(1-e^{-\sqrt{6}\frac{r}{r_0}}\right).(9)$

The right diagram in Figure 1 shows the truncated Coulomb potential $u(r)$ together with the increasingly more accurate approximations $u_1(r)$, $u_2(r)$, $u_3(r)$, and $u_4(r)$. In the present work, we will focus on the first-order approximation $u_1(r)$ only, leading to a fourth-order Poisson equation ${\mathcal{A}}_1\mathrm{\Psi }=-\eta /\nu$ that utilizes the differential operator ${\mathcal{A}}_1={\nabla }^2[1-{(r_0\nabla )}^2/6]/(4\pi l_B)$ similar to the BSK model (Bazant et al., 2011). However, Equations 7, 8 provide a method for the incorporation of higher-order approximations in a straightforward manner.

Figure 1

Figure 1. Left diagram: plot of the truncated Coulomb potential $u(r)$ as a function of the radial distance $r$ according to Equation 5. The truncation distance is $r_0$, and the Bjerrum length is $l_B$. Right diagram: truncated Coulomb potential $u(r)$ (thick solid line) together with the increasingly more accurate approximations $u_1(r)$ (first order, shown in red color), $u_2(r)$ (second order, orange), $u_3(r)$ (third order, green), and $u_4(r)$ (fourth order, blue), as defined in Equation 8. The dotted line marks the Coulomb potential $l_B/r$.

We point out that the BSK model (Bazant et al., 2011) is an effective mean-field approach based on the pair potential $u_1(r)$. The emergence of $u_1(r)$ from the bare Coulomb potential can be the result of charge discreteness, which we consider in the present work, or of short-range correlations. Our present work considers electrostatic correlations beyond the short-range correlations of the BSK model. Electrostatic correlations between ion pairs are not accounted for by the BSK model.

2.2 Free energy minimization

Before formulating and minimizing a free energy that accounts for charge discreteness and electrostatic correlations between ion pairs, we start with three comments: first, choosing the truncation radius $r_0$ to be equal or slightly larger than the lattice spacing $b$ enables us to explicitly include the electrostatic interactions between neighboring ions: ${\omega }_{11}={\omega }_{22}=\omega =l_B/b$ between cation–cation and anion–anion pairs and ${\omega }_{12}={\omega }_{21}=-\omega =-l_B/b$ between cation–anion and anion–cation pairs. This introduces the definition $\omega =l_B/b$ and the notation ${\omega }_{ij}$ that we use in Equation 10. Second, we account for electrostatic correlations between ions located at neighboring lattice sites through QCA (Davis, 1996). Correlations involving more than two ions can, in principle, be incorporated into QCA (Bossa et al., 2015), but they are ignored in the present work. To account for correlated ion pairs, we consider the local fraction ${\varphi }_{ij}={\varphi }_{ij}(\mathbf{r})$ of $i$-$j$ pairs, where the indices $i$ and $j$ adopt the value 1 for a cation and 2 for an anion. Third, we introduce the chemical potentials ${\mu }_1$ and ${\mu }_2$ of the cations and anions, respectively, and choose them to ensure coexistence with a bulk phase of the ionic liquid, where ${\varphi }_1={\varphi }_2=1/2$. We also recall the definition of the lattice coordination number $z$. With that, we express the free energy of a lattice-based ionic liquid as follows (Downing et al., 2018b):

$\begin{array}{ccc}\hfill F& =\hfill & \frac{1}{\nu }\int d^3\mathbf{r}\left[-\frac{\nu }{2}\mathrm{\Psi }\mathcal{A}\mathrm{\Psi }+\left(1-z\right)\left({\varphi }_1\ln {\varphi }_1+{\varphi }_2\ln {\varphi }_2\right)\right.\hfill \\ \hfill & \left.+\frac{z}{2}{\displaystyle \sum _{ij=1}^2}\left({\varphi }_{ij}\ln {\varphi }_{ij}+{\omega }_{ij}{\varphi }_{ij}\right)-{\mu }_1{\varphi }_1-{\mu }_2{\varphi }_2\right].\hfill \end{array}(10)$

The four contributions to this free energy account for the electrostatic energy due to the long-range part of the Coulomb potential, the mixing entropy of ion pairs according to QCA, the electrostatic energy associated with neighboring ion pairs, and a Legendre transformation that fixes the chemical potentials of the cations and anions in the ionic liquid. The ${\varphi }_{ij}$ in Equation 10 must fulfill the three conservation relations ${\varphi }_1={\varphi }_{11}+{\varphi }_{12}$, ${\varphi }_2={\varphi }_{21}+{\varphi }_{22}$, and ${\varphi }_{12}={\varphi }_{21}$. This leaves one degree of freedom for the local distribution of ion pairs, which we define conveniently as $\bar{\varphi }={\varphi }_{12}+{\varphi }_{21}$. The difference in mole fractions $\eta ={\varphi }_1-{\varphi }_2$ constitutes another degree of freedom. We thus have the following equation:

${\varphi }_1=\frac{1+\eta }{2},{\varphi }_2=\frac{1-\eta }{2},{\varphi }_{11}={\varphi }_1-\frac{\bar{\varphi }}{2},{\varphi }_{22}={\varphi }_2-\frac{\bar{\varphi }}{2},{\varphi }_{12}={\varphi }_{21}=\frac{\bar{\varphi }}{2}.(11)$

Note that both $\bar{\varphi }=\bar{\varphi }(\mathbf{r})$ and $\eta =\eta (\mathbf{r})$ depend on the position $\mathbf{r}$ inside the ionic liquid. Calculation of the first variation in the free energy $F=F(\eta ,\bar{\varphi })$ leads to

$\begin{array}{ccc}\hfill \delta F& =\hfill & \frac{1}{\nu }\int d^3\mathbf{r}\left\{\delta \eta \left[\mathrm{\Psi }+\frac{1-z}{2}\ln \frac{1+\eta }{1-\eta }+\frac{z}{4}\ln \frac{1+\eta -\bar{\varphi }}{1-\eta -\bar{\varphi }}-\frac{1}{2}\left({\mu }_1-{\mu }_2\right)\right]\right.\hfill \\ \hfill & +\hfill & \left.\delta \bar{\varphi }\left[\frac{z}{4}\ln \frac{{\bar{\varphi }}^2}{\left(1+\eta -\bar{\varphi }\right)\left(1-\eta -\bar{\varphi }\right)}-z\omega \right]\right\}.\hfill \end{array}(12)$

We have not included boundary terms in Equation 12 because they will be discussed separately below. By symmetry, the chemical potentials ${\mu }_1={\mu }_2$ are equal. Thermal equilibrium demands $\delta F=0$, thus leading to the following two equations:

$e^{4\omega }=\frac{{\bar{\varphi }}^2}{{\left(1-\bar{\varphi }\right)}^2-{\eta }^2},0=2\mathrm{\Psi }+\left(1-z\right)\ln \frac{1+\eta }{1-\eta }+\frac{z}{2}\ln \frac{1+\eta -\bar{\varphi }}{1-\eta -\bar{\varphi }},(13)$

which must be fulfilled at each position $\mathbf{r}$. Equation 13 defines the two relations $\eta =\eta (\mathrm{\Psi })$ and $\bar{\varphi }=\bar{\varphi }(\mathrm{\Psi })$, with $\omega$ and $z$ as additional parameters. Spatial compositional variations arise from the solutions of Equation 3. However, to simplify calculations, we replace the operator $\mathcal{A}$ by its first-order approximation ${\mathcal{A}}_1$, as defined in Equation 7. This renders the level of modeling the long-range part of the electrostatic interactions similar to the BSK model (Bazant et al., 2011) yet with different boundary conditions as discussed below. Equation 3 then reads as follows:

$l^2{\nabla }^2\left(1-\frac{r_0^2}{6}{\nabla }^2\right)\mathrm{\Psi }=-\eta \left(\mathrm{\Psi }\right),(14)$

where we have defined the length $l=\sqrt{\nu /(4\pi l_B)}$. For example, $l_B=1\text{nm}$ and $\nu =1{\text{nm}}^3$ imply $l=0.3\text{nm}$. We will use $l$ as unit length throughout this work. Equation 14 is a self-consistency relationship that takes the form of a fourth-order differential equation with the function $\eta (\mathrm{\Psi })$ defined through a set of two algebraic equations. As pointed out above, electrostatic interactions are represented by the potential $u_1(r)$ in Equation 9. Unlike interpretations in terms of short-range ion correlations (Storey and Bazant, 2012; Moon et al., 2015; Blossey et al., 2017; Shalabi et al., 2019; de Souza and Bazant, 2020), we associate the use of $u_1(r)$ with an approximation of the truncated Coulomb potential $u(r)$ in Equation 5, which originates in our goal to model charge discreteness rather than correlations. A complementary interpretation of $u_1(r)$ would be in terms of excluded volume interactions between adjacent ions, and such an interpretation has been proposed recently (Gupta et al., 2020a).

The function $\eta =\eta (\mathrm{\Psi })$ that is needed to solve Equation 14 satisfies Equation 13 and is displayed together with $\bar{\varphi }(\mathrm{\Psi })$ in Figure 2 for different parameters $z$ and $\omega$. The solution of Equation 13 for $\omega =0$ is $\eta =-\tanh (\mathrm{\Psi })$, irrespective of $z$, and is shown as a black line in the left diagram of Figure 2. Curves of the same color refer to $z=2$ (red), $z=4$ (blue), and $z=6$ (green) and are calculated for different values of $\omega$, as indicated in the legend. Note that for $z=2$ (the red lines in Figure 2), Equation 13 yields the following simple analytic result:

$\eta \left(\mathrm{\Psi }\right)=-\frac{\sinh \mathrm{\Psi }}{\sqrt{e^{4\omega }+{\sinh }^2\mathrm{\Psi }}},\bar{\varphi }\left(\mathrm{\Psi }\right)=\frac{e^{2\omega }-\sqrt{1+\left(e^{4\omega }-1\right)\eta {\left(\mathrm{\Psi }\right)}^2}}{2\sinh \left(2\omega \right)},(15)$

which corresponds to the solution of the one-dimensional Ising model (Davis, 1996) in an external field $\mathrm{\Psi }$.

Figure 2

Figure 2. Relationships $\eta =\eta (\mathrm{\Psi })$ (left diagram) and $\bar{\varphi }(\mathrm{\Psi })$ (right diagram) according to Equation 13. Curves of the same color refer to $z=2$ (red), $z=4$ (blue), and $z=6$ (green) and are calculated for $\omega =1/2$ (solid colored lines), $\omega =1$ (broken lines), and $\omega =1.5$ (dotted lines). For $\omega =0$, Equation 13 predicts $\eta =-\tanh (\mathrm{\Psi })$ irrespective of $z$ (the black curve on the left diagram). All curves colored red $(z=2)$ follow the explicit expressions in Equation 15. The thick transparent lines in the left diagram represent third-order approximations of $\eta (\mathrm{\Psi })$, as specified in Equations 29, 30. They will be utilized in Sections 3.1, 3.2 to calculate the differential capacitance for the linearized model and beyond that using a perturbation approach.

2.3 Boundary conditions

From this point forward, we consider a single charged planar electrode that carries a uniform surface charge density $\sigma$ and is in contact with the ionic liquid. We locate the origin of a Cartesian coordinate system at the electrode surface, with its $x$-axis pointing normal into the ionic liquid. The electrostatic potential $\mathrm{\Psi }=\mathrm{\Psi }(x)$ then depends only on the distance $x$ to the electrode, and Equation 14 reads as follows:

$l^2{\mathrm{\Psi }}^{″}\left(x\right)-l^2\frac{r_0^2}{6}\mathrm{\Psi }″″\left(x\right)=-\eta \left(\mathrm{\Psi }\left(x\right)\right),(16)$

where a prime denotes the derivative with respect to the position $x$. We derive solutions of this fourth-order non-linear differential equation for $x\ge 0$ subject to the boundary conditions $\mathrm{\Psi }(x\to \infty )={\mathrm{\Psi }}^{\prime }(x\to \infty )=0$ as well as

${\mathrm{\Psi }}^{\prime }\left(0\right)=0,\mathrm{\Psi }‴\left(0\right)=\frac{6s}{l{r}_0^2},(17)$

where $s=\nu \sigma /(le)=4\pi l_{B}l\sigma /e$ is a scaled (dimensionless) surface charge density. Different sets of boundary conditions at the electrode surface $x=0$ have been proposed in the past for fourth-order differential equations of the type in Equation 16: ${\mathrm{\Psi }}^{\prime }(0)=-s/l$ and $\mathrm{\Psi }‴(0)=0$ in the original BSK model (Bazant et al., 2011) and ${\mathrm{\Psi }}^{\prime }(0)=-s/l$ and ${\mathrm{\Psi }}^{″}(0)=0$ based on variational free energy minimization (Gupta et al., 2020b). In addition, ${\mathrm{\Psi }}^{\prime }(0)=-s/l$ and $(r_0/\sqrt{6})\times \mathrm{\Psi }‴(0)={\mathrm{\Psi }}^{″}(0)$ have been proposed based on the continuity of the Maxwell stress tensor at the charged electrode (de Souza and Bazant, 2020) and on interpreting $u_1(r)$ in Equation 9 as a specific combination of Coulomb and Yukawa potentials (Bossa and May, 2020). We take the point of view that the truncated Coulomb potential applies to all charges in the system, including those within the ionic liquid and those at the electrode surface. In this case, the energy (per unit area $A$) of the system associated with the long-range part of the electrostatic interactions is as follows:

$\begin{array}{cccc}\frac{U}{A}\hfill & =\hfill & \frac{l^2}{2\nu }\underset{0}{\overset{\infty }{\int }}dx\left({{\mathrm{\Psi }}^{\prime }}^2+\frac{r_0^2}{6}{{\mathrm{\Psi }}^{″}}^2\right)\hfill & =-\frac{l^2}{2\nu }\left.\{,\mathrm{\Psi }\left(0\right)\left({\mathrm{\Psi }}^{\prime }\left(0\right)-\frac{r_0^2}{6}\mathrm{\Psi }‴\left(0\right)\right)+\frac{r_0^2}{6}{\mathrm{\Psi }}^{\prime }\left(0\right){\mathrm{\Psi }}^{″}\left(0\right)+\underset{0}{\overset{\infty }{\int }}dx\mathrm{\Psi }\left({\mathrm{\Psi }}^{″}-\frac{r_0^2}{6}\mathrm{\Psi }″″\right)\right\}\}.\hfill \end{array}(18)$

Its first variation is as follows:

$\frac{\delta U}{A}=-\frac{l^2}{\nu }\left.\{,\mathrm{\Psi }\left(0\right)\left(\delta {\mathrm{\Psi }}^{\prime }\left(0\right)-\frac{r_0^2}{6}\delta \mathrm{\Psi }‴\left(0\right)\right)+\frac{r_0^2}{6}{\mathrm{\Psi }}^{\prime }\left(0\right)\delta {\mathrm{\Psi }}^{″}\left(0\right)+\underset{0}{\overset{\infty }{\int }}dx\mathrm{\Psi }\left(\delta {\mathrm{\Psi }}^{″}-\frac{r_0^2}{6}\delta \mathrm{\Psi }″″\right)\right\}\}.(19)$

The integral in Equation 19 is $(1/\nu ){\int }_0^{\infty }dx\mathrm{\Psi }\delta \eta$, which, when combined with the variation in the non-electrostatic contributions to the free energy, will vanish, as demonstrated in Equation 12. Poisson’s equation for the truncated Coulomb potential, $l^2/\nu [{\mathrm{\Psi }}^{″}-(r_0^2/6)\mathrm{\Psi }″″]=-\eta /\nu$, relates the potential to the local volume charge density $e\eta /\nu$. Integrating that equation across the surface of the electrode yields $l[{\mathrm{\Psi }}^{\prime }(0)-(r_0^2/6)\mathrm{\Psi }‴(0)]=-s$. Because the scaled surface charge density $s$ is fixed, the ensuing relation $\delta {\mathrm{\Psi }}^{\prime }(0)-(r_0^2/6)\times \delta \mathrm{\Psi }‴(0)=0$ causes the first surface term in Equation 19 to vanish. Vanishing of the second surface term demands ${\mathrm{\Psi }}^{\prime }(0)=0$. Because of ${\mathrm{\Psi }}^{\prime }(0)=0$, the term $(r_0^2/6)\times {\mathrm{\Psi }}^{\prime }(0){\mathrm{\Psi }}^{″}(0)$ in the third line of Equation 18 vanishes. The remaining two terms, the electrostatic energy due to the surface charges on the electrode and the electrostatic energy due the ions in the ionic liquid, then render Equation 18 identical to Equation 4. Hence, the boundary conditions at the electrode, ${\mathrm{\Psi }}^{\prime }(0)=0$ and $\mathrm{\Psi }‴(0)=6s/(l{r}_0^2)$, lead to the consistent vanishing of $\delta F$ with $F$ defined in Equation 10, including all surface terms.

The following three comments complete our discussion of the boundary conditions. First, for $x<0$, the potential is constant, $\mathrm{\Psi }(x)=\mathrm{\Psi }(0)$. Hence, because the region $x<0$ does not contribute to $U$, the integration in Equation 18 may start from $x=0$ instead of also including negative values of $x$. Second, short-range electrostatic interactions of the ionic liquid with the electrode surface contribute another surface term $-\eta (0)sl/b$ to $U/A$ in Equation 18. Because the variation of this term vanishes, it does not add to the boundary conditions. Third, our approach recovers the boundary condition ${\mathrm{\Psi }}^{\prime }(0)=-s/l$, used previously to solve the second-order self-consistency differential equation $l^2{\mathrm{\Psi }}^{″}(x)=-\eta (\mathrm{\Psi }(x))$ that results from our model for $r_0=0$ (Downing et al., 2018b).

The solution $\mathrm{\Psi }(x)$ of Equation 16 reveals how the surface potential ${\mathrm{\Psi }}_0=\mathrm{\Psi }(x=0)$ depends on the scaled surface charge density $s$, allowing us to compute the differential capacitance as follows:

$C_{\mathit{diff}}=\frac{ϵ{ϵ}_0}{l}\frac{ds}{d{\mathrm{\Psi }}_0}=\frac{ϵ{ϵ}_0}{l}{\bar{C}}_{\mathit{diff}}.$

Calculation of the scaled differential capacitance ${\bar{C}}_{\mathit{diff}}=ds/d{\mathrm{\Psi }}_0$ is a major focus of the present work.

3 Results and discussion

Before computing general results for ${\bar{C}}_{\mathit{diff}}(s)$ numerically, we analytically derive quadratic expressions of the form ${\bar{C}}_{\mathit{diff}}={\bar{C}}_{\mathit{diff}}^{\mathit{lin}}+s^2\times \mathrm{\Omega }/2$ through linearizing our model followed by a perturbation approach to access the nonlinear regime. The determination of the two constants ${\bar{C}}_{\mathit{diff}}^{\mathit{lin}}$ and $\mathrm{\Omega }$ as functions of $r_0$, $\omega$, and $z$ provides insights about the role of charge discreteness and electrostatic correlations between ion pairs.

3.1 Linearized theory

For small potentials, $\mathrm{\Psi }\ll 1$, only the linear part of the relationship $\eta (\mathrm{\Psi })$ is significant. We use Equation 13 to perform a series expansion of $\eta (\mathrm{\Psi })$ up to the linear order. The result, $\eta =-\mathrm{\Psi }/\left[1+(z/2)\times \left(e^{2\omega }-1\right)\right]$, can conveniently be expressed as $\eta =-(3/2)\times {(l/r_c)}^2\times \mathrm{\Psi }$, where we have defined the length:

$r_c=l\sqrt{\frac{3}{2}\left[1+\frac{z}{2}\left(e^{2\omega }-1\right)\right]}.(20)$

This length will be shown to play a role in separating different regimes of the solution $\mathrm{\Psi }(x)$ for the linear problem. We thus solve the equation

${\mathrm{\Psi }}^{″}\left(x\right)-\frac{r_0^2}{6}\mathrm{\Psi }″″\left(x\right)=\frac{3}{2r_c^2}\mathrm{\Psi }\left(x\right),$

subject to ${\mathrm{\Psi }}^{\prime }(0)=0$, $\mathrm{\Psi }‴(0)=(6/r_0^2)\times (s/l)$, $\mathrm{\Psi }(x\to \infty )={\mathrm{\Psi }}^{\prime }(x\to \infty )=0$. The solution of the linear problem, denoted in the following by ${\mathrm{\Psi }}_{\mathit{lin}}(x)$, can generally be expressed as the sum of two contributions:

${\mathrm{\Psi }}_{\mathit{lin}}\left(x\right)=s\left[a_1{e}^{-{\lambda }_1\frac{x}{r_0}}+a_2{e}^{-{\lambda }_2\frac{x}{r_0}}\right],(21)$

with

$\frac{{\lambda }_1}{\sqrt{3}}=\sqrt{1+\sqrt{1-{\left(\frac{r_0}{r_c}\right)}^2}},\frac{{\lambda }_2}{\sqrt{3}}=\sqrt{1-\sqrt{1-{\left(\frac{r_0}{r_c}\right)}^2}}.(22)$

The two constants $a_1$ and $a_2$ follow from the boundary conditions in Equation 17,

$a_1=-\frac{r_0}{\sqrt{3}l}\frac{1}{\sqrt{1-{\left(\frac{r_0}{r_c}\right)}^2}\sqrt{1+\sqrt{1-{\left(\frac{r_0}{r_c}\right)}^2}}},a_2=\frac{r_0}{\sqrt{3}l}\frac{1}{\sqrt{1-{\left(\frac{r_0}{r_c}\right)}^2}\sqrt{1-\sqrt{1-{\left(\frac{r_0}{r_c}\right)}^2}}}.(23)$

In the following, we discuss how the structure of the solution changes as $r_0$ is varied. For $r_0=0$, the corresponding potential,

${\mathrm{\Psi }}_{\mathit{lin}}\left(x\right)=\sqrt{\frac{2}{3}}\frac{r_c}{l}s{e}^{-\left(\sqrt{\frac{3}{2}}\frac{x}{r_c}\right)},(24)$

is characterized by a single decay length. For $0<r_0<r_c$, there is a double-exponential decay with the two decay lengths $r_0/{\lambda }_1$ and $r_0/{\lambda }_2$. At $r_0=r_c$, the decay length $r_0/{\lambda }_2$ diverges, and the potential becomes as follows:

${\mathrm{\Psi }}_{\mathit{lin}}\left(x\right)=s\frac{r_c}{l}\left(\frac{1}{\sqrt{3}}+\frac{x}{r_c}\right)e^{-\left(\sqrt{3}\frac{x}{r_c}\right)}.(25)$

Finally, for $r_0>r_c$, the quantities ${\lambda }_1={\lambda }_3+i{\lambda }_4$ and ${\lambda }_2={\lambda }_3-i{\lambda }_4$, as well as $a_1=a_3-i{a}_4$ and $a_2=a_3+i{a}_4$, adopt complex conjugate values with

${\lambda }_3=\sqrt{\frac{3}{2}\left(\frac{r_0}{r_c}+1\right)}, {\lambda }_4=\sqrt{\frac{3}{2}\left(\frac{r_0}{r_c}-1\right)}$

and

$a_3=\frac{r_c}{\sqrt{6}l}\frac{1}{\sqrt{\frac{r_0}{r_c}+1}}, a_4=\frac{r_c}{\sqrt{6}l}\frac{1}{\sqrt{\frac{r_0}{r_c}-1}}.$

The potential

${\mathrm{\Psi }}_{\mathit{lin}}\left(x\right)=2s{e}^{-\left(\frac{{\lambda }_{3}x}{r_0}\right)}\left[a_3\cos \left(\frac{{\lambda }_{4}x}{r_0}\right)+a_4\sin \left(\frac{{\lambda }_{4}x}{r_0}\right)\right](26)$

then exhibits exponentially decaying oscillations. Figure 3 shows ${\mathrm{\Psi }}_{\mathit{lin}}(x)/s$ for $\omega =1$ (upper three curves) and $\omega =0$ (lower three curves) and three different ratios of $r_0/r_c$, namely $r_0=0$ (black), $r_0=r_c$ (red), and $r_0=2r_c$ (blue). Figure 3 suggests the tendency of $\omega$ to increase the magnitude and decay length of the potential. The truncation length $r_0$ only moderately modifies this tendency. Irrespective of the ratio $r_0/r_c$, the same linear relationship,

${\mathrm{\Psi }}_0^{\mathit{lin}}=s\sqrt{\frac{2}{3}}\frac{r_c}{l}\frac{1}{\sqrt{1+\frac{r_0}{r_c}}},(27)$

between surface potential ${\mathrm{\Psi }}_0^{\mathit{lin}}={\mathrm{\Psi }}_{\mathit{lin}}(x=0)$ and scaled surface charge density $s$ results. Hence, Equation 27 is valid within the linearized model for any choice of $z$, $\omega$, and $r_0$. This leads to a prediction for the differential capacitance ${\bar{C}}_{\mathit{diff}}(s=0)={\bar{C}}_{\mathit{diff}}^{\mathit{lin}}$ in the limit of an uncharged electrode, where both $s=0$ and ${\mathrm{\Psi }}_0=0$,

${\bar{C}}_{\mathit{diff}}^{\mathit{lin}}=\frac{ds}{d{\mathrm{\Psi }}_0^{\mathit{lin}}}=\frac{\sqrt{1+\frac{r_0}{r_c}}}{\sqrt{1+\frac{z}{2}\left(e^{2\omega }-1\right)}}.(28)$

Recall that $r_c$ depends on $z$ and $\omega$ through Equation 20. Clearly, increasing $\omega$ lowers ${\bar{C}}_{\mathit{diff}}^{\mathit{lin}}$, whereas larger $r_0$ increases ${\bar{C}}_{\mathit{diff}}^{\mathit{lin}}$.

Figure 3

Figure 3. ${\mathrm{\Psi }}_{\mathit{lin}}(x)/s$ for $\omega =1$ (upper three curves) and $\omega =0$ (lower three curves, magnified in the inset), with $r_0=0$ (black), $r_0=r_c$ (red), and $r_0=2r_c$ (blue). All curves are calculated for $z=6$. The black curves $(r_0=0)$ are single exponentials according to Equation 24 with a decay length $l$ for $\omega =0$ and $4.5\times l$ for $\omega =1$. The red curves $(r_0=r_c)$ are specified by Equation 25. The blue curves, which follow Equation 26, exhibit exponentially decaying oscillations.

3.2 Non-linear model: a perturbation approach

The linearized model is based on the relationship $\eta \sim \mathrm{\Psi }$. In this section, we employ a perturbation approach to analytically model the nonlinear region, where the potential deviates only slightly from ${\mathrm{\Psi }}_{\mathit{lin}}$. To this end, we expand the relationship $\eta (\mathrm{\Psi })$ up to the third order,

$-\eta =\frac{3}{2}{\left(\frac{l}{r_c}\right)}^2\mathrm{\Psi }+b_3{\mathrm{\Psi }}^3.(29)$

From Equation 13, we find the third-order coefficient as follows:

$b_3=\frac{\frac{z}{2}\left(e^{2\omega }+2\right){\left(e^{2\omega }-1\right)}^2-2}{6{\left[1+\frac{z}{2}\left(e^{2\omega }-1\right)\right]}^4}.(30)$

The left diagram in Figure 2 displays the third-order approximation, as specified in Equations 29, 30 (color-matching thick transparent lines). Hence, we aim to solve the nonlinear differential equation

${\mathrm{\Psi }}^{″}\left(x\right)-\frac{r_0^2}{6}\mathrm{\Psi }″″\left(x\right)=\frac{3}{2r_c^2}\mathrm{\Psi }\left(x\right)+\frac{b_3}{l^2}{\mathrm{\Psi }}^3\left(x\right),(31)$

subject to ${\mathrm{\Psi }}^{\prime }(0)=0$, $\mathrm{\Psi }‴(0)=(6/r_0^2)\times (s/l)$, and $\mathrm{\Psi }(x\to \infty )={\mathrm{\Psi }}^{\prime }(x\to \infty )=0$. We express the potential $\mathrm{\Psi }(x)={\mathrm{\Psi }}_{\mathit{lin}}(x)+b_3{\mathrm{\Psi }}_{\mathit{per}}(x)$ as the sum of the solution for the linear problem $({\mathrm{\Psi }}_{\mathit{lin}})$ and a small perturbation $(b_3{\mathrm{\Psi }}_{\mathit{per}})$. Substituting $\mathrm{\Psi }(x)$ into Equation 31 leads to the linear, inhomogeneous differential equation

${\mathrm{\Psi }}_{per}^{″}\left(x\right)-\frac{r_0^2}{6}\mathrm{\Psi }″{″}_{per}\left(x\right)=\frac{3}{2r_c^2}{\mathrm{\Psi }}_{\mathit{per}}\left(x\right)+\frac{1}{l^2}{\mathrm{\Psi }}_{\mathit{lin}}^3\left(x\right)(32)$

for the perturbation contribution of the potential ${\mathrm{\Psi }}_{\mathit{per}}(x)$, where ${\mathrm{\Psi }}_{\mathit{lin}}(x)$ is given by Equation 21. The solution ${\mathrm{\Psi }}_{\mathit{per}}(x)={\mathrm{\Psi }}_{\mathit{inh}}(x)+{\mathrm{\Psi }}_{\mathit{hom}}(x)$ of Equation 32 can be expressed as the sum of a specific solution for the inhomogeneous equation (index “inh”):

$\begin{array}{ccc}\hfill {\mathrm{\Psi }}_{\mathit{inh}}\left(x\right)& =\hfill & s^3{\left(\frac{r_0}{l}\right)}^2\left[\frac{a_1^3{e}^{-3{\lambda }_1\frac{x}{r_0}}}{{\left(3{\lambda }_1\right)}^2-\frac{1}{6}{\left(3{\lambda }_1\right)}^4-\frac{3}{2}{\left(\frac{r_0}{r_c}\right)}^2}+\frac{3a_1^2{a}_2{e}^{-\left(2{\lambda }_1+{\lambda }_2\right)\frac{x}{r_0}}}{{\left(2{\lambda }_1+{\lambda }_2\right)}^2-\frac{1}{6}{\left(2{\lambda }_1+{\lambda }_2\right)}^4-\frac{3}{2}{\left(\frac{r_0}{r_c}\right)}^2}\right.\hfill \\ \hfill & \left.+\frac{3a_1{a}_2^2{e}^{-\left({\lambda }_1+2{\lambda }_2\right)\frac{x}{r_0}}}{{\left({\lambda }_1+2{\lambda }_2\right)}^2-\frac{1}{6}{\left({\lambda }_1+2{\lambda }_2\right)}^4-\frac{3}{2}{\left(\frac{r_0}{r_c}\right)}^2}+\frac{a_2^3{e}^{-3{\lambda }_2\frac{x}{r_0}}}{{\left(3{\lambda }_2\right)}^2-\frac{1}{6}{\left(3{\lambda }_2\right)}^4-\frac{3}{2}{\left(\frac{r_0}{r_c}\right)}^2}\right],\hfill \end{array}$

where ${\lambda }_1$, ${\lambda }_2$, $a_1$, and $a_2$ are defined in Equations 22, 23, and a solution for the homogeneous equation (index “hom”):

${\mathrm{\Psi }}_{\mathit{hom}}\left(x\right)=s^3\left[c_1{e}^{-{\lambda }_1\frac{x}{r_0}}+c_2{e}^{-{\lambda }_2\frac{x}{r_0}}\right].(33)$

The two constants $c_1$ and $c_2$ in Equation 33 can be determined such that besides ${\mathrm{\Psi }}_{\mathit{per}}(x\to \infty )={\mathrm{\Psi }}_{\mathit{per}}^{\prime }(x\to \infty )=0$, Equation 32 also satisfies the boundary conditions ${\mathrm{\Psi }}_{\mathit{per}}^{\prime }(0)=0$ and ${\mathrm{\Psi }}_{\mathit{per}}^{‴}(0)=0$. This results in an expression for the surface contribution of the perturbation potential ${\mathrm{\Psi }}_{\mathit{per}}(0)=s^{3}B$ with

$B=-\frac{1}{4}{\left[1+\frac{z}{2}\left(e^{2\omega }-1\right)\right]}^{\frac{5}{2}}\times \frac{1+\frac{25}{6}\frac{r_0}{r_c}+\frac{41}{12}{\left(\frac{r_0}{r_c}\right)}^2}{{\left(1+\frac{r_0}{r_c}\right)}^{\frac{5}{2}}\left(1+\frac{5}{3}\frac{r_0}{r_c}\right)}.(34)$

The total surface potential ${\mathrm{\Psi }}_0={\mathrm{\Psi }}_{\mathit{lin}}(0)+b_3{\mathrm{\Psi }}_{\mathit{per}}(0)=s/{\bar{C}}_{\mathit{diff}}^{\mathit{lin}}+s^{3}B{b}_3$ can be used to calculate the scaled differential capacitance,

${\bar{C}}_{\mathit{diff}}=\frac{1}{\frac{d{\mathrm{\Psi }}_0}{ds}}={\bar{C}}_{\mathit{diff}}^{\mathit{lin}}-3B{b}_3{\left({\bar{C}}_{\mathit{diff}}^{\mathit{lin}}\right)}^2{s}^2={\bar{C}}_{\mathit{diff}}^{\mathit{lin}}+\frac{\mathrm{\Omega }}{2}{s}^2,(35)$

analytically up to quadratic order in $s$. Using the expressions for ${\bar{C}}_{\mathit{diff}}^{\mathit{lin}}$, $b_3$, and $B$ in Equations 28, 30, 34, respectively, to calculate the quadratic-order coefficient $\mathrm{\Omega }=-6B{b}_3{({\bar{C}}_{\mathit{diff}}^{\mathit{lin}})}^2$ yields the following:

$\mathrm{\Omega }=\frac{\frac{z}{2}\left(e^{2\omega }+2\right){\left(e^{2\omega }-1\right)}^2-2}{4{\left[1+\frac{z}{2}\left(e^{2\omega }-1\right)\right]}^{\frac{5}{2}}}\times \frac{1+\frac{25}{6}\frac{r_0}{r_c}+\frac{41}{12}{\left(\frac{r_0}{r_c}\right)}^2}{{\left(1+\frac{r_0}{r_c}\right)}^{\frac{3}{2}}\left(1+\frac{5}{3}\frac{r_0}{r_c}\right)}.(36)$

The quadratic expression for ${\bar{C}}_{\mathit{diff}}(s)$ in Equation 35 with the analytic results for ${\bar{C}}_{\mathit{diff}}^{\mathit{lin}}$ and $\mathrm{\Omega }$ given in Equations 28, 36 is the central outcome of this work, accounting (albeit on an approximate level) for both charge discreteness and electrostatic correlations.

3.3 Limiting cases

3.3.1 Continuum mean-field model, $r_0=0$ and $\omega =0$

For $r_0=0$ and $\omega =0$, Equation 16 reduces to the well-known classical mean-field equation $l^2{\mathrm{\Psi }}^{″}(x)=\tanh \mathrm{\Psi }$ that emerges from the lattice-based Poisson–Boltzmann framework of a solvent-free ionic liquid (Borukhov et al., 1997; Kornyshev, 2007). Neither charge discreteness nor ion correlations are accounted for. Subject to the boundary conditions ${\mathrm{\Psi }}^{\prime }(0)=-s/l$ and $\mathrm{\Psi }(x\to \infty )=0$, we obtain a surface potential ${\mathrm{\Psi }}_0=\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{h}(e^{s^2/2})$ for $x\ge 0$ and ${\mathrm{\Psi }}_0=-\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{h}(e^{s^2/2})$ for $x\le 0$, and thus, a scaled differential capacitance ${\bar{C}}_{\mathit{diff}}={(d{\mathrm{\Psi }}_0/ds)}^{-1}=\sqrt{1-e^{-s^2}}/|s|$. For small $|s|$, this can be represented by ${\bar{C}}_{\mathit{diff}}={\bar{C}}_{\mathit{diff}}^{\mathit{lin}}+s^2\mathrm{\Omega }/2$ with ${\bar{C}}_{\mathit{diff}}^{\mathit{lin}}=1$ and $\mathrm{\Omega }=-1/2$. All models discussed in the present work recover this limit if both $r_0=0$ and $\omega =0$.

3.3.2 Vanishing truncation length, $r_0=0$

For $r_0=0$, the expressions for ${\bar{C}}_{\mathit{diff}}^{\mathit{lin}}$ and $\mathrm{\Omega }$ in Equations 28, 36 simplify to the following:

${\bar{C}}_{\mathit{diff}}^{\mathit{lin}}=\frac{1}{\sqrt{1+\frac{z}{2}\left(e^{2\omega }-1\right)}}, \mathrm{\Omega }=\frac{1}{4}\times \frac{\frac{z}{2}\left(e^{2\omega }+2\right){\left(e^{2\omega }-1\right)}^2-2}{{\left[1+\frac{z}{2}\left(e^{2\omega }-1\right)\right]}^{5/2}}.$

This is identical to the results derived previously by Downing et al. (2018b), where an approach analogous to the present one yet without truncating the Coulomb potential was proposed. Not truncating the Coulomb potential while still adding nearest-neighbor electrostatic interactions separately into the model implies a double counting of nearest neighbor ion–ion interactions. The truncation of the Coulomb potential introduced in the present work eliminates this inconsistency.

As pointed out by Downing et al. (2018b), $\mathrm{\Omega }$ changes sign upon increasing $\omega$, indicating a transition from bell shape to camel shape of the differential capacitance ${\bar{C}}_{\mathit{diff}}$. This transition is caused by electrostatic correlations between ion pairs, as shown below by comparing the predictions for ${\bar{C}}_{\mathit{diff}}$ in the presence and absence of electrostatic correlations.

3.3.3 Vanishing central charge, $\omega =0$

The choice $\omega =0$ eliminates the explicit account of electrostatic interactions between nearest neighbors. Without these interactions, ion correlations between neighboring ions are no longer present. The solution of Equation 13 is $\bar{\varphi }=2{\varphi }_1{\varphi }_2$, and the corresponding relation $\eta =-\tanh \mathrm{\Psi }$ implies mean-field electrostatics subject to a truncated Coulomb potential. The resulting fourth-order differential equation $l^2{\mathrm{\Psi }}^{″}(x)-l^2{r}_0^2\mathrm{\Psi }″″(x)/6=\tanh \mathrm{\Psi }(x)$ is equivalent to the modified Poisson–Boltzmann equation of the solvent-free BSK model (Bazant et al., 2011) yet with different boundary conditions. For $\omega =0$, we obtain the following from Equations 28, 36:

${\bar{C}}_{\mathit{diff}}^{\mathit{lin}}=\sqrt{1+\sqrt{\frac{2}{3}}\frac{r_0}{l}}, \mathrm{\Omega }=-\frac{1}{2}\times \frac{1+\frac{25}{6}\sqrt{\frac{2}{3}}\frac{r_0}{l}+\frac{41}{12}{\left(\sqrt{\frac{2}{3}}\frac{r_0}{l}\right)}^2}{{\left(1+\sqrt{\frac{2}{3}}\frac{r_0}{l}\right)}^{\frac{3}{2}}\left(1+\frac{5}{3}\sqrt{\frac{2}{3}}\frac{r_0}{l}\right)}.(37)$

The negative sign of $\mathrm{\Omega }$ signifies a bell shape of the differential capacitance ${\bar{C}}_{\mathit{diff}}(s)$. Hence, for $\omega =0$, ${\bar{C}}_{\mathit{diff}}(s)$ cannot exhibit a camel shape.

3.3.4 The influence of electrostatic correlations between ion pairs

The usage of the QCA approximation accounts for electrostatic correlations between neighboring ion pairs, whereas three-body and higher-order correlations are neglected. To investigate the role of electrostatic ion pair correlations, we compare QCA to a mean-field model, which is based on a random mixing approximation (RMA) (Davis, 1996). RMA ignores ion pair correlations by modeling the entropy contribution of the free energy through an ideal lattice gas. This is accomplished by imposing $\bar{\varphi }=2{\varphi }_1{\varphi }_2$, which, when used in Equation 11, leads to the free energy

$F=\frac{1}{\nu }\int d^3\mathbf{r}\left[-\frac{\nu }{2}\mathrm{\Psi }\mathcal{A}\mathrm{\Psi }+{\varphi }_1\ln {\varphi }_1+{\varphi }_2\ln {\varphi }_2+\omega {\left({\varphi }_1-{\varphi }_2\right)}^2-{\mu }_1{\varphi }_1-{\mu }_2{\varphi }_2\right].$

Electrostatic interactions between neighboring ions are thus described on the level of the familiar Bragg–Williams model (Davis, 1996), whereas the long-range components are accounted for by the truncated Coulomb potential. As above, we use ${\varphi }_1=(1+\eta )/2$ and ${\varphi }_2=(1-\eta )/2$ and note that ${\mu }_1={\mu }_2$. Vanishing first variation of the free energy

$\delta F=\frac{1}{\nu }\int d^3\mathbf{r}\delta \eta \left[\mathrm{\Psi }-\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{t}\mathrm{a}\mathrm{n}\mathrm{h}\left(\eta \right)-2\omega \eta \right]$

then implies the relation

$\mathrm{\Psi }=-\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{t}\mathrm{a}\mathrm{n}\mathrm{h}\left(\eta \right)-2\omega \eta .(38)$

This implicitly defines the function $\eta (\mathrm{\Psi })$ to be used in the fourth-order self-consistency differential equation, Equation 14, namely, $l^2{\mathrm{\Psi }}^{″}(x)-l^2{r}_0^2\mathrm{\Psi }″″(x)/6=-\eta (\mathrm{\Psi }(x))$. Up to the third order in $\mathrm{\Psi }$, this function is given by $\eta =-\mathrm{\Psi }/(1+2\omega )+{\mathrm{\Psi }}^3/[3{(1+2\omega )}^4]$. As above, we solve the linearized model and perform a perturbation approach, yielding again an analytic expression for the scaled differential capacitance ${\bar{C}}_{\mathit{diff}}={\bar{C}}_{\mathit{diff}}^{\mathit{lin}}+s^2\mathrm{\Omega }/2$ up to quadratic order in the scaled surface charge density $s$. We find that

${\bar{C}}_{\mathit{diff}}^{\mathit{lin}}=\frac{\sqrt{1+\frac{r_0}{{\bar{r}}_c}}}{\sqrt{1+2\omega }}, \mathrm{\Omega }=-\frac{1}{2{\left(1+2\omega \right)}^{5/2}}\frac{1+\frac{25}{6}\frac{r_0}{{\bar{r}}_c}+\frac{41}{12}{\left(\frac{r_0}{{\bar{r}}_c}\right)}^2}{{\left(1+\frac{r_0}{{\bar{r}}_c}\right)}^{3/2}\left(1+\frac{5}{3}\frac{r_0}{{\bar{r}}_c}\right)},(39)$

where we have defined ${\bar{r}}_c=l\sqrt{3(1+2\omega )/2}$. Equation 39 accounts for the discrete nature of the ionic charges, but neglects electrostatic correlations between ion pairs. As for Equation 37, the negative sign of $\mathrm{\Omega }$ implies a bell shape of ${\bar{C}}_{\mathit{diff}}(s)$. Hence, in the absence of electrostatic correlations between ion pairs, our model does not predict a transition to a camel shape of ${\bar{C}}_{\mathit{diff}}(s)$.

3.4 Numerical results

In this section, we present numerical results for ${\bar{C}}_{\mathit{diff}}(s)$, with various parameter choices for $\omega$, $r_0$, and $z$, and then contrast the predictions of the model in the presence and absence of electrostatic ion–ion correlations. Solutions of Equation 14 employ the function $\eta (\mathrm{\Psi })$ according to Equation 13 in the presence of electrostatic ion correlations (modeled according to QCA) and according to Equation 38 in the absence of electrostatic ion correlations (modeled according to RMA).

Figure 4 summarizes the predictions of our model according to Equation 16 with $\eta (\mathrm{\Psi })$ specified in Equation 13 and the boundary conditions in Equation 17. Recall that this model incorporates both the truncated Coulomb potential and QCA, thus accounting, approximatively, for both charge discreteness and electrostatic correlations. The scaled differential capacitance ${\bar{C}}_{\mathit{diff}}$ is displayed as a function of the scaled surface charge density $s$ with $z=2$ (upper three diagrams) and $z=6$ (lower three diagrams), as well as $r_0=0$ (left, colored purple), $r_0=l$ (middle, colored red), and $r_0=3l$ (right, colored green). Each graph contains 11 curves, varying from $\omega =0$ (the light colored curve on the top of each diagram) to $\omega =1$ (the black curve on the bottom of each diagram) in steps of 0.1.

Figure 4

Figure 4. Scaled differential capacitance ${\bar{C}}_{\mathit{diff}}$ as a function of the scaled surface charge density $s$ based on solving Equations 13, 16, 17. Each diagram shows 11 curves for $\omega =0$ (uppermost curve) to $\omega =1$ (black curve on the bottom), with $\omega$ changing in steps of 0.1. Diagrams refer to $z=2$ (upper three diagrams) and $z=6$ (lower three diagrams), as well as $r_0=0$ (left column, colored purple), $r_0=l$ (middle column, colored red), and $r_0=3l$ (right column, colored green).

We point out that in the vicinity of $s=0$, where ${\bar{C}}_{\mathit{diff}}(s)={\bar{C}}_{\mathit{diff}}^{\mathit{lin}}+s^2\mathrm{\Omega }/2$ can be represented by a quadratic function, each curve reproduces our analytic expression, as specified in Equations 28, 36. Hence, Figure 4 reinforces the validity of our analytic results and provides an extension of ${\bar{C}}_{\mathit{diff}}(s)$ beyond the quadratic regime. As already predicted analytically, we observe the transition from a bell-shaped to a camel-shaped profile of the differential capacitance as the strength of the electrostatic interaction $\omega =l_B/b$ is increased. The transition occurs at a value of $\omega$ that satisfies the equation $2-3e^{2\omega }+e^{6\omega }=4/z$. For $z=2$ and $z=6$, this amounts to $\omega =1/4\times \ln 3\approx 0.275$ and $\omega =0.182$, respectively, independently of $r_0$. The truncation length $r_0$ does affect the value of ${\bar{C}}_{\mathit{diff}}(s=0)$ where the bell-to-camel shape transition occurs, but the value of $\omega$ where the transition happens is independent of $r_0$. The role of $r_0$ is to amplify the magnitudes of ${\bar{C}}_{\mathit{diff}}$. This is most clearly evidenced by the explicit relationship for the dependence of ${\bar{C}}_{\mathit{diff}}^{\mathit{lin}}$ on $r_0$ in Equation 28.

Figure 4 displays the dependence of ${\bar{C}}_{\mathit{diff}}(s)$ for the two parameters $r_0/l$ and $\omega$ (given $z$ is fixed). Yet, if we identify the truncation radius $r_0\sim b$ with the lattice spacing, there should be a relationship between $r_0/l$ and $\omega$. Recall the definitions $\omega =l_B/b$ and $l^2=\nu /(4\pi l_B)$. The lattice spacing $b\approx {\nu }^{1/3}$ represents the size of the ions in the ionic liquid, typically on the order of a nanometer. If we identify the truncation radius $r_0=b$ with the lattice spacing, we obtain $r_0/l\approx \sqrt{4\pi }\times \sqrt{\omega }$. Hence, the choice $\omega =1$ (implying $l_B=b$) corresponds to $r_0/l\approx 3$, which is displayed by the black curve in the right column of Figure 4. Further increasing the Bjerrum length (implying $\omega >1$) yields larger values for $r_0/l$. These cases are not shown in Figure 4 but are easily accessible through our analytic expression for ${\bar{C}}_{\mathit{diff}}(s)={\bar{C}}_{\mathit{diff}}^{\mathit{lin}}+s^2\mathrm{\Omega }/2$.

Figure 5 shows a comparison of our model with and without accounting for electrostatic correlations between ion pairs. The lower three diagrams in Figure 5 reproduce the lower three diagrams of Figure 4, calculated for $z=6$ and using QCA where ion pair correlations are accounted for. The upper three diagrams show analogous results, yet in the absence of ion pair correlations, employing RMA, as introduced in Section 3.3.4. We have verified the agreement of the numerical results for RMA (upper three diagrams) in the region of small $|s|$ with the analytic expression ${\bar{C}}_{\mathit{diff}}(s)={\bar{C}}_{\mathit{diff}}^{\mathit{lin}}+s^2\mathrm{\Omega }/2$ specified in Equation 39. As already pointed out, RMA does not predict a transition from a bell-shaped to a camel-shaped profile of the differential capacitance as the strength of the electrostatic interaction $\omega =l_B/b$ is increased. Hence, we conclude that electrostatic correlations rather than charge discreteness or excluded volume interactions cause camel-shape profiles of ${\bar{C}}_{\mathit{diff}}(s)$. Starting from the classical mean-field approach, ${\bar{C}}_{\mathit{diff}}=\sqrt{1-e^{-s^2}}/|s|$, displayed by the uppermost purple line on the left two diagrams in Figure 5, and introducing charge discreteness without accounting for electrostatic correlations corresponds to transitioning to the black curve on the right top diagram of Figure 5: the curve widens but without qualitatively changing its shape.

Figure 5

Figure 5. Scaled differential capacitance ${\bar{C}}_{\mathit{diff}}$ as a function of the scaled surface charge density $s$ based on RMA (upper three diagrams) and QCA (lower three diagrams). The lower three diagrams of Figures 4, 5 are identical. Each diagram shows 11 curves for $\omega =0$ (uppermost curve) to $\omega =1$ (black curve on the bottom), with $\omega$ changing in steps of 0.1. Diagrams refer to $r_0=0$ (left column, colored purple), $r_0=l$ (middle column, colored red), and $r_0=3l$ (right column, colored green). All curves are calculated for $z=6$. All models are based on solving Equations 16, 17, yet with the relationship $\eta (\mathrm{\Psi })$ emerging from Equation 38 for the upper three diagrams and from Equation 13 for the lower three diagrams.

Because the presence of the bell-to-camel shape transition is independent of $r_0$, the reason of its existence for a compact, solvent-free ionic liquid reflects the role of electrostatic ion–ion correlations. Correlated clusters of ions increase the decay length of the potential ${\mathrm{\Psi }}_{\mathit{lin}}(x)$, as shown in Figure 3. Hence, if the system of electrode and diffuse layer of counterions is represented by an effective parallel-plate capacitor, the distance between the capacitor plates is larger in the presence of electrostatic ion–ion correlations. The correspondingly smaller value of ${\bar{C}}_{\mathit{diff}}^{\mathit{lin}}$ is consistent with the existence of a transition from the bell-shaped to a camel-shaped profile of ${\bar{C}}_{\mathit{diff}}$ when electrostatic ion–ion correlations are accounted for.

4 Conclusion

The classical lattice-based mean-field theory for a densely packed, solvent-free ionic liquid employs continuum electrostatics and ignores ion–ion correlations, leading to a bell-shaped profile according to ${\bar{C}}_{\mathit{diff}}=\sqrt{1-e^{-s^2}}/|s|$ for the scaled differential capacitance as a function of the scaled charge density $s$. The present work is an attempt to analyze how ${\bar{C}}_{\mathit{diff}}$ changes if charge discreteness and electrostatic correlations between ion pairs are accounted for. To account for charge discreteness, we employ a truncated Coulomb potential (shown in Figure 1) to model the long-range part of the electrostatic interactions but exclude those of a given ion with its nearest neighbors. These short-range interactions are modeled explicitly as interactions between individual point charges. Electrostatic correlations between pairs of neighboring ions are accounted for using the quasi-chemical approximation (QCA) approach. Downing et al. (2018b) included QCA similarly to the present approach but without truncating the Coulomb potential, which led to a double counting of electrostatic interactions between nearest neighbor ions. The present work is no longer subject to this inconsistency. We emphasize that our approach is still subject to approximations: the truncation of the Coulomb potential is modeled only on the level of fourth-order electrostatics (${\mathcal{A}}_1$ in Equation 7), and QCA ignores correlations of higher order and between pairs of ions separated by distances larger than one lattice spacing. Yet, these approximations can, in principle, be overcome by considering sixth- or higher-order electrostatics (${\mathcal{A}}_2$ or higher in Equation 7) and by incorporating larger clusters into QCA (Bossa et al., 2015).

Our approximation of the full self-consistency relationship, $\mathcal{A}\mathrm{\Psi }=-\eta (\mathrm{\Psi })/\nu$, by the fourth-order differential equation, ${\mathcal{A}}_1\mathrm{\Psi }=-\eta (\mathrm{\Psi })/\nu$, renders our model structurally equivalent to the BSK model (Bazant et al., 2011). An important aspect of our study includes the boundary conditions associated with that equation (Equation 17), which emerge as part of the functional minimization of the free energy (Equations 12, 19) and differ from previously proposed sets of boundary conditions (Bazant et al., 2011; Gupta et al., 2020b; de Souza and Bazant, 2020; Bossa and May, 2020).

A major result of our study is analytic expressions for the coefficients ${\bar{C}}_{\mathit{diff}}^{\mathit{lin}}$ and $\mathrm{\Omega }$ in an expansion of ${\bar{C}}_{\mathit{diff}}(s)={\bar{C}}_{\mathit{diff}}^{\mathit{lin}}+s^2\mathrm{\Omega }/2$ up to quadratic order in $s$. We have obtained ${\bar{C}}_{\mathit{diff}}^{\mathit{lin}}$ from solving the linearized model and $\mathrm{\Omega }$ from performing a perturbation approach into the non-linear regime. The quadratic expression of ${\bar{C}}_{\mathit{diff}}(s)$ reveals the existence of a transition from a bell-shaped to a camel-shaped profile of the differential capacitance as a function of the electrostatic interaction strength $\omega$ between neighboring ions. However, such a transition exists only in the presence of electrostatic correlations, not in their absence. Hence, we conclude that electrostatic correlations between ion pairs are able to turn the bell-shaped profile predicted by mean-field theory, ${\bar{C}}_{\mathit{diff}}=\sqrt{1-e^{-s^2}}/|s|$, into a camel-shaped profile. Ion discreteness, on the other hand, enhances the magnitude but does not qualitatively alter the profile of ${\bar{C}}_{\mathit{diff}}(s)$. Our analytic predictions agree with numerical calculations of the full profiles for ${\bar{C}}_{\mathit{diff}}(s)$, i.e., beyond the quadratic regime, that we have presented in Figures 4, 5.

Data availability statement

The original contributions presented in the study are included in the article/supplementary material; further inquiries can be directed to the corresponding author.

Author contributions

GB: conceptualization, formal analysis, writing–original draft, and writing–review and editing. SM: conceptualization, formal analysis, writing–original draft, and writing–review and editing.

Funding

The author(s) declare that financial support was received for the research, authorship, and/or publication of this article. Author G.V.B. thanks Fondecyt (grant no. 11240064) and VIDCA/UACh (grant INS-INV 2022-09) for the support.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Generative AI statement

The author(s) declare that no Generative AI was used in the creation of this manuscript.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

Armand, M., Endres, F., MacFarlane, D. R., Ohno, H., and Scrosati, B. (2009). Ionic-liquid materials for the electrochemical challenges of the future. Nat. Mater. 8, 621–629. doi:10.1038/nmat2448

PubMed Abstract | CrossRef Full Text | Google Scholar

Avni, Y., Adar, R. M., and Andelman, D. (2020). Charge oscillations in ionic liquids: a microscopic cluster model. Phys. Rev. E 101, 010601. doi:10.1103/physreve.101.010601

PubMed Abstract | CrossRef Full Text | Google Scholar

Baker, N. A., Hünenberger, P. H., and McCammon, J. A. (1999). Polarization around an ion in a dielectric continuum with truncated electrostatic interactions. J. Chem. Phys. 110, 10679–10692. doi:10.1063/1.479013

CrossRef Full Text | Google Scholar

Bazant, M. Z., Storey, B. D., and Kornyshev, A. A. (2011). Double layer in ionic liquids: overscreening versus crowding. Phys. Rev. Lett. 106, 046102. doi:10.1103/physrevlett.106.046102

PubMed Abstract | CrossRef Full Text | Google Scholar

Bhuiyan, L. B., Lamperski, S., Wu, J., and Henderson, D. (2012). Monte Carlo simulation for the double layer structure of an ionic liquid using a dimer model: a comparison with the density functional theory. J. Phys. Chem. B 116, 10364–10370. doi:10.1021/jp304362y

PubMed Abstract | CrossRef Full Text | Google Scholar

Blossey, R., Maggs, A., and Podgornik, R. (2017). Structural interactions in ionic liquids linked to higher-order Poisson-Boltzmann equations. Phys. Rev. E 95, 060602. doi:10.1103/physreve.95.060602

PubMed Abstract | CrossRef Full Text | Google Scholar

Borukhov, I., Andelman, D., and Orland, H. (1997). Steric effects in electrolytes: a modified Poisson-Boltzmann equation. Phys. Rev. Lett. 79, 435–438. doi:10.1103/physrevlett.79.435

CrossRef Full Text | Google Scholar

Bossa, G. V., and May, S. (2020). Stability of ionic liquid modeled by composite Coulomb-Yukawa potentials. Phys. Rev. Res. 2, 032040. doi:10.1103/physrevresearch.2.032040

CrossRef Full Text | Google Scholar

Bossa, G. V., Roth, J., and May, S. (2015). Modeling lipid–lipid correlations across a bilayer membrane using the quasi-chemical approximation. Langmuir 31, 9924–9932. doi:10.1021/acs.langmuir.5b01719

PubMed Abstract | CrossRef Full Text | Google Scholar

Budkov, Y. A. (2020). Statistical field theory of ion–molecular solutions. Phys. Chem. Chem. Phys. 22, 14756–14772. doi:10.1039/d0cp02432e

PubMed Abstract | CrossRef Full Text | Google Scholar

Chao, H., and Wang, Z.-G. (2020). Effects of surface transition and adsorption on ionic liquid capacitors. J. Phys. Chem. Lett. 11, 1767–1772. doi:10.1021/acs.jpclett.0c00023

PubMed Abstract | CrossRef Full Text | Google Scholar

Coles, S. W., Park, C., Nikam, R., Kanduc, M., Dzubiella, J., and Rotenberg, B. (2020). Correlation length in concentrated electrolytes: insights from all-atom molecular dynamics simulations. J. Phys. Chem. B 124, 1778–1786. doi:10.1021/acs.jpcb.9b10542

PubMed Abstract | CrossRef Full Text | Google Scholar

Davis, H. T. (1996). “Statistical mechanics of phases, interfaces, and thin films,” in Advances in interfacial engineering series (VCH).

Google Scholar

Del Pópolo, M. G., and Voth, G. A. (2004). On the structure and dynamics of ionic liquids. J. Phys. Chem. B 108, 1744–1752. doi:10.1021/jp0364699

CrossRef Full Text | Google Scholar

Démery, V., Dean, D. S., Hammant, T. C., Horgan, R. R., and Podgornik, R. (2012). The one-dimensional Coulomb lattice fluid capacitor. J. Chem. Phys. 137, 064901. doi:10.1063/1.4740233

PubMed Abstract | CrossRef Full Text | Google Scholar

de Souza, J. P., and Bazant, M. Z. (2020). Continuum theory of electrostatic correlations at charged surfaces. J. Phys. Chem. C 124, 11414–11421. doi:10.1021/acs.jpcc.0c01261

CrossRef Full Text | Google Scholar

Downing, R., Berntson, B. K., Bossa, G. V., and May, S. (2018a). Differential capacitance of ionic liquids according to lattice-gas mean-field model with nearest-neighbor interactions. J. Chem. Phys. 149, 204703. doi:10.1063/1.5047490

PubMed Abstract | CrossRef Full Text | Google Scholar

Downing, R., Bossa, G. V., and May, S. (2018b). The role of ion–ion correlations for the differential capacitance of ionic liquids. J. Phys. Chem. C 122, 28537–28544. doi:10.1021/acs.jpcc.8b09756

CrossRef Full Text | Google Scholar

Eyckens, D. J., and Henderson, L. C. (2019). A review of solvate ionic liquids: physical parameters and synthetic applications. Front. Chem. 7, 263. doi:10.3389/fchem.2019.00263

PubMed Abstract | CrossRef Full Text | Google Scholar

Fedorov, M. V., and Kornyshev, A. A. (2008). Towards understanding the structure and capacitance of electrical double layer in ionic liquids. Electroch. Acta 53, 6835–6840. doi:10.1016/j.electacta.2008.02.065

CrossRef Full Text | Google Scholar

Frydel, D. (2016). Mean field electrostatics beyond the point charge description. Adv. Chem. Phys. 160, 209–260. doi:10.1002/9781119165156.ch4

CrossRef Full Text | Google Scholar

Goodwin, Z. A., de Souza, J. P., Bazant, M. Z., and Kornyshev, A. A. (2023). “Mean-field theory of the electrical double layer in ionic liquids,” in Encyclopedia of ionic liquids (Springer), 837–850.

CrossRef Full Text | Google Scholar

Goodwin, Z. A., Feng, G., and Kornyshev, A. A. (2017). Mean-field theory of electrical double layer in ionic liquids with account of short-range correlations. Electroch. Acta 225, 190–197. doi:10.1016/j.electacta.2016.12.092

CrossRef Full Text | Google Scholar

Gupta, A., Govind Rajan, A., Carter, E. A., and Stone, H. A. (2020a). Ionic layering and overcharging in electrical double layers in a Poisson-Boltzmann model. Phys. Rev. Lett. 125, 188004. doi:10.1103/physrevlett.125.188004

PubMed Abstract | CrossRef Full Text | Google Scholar

Gupta, A., Govind Rajan, A., Carter, E. A., and Stone, H. A. (2020b). Thermodynamics of electrical double layers with electrostatic correlations. J. Phys. Chem. C 124, 26830–26842. doi:10.1021/acs.jpcc.0c08554

CrossRef Full Text | Google Scholar

Hayes, R., Warr, G. G., and Atkin, R. (2015). Structure and nanostructure in ionic liquids. Chem. Rev. 115, 6357–6426. doi:10.1021/cr500411q

PubMed Abstract | CrossRef Full Text | Google Scholar

Henderson, D., Lamperski, S., Jin, Z., and Wu, J. (2011). Density functional study of the electric double layer formed by a high density electrolyte. J. Phys. Chem. B 115, 12911–12914. doi:10.1021/jp2078105

PubMed Abstract | CrossRef Full Text | Google Scholar

Khlyupin, A., Nesterova, I., and Gerke, K. (2023). Molecular scale roughness effects on electric double layer structure in asymmetric ionic liquids. Electrochim. Acta 450, 142261. doi:10.1016/j.electacta.2023.142261

CrossRef Full Text | Google Scholar

Kondrat, S., Feng, G., Bresme, F., Urbakh, M., and Kornyshev, A. A. (2023). Theory and simulations of ionic liquids in nanoconfinement. Chem. Rev. 123, 6668–6715. doi:10.1021/acs.chemrev.2c00728

PubMed Abstract | CrossRef Full Text | Google Scholar

Kondrat, S., Georgi, N., Fedorov, M. V., and Kornyshev, A. A. (2011). A superionic state in nano-porous double-layer capacitors: insights from Monte Carlo simulations. Phys. Chem. Chem. Phys. 13, 11359–11366. doi:10.1039/c1cp20798a

PubMed Abstract | CrossRef Full Text | Google Scholar

Kornyshev, A. A. (2007). Double-layer in ionic liquids: paradigm change? J. Phys. Chem. B 111, 5545–5557. doi:10.1021/jp067857o

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, A. A., Kondrat, S., Vella, D., and Goriely, A. (2015). Dynamics of ion transport in ionic liquids. Phys. Rev. Lett. 115, 106101. doi:10.1103/physrevlett.115.106101

PubMed Abstract | CrossRef Full Text | Google Scholar

Linse, P., and Andersen, H. C. (1986). Truncation of Coulombic interactions in computer simulations of liquids. J. Chem. Phys. 85, 3027–3041. doi:10.1063/1.451011

CrossRef Full Text | Google Scholar

Ma, K., Woodward, C. E., and Forsman, J. (2014). Classical density functional study on interfacial structure and differential capacitance of ionic liquids near charged surfaces. J. Phys. Chem. C 118, 15825–15834. doi:10.1021/jp504001u

CrossRef Full Text | Google Scholar

May, S. (2019). Differential capacitance of the electric double layer: mean-field modeling approaches. Curr. Opin. Electrochem. 13, 125–131. doi:10.1016/j.coelec.2018.12.002

CrossRef Full Text | Google Scholar

Moon, G. J., Ahn, M. M., and Kang, I. S. (2015). Osmotic pressure of ionic liquids in an electric double layer: prediction based on a continuum model. Phys. Rev. E 92, 063020. doi:10.1103/physreve.92.063020

PubMed Abstract | CrossRef Full Text | Google Scholar

Nakamura, I. (2015). Dipolar self-consistent field theory for ionic liquids: effects of dielectric inhomogeneity in ionic liquids between charged plates. J. Phys. Chem. C 119, 7086–7094. doi:10.1021/jp511770r

CrossRef Full Text | Google Scholar

Nesterova, I., Evstigneev, N. M., Ryabkov, O. I., Gerke, K. M., and Khlyupin, A. (2025). Mechanism of overscreening breakdown by molecular-scale electrode surface morphology in asymmetric ionic liquids. J. Colloid Interface Sci. 677, 396–405. doi:10.1016/j.jcis.2024.08.040

PubMed Abstract | CrossRef Full Text | Google Scholar

Nordness, O., and Brennecke, J. F. (2020). Ion dissociation in ionic liquids and ionic liquid solutions. Chem. Rev. 120, 12873–12902. doi:10.1021/acs.chemrev.0c00373

PubMed Abstract | CrossRef Full Text | Google Scholar

Shalabi, A., Daniels, L., Scott, M., and Mišković, Z. (2019). Differential capacitance of ionic liquid interface with graphene: the effects of correlation and finite size of ions. Electrochim. Acta 319, 423–434. doi:10.1016/j.electacta.2019.06.171

CrossRef Full Text | Google Scholar

Sher, A., van Schilfgaarde, M., Chen, A.-B., and Chen, W. (1987). Quasichemical approximation in binary alloys. Phys. Rev. B 36, 4279–4295. doi:10.1103/physrevb.36.4279

PubMed Abstract | CrossRef Full Text | Google Scholar

Silva, W., Zanatta, M., Ferreira, A. S., Corvo, M. C., and Cabrita, E. J. (2020). Revisiting ionic liquid structure-property relationship: a critical analysis. Int. J. Mol. Sci. 21, 7745. doi:10.3390/ijms21207745

PubMed Abstract | CrossRef Full Text | Google Scholar

Sowmiah, S., Srinivasadesikan, V., Tseng, M.-C., and Chu, Y.-H. (2009). On the chemical stabilities of ionic liquids. Molecules 14, 3780–3813. doi:10.3390/molecules14093780

PubMed Abstract | CrossRef Full Text | Google Scholar

Storey, B. D., and Bazant, M. Z. (2012). Effects of electrostatic correlations on electrokinetic phenomena. Phys. Rev. E 86, 056303. doi:10.1103/physreve.86.056303

PubMed Abstract | CrossRef Full Text | Google Scholar

Torimoto, T., Tsuda, T., Okazaki, K.-i., and Kuwabata, S. (2010). New frontiers in materials science opened by ionic liquids. Adv. Mater. 22, 1196–1221. doi:10.1002/adma.200902184

PubMed Abstract | CrossRef Full Text | Google Scholar

Varner, S., and Wang, Z.-G. (2022). Effects of dilution in ionic liquid supercapacitors. Phys. Chem. Chem. Phys. 24, 27362–27374. doi:10.1039/d2cp03398d

PubMed Abstract | CrossRef Full Text | Google Scholar

Walden, P. (1914). Molecular weights and electrical conductivity of several fused salts. Bull. Acad. Imper. Sci. St. Petersbg. 1800.

Google Scholar

Wang, Y., Jiang, W., Yan, T., and Voth, G. A. (2007). Understanding ionic liquids through atomistic and coarse-grained molecular dynamics simulations. Acc. Chem. Res. 40, 1193–1199. doi:10.1021/ar700160p

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, Y.-L., Li, B., Sarman, S., Mocci, F., Lu, Z.-Y., Yuan, J., et al. (2020). Microstructural and dynamical heterogeneities in ionic liquids. Chem. Rev. 120, 5798–5877. doi:10.1021/acs.chemrev.9b00693

PubMed Abstract | CrossRef Full Text | Google Scholar

Watanabe, M., Thomas, M. L., Zhang, S., Ueno, K., Yasuda, T., and Dokko, K. (2017). Application of ionic liquids to energy storage and conversion materials and devices. Chem. Rev. 117, 7190–7239. doi:10.1021/acs.chemrev.6b00504

PubMed Abstract | CrossRef Full Text | Google Scholar

Weingärtner, H. (2008). Understanding ionic liquids at the molecular level: facts, problems, and controversies. Angew. Chem. Int. Ed. 47, 654–670. doi:10.1002/anie.200604951

CrossRef Full Text | Google Scholar

Wishart, J. F. (2009). Energy applications of ionic liquids. Energy Environ. Sci. 2, 956–961. doi:10.1039/b906273d

CrossRef Full Text | Google Scholar

Wu, J., Jiang, T., Jiang, D., Jin, Z., and Henderson, D. (2011). A classical density functional theory for interfacial layering of ionic liquids. Soft Matter 7, 11222–11231. doi:10.1039/c1sm06089a

CrossRef Full Text | Google Scholar

Ye, B. B., and Wang, Z.-G. (2022). A coarse-grained model of room-temperature ionic liquids between metal electrodes: a molecular dynamics study. Phys. Chem. Chem. Phys. 24, 11573–11584. doi:10.1039/d2cp00166g

PubMed Abstract | CrossRef Full Text | Google Scholar

Yochelis, A. (2014). Transition from non-monotonic to monotonic electrical diffuse layers: impact of confinement on ionic liquids. Phys. Chem. Chem. Phys. 16, 2836–2841. doi:10.1039/c3cp55002h

PubMed Abstract | CrossRef Full Text | Google Scholar

Zeman, J., Kondrat, S., and Holm, C. (2021). Ionic screening in bulk and under confinement. J. Chem. Phys. 155, 204501. doi:10.1063/5.0069340

PubMed Abstract | CrossRef Full Text | Google Scholar