# Image charge effects under metal and dielectric boundary conditions

Tingtao Zhou<EMAIL_ADDRESS>Divisions of Engineering and Applied
Sciences and Chemistry and Chemical Engineering, California Institute of
Technology, Pasadena, California 91125, USA  Dorian Bruch  Division of
Chemistry and Chemical Engineering, California Institute of Technology,
Pasadena, California 91125, USA  Zhen-Gang Wang<EMAIL_ADDRESS>Division of Chemistry and Chemical Engineering, California Institute of
Technology, Pasadena, California 91125, USA

###### Abstract

Image charge effect is a fundamental problem in electrostatics. However, a
proper treatment at the continuum level for many-ion systems, such as
electrolyte solutions or ionic liquids, remains an open theoretical question.
Here, we demonstrate and systematically compare the image charge effects under
metal and dielectric boundary conditions (BCs), based on a renormalized
Gaussian-fluctuation theory. Our calculations for a simple 1:1 symmetric
electrolyte in the point-charge approximation show that the double-layer
structure, capacitance, and interaction forces between like-charged plates
depend strongly on the types of boundaries, even in the weak-coupling regime.
Like-charge attraction is predicted for both metal and dielectric BCs.
Finally, we comment on the effects of a dielectrically-saturated solvent layer
on the metal surface. We provide these results to serve as a baseline for
comparison with more realistic molecular dynamics simulations and experiments.

##  I Introduction

Image charge (IC) interaction has wide implications in electrostatics problems
with interfaces. A classic example is the study on surface tension of
electrolyte solutions pioneered by Wagner [1] and Onsager and Samaras [2],
where the IC repulsion at air/water interface was considered responsible for
the increase in the surface tension due to added salt ions. Monte Carlo
simulations [3, 4] and an integral equation theory based on hypernetted chain
approximation [5] have shown ion depletion near a weakly charged dielectric
interface, while ion accumulation near a metal interface is expected.
Nanoscale supercapacitors have received recent interest for nanoscale device
applications [6, 7, 8] and may be used to probe material characteristics at
the molecular level [9]. These nanoscale devices require proper treatment of
IC effects. IC interaction also affects the study of capacitive desalination
[10], capillary freezing of ionic liquids [11], surface adsorption of ions
[12], ion transport in charged nano-porous materials [13], and may even play a
role in planet formation [14].

The Poisson–Boltzmann (PB) theory has been widely accepted as the leading
order theory for the weak-coupling regime of electrostatics [3, 15, 16, 17,
18], but it does not include ion-correlation or IC effect. Theoretically, a
primitive treatment is to explicitly construct the image charge of a single
ion and modify the free energy at the interface [2]. In field-theoretic
formulations, perturbative loop-expansion treatments have been used to derive
corrections to the PB theory [16, 19, 20] to include IC interactions. However,
by analysis of a single point charge approaching a bounding interface, one can
show that the image charge interaction diverges as the charge approaches the
interface [21]. For a dielectric interface, this singular behavior results in
a singular boundary layer of the electrolyte close to the interface, eluding
regular perturbation methods [15]. Ref. 22 derived a renormalized field-theory
framework with implicit solvent. There, the solvent is modeled as a continuum
with an effective dielectric constant. Sharp dielectric interfaces have been
studied using this framework [23, 24] and lattice Monte Carlo simulations [25]
have shown good agreement with the predictions from Ref. 23. Other field-
theoretic treatments including IC interactions usually introduce an arbitrary
cutoff to avoid divergence of self-energy [26, 27]. For the metal boundary,
several techniques for Molecular Dynamics (MD) simulations have been devised
to account for the IC interactions, such as iterative methods [28, 29, 30, 31]
or explicit construction of image charges [32, 33, 34, 35, 36], and ion
accumulation on the metal surfaces has been observed [37]. Simulations using
the method of periodic Green’s functions have also been constructed for both
dielectric and metal boundaries [37]. In all of these simulations, the solvent
is treated implicitly where the polarization of the solvent is modeled by a
bulk dielectric constant. However, as found in Ref. 38, IC interaction is
strongly affected by the polarizability of both solvent molecules and ions — a
polarizable solvent even cancels out some effects of the IC, which is
neglected in models with implicit solvent background.

IC effects are known to be important for electric double layer capacitors. For
example, experimental evidence shows that surface polarizability leads to a
much higher differential capacitance than expected [39], and IC interaction
may even lead to a phase transition at the electrolyte-metal interface [40,
41]. For room-temperature ionic liquids, MD simulations and 1D model analysis
[42, 43] showed importance of IC effects on capacitance curves. For lower bulk
salt concentrations, recent MD simulations [38] have shown significant
enhancement of charge separation in a nano-supercapacitor with conducting
electrodes. MD simulations also showed [44] non-trivial properties of double-
layer structures of polyelectrolytes due to surface polarization for both
dielectric and metal boundaries.

IC effects also influence colloidal interactions [45, 46], which play a
critical role in many soft matter systems. The classic DLVO theory [47]
predicts only screened repulsive electrostatic interactions. More recent
efforts extending the PB theory, such as the renormalized jellium model [48,
49, 50, 51, 52], result in a renormalized charge but qualitatively same
repulsive tendency. In fact, Trizac has shown that generally local density
approximations lead to repulsive pair potentials [53]. Nevertheless, like-
charged colloidal attraction have been reported [54, 55, 56] and has wide
implications. For example, during cement setting, attraction between like-
charged surfaces may allow the cohesion strength to increase beyond the order
of magnitude of capillary stress [57, 58, 59] or van der Waals interaction.
Without systematic considerations of IC effects, the like-charge attraction
has been mostly attributed to ion-ion correlations and fluctuations [5, 60,
61, 62, 63, 64, 15, 65, 66]. These correlation effects are generally thought
to become significant at the strong-coupling limit [67, 68, 69, 70], typically
from high surface charge density or multi-valent ions. In the weak-coupling
regime, Ref. [23] has shown like-charge attraction due to IC-induced ion
depletion between dielectric plates. For 1:1 electrolyte with metal boundary
conditions and fixed surface charge density, dos Santos & Levin [71] concluded
that attraction between like-charged particles can happen for spheres even at
the level of the PB theory.

In reality, IC interaction can be further complicated by various other
factors. These interfacial effects are often ion-specific [72], such as the
famous Hofmeister series [73] for protein salt-in/salt-out. The hydration
shell structure of the ions also changes when they are adsorbed onto a metal
surface, which affects their hydration energy. In addition, more complex and
unexpected behaviors of the differential capacitance have been reported in MD
simulations with polyelectrolytes [44]. In the case of a free interface, the
IC interaction is also coupled with surface capillary waves [74, 75, 76].
These complications are beyond the scope of this work.

In this paper, we provide a simple theoretical picture of the IC effects as a
reference for further studies with more realistic molecular models for
electrolytes. To this end, we treat all ions as point charges and assume an
implicit solvent background. We present a systematic comparison between the
metal BC, dielectric BC, and PB theory in terms of the differential
capacitance of nano-capacitors and forces between like-charged plates by
numerically solving the variational field theory for an incompressible
electrolyte. Our calculations show that the differential capacitance depends
strongly on the boundary condition (metal or dielectric) as well as the bulk
salt concentration. For example, even at low bulk concentration but if the
slit width is small, the metal BC leads to a “bird-shaped" capacitance curve
[77]. For a simple 1:1 electrolyte, we find the interplate forces change from
repulsion to attraction as the separation distance decreases between
dielectric plates, and stronger and long-ranged attractions between like-
charged metal plates, due to the IC interaction. Moreover, to account for the
distance of closest approach of real ions to a surface, we consider a thin
layer of polarizable solvent on the metal surface, preventing direct contact
of ions with the boundary. We show that even a very thin surface dielectric
layer cancels out some IC effects, consistent with the findings in Ref. 38.

##  II Theory

###  II.1 General Gaussian variational approach

We start by recapitulating the field-theoretic framework as presented in Ref.
22. Consider a system with cations (of valence $z_{+}$) and anions (of valence
$-z_{-}$) in an electrolyte solution, and external charges on the boundary
surfaces $e\rho_{ex}$. The solution is connected to a reservoir with bulk
cation and anion concentrations $c_{+,0}$ and $c_{-,0}$, respectively. The
total microscopic charge density is

$e\rho(\mathbf{r})=e\left(\rho_{ex}(\mathbf{r})+z_{+}\sum_{i+}h_{+}(\mathbf{r}-\mathbf{r}_{+}^{i})-z_{-}\sum_{j-}h_{-}(\mathbf{r}-\mathbf{r}_{-}^{j})\right)$
(1)

where the salt ions have a charge spread kernel, which is taken to the point
charge limit
$h_{\pm}(\mathbf{r}-\mathbf{r}_{\pm}^{i})=\delta(\mathbf{r}-\mathbf{r}_{\pm}^{i})$
throughout this work. The Coulomb energy is

$H=\frac{e^{2}}{2}\int
d\mathbf{r}d\mathbf{r}^{\prime}\rho(\mathbf{r})G_{0}(\mathbf{r},\mathbf{r}^{\prime})\rho(\mathbf{r}^{\prime})$
(2)

where the Coulomb operator $G_{0}$ is defined by

$-\nabla\cdot\varepsilon(\mathbf{r})\nabla
G_{0}(\mathbf{r},\mathbf{r}^{\prime})=\delta(\mathbf{r}-\mathbf{r}^{\prime})$
(3)

The grand canonical partition function is

$\begin{split}\Xi=&\sum_{n_{+}}\sum_{n_{-}}\frac{e^{n_{+}\mu_{+}}e^{n_{-}\mu_{-}}}{n_{+}!n_{-}!v_{+}^{n_{+}}v_{-}^{n_{-}}}\int
d\mathbf{r}_{i+}d\mathbf{r}_{j-}e^{-\beta H}\\\
=&\frac{1}{\sqrt{\det(G_{0})}}\int D[\phi]e^{-\int
d\mathbf{r}\mathcal{L}[\phi(\mathbf{r})]}\end{split}$ (4)

where the action density is

$\mathcal{L}[\phi]=\frac{1}{8\pi\ell_{B}}\left(\nabla\phi(\mathbf{r})\right)^{2}+i\rho_{ex}(\mathbf{r})\phi(\mathbf{r})-\lambda_{+}e^{-iz_{+}\phi(\mathbf{r})}+\lambda_{-}e^{iz_{-}\phi(\mathbf{r})}$
(5)

For convenience Equation. (5) and subsequent ones throughout the paper are
written in dimensionless variables. The Bjerrum length is
$\ell_{B}=e^{2}/(4\pi\epsilon k_{B}T)$. The complex-valued auxiliary field
$\phi(\mathbf{r})$ conjugate to the charge density $\rho(\mathbf{r})$ is
introduced through the standard Hubbard–Stratonovich transformation [78, 79]
and then non-dimensionalized by $k_{B}T/e$.
$\lambda_{\pm}=\frac{e^{\mu_{\pm}}}{v_{\pm}}$ are the bulk activities of salt
ions, where $\mu_{\pm}$ are the chemical potentials and $v_{\pm}$ are the
volume scales which have no thermodynamic consequences other than shifting the
reference of the chemical potential; therefore, for convenience, we choose
$v_{\pm}$ to be simply the volume of the ions. A renormalization of the field
theory then can be derived through a variational approach by extremizing [22,
80] the grand free energy

$W\leq W_{ref}+\left<L[\phi]-L_{ref}[\phi]\right>_{ref}$ (6)

The auxiliary field is decomposed into

$\phi(\mathbf{r})=-i\psi(\mathbf{r})+\chi(\mathbf{r})$ (7)

where the mean-field part $-i\psi(\mathbf{r})$ and fluctuation parts
$\chi(\mathbf{r})$ are purely imaginary and real fields, respectively. The
reference action is chosen as a Gaussian form for the fluctuation part for the
auxiliary field

$L_{ref}=\frac{1}{2}\int
d\mathbf{r}d\mathbf{r}^{\prime}\chi(\mathbf{r})G^{-1}(\mathbf{r},\mathbf{r}^{\prime})\chi(\mathbf{r}^{\prime})$
(8)

The Green’s function $G(\mathbf{r},\mathbf{r}^{\prime})$ is unknown and will
be determined together with the mean-field $\psi(\mathbf{r})$ by the
extremization of $W$, which results in

$\displaystyle-\nabla\cdot\left[\frac{1}{4\pi\ell_{B}(\textbf{r})}\nabla\psi(\textbf{r})\right]=\rho_{ex}(\textbf{r})+z_{+}c_{+}(\textbf{r})-z_{-}c_{-}(\textbf{r})$
(9) $\displaystyle-\nabla\cdot\left[\frac{1}{4\pi\ell_{B}(\textbf{r})}\nabla
G(\textbf{r},\textbf{r}^{\prime})\right]+2I(\mathbf{r})G(\textbf{r},\textbf{r}^{\prime})=\delta(\textbf{r}-\textbf{r}^{\prime})$
(10)

where the ionic strength is
$I(\mathbf{r})=\frac{1}{2}\left[z_{+}^{2}c_{+}(\mathbf{r})+z_{-}^{2}c_{-}(\mathbf{r})\right]$,
and ion concentrations are

$c_{\pm}(\textbf{r})=\lambda_{\pm}\exp\left\\{\mp
z_{\pm}\psi(\textbf{r})-u_{\pm}(\mathbf{r})\right\\}$ (11)

with the self energy

$u_{\pm}(\mathbf{r})=\frac{1}{2}z_{\pm}^{2}G(\textbf{r},\textbf{r}^{\prime}\rightarrow\mathbf{r})$
(12)

Here the self-energy appears in the exponential of the Boltzmann factor of the
charge distributions and it is the key feature of this framework. It includes
a singular part, which is present even in a homogeneous solution without
boundaries; this singular part can be regularized by the introduction of a
smearing function [22], with a spread chosen to reproduce the Born solvation
energy in a bulk electrolyte solution and can be absorbed into the bulk ion
chemical potentials. The remaining part of the self-energy comes from spatial
variations of charge distribution, correlation energy, and the existence of
boundaries and image charges. This part of the self-energy is finite inside a
continuum electrolyte, except at dielectric discontinuities and boundaries.

###  II.2 Incompressible electrolyte

Even within the PB framework, ion densities near electrolyte–metal interfaces
can be overestimated and exceeds the physically maximum value of close packing
[81], due to the unconstrained exponential dependence of ion concentration on
the electric potential. The point-charge approximation in the variational
field theory presented in the last section still exhibits divergent self
energy near an interface [21], resulting in complete ion depletion (divergent
accumulation) near a dielectric (metal) wall. Here we incorporate the
incompressibility constraint [82, 83] into the field-theoretic framework to
regulate ion concentration near the boundary surfaces. This way, the IC
interaction is still calculated for point charges that can approach the
interface indefinitely close, and hence the IC force may diverge as expected,
but the ion concentration is regulated by a saturation value. The
incompressibility condition is represented by a product of delta functions
that enforces the solution density $\frac{1}{v_{0}}$ at every point in space,
$\prod_{\textbf{r}}\delta\left(\hat{c}_{+}(\mathbf{r})+\hat{c}_{-}(\mathbf{r})+\hat{c}_{s}(\mathbf{r})-\frac{1}{v_{0}}\right)$,
where the number densities of cation, anion and solvent molecules are

$\begin{split}\hat{c}_{+}(\mathbf{r})=&\sum_{i+}\delta(\mathbf{r}-\mathbf{r}_{+}^{i})\\\
\hat{c}_{-}(\mathbf{r})=&\sum_{i-}\delta(\mathbf{r}-\mathbf{r}_{-}^{i})\\\
\hat{c}_{s}(\mathbf{r})=&\sum_{s-}\delta(\mathbf{r}-\mathbf{r}_{s}^{i})\end{split}$
(13)

For simplicity we have assumed that the volumes of all species (cation, anion,
and solvent) are equal $v_{+}=v_{-}=v_{s}=v_{0}=(0.5~{}\text{nm})^{3}$, so
that the total concentration of all species is $c_{0}\approx 13.33$ M
everywhere.

To deal with the incompressibility constraint, we represent the product of
delta functions as its Fourier transform by introducing an auxiliary field
$\eta(\mathbf{r})$

$\begin{split}\prod_{\textbf{r}}\delta\left(\hat{c}_{+}(\mathbf{r})+\hat{c}_{-}(\mathbf{r})+\hat{c}_{s}(\mathbf{r})-\frac{1}{v_{0}}\right)&=\\\
\int D[\eta(\mathbf{r})]&e^{i\int
d\mathbf{r}\eta(\hat{c}_{+}+\hat{c}_{-}+\hat{c}_{s}-\frac{1}{v_{0}})}\end{split}$
(14)

where the notation $\int D[\eta(\mathbf{r})]$ is the functional integral with
respect to the field $\eta(\mathbf{r})$. Now $\hat{c}_{\pm,s}(\mathbf{r})$
still depend on the ion/solvent positions $\\{\mathbf{r}^{i}\\}$, which are
integrated over in the grand partition function.

With this constraint incorporated into the grand partition function,
extremizing the grand free energy with respect to the new auxiliary fields
gives modified equations for concentration fields

$\displaystyle c_{\pm}(\textbf{r})=$
$\displaystyle\lambda_{\pm}e^{\mp(z_{\pm}\psi(\textbf{r}))+\eta(\textbf{r})-u_{\pm}(\textbf{r})}$
(15) $\displaystyle c_{s}(\textbf{r})=$
$\displaystyle\frac{e^{\mu_{s}}}{v_{0}}e^{\eta(\textbf{r})}$ (16)
$\displaystyle\frac{1}{v_{0}}=$ $\displaystyle
c_{+}(\textbf{r})+c_{-}(\textbf{r})+c_{s}(\textbf{r})$ (17)

where these Equations (15) and (16) are obtained from the integrand of
chemical potential derivatives of the grand free energy. Note that the
differential equations (9) – (10) for the electric potential $\psi$ and
Green’s function $G$ remain the same. Hence, the treatment for
incompressibility is at the mean-field level. For consistency in the
comparisons throughout this paper, we impose the incompressibility condition
for all cases—PB theory, metal and dielectric BCs, and for a thin dielectric
surface layer on a metal wall.

###  II.3 Force and capacitance between two plates

Consider an electrolyte solution confined between two infinite parallel
plates. In this work, we focus on comparing four cases of boundary conditions
(BCs): (1) both plates are perfect metal; (2) both plates are dielectric
medium with a dielectric constant $\epsilon_{p}=2$, lower than that of the
solution, assumed to be $\epsilon_{r}=80$; (3) PB theory without image charge
corrections; and (4) both plates are metal with a thin dielectric surface
layer with a dielectric constant $\epsilon_{\text{L}}$, lower than that of
solution. For simplicity, cases (2) and (4) will be referred to as dielectric
and layer BC from now on, respectively. The electrolyte is connected to a bulk
1:1 solution of salt ion concentrations $c_{+,0}=c_{-,0}=c_{0}$; hence, the
calculations are in the grand canonical ensemble.

Due to the translational and rotational symmetry in the directions parallel to
the plates, it is convenient to use cylindrical coordinates ($\rho,\theta,z$).
The Equations (9), (10) and (15) – (17) are simplified to a set of 1D ordinary
differential equations (ODEs) along the $z$-axis perpendicular to the plates

$\displaystyle-\frac{1}{4\pi\ell_{B}}\frac{\partial^{2}}{\partial
z^{2}}\psi(z)=\rho_{ex}(z)+c_{+}(z)-c_{-}(z)$ (18)
$\displaystyle\delta(z-z^{\prime})=\frac{1}{4\pi\ell_{B}}\left\\{\left[k^{2}+\kappa(z)^{2}\right]-\frac{\partial^{2}}{\partial
z^{2}}\right\\}G(k,z,z^{\prime})$ (19) $\displaystyle
c_{\pm}(z)=\lambda_{\pm}e^{\eta(z)-u(z)\mp\psi(z)}$ (20) $\displaystyle
u(z)=u_{Born}+u_{DH}+\frac{1}{4\pi}\int_{0}^{\infty}dk\left(kG(k,z,z)-2\pi\ell_{B}\right)$
(21) $\displaystyle\kappa(z)^{2}=4\pi\ell_{B}\left[c_{+}(z)+c_{-}(z)\right]$
(22) $\displaystyle c_{s}(z)=\lambda_{s}e^{\eta(z)}$ (23)
$\displaystyle\frac{1}{v_{0}}=c_{+}(z)+c_{-}(z)+c_{s}(z)$ (24)

where we use the partial Fourier transform $G(k,z,z^{\prime})$ of the Green
function $G(\mathbf{r},\mathbf{r}^{\prime})$: by symmetry of the geometry
$G(\textbf{r},\textbf{r}^{\prime})=G(|\rho-\rho^{\prime}|,z,z^{\prime})$. One
then only Fourier transforms the separation $s=|\rho-\rho^{\prime}|$ in the
parallel direction

$G(|\rho-\rho^{\prime}|,z,z^{\prime})=\frac{1}{2\pi}\int_{0}^{\infty}kdkJ_{0}(ks)G(k,z,z^{\prime})$
(25)

where $J_{0}(x)$ is the zeroth-order Bessel function. After the partial
Fourier transform, Equation (10) becomes Equation (19). Inverse Fourier
transforming $G(k,z,z^{\prime})$ back and taking the limit of
$\rho\rightarrow\rho^{\prime}$, $z\rightarrow z^{\prime}$ according to
Equation (12) gives Equation (21). In Equation (21), the singular part of the
self-energy is replaced with a constant reference energy, including the Born
energy $u_{Born}=z^{2}\ell_{B}/(2a)$ and a Debye-Hückel correlation term
$u_{DH}=-z^{2}\ell_{B}\kappa_{B}/2$, which is assumed homogeneous inside the
domain and in the reservoir, and hence can be absorbed into the bulk fugacity.
$\kappa_{B}=\sqrt{8\pi c_{0}\ell_{B}}$ is the inverse Debye length in the bulk
solution. Only the excess part (second term of RHS in Equation (21)) will be
relevant for our purpose here. By setting $\eta=0$ for the reservoir, the
chemical potentials are determined by the bulk concentrations

$\displaystyle\lambda_{s}=\frac{1}{v_{0}}-2c_{0}$ (26)
$\displaystyle\lambda_{\pm}e^{-u_{\infty}}=c_{\pm,0}=c_{0}$ (27)

The BCs for the ODE Eq. (18) can be set by fixing the surface potential
(Dirichlet) or fixing the surface charge density (Neumann), depending on the
specific application. We now specify the BCs for the ODE Eq. (19) for the
Green’s function at the surfaces $z=-L/2$ and $z=L/2$. With metal plates,
there is no fluctuation inside the metal, so
$G(z=-L/2,z^{\prime})=G(z=L/2,z^{\prime})=0$ (noting that
$G(\mathbf{r},\mathbf{r}^{\prime})=\left<\chi(\mathbf{r})\chi(\mathbf{r}^{\prime})\right>$
is the variance of the fluctuation field $\chi(\mathbf{r})$ from its defining
Equation. (8)). With dielectric plates of relative permitivity $\epsilon_{p}$,
different from the solvent relative permitivity $\epsilon_{r}$, the equation
can be extended to the whole space and the dielectric constant is regarded as
having a discontinuity at $z=\pm L/2$ so that $\epsilon(z)=\epsilon_{p}$
outside $|z|>L/2$ and $\epsilon(z)=\epsilon_{r}$ inside $|z|<L/2$. Accounting
for this jump one arrives at [24] a Robin BC for $G(k,z,z^{\prime})$

$\displaystyle\epsilon_{r}\partial_{z}G(k,z=\pm L/2,z^{\prime})$
$\displaystyle=\mp k\epsilon_{p}G(k,z=\pm L/2,z^{\prime})$ (28)

For the case of a thin dielectric layer on the metal surface, the metal
surface is at fixed potential, and Equations (18) and (19) have Robin boundary
conditions, derived in Appendix C.

The differential capacitance $C_{d}=d\sigma_{s}/dV$ can be readily obtained by
sweeping the voltage $V$, which is symmetrically applied on both plates such
that $\psi(z=\pm L/2)=V$. $\sigma_{s}$ is the surface charge density on a
plate. We also examine forces between like-charged plates by continuously
varying their separation distance and computing the total grand free energy
$W$. In any case, the free energy per unit area of the system is [24]

$\begin{split}W&=\frac{1}{2}\int
dz\psi(z)\left\\{\rho_{ex}(z)-c_{+}(z)+c_{-}(z)\right\\}\\\ &+\int
dz(\eta(z)-1)\left\\{c_{+}(z)+c_{-}(z)+c_{s}(z)\right\\}\\\ &+\int
dzI(z)\int_{0}^{1}[G(z,z;\xi)-G(z,z)]d\xi\end{split}$ (29)

for fixed surface charges $\sigma_{s}$, and performing a Legendre transform
gives the free energy $Y=W-2\sigma_{s}V$ for fixed surface potentials at
$\psi(z=\pm L/2)=V$. The “charging” Green’s function $G(z,z^{\prime};\xi)$ is
obtained by solving Eq. (19) but with $\kappa(z)^{2}$ replaced by
$\xi\kappa(z)^{2}$. The force per unit area is

$\Pi=-\frac{\partial\left(\mathcal{W}(L)-\mathcal{W}(\infty)\right)}{\partial
L}$ (30)

and $\mathcal{W}$ is $W$ or $Y$ depending on whether surface charge or surface
potential is specified [84].

##  III Numerical results and discussions

We numerically evaluate the Equations (18) – (24) iteratively. For the fixed
voltage (Dirichlet) boundary conditions, the Poisson equation (18) is solved
by the Chebyshev spectral method, and we combine the Chebyshev and shooting
methods to solve the fixed charge (Neumann) boundary condition. For a given
pair of $k$ and $z^{\prime}$ values, the Green’s function is solved by a
finite difference method with non-uniform grid size for accuracy. For the
inverse Fourier transform, the integration over $k$ in Equation (21) is
performed with Gauss–Laguerre quadrature.

Figure 1:  Charge accumulation/depletion at the planar boundary due to image
charge interactions. No applied voltage or fixed surface charge on the wall.
The bulk salt concentration is $c_{0}=0.1$ M for each case. The dielectric
constant of the implicit solvent is $\epsilon_{r}=80$; for the dielectric
plates, $\epsilon_{p}=2$. (a) Self-energy profiles corresponding to metal BC
(blue) and dielectric BC (red). (b) Concentration profiles: the blue line
represents the metal BC, corresponding to the left y-axis. The red line
represents the dielectric BC, corresponding to the right y-axis. The inset
shows the deviations from bulk concentration $c_{0}$ (absolute value
$|c/c_{0}-1|$) on the same log-scale y-axis.  Figure 2:  Finite-size effect on
the charge accumulation/depletion at the boundaries due to image charge
interactions without applied voltage or fixed charge. The slit pore width is
$L=1$ nm and the bulk salt concentration is $c_{0}=0.1$ M for each case. The
dielectric constant of the implicit solvent is $\epsilon_{r}=80$; for the
dielectric plates, $\epsilon_{p}=2$. (a) Self-energy profiles corresponding to
metal BC (blue) and dielectric BC (red). (b) Concentration profiles: the
cation and anion profiles overlap for either metal or dielectric BC since
there is no voltage applied. The blue line shows the metal BC and red line
shows the dielectric BC. The inset shows both profiles in semilog-scale.
Figure 3:  Single wall fixed at surface potential $V=1$. Bulk ion
concentrations are $c_{0,+}=c_{0,-}=0.1$ M. To amplify the boundary layer, we
zoom into the region $0<x<1$ nm. (a) Electric potential profiles for the metal
(blue solid), dielectric (red solid) and PB (black dashed) BCs. (b) Cation
profiles with metal BCs (blue, left-y axis), dielectric (red, right y-0axis)
BCs, and PB theory (black, right y-axis). (c) Anion profiles: legends and axes
are the same as in (b).  Figure 4:  Symmetrically charged slit pore of 1 nm
with both plates at fixed surface potentials $V_{L}=V_{R}=1$. Bulk ion
concentrations are $c_{0,+}=c_{0,-}=0.1$ M. (a) Electric potential profiles
for the metal BCs (blue solid), dielectric (red solid) BCs, and PB theory
(black dashed). (b) Cation profiles with metal BCs (blue, left y-axis),
dielectric BCs (red, right y-axis) and PB theory (black, right y-axis). (c)
Anion profiles: legends and axes are the same as in (b).  Figure 5:  Ion
concentration versus distance from a single metal wall at fixed surface
potential $V=0$ with a dielectric surface layer. The dielectric layer has size
$a=0.1$ nm and dielectric constant $\epsilon_{\text{L}}=60$ (blue) or
$\epsilon_{\text{L}}=2$ (red). The bulk ion concentrations are
$c_{0,+}=c_{0,-}=0.1$ M. For comparison, the ion concentration near a
dielectric wall ($\epsilon_{\text{p}}=2$) with no layer is shown in black. The
implicit solvent has dielectric constant $\epsilon_{\text{r}}=80.$ Figure 6:
Ion concentration versus distance for a symmetrically charged slit pore of
$L=1$ nm. Both walls are metal with a dielectric layer of size $a=0.1$ nm on
each metal surface. The dielectric layer has dielectric constant
$\epsilon_{\text{L}}=60$ (blue) or $\epsilon_{\text{L}}=2$ (red). The bulk ion
concentrations are $c_{0,+}=c_{0,-}=0.1$ M. The implicit solvent has
dielectric constant $\epsilon_{\text{r}}=80$. (a) Applied potential $V=0$. (b)
Positive ion densities for applied potential $V=1$. (c) Negative ion densities
for applied potential $V=1$.

###  III.1 Charge accumulation/depletion at the surface

For a single point charge near a planar surface, its image charge is known to
be attractive for a metal boundary and repulsive for a dielectric medium with
a lower dielectric constant than that of the electrolyte. As the point charge
approaches the boundary interface, its distance to the image charge decreases
to 0, resulting in a divergent image charge interaction. This is reflected in
the self-energy diverging at the domain boundaries as shown in Fig. 1(a). The
wall surface for both curves in Fig. 1 are fixed at $V=0$, and bulk ion
concentrations are $c_{+,0}=c_{-,0}=0.1$ M. The self-energy near the surface
has two contributions: (1) there are no other salt ions inside the plate,
hence less correlation energy. (2) the dielectric property of the plate
differs from the solution, leading to IC interactions. For dielectric BC, salt
ions are completely depleted from the surface; for metal BC, the ion density
saturates at the surface as shown in Fig. 1(b). The cation and anion profiles
overlap for either metal or dielectric BC since there is no applied voltage or
fixed charges. Although in both cases the electric potential is constant 0
everywhere, the concentration profiles are nontrivial and nonuniform,
different from the PB solution that gives trivial constant concentrations
equal the bulk value of 0.1 M. Away from the boundary-layer zone, both
dielectric and metal BCs show exponentially decaying deviations of ion
densities from the bulk with a similar screening length, while the metal curve
shows a significantly different screening length for the boundary layer zone,
as shown in the inset of Fig. 1(b).

For two parallel plates (a slit pore) with $V=0$ on both, the double layers
from the two surfaces start to overlap as the separation between the plates
decreases. Hence, the middle-point self-energy and ion concentrations deviate
from the bulk values as shown in Fig. 2. For metal BC (blue line), the image
charge attraction is not fully screened at the middle point so that
$c(x=0)/c_{0}\approx 1.27>1$. This is more clearly shown in the inset where
the blue line does not reach $10^{0}$ and is not flat in the middle. For
dielectric BC (red line), the image charge repulsion leads to ion depletion in
the entire slit, as shown by the right y-axis of Fig. 2(b), and the middle
point concentration is significantly reduced from the bulk value.

We now examine the double layer structure when a weak but non-zero electric
potential is maintained at the boundary surface(s). Fig. 3 shows the ion
concentrations and electric potential profiles near a single wall fixed at
$V=1$. The bulk ion concentrations are $c_{+,0}=c_{-,0}=0.1$ M. The electric
potential profiles for all three boundary conditions are deceptively similar
everywhere, as shown by the solid lines corresponding to the left y-axis of
Fig. 3(a). However, the charge densities are different at a very short
distance from the the wall due to the drastic differences in self-energy, as
shown by the dashed lines corresponding to the right y-axis of Fig. 3(a).
Beyond the boundary layer, the three charge densities converge to the same
outer solution. The IC interaction also leads to non-monotonic concentrations
of anions for dielectric BC near the wall, as shown by the red line in Fig.
3(c).

For a slit pore, as its width decreases with the surface potential of both
plates fixed at $V=1$, the double layers from the two surfaces start to
overlap. For this case, one can see the difference between the electric
potentials of the different BCs as shown in Fig. 4(a). Overall, the
concentration inside the slit pore is largest for metal BC and weakest for
dielectric BC. This leads to a different magnitude of the screening effect,
which mostly comes from the boundary layers. Going from the surface into the
solution, for metal BC the electrostatic potential drops faster than the PB
theory, whereas for dielectric BC, it drops slower than the PB theory.
Compared to the single-wall case (Fig. 3(b) and (c)), the overlapping double
layers and stronger IC effects in the 1 nm pore lead to two main differences:
(1) the boundary layers for both the metal and dielectric BCs become thicker;
(2) the anion concentration for dielectric BC increases monotonically from the
boundary until the middle-point as shown in Fig. 4(b) and (c). For metal BC,
the stronger IC attraction leads to a peak of cation concentration close to
the wall, as shown in Fig. 3(b). Moreover, both cations and anions enrich near
the surface for metal BC despite the positive applied potential.

Now, we turn to the case of a dielectric layer with dielectric constant
$\epsilon_{\text{L}}$ on the metal surface that is lower than the dielectric
constant of the solution. The dielectric layer generates a repulsive image
charge that competes with the attractive image charge of the metal. This
repulsive image charge is stronger with greater disparity between the layer
and solution dielectric constant, and this effect is shown for a single planar
boundary in Fig. 5. We see that the ion density completely depletes from the
surface, meaning the repulsive IC from the layer overpowers the metal IC at
all distances. For a layer with a low dielectric constant
$\epsilon_{\text{L}}=2$, the repulsive IC is so strong that the ion profile
nearly matches that of a dielectric wall with the same dielectric constant,
shown by the dashed line in Fig. 5.

Next, we examine the double-layer structure between two parallel metal plates
with a dielectric layer on each surface, shown in Fig. 6. For no applied
potential and $\epsilon_{\text{L}}=60$ in Fig. 6(a), the ions adsorb in the
pore due to the metal’s attractive IC. However, the ions deplete close to the
wall due to the layer’s repulsive IC, giving a peak in the ion density
profile. The layer prevents ion saturation at the metal wall previously seen
in Fig. 2, yielding density profiles that qualitatively resemble those of Ref.
38. If the layer dielectric constant is low ($\epsilon_{\text{L}}=2$), the
ions completely deplete in the pore, similar to the profiles of the dielectric
wall in Fig. 2, albeit with a higher maximum density due to the metal’s
attractive IC. For applied potential $V=1$ in Fig. 6(b) and (c) at
$\epsilon_{\text{L}}=60$, the positive and negative ion density profiles are
nonmonotonic for the same reasons as the zero-potential case.

###  III.2 Differential capacitance

We next examine the differential capacitance $C_{d}=d\sigma_{s}/dV$ for
nanometer-sized slit pores. The two plates of the slit pore are at the same
surface potential $\psi(z=\pm L/2)=V$ to simulate a nanoscale pore inside a
porous electrode that is maintained at a constant potential.

Figure 7:  Differential capacitance curves for four cases: metal (blue),
dielectric (red), and layer BCs (yellow), as well as PB theory (black). The
capacitance curves are shown under conditions: (a) $c_{0}=0.1$ M, $L=4$ nm.
(b) $c_{0}=5$ M, $L=4$ nm. (c) $c_{0}=0.1$ M, $L=1$ nm.  Figure 8:
Differential capacitance curves for four cases: metal (blue), dielectric
(red), and layer BCs (yellow), as well as PB theory (black) as a function of
plate separation for $V=0$. (a) Bulk ion concentration is $c_{0}=5$ M. (b)
Bulk ion concentration is $c_{0}=0.1$ M.

In all three system conditions explored in Fig. 7, the initial capacitance (at
very low voltage) of the metal BC is higher than the PB theory, dielectric BC,
or layer BC. This is consistent with the qualitative understanding that metal
BC image charge attraction significantly enhances charge separation, where
charge separation refers to the difference in anion and cation densities [38].
On the other hand, at weak electric potentials, the image charge repulsion of
dielectric and layer BCs strongly depletes salt ions and reduces charge
separation. As a result, the initial $C_{d}$ is low. Moreover, since ions
cannot closely approach the metal surface due to exclusion by the layer, the
metal’s surface charge is less sensitive to the applied potential. This
effect, combined with the layer’s repulsive image charge, gives a lower
capacitance than the other three cases for nearly all conditions. The only
exception is at low salt concentration, small separation, and low applied
potential in Fig. 7(c), where the layer capacitance slightly exceeds
dielectric BC and PB theory. The reason is that, for very small separations,
the attractive IC from the metal dominates the layer’s repulsive IC, causing
greater accumulation of ions in the pore compared to dielectric BC and PB
theory.

At low bulk salt concentration $c_{0}=0.1$ M, and large separation between the
plates ($L=4$ nm, compared to the double layer thickness $<1$ nm), the
differential capacitance displays the “camel-shape” for dielectric and layer
BCs, as well as PB theory [85], whereas that of the metal BC has a slight
“bird shape” [77], as shown in Fig. 7(a). For the PB theory, this is due to
the nonlinear dependence of charge density on the electric potential. For
dielectric BC, as surface potential increases, the attraction of the surface
for counter-ions overcomes the image charge repulsion in the middle of the
slit; therefore, the capacitance is enhanced. The same is true for the layer
BC; however, since the metal’s surface charge is less sensitive to the applied
voltage for reasons mentioned earlier, the capacitance never exceeds the PB
result except for small separations with low applied potential and low bulk
concentration.

It is known that the overestimate of surface ion density in the Gouy–Chapman
theory leads to a differential capacitance curve that is always “camel-
shaped”, with a minimum at 0 applied voltage. Kornyshev [85] pointed out that
the “lattice saturation” effect [86, 87], which accounts for finite ion sizes,
can explain the observed deviations from the Gouy–Chapman theory, i.e. a
“bell” shaped capacitance curve for room temperature ionic liquids [88]. Fig.
7(b) shows this scenario with a high bulk salt concentration $c_{0}=5$ M, and
large separation $L=4$ nm. The solvent concentration is only $c_{s}\approx
3.33$ M in the bulk reservoir, lower than the ion concentrations, as in a
‘water-in-salt’ electrolyte. At small voltage, initial charge separation is
enhanced due to the large number of salt ions in the slit, giving a large
$C_{d}$. Then, as the electric potential increases, the double layer saturates
quickly, regardless of the BC. Thus, the overall $C_{d}$ displays a bell shape
for all four cases. These results are consistent with previous studies [89,
90, 85, 91]. Moreover, all curves converge at high applied potential due to
double-layer saturation. We expect this convergence to occur for sufficiently
high potential under the conditions for Fig. 7(a) and (c) as well. Similar to
the modified PB theory, the ion-saturation effect plays a key role to render
the bell shaped curves.

The more interesting regime is at low bulk salt concentration $c_{0}=0.1$ M
and very small separation $L=1.0$ nm. In this case, the boundary accumulation
/ depletion layers are very strong and dictate the charging behaviors at low
surface potentials. Specifically, the IC attraction of metal BC strongly
enhances the charge separation, resulting in a high initial capacitance. But,
as the voltage increases, the saturated boundary layers screen out the
influence of the surface potential in the middle of the slit, as seen in Fig.
4. This leads to a small second peak of the blue solid line in Fig. 7(c),
giving a “bird-shaped” capacitance curve. On the other hand, for the
dielectric and layer BCs, the IC repulsion at small distance is too strong at
weak electric potential, so very few counter-ions migrate into the slit when
the voltage increases initially. Hence, the low-voltage capacitance is quite
low for the dielectric and layer BCs. Additionally, the linear response regime
is much smaller for metal BC under low bulk salt and small separation compared
to dielectric and layer BCs. The reason for this is that the charge
fluctuation in the metal is much higher under these conditions, causing
significant deviations from the linear response capacitance even at low
applied potentials.

Lastly, we investigate the capacitance as a function of plate separation at
$V=0$ for two bulk ion concentrations $c_{0}=5$ and $0.1$ M, shown in Fig. 8
(a) and (b), respectively. At both concentrations, we see that the metal
capacitance exceeds PB theory due to the attractive metal IC. Moreover, the
metal capacitance increases as the pore size decreases, with an anomalously
large increase for $c_{0}=0.1$ M compared to PB theory, similar to the results
shown in Ref. [42]. For a dielectric boundary, the capacitance is lower than
PB theory for both concentrations due to ion depletion in the pore caused by
the dielectric wall’s repulsive IC. For the layer BC in both cases, the
capacitance is initially high due to the metal’s attractive IC dominating at
strong confinement. However, as the plates move apart, the capacitance
decreases due to the layer’s repulsive IC overtaking the attractive IC from
the metal, causing ions to deplete in the pore. This effect creates a minimum
in the capacitance. For large separations ($>1.5$) nm, the layer’s repulsive
IC has a weaker effect since it is a boundary layer effect, so the capacitance
gradually increases.

Figure 9:  Force between neutral and like-charged plates with different
material boundaries. In all plots the bulk concentrations are colored by:
$c_{0}=0.1$ M (blue, red and black) and $c_{0}=0.2$ M (green). (a) Two
dielectric plates with same surface charge densities. Surface charge densities
are $\sigma=0.01,0.02,0.05$ /nm2 for the blue, red, and black line,
respectively, and $\sigma=0.02$ /nm2 for the green line. (b) Two dielectric
plates with the same fixed surface potentials. The surface potentials are
$V=0,1,2$ for the blue, red, and black line, respectively, and $V=0$ for the
green line. (c) Two metal plates with the same fixed surface potentials. The
surface potentials are $V=0,1,2$ for the blue, red, and black line,
respectively, and $V=0$ for the green line.  Figure 10:  Force between neutral
and like-charged metal plates with a $0.1$ nm thick dielectric layer with
dielectric constant $\epsilon_{\text{L}}=60$ on each plate. The metal plates
have fixed surface potentials of either $V=0,1,2$. The blue, red, and black
curves have bulk ion concentration $c_{0}=0.1$ M, and the green curve has bulk
ion concentration $c_{0}=0.2$ M.

###  III.3 Force between like-charged plates

Finally, we investigate the forces between two symmetrically charged plates.
In reality, it is more natural to control surface charge density for
dielectric plates, or to maintain constant surface electric potential for
metal plates, especially if the surface geometry is not simple. For ease of
conceptual comparison, we also present results for dielectric BC with fixed
surface potentials. In Fig. 9 (a), we investigate three surface charges
$\sigma=0.01,0.02,0.05$ e/nm2 at bulk concentration $c_{0}=0.1$ M (blue,
black, and red curves, respectively). For lower surface charges
$\sigma=0.01,0.02$ e/nm2, the ion depletion due to IC of dielectric BC results
in an attractive potential well. There is also a small but noticeable
repulsive regime for plate separations between 1$\sim$2 nm, which is
consistent with the findings of Ref. [23], although there it was computed
using a WKB approximation. At higher surface charge $\sigma=0.05$ e/nm2, the
same-charge repulsion overcomes the IC effect. Moreover, for $c_{0}=0.1$ M,
the attractive well gradually disappears and completely vanishes when
increasing the surface charge from $0.01$ to $0.05$ e/nm2 in Fig. 9 (a).
Moving to the effect of bulk concentration, the green curve is for a fixed
surface charge $\sigma=0.02$ e/nm2 with bulk concentration $c_{0}=0.2$ M.
Comparing the green curve to the red curve, it is obvious that the higher bulk
concentration leads to a stronger depletion force at the same surface charge
density. This is also qualitatively consistent with Ref. [23]. However, when
the boundary condition is changed to fixed surface potentials, the purely
repulsive regime is no longer observed for reasonably weak applied potentials,
as shown in Fig. 9(b). The repulsion at short distance in Fig. 9(a) comes from
the strong repulsion between the fixed surface charges, while in Fig. 9(b) the
surface charges adjust to avoid strong repulsion and to lower the total free
energy. At fixed surface potential, both higher surface potential (black line)
and higher bulk salt concentration (green line) amplify the variation in the
force, as the IC-induced depletion effect is enhanced.

In the case of metal BC, the IC attraction is very strong and dominates the
interplate force. This is shown by Fig. 9(c), where the magnitude of the
attraction is much larger than in Fig. 9(b). In addition, no repulsive regime
is observed and the attractive force decays over a longer lengthscale. Recall
that we use point-charge model throughout the calculations here for
consistency and fair comparisons. This means that ions can approach the
boundaries infinitely close. For dielectric BC, the IC repulsion self-
regulates this effect by suppressing the ion distributions at the surface.
However, for metal BC, the IC attraction results in high concentration of
counter-ions at the surface close to the saturation concentration, with a very
large self-energy. As a result the attraction force computed here at short
separations between the two electrodes may be overestimated compared with
experimental or simulation results with finite ion sizes. For this reason, we
introduce the thin dielectric layer to the metal surfaces to account for the
distance of closest approach for ions, which we discuss next.

To finish, we investigate the effect of a $0.1$ nm thick dielectric layer
(dielectric constant $\epsilon_{\text{L}}=60$) on metal plates, where the
metal plates are at the same fixed potential. In Fig. 10, we see that the
layer’s repulsive IC cancels the metal’s attractive IC for separations $>1$
nm, preventing the plates from interacting until close separation ($<1$) nm
for all conditions considered. Consequently, the lengthscale of interaction
between plates with layer BC is much smaller than that of metal and dielectric
BCs. Moreover, at separations of $1$ nm and above, the ions are depleted from
the walls. As the plates move from 1 nm to 0.9 nm spacing, the metal’s
attractive IC overcomes the repulsive IC from the layer, causing ions to flood
into the pore that generates an attractive force between the plates.
Furthermore, the peak in the $V=2$ curve stems from a stronger flooding of
ions into the pore compared to the lower potentials. The slight bump in the
force curves near 0.85 nm separation is due to the ions saturating the pore
walls.

##  IV Conclusions

In this work, we discuss the image charge effect and its consequences in nano-
capacitors and forces between like-charged plates. We treat the IC effects
non-perturbatively for a point-charge model, based on a renormalized Gaussian
fluctuation theory. We apply the incompressibility condition to account for
the excluded volume interactions between the ions and solvent and to prevent
the ion density from diverging at the interfaces. We obtain full numerical
solutions to the renormalized field theory subject to different boundary
conditions. We focus on a comparison of these calculations under the (1)
perfectly conducting metal plates, (2) dielectric plates, (3) metal plates
with a thin dielectric layer, and (4) the PB theory. We demonstrate these
comparisons using a simple 1:1 electrolyte.

We first show that the double-layer structures are significantly affected by
the IC interaction. For a single surface without applied voltage or
externally-imposed surface charge, IC leads to saturation of ions at the
surface for the metal BC, and strong ion depletion for the dielectric and
layer BCs. For the layer BC, the depletion strengthens with lower dielectric
constant relative to the solvent. For a narrow slit pore, the middle-point
concentrations significantly deviate from the bulk values due to IC. When a
single surface is fixed at a weak applied potential, the electrostatic
potential profiles are not sensitive to the BCs. However, ion densities near
the surface depends strongly on the type of BC applied. For a narrow slit pore
with both surfaces at a weak fixed potential, both the ion density and the
electrostatic potential profiles are strongly affected by the BCs. Moreover,
the layer BC gives nonmonotonic ion density profiles that are qualitatively
similar to those observed in MD simulations [38].

Based on the understanding about the double-layer structures, our calculations
further show that the boundary layers due to IC interaction lead to
significant differences in the differential capacitance of a nano-pore under
different boundary conditions. The dielectric and layer BCs show a
qualitatively similar trend with the PB theory, where the capacitance curve
changes from a camel shape to a bell shape as bulk salt concentration
increases. Moreover, the metal BC gives a “bird-shaped” capacitance at low
bulk ion concentrations and small pore width due to IC-induced boundary layer
accumulation and ion saturation. The layer BC shows that the ion exclusion
near the metal surface combined with the repulsive IC yields a lower
capacitance than PB, metal, or dielectric BCs at nearly all conditions. Notice
that modified PB treatment [42] only accounts for the IC from the mean-charge
density $\rho=c_{+}-c_{-}$ in the Poisson equation, while our treatment
captures the IC effect even at $\rho=0$, and the ion density profiles $c_{+}$
and $c_{-}$ are non-uniform. As a result, our predictions of the capacitance
are different from the modified PB treatment, especially at weak applied
potential.

The differences between the BCs are also manifested in the forces between
like-charged planar walls. For dielectric BC, at fixed surface charge density,
the attraction well gradually diminishes as surface charge increases; at fixed
surface potential, the attraction well is stronger at higher potential values.
For metal BC with fixed surface potentials on both plates, which is more
directly relevant to nano-pores inside porous electrodes, the attraction is
much stronger than the dielectric BC and decays slower; no repulsion is
observed in this ensemble in the weak-coupling regime. Increasing either
surface potential or bulk concentration slightly enhances the attraction
force. The prediction of pure attraction for planar metal plates is different
from the nonlinear PB theory, which only gives repulsion between like-charged
metal plates—Ref. 71 showed that in the PB theory, one has to invoke curvature
effects to induce like-charge attraction with metal BC. For layer BC at fixed
surface potentials, the interaction lengthscale is much smaller than that of
metal or dielectric BCs due to the layer’s repulsive IC canceling the metal’s
attractive IC at separations $>1$ nm. For separations $<1$ nm, the plates are
strongly attractive, similar to the metal BC, due to ions flooding the pore at
close separation when the metal IC overcomes the layer’s repulsive IC.

Our work presents a fundamental perspective on the IC interaction,
highlighting its non-negligible influence on double layer structures, and
consequently, nanoscale capacitors and forces between like-charged plates.
These effects are often omitted or treated in a mean-field fashion. However,
we show that even at zero mean-charge density, where mean-field predictions
will neglect IC, it still has a strong effect. Our results provide a starting
point to include IC in various realistic applications.

###### Acknowledgements.

T.Z. acknowledges the support of the Cecil and Sally Drinkward Postdoc
Fellowship. He also thanks Z. Peng for helpful discussions on numerical
methods. Z.-G.W. acknowledges financial support from the Hong Kong Quantum AI
Lab Ltd.

## Data Availability Statement

All data that support the findings of this study are available within the
article.

##  Appendix A Brief summary of the key equations

Most of the derivations were already presented in Ref. 22 and earlier
references cited there. Here we sketch the outline of the derivations,
including the incompressibility constraint, for the reader’s convenience. We
start with the Hamiltonian in SI unit

$H=\int
d\mathbf{r}d\mathbf{r}^{\prime}\frac{e^{2}}{2}\rho(\mathbf{r})G_{0}(\mathbf{r},\mathbf{r}^{\prime})\rho(\mathbf{r}^{\prime})$
(31)

where the Coulomb operator is defined by

$-\nabla\cdot\varepsilon\nabla^{\prime}G_{0}(\mathbf{r},\mathbf{r}^{\prime})=\delta(\mathbf{r}-\mathbf{r}^{\prime})$
(32)

The total charge density can be decomposed into external charge, cations and
anions

$\rho(\mathbf{r})=\rho_{ex}(\mathbf{r})+z_{+}\sum_{i+}h_{+}(\mathbf{r}-\mathbf{r}_{+}^{i})-z_{-}\sum_{j-}h_{-}(\mathbf{r}-\mathbf{r}_{-}^{j})$
(33)

where $h_{\pm}$ represents the spread shape of ion charge. The
incompressibility of the liquid requires

$\hat{c}_{tot}=\hat{c}_{+}(\mathbf{r})+\hat{c}_{-}(\mathbf{r})+\hat{c}_{s}(\mathbf{r})=1/v_{0}$
(34)

where the density of species

$\begin{split}\hat{c}_{+}(\mathbf{r})=&\sum_{i+}h_{+}(\mathbf{r}-\mathbf{r}_{+}^{i})\\\
\hat{c}_{-}(\mathbf{r})=&\sum_{i-}h_{-}(\mathbf{r}-\mathbf{r}_{-}^{i})\\\
\hat{c}_{s}(\mathbf{r})=&\sum_{s-}h_{s}(\mathbf{r}-\mathbf{r}_{s}^{i})\end{split}$
(35)

subscript $s$ stands for solvent.

The canonical partition function is

$\begin{split}Q=&\frac{1}{n_{+}!n_{-}!n_{s}!v_{+}^{n_{+}}v_{-}^{n_{-}}v_{s}^{n_{s}}}\int
d\mathbf{r}_{i+}d\mathbf{r}_{j-}d\mathbf{r}_{s}e^{-\beta
H}\prod_{r}\delta\left(\hat{c}_{tot}-1/v_{0}\right)\\\
=&\frac{1}{n_{+}!n_{-}!n_{s}!v_{+}^{n_{+}}v_{-}^{n_{-}}v_{s}^{n_{s}}}\int
d\mathbf{r}_{i+}d\mathbf{r}_{j-}d\mathbf{r}_{s}\\\ &\times
e^{-\beta\frac{e^{2}}{2}\int
d\mathbf{r}d\mathbf{r}^{\prime}(\rho_{ex}+\rho_{ion})(\mathbf{r})G_{0}(\mathbf{r},\mathbf{r}^{\prime})(\rho_{ex}+\rho_{ion})(\mathbf{r}^{\prime})}\\\
&\times\int dc_{\pm,s}\int D[\eta]e^{i\int
d\mathbf{r}\eta(\mathbf{r})({c}_{+}+{c}_{-}+{c}_{s}-1/v_{0})(\mathbf{r})}\\\
&\times\int dw_{\pm,s}e^{i\int
d\mathbf{r}\left[w_{+}(\hat{c}_{+}-c_{+})+w_{-}(\hat{c}_{-}-c_{-})+w_{s}(\hat{c}_{s}-c_{s})\right]}\end{split}$
(36)

where auxiliary fields $w_{\pm,s}$ are also introduced. In general, if there
are other quadratic or higher-order interaction terms involving the
concentrations, it will be useful to deal with the $w$ fields instead of
directly manipulating the concentration operators. Using the Hubbard-
Stratonovich transform

$\int dx\exp\left(-\frac{1}{2}x^{T}A^{-1}x\pm
ixy\right)=\sqrt{\det(A)}\exp\left(-\frac{1}{2}y^{T}Ay\right)$ (37)

the exponential factor is

$\begin{split}&e^{-\beta\frac{e^{2}}{2}\int
d\mathbf{r}d\mathbf{r}^{\prime}\rho(\mathbf{r})G_{0}(\mathbf{r},\mathbf{r}^{\prime})\rho(\mathbf{r}^{\prime})}\\\
&=\frac{1}{\sqrt{\det(G_{0})}}\int D[\phi]e^{-\int
d\mathbf{r}d\mathbf{r}^{\prime}\left(\frac{1}{8\pi\ell_{B}}\nabla\phi(\mathbf{r})\cdot\nabla^{\prime}\phi(\mathbf{r}^{\prime})+i\rho(\mathbf{r})\phi(\mathbf{r}^{\prime})\right)\delta(\mathbf{r}-\mathbf{r}^{\prime})}\\\
\end{split}$ (38)

here $\ell_{B}=e^{2}/(4\pi\epsilon k_{B}T)$ is the Bjerrum length and $\phi$
is the dimensionless potential scaled by $\frac{k_{B}T}{e}$.

Now the grand partition function is

$\begin{split}\Xi=&\sum_{n_{+}}\sum_{n_{-}}\sum_{n_{s}}\frac{Qe^{n_{+}\mu_{+}}e^{n_{-}\mu_{-}}e^{n_{s}\mu_{s}}}{\sqrt{\det(G_{0})}}\int
D[\phi]D[\eta]dw_{\pm,s}dc_{\pm,s}e^{-\mathcal{L}}\end{split}$ (39)

where the action is

$\begin{split}\mathcal{L}=&\int
d\mathbf{r}d\mathbf{r}^{\prime}\left[\frac{1}{8\pi\ell_{B}}\nabla\phi(\mathbf{r})\cdot\nabla^{\prime}\phi(\mathbf{r}^{\prime})+i\rho_{ex}(\mathbf{r})\phi\mathbf{r}^{\prime})\right]\delta(\mathbf{r}-\mathbf{r}^{\prime})\\\
&-\int
d\mathbf{r}\left(\lambda_{+}e^{-i\left(h_{+}*[z_{+}\phi+iw_{+}]\right)(\mathbf{r})}+\lambda_{-}e^{i\left(h_{-}*[z_{-}\phi-
iw_{-}]\right)(\mathbf{r})}\right)\\\ &+\int
d\mathbf{r}\left(w_{+}c_{+}+w_{-}c_{-}+w_{s}c_{s}-\eta(c_{+}+c_{-}+c_{s}-1/v_{0})\right)\\\
&-\int d\mathbf{r}\lambda_{s}e^{h_{s}*w_{s}}\end{split}$ (40)

$*$ stands for the convolution operator

$(h*\phi)(\mathbf{r})=\int
d\mathbf{r}^{\prime}h(\mathbf{r}^{\prime}-\mathbf{r})\phi(\mathbf{r}^{\prime})$
(41)

and $\lambda_{\pm}=\frac{e^{\mu_{\pm}}}{v_{\pm}}$ is the fugacity.

The variational approach decomposes the field $\phi=-i\psi+\chi$, and the
action becomes

$\begin{split}\mathcal{L}[\phi]&=\int
d\mathbf{r}\left\\{\frac{1}{8\pi\ell_{B}}\left(\nabla\left(-i\psi(\mathbf{r})+\chi(\mathbf{r})\right)\right)^{2}+\rho_{ex}(\mathbf{r})(-i\psi(\mathbf{r})+\chi(\mathbf{r}))\right\\}\\\
&-\int
d\mathbf{r}\lambda_{+}e^{-z_{+}(h_{+}*\psi)(\mathbf{r})+(h_{+}*w_{+})(\mathbf{r})-iz_{+}(h_{+}*\chi)(\mathbf{r})}\\\
&-\int
d\mathbf{r}\lambda_{-}e^{z_{-}(h_{-}*\psi)(\mathbf{r})+(h_{-}*w_{-})(\mathbf{r})+iz_{-}(h_{-}*\chi)(\mathbf{r})}\\\
&+\int
d\mathbf{r}\left(w_{+}c_{+}+w_{-}c_{-}+w_{s}c_{s}-\eta(c_{+}+c_{-}+c_{s}-1/v_{0})\right)\\\
&-\int d\mathbf{r}\lambda_{s}e^{h_{s}*w_{s}}\end{split}$ (42)

The average number (density) of ions is

$\begin{split}\left<n_{\pm}\right>=\frac{\partial\ln\Xi}{\partial\mu_{\pm}}=&\int
d\mathbf{r}\lambda_{\pm}\left<e^{\mp
z_{\pm}(h_{\pm}*(\psi+i\chi))(\mathbf{r})+(h_{\pm}*w_{\pm})(\mathbf{r})}\right>\\\
c_{\pm}(\mathbf{r})=&\lambda_{\pm}\left<e^{\mp
z_{\pm}(h_{\pm}*(\psi+i\chi))(\mathbf{r})+(h_{\pm}*w_{\pm})(\mathbf{r})}\right>\end{split}$
(43)

For the variatioal renormalization, a Gaussian reference action is chosen

$\mathcal{L}_{ref}=\int
d\mathbf{r}d\mathbf{r}^{\prime}\frac{1}{2}\chi(\mathbf{r})G^{-1}(\mathbf{r},\mathbf{r}^{\prime})\chi(\mathbf{r}^{\prime})$
(44)

and the extremized free energy is

$\begin{split}W=&W_{ref}+\left<\mathcal{L}[\phi]-\mathcal{L}_{ref}[\phi]\right>_{ref}\\\
=&-\frac{1}{2}\ln\left(\frac{\det(G)}{\det(G_{0})}\right)+\left<L[\phi]\right>_{ref}-\frac{1}{2}\end{split}$
(45)

To evaluate $\left<\mathcal{L}[\phi]\right>_{ref}$, we note that

$\begin{split}&\left<\int
d\mathbf{r}d\mathbf{r}^{\prime}\frac{1}{4\pi\ell_{B}(\textbf{r})}\delta(\mathbf{r}-\mathbf{r}^{\prime})\nabla\chi(\mathbf{r})\cdot\nabla^{\prime}\chi(\mathbf{r}^{\prime})\right>_{ref}\\\
=&\int
d\mathbf{r}d\mathbf{r}^{\prime}\left\\{\left<\chi(\mathbf{r})\chi(\mathbf{r}^{\prime})\right>_{ref}\nabla\left[\frac{1}{4\pi\ell_{B}(\textbf{r})}\right]\cdot\nabla^{\prime}\delta(\mathbf{r}-\mathbf{r}^{\prime})\right\\}\\\
=&\int
d\mathbf{r}d\mathbf{r}^{\prime}\left\\{G(\mathbf{r},\mathbf{r}^{\prime})\nabla\left[\frac{1}{4\pi\ell_{B}(\textbf{r})}\right]\cdot\nabla^{\prime}\delta(\mathbf{r}-\mathbf{r}^{\prime})\right\\}\\\
\end{split}$ (46)

and

$\begin{split}\left<e^{\mp
iz_{\pm}(h_{\pm}*\chi)(\mathbf{r})}\right>_{ref}=&e^{-\frac{1}{2}z_{\pm}^{2}\int
d\mathbf{y}d\mathbf{y}^{\prime}h_{\pm}(\mathbf{y}-\mathbf{r})G(\mathbf{y},\mathbf{y}^{\prime})h_{\pm}(\mathbf{y}^{\prime}-\mathbf{r})}\end{split}$
(47)

Now we can simplify the averaged action as

$\left<\mathcal{L}[\psi]\right>_{ref}=\int
dr\left\\{-\frac{1}{8\pi\ell_{B}}(\nabla\psi(\mathbf{r}))^{2}+\rho_{ex}(\mathbf{r})\psi(\mathbf{r})-c_{\pm}(\mathbf{r})\right\\}\\\
+\frac{1}{2}\int
drdr^{\prime}G(\mathbf{r},\mathbf{r}^{\prime})\nabla\cdot\left[\frac{1}{4\pi\ell_{B}(\textbf{r})}\right]\nabla^{\prime}\delta(\mathbf{r}-\mathbf{r}^{\prime})\\\
+\int
d\mathbf{r}\left(w_{+}c_{+}+w_{-}c_{-}+w_{s}c_{s}-\eta(c_{+}+c_{-}+c_{s}-1/v_{0})\right)\\\
-\int d\mathbf{r}\lambda_{s}e^{h_{s}*w_{s}}$ (48)

The Euler-Lagrange equations in the main text for this field theory then
result from variational derivatives of $W$ with respect to $\psi$, $G$,
$w_{\pm,s}$, $c_{\pm,s}$ and $\eta$.

##  Appendix B Image charge of a single point charge near layered interfaces

###  B.1 Metal-dielectric interface

Suppose for half space $z<0$, it’s a metal plate. For half space $z>a$, it’s a
dielectric medium with relative permitivity $\epsilon_{2}$. In between $0<z<a$
it’s another dielectric with $\epsilon_{1}$. A point charge $q$ is placed at
$(0,0,a+b)$. We now solve for the potential inside the two dielectric media
$\phi_{1}$ and $\phi_{2}+\phi_{3}$. Due to cylindrical symmetry, the Poisson
equation can be written as

$\begin{split}&\left(\partial_{z}^{2}+\frac{1}{r}\partial_{r}r\partial_{r}+\frac{1}{r^{2}}\partial_{\theta}^{2}\right)\phi_{i}=0\qquad
i=1,2\\\
&\phi_{3}=\frac{eq}{4\pi\epsilon_{0}\epsilon_{2}\sqrt{r^{2}+[z-(a+b)]^{2}}}\end{split}$
(49)

Nondimensionalize by the Bjerrum length in the liquid

$\ell_{B}=\frac{e^{2}}{4\pi\epsilon_{0}\epsilon_{2}k_{B}T}$ (50)
$\psi=\frac{e\phi}{k_{B}T}=\frac{q\ell_{B}}{\sqrt{r^{2}+[z-(a+b)]^{2}}}$ (51)

For water $\epsilon_{2}=80$, $\ell_{B}\sim 0.7$ nm.

On the metal-dielectric interface $z=a$, the boundary conditions are

$\begin{split}\epsilon_{1}\partial_{z}\phi_{1}(r,\theta,a)&=\epsilon_{2}\partial_{z}\left(\phi_{2}(r,\theta,a)+\phi_{3}(r,\theta,a)\right)\\\
\phi_{1}(r,\theta,a)&=\phi_{2}(r,\theta,a)+\phi_{3}(r,\theta,a)\\\
\end{split}$ (52)

the second condition implies already the continuity of the tangential $E$
field

$\partial_{r}\phi_{1}(r,\theta,a)=\partial_{r}\left(\phi_{2}(r,\theta,a)+\phi_{3}(r,\theta,a)\right)$
(53)

At $z=0$,

$\phi_{1}(r,\theta,0)=0$ (54)

At $z\rightarrow\infty$,

$\phi_{2}(r,\theta,\infty)=0$ (55)

The general solution satifying BCs Eqn. 54 and Eqn. 55

$\begin{split}\phi_{1}&=\int_{0}^{\infty}d\lambda\sinh(\lambda z)J_{0}(\lambda
r)f_{1}(\lambda)\\\ \phi_{2}&=\int_{0}^{\infty}d\lambda e^{-\lambda
z}J_{0}(\lambda r)f_{2}(\lambda)\end{split}$ (56)

Then one use the fact

$\frac{1}{\sqrt{r^{2}+(z-a-b)^{2}}}=\int_{0}^{\infty}d\lambda J_{0}(\lambda
r)e^{-\lambda|z-a-b|}$ (57)

To match the BC at $z=a$

$\begin{split}q\epsilon_{2}e^{-\lambda
b}&=f_{1}(\lambda)\epsilon_{1}\cosh(a\lambda)+f_{2}(\lambda)\epsilon_{2}e^{-\lambda
a}\\\ qe^{-\lambda b}&=f_{1}(\lambda)\sinh(a\lambda)-f_{2}(\lambda)e^{-\lambda
a}\end{split}$ (58)

As a result,

$\begin{split}f_{1}(\lambda)&=\frac{2\epsilon_{2}qe^{-\lambda
b}}{\epsilon_{1}\cosh(\lambda a)+\epsilon_{2}\sinh(\lambda a)}\\\
f_{2}(\lambda)&=q\frac{\epsilon_{2}\sinh(\lambda a)-\epsilon_{1}\cosh(\lambda
a)}{\epsilon_{1}\cosh(\lambda a)+\epsilon_{2}\sinh(\lambda
a)}e^{\lambda(a-b)}\\\ \end{split}$ (59)

For the IC force, at $z=a_{+0}$, the excess E field is

$\begin{split}E_{r}&=-\int_{0}^{\infty}d\lambda e^{-\lambda a}\lambda
J_{0}^{\prime}(\lambda r)f_{2}(\lambda)\\\ &=\int_{0}^{\infty}d\lambda
e^{-\lambda a}\lambda J_{1}(\lambda r)f_{2}(\lambda)\\\
E_{z}&=\int_{0}^{\infty}d\lambda e^{-\lambda a}\lambda J_{0}(\lambda
r)f_{2}(\lambda)\end{split}$ (60)

in addition to the E field generated by the point charge potential $\phi_{3}$.

The surface charge density $\sigma(z=a)$ on the interface at $z=a$ can be
extracted by the difference of the E field

$\begin{split}E_{z}(z=a+0)-E_{z}(z=a-0)=\sigma(z=a)\left(\frac{1}{\epsilon_{1}}+\frac{1}{\epsilon_{2}}\right)\end{split}$
(61)

Similarly, for the surface charge density $\sigma(z=0)$, we have

$\sigma(z=0)=\epsilon_{1}E_{z}(z=0+0)$ (62)

##  Appendix C  IC with a dielectric surface layer

Suppose $z>0$ half space is filled by the electrolyte solution with dielectric
constant $\epsilon_{r}$, while $-a<z<0$ is a dielectric layer with
$\epsilon_{\text{L}}$.

The BC for Poisson equation with fixed surface potentials is a Robin type

$\epsilon_{\text{L}}\phi(0)-\epsilon_{r}a\phi^{\prime}(0)=\epsilon_{\text{L}}V_{1}$
(63)

For another layer at $L<z<L+a$, we have similarly

$\epsilon_{\text{L}}\phi(L)+\epsilon_{r}a\phi^{\prime}(L)=\epsilon_{\text{L}}V_{2}$
(64)

For fixed surface charge $e\sigma_{S}$ at $z=-a$, the BC at $z=0$ is still a
Neumann BC

$\begin{split}e\sigma_{S}=-\epsilon_{\text{L}}\phi_{\text{L}}^{\prime}(0)=-\epsilon_{r}\phi_{r}^{\prime}(0)\end{split}$
(65)

For the self-energy Green’s function, a Robin type BC now applies to both
metal and dielectric plates at $z<-a$. For a dielectric plate with
$\epsilon_{p}$

$G^{\prime}(0)-k\frac{\epsilon_{\text{L}}}{\epsilon_{r}}\frac{\epsilon_{\text{L}}\sinh(ka)+\epsilon_{p}\cosh(ka)}{\epsilon_{p}\sinh(ka)+\epsilon_{\text{L}}\cosh(ka)}G(0)=0$
(66)

Taking the limit of $\epsilon_{p}\rightarrow\infty$, we have the BC for the
metal plate

$G^{\prime}(0)-k\frac{\epsilon_{\text{L}}}{\epsilon_{r}}\coth(ka)G(0)=0$ (67)

## References

  * Wagner [1924] C. Wagner, Phys. Z 25, 474 (1924).
  * Onsager and Samaras [1934] L. Onsager and N. N. Samaras, The Journal of Chemical Physics 2, 528 (1934).
  * Naji  _et al._ [2005] A. Naji, S. Jungblut, A. G. Moreira, and R. R. Netz, Physica A: Statistical Mechanics and its Applications 352, 131 (2005).
  * Bakhshandeh  _et al._ [2011] A. Bakhshandeh, A. P. Dos Santos, and Y. Levin, Physical review letters 107, 107801 (2011).
  * Kjellander and Marcělja [1984] R. Kjellander and S. Marcělja, Chemical physics letters 112, 49 (1984).
  * Lemay  _et al._ [2013] S. G. Lemay, S. Kang, K. Mathwig, and P. S. Singh, Accounts of chemical research 46, 369 (2013).
  * Babel  _et al._ [2018] S. Babel, M. Eikerling, and H. Löwen, The Journal of Physical Chemistry C 122, 21724 (2018).
  * Han  _et al._ [2019] J.-W. Han, M.-L. Seol, D.-I. Moon, G. Hunter, and M. Meyyappan, Nature Electronics 2, 405 (2019).
  * Li  _et al._ [2010] T. Li, W. Hu, and D. Zhu, Advanced Materials 22, 286 (2010).
  * Biesheuvel  _et al._ [2014] P. Biesheuvel, S. Porada, M. Levi, and M. Z. Bazant, Journal of solid state electrochemistry 18, 1365 (2014).
  * Comtet  _et al._ [2017] J. Comtet, A. Nigués, V. Kaiser, B. Coasne, L. Bocquet, and A. Siria, Nature materials 16, 634 (2017).
  * Chao and Wang [2020] H. Chao and Z.-G. Wang, The journal of physical chemistry letters 11, 1767 (2020).
  * Zhou  _et al._ [2020] T. Zhou, M. Mirzadeh, R. J.-M. Pellenq, and M. Z. Bazant, Physical Review Fluids 5, 124201 (2020).
  * Steinpilz  _et al._ [2020] T. Steinpilz, K. Joeris, F. Jungmann, D. Wolf, L. Brendel, J. Teiser, T. Shinbrot, and G. Wurm, Nature Physics 16, 225 (2020).
  * Netz and Orland [2000] R. R. Netz and H. Orland, The European Physical Journal E 1, 203 (2000).
  * Kanduč and Podgornik [2007] M. Kanduč and R. Podgornik, The European Physical Journal E 23, 265 (2007).
  * Podgornik [1989] R. Podgornik, The Journal of chemical physics 91, 5840 (1989).
  * Naji  _et al._ [2013] A. Naji, M. Kanduč, J. Forsman, and R. Podgornik, The Journal of chemical physics 139, 150901 (2013).
  * Lau and Sokoloff [2020] A. Lau and J. Sokoloff, Physical Review E 102, 052606 (2020).
  * Solis and de la Cruz [2021] F. Solis and M. O. de la Cruz, The Journal of Chemical Physics 155, 104703 (2021).
  * Wang and Wang [2014] R. Wang and Z.-G. Wang, Physical review letters 112, 136101 (2014).
  * Wang [2010] Z.-G. Wang, Physical Review E 81, 021501 (2010).
  * Wang and Wang [2013] R. Wang and Z.-G. Wang, The Journal of chemical physics 139, 124702 (2013).
  * Wang and Wang [2015] R. Wang and Z.-G. Wang, The Journal of chemical physics 142, 104705 (2015).
  * Jiang and Wang [2018] J. Jiang and Z.-G. Wang, The Journal of chemical physics 148, 114105 (2018).
  * Podgornik and Žekš [1988] R. Podgornik and B. Žekš, Journal of the Chemical Society, Faraday Transactions 2: Molecular and Chemical Physics 84, 611 (1988).
  * Dean and Horgan [2004] D. S. Dean and R. R. Horgan, Physical Review E 69, 061603 (2004).
  * Siepmann and Sprik [1998] J. I. Siepmann and M. Sprik, The Journal of chemical physics 102, 511 (1998).
  * Reed  _et al._ [2007] S. K. Reed, O. J. Lanning, and P. A. Madden, The Journal of chemical physics 126, 084704 (2007).
  * Limmer  _et al._ [2013] D. T. Limmer, C. Merlet, M. Salanne, D. Chandler, P. A. Madden, R. Van Roij, and B. Rotenberg, Physical review letters 111, 106102 (2013).
  * Nguyen  _et al._ [2019] T. D. Nguyen, H. Li, D. Bagchi, F. J. Solis, and M. O. de la Cruz, Computer Physics Communications 241, 80 (2019).
  * Lin  _et al._ [2009] Y. Lin, A. Baumketner, S. Deng, Z. Xu, D. Jacobs, and W. Cai, The Journal of chemical physics 131, 10B608 (2009).
  * Gan and Xu [2011] Z. Gan and Z. Xu, Physical Review E 84, 016705 (2011).
  * Xu [2013] Z. Xu, Physical Review E 87, 013307 (2013).
  * dos Santos  _et al._ [2016] A. P. dos Santos, M. Girotto, and Y. Levin, The Journal of chemical physics 144, 144103 (2016).
  * Dwelle and Willard [2019] K. A. Dwelle and A. P. Willard, The Journal of Physical Chemistry C 123, 24095 (2019).
  * Dos Santos  _et al._ [2017] A. P. Dos Santos, M. Girotto, and Y. Levin, The Journal of chemical physics 147, 184105 (2017).
  * Son and Wang [2021] C. Y. Son and Z.-G. Wang, Proceedings of the National Academy of Sciences 118 (2021).
  * Mariappan  _et al._ [2010] C. Mariappan, T. Heins, and B. Roling, Solid State Ionics 181, 859 (2010).
  * Rotenberg and Salanne [2015] B. Rotenberg and M. Salanne, The journal of physical chemistry letters 6, 4978 (2015).
  * Lee and Perkin [2016] A. A. Lee and S. Perkin, The journal of physical chemistry letters 7, 2753 (2016).
  * Kondrat  _et al._ [2011] S. Kondrat, N. Georgi, M. V. Fedorov, and A. A. Kornyshev, Physical Chemistry Chemical Physics 13, 11359 (2011).
  * Kondrat  _et al._ [2014] S. Kondrat, A. A. Kornyshev, _et al._ , Physical review letters 113, 048701 (2014).
  * Bagchi  _et al._ [2020] D. Bagchi, T. D. Nguyen, and M. O. de la Cruz, Proceedings of the National Academy of Sciences 117, 19677 (2020).
  * Israelachvili [2015] J. N. Israelachvili,  _Intermolecular and surface forces_ (Academic press, 2015).
  * Belloni [2000] L. Belloni, Journal of Physics: Condensed Matter 12, R549 (2000).
  * Verwey  _et al._ [1948] E. J. W. Verwey, J. T. G. Overbeek, and K. Van Nes,  _Theory of the stability of lyophobic colloids: the interaction of sol particles having an electric double layer_ (Elsevier Publishing Company, 1948).
  * Trizac and Levin [2004] E. Trizac and Y. Levin, Physical Review E 69, 031403 (2004).
  * Pianegonda  _et al._ [2007] S. Pianegonda, E. Trizac, and Y. Levin, The Journal of chemical physics 126, 014702 (2007).
  * Colla and Levin [2010] T. E. Colla and Y. Levin, The Journal of chemical physics 133, 234105 (2010).
  * Colla  _et al._ [2012] T. E. Colla, A. P. dos Santos, and Y. Levin, The Journal of chemical physics 136, 194103 (2012).
  * De Soria  _et al._ [2016] M. I. G. De Soria, C. E. Álvarez, and E. Trizac, Physical Review E 94, 042609 (2016).
  * Trizac [2000] E. Trizac, Physical Review E 62, R1465 (2000).
  * Kjellander  _et al._ [1990] R. Kjellander, S. Marčelja, R. Pashley, and J. Quirk, The Journal of chemical physics 92, 4399 (1990).
  * Bowen and Sharif [1998] W. R. Bowen and A. O. Sharif, Nature 393, 663 (1998).
  * Gelbart  _et al._ [2000] W. M. Gelbart, R. F. Bruinsma, P. A. Pincus, and V. A. Parsegian, Physics today 53, 38 (2000).
  * Zhou  _et al._ [2019a] T. Zhou, K. Ioannidou, E. Masoero, M. Mirzadeh, R. J.-M. Pellenq, and M. Z. Bazant, Langmuir 35, 4397 (2019a).
  * Zhou  _et al._ [2019b] T. Zhou, K. Ioannidou, F.-J. Ulm, M. Z. Bazant, and R.-M. Pellenq, Proceedings of the National Academy of Sciences 116, 10652 (2019b).
  * Monfared  _et al._ [2020] S. Monfared, T. Zhou, J. E. Andrade, K. Ioannidou, F. Radjaï, F.-J. Ulm, and R. J.-M. Pellenq, Physical Review Letters 125, 255501 (2020).
  * Kjellander  _et al._ [1988] R. Kjellander, S. Marcelja, R. Pashley, and J. Quirk, The Journal of Physical Chemistry 92, 6489 (1988).
  * Pellenq  _et al._ [1997] R.-M. Pellenq, J. Caillol, and A. Delville, The Journal of Physical Chemistry B 101, 8584 (1997).
  * Pincus and Safran [1998] P. Pincus and S. Safran, EPL (Europhysics Letters) 42, 103 (1998).
  * Levin [1999] Y. Levin, Physica A: Statistical Mechanics and its Applications 265, 432 (1999).
  * Linse and Lobaskin [1999] P. Linse and V. Lobaskin, Physical review letters 83, 4208 (1999).
  * Netz [2001] R. R. Netz, The European Physical Journal E 5, 557 (2001).
  * Pellenq and Van Damme [2004] R. J.-M. Pellenq and H. Van Damme, Mrs Bulletin 29, 319 (2004).
  * Moreira and Netz [2000] A. G. Moreira and R. R. Netz, EPL (Europhysics Letters) 52, 705 (2000).
  * Moreira and Netz [2001] A. G. Moreira and R. R. Netz, Physical review letters 87, 078301 (2001).
  * Šamaj and Trizac [2011a] L. Šamaj and E. Trizac, Physical review letters 106, 078301 (2011a).
  * Šamaj and Trizac [2011b] L. Šamaj and E. Trizac, Physical Review E 84, 041401 (2011b).
  * Dos Santos and Levin [2019] A. P. Dos Santos and Y. Levin, Physical review letters 122, 248005 (2019).
  * Jungwirth and Tobias [2006] P. Jungwirth and D. J. Tobias, Chemical reviews 106, 1259 (2006).
  * Hofmeister [1888] F. Hofmeister, Pathol. Pharmakol 24, 247 (1888).
  * Otten  _et al._ [2012] D. E. Otten, P. R. Shaffer, P. L. Geissler, and R. J. Saykally, Proceedings of the National Academy of Sciences 109, 701 (2012).
  * Stern  _et al._ [2013] A. C. Stern, M. D. Baer, C. J. Mundy, and D. J. Tobias, The Journal of chemical physics 138, 114709 (2013).
  * Vaikuntanathan and Geissler [2014] S. Vaikuntanathan and P. L. Geissler, Physical review letters 112, 020603 (2014).
  * Cruz  _et al._ [2018] C. Cruz, A. Ciach, E. Lomba, and S. Kondrat, The Journal Physical Chemistry C 123, 1596–1601 (2018).
  * Hubbard [1959] J. Hubbard, Physical Review Letters 3, 77 (1959).
  * Stratonovich [1957] R. Stratonovich, in  _Soviet Physics Doklady_ , Vol. 2 (1957) p. 416.
  * Frydel [2015] D. Frydel, European Journal of Physics 36, 065050 (2015).
  * Cuvillier  _et al._ [1997] N. Cuvillier, M. Bonnier, F. Rondelez, D. Paranjape, M. Sastry, and P. Ganguly, Trends in Colloid and Interface Science XI , 118 (1997).
  * Borukhov  _et al._ [1997] I. Borukhov, D. Andelman, and H. Orland, Physical review letters 79, 435 (1997).
  * Borukhov  _et al._ [2000] I. Borukhov, D. Andelman, and H. Orland, Electrochimica Acta 46, 221 (2000).
  * Bruch  _et al._ [2022] D. Bruch, C. Balzer, and Z.-G. Wang, Journal of Chemical Physics 156 (2022).
  * Kornyshev [2007] A. A. Kornyshev, Double-layer in ionic liquids: paradigm change? (2007).
  * Eigen and Wicke [1954] M. Eigen and E. Wicke, The Journal of Physical Chemistry 58, 702 (1954).
  * Kornyshev and Vorotyntsev [1981] A. A. Kornyshev and M. A. Vorotyntsev, Electrochimica Acta 26, 303 (1981).
  * Nanjundiah  _et al._ [1997] C. Nanjundiah, S. McDevitt, and V. Koch, Journal of The Electrochemical Society 144, 3392 (1997).
  * Lockett  _et al._ [2008] V. Lockett, R. Sedev, J. Ralston, M. Horne, and T. Rodopoulos, The Journal of Physical Chemistry C 112, 7486 (2008).
  * Fedorov and Kornyshev [2008] M. V. Fedorov and A. A. Kornyshev, The Journal of Physical Chemistry B 112, 11868 (2008).
  * Qing  _et al._ [2021] L. Qing, S. Zhao, and Z.-G. Wang, The Journal of Physical Chemistry B 125, 625 (2021).