# The Differential Capacitance as a probe for the electric double layer
structure and the electrolyte bulk composition

Peter Cats<EMAIL_ADDRESS>Institute for Theoretical Physics, Center for Extreme
Matter and Emergent Phenomena, Utrecht University, Princetonplein 5, 3584 CC,
the Netherlands  René van Roij<EMAIL_ADDRESS>Institute for Theoretical
Physics, Center for Extreme Matter and Emergent Phenomena, Utrecht University,
Princetonplein 5, 3584 CC, the Netherlands

###### Abstract

In this work we theoretically study the differential capacitance of an aqueous
electrolyte in contact with a planar electrode, using classical Density
Functional Theory, and show how this measurable quantity can be used as a
probe to better understand the structure and composition of the electric
double layer at play. Specifically, we show how small trace amounts of
divalent ions can influence the differential capacitance greatly, and also how
small ions dominate its behaviour for high electrode potentials. In this
study, we consider primitive model electrolytes, and not only use the standard
definition of the differential capacitance, but also derive a new expression
from mechanical equilibrium in a planar geometry. This expression reveals
explicitly that the first layer of ions near the charged surface is key to its
understanding. Our insights might be used as a guide in experiments to better
understand the electrolyte-electrode interface as well as the (composition of
the) bulk electrolyte.

††preprint: AIP/123-QED

##  I Introduction

At the moment, energy storage is one of the major challenges in the battle
against climate change. There are many “green" methods nowadays to transform
one energy form into another [1], with electricity as the favoured end
product, but the options to store it efficiently on a large scale are limited.
With solar and wind energy on the rise, the era of fossil energy sources will
come to and end. However, sunshine and wind are intermittent and therefore the
amount of harvested energy is not constant and can fluctuate quite drastically
over time and geographic position, which requires an increasing flexibility of
the power system [2, 3]. Efficient ways of storing energy can increase the
flexibility of the power system [4, 5, 6, 7], allowing us to successfully make
the transition towards renewable energy.

Technologies for energy storage vary immensely, and several classifications
can be considered, for instance the actual technology used (e.g. chemical,
electrochemical, electrical, mechanical and thermal) or the scale of storage
(e.g. small scale, short term, large scale, long term); each of these
subcategories has its own strengths, weaknesses and range of applicability [8,
6]. One of these technologies employs electric double layers (EDLs), in which
the spatial distribution of mobile charges of electrolytes in contact with
electrodes allows for storage of electric energy. Basically, such a device is
comprised of a charged conductor (an electrode) immersed in an electrolyte,
where the electron charge on the electrode is screened by the ions in the
electrolyte. Such a system behaves in some sense similar to a conventional
capacitor of two electrodes, however with a much smaller distance between the
oppositely charged entities, which in this case consists of an immobile
electrode with, say, a positive charge in close contact with a layer of mobile
ions in the electrolyte with a net negative charge. Generally, the electrodes
of EDL capacitors have a porous structure with a gigantic surface-to-volume
ratio (up to $2000$ m2 g-1) [9], and no chemical reactions take place at the
surface. The large area leads to a large capacitance and therefore a large
amount of energy stored compared to conventional capacitors and the absence of
reactions leads to long lifetimes compared to conventional batteries, which
makes them “greener" as well. One of the disadvantages of EDL capacitors is
the rather small energy density compared to common Lithium-ion batteries (by
roughly a factor of 20) [4].

Within a density functional theory that captures Coulomb and hard-sphere
interactions accurately [10], we calculate the differential capacitance for
several electrolytes, changing the size and valency of the ions similar to
work done in Refs. 11, 12, 13, 14, 15, 16, 17, 18, 19, 20. However, we derive
a new relation for the differential capacitance, which explicitly shows that
the response of the first layer of ions near the electrode to the applied
potential on the electrode is key to its understanding. We utilize this
insight to understand two- and three-component electrolytes and explain how
impurities can have a large effect on the differential capacitance.

##  II Differential Capacitance

###  II.1 The System

As mentioned in the introduction, we wish to focus on EDL capacitors, of which
there is a variety of different types, each having a different electrode-
electrolyte combination [21]. For now, we consider two porous (carbon)
electrodes held at a potential difference $\Psi$ immersed in an aqueous
electrolyte at room temperature $T$ and dielectric constant $\epsilon$. We
will focus on a single pore inside one electrode, and model it as two planar
surfaces with the same static surface potential $\Phi_{0}$, separated by a
distance $H$ [22, 23], as depicted in Fig. 1. We stress that the surface
potential $\Phi_{0}$ on the surfaces of the pore is w.r.t. a potential in a
macroscopic electrolyte reservoir (the space between the two electrodes),
while $\Psi$ is the potential difference between the two electrodes in the EDL
capacitor.

Figure 1: An illustration of an EDL capacitor and its structure. The EDL
capacitor within the dashed box contains two porous electrodes (black dots),
held at a potential difference $\Psi$, through which the electrolyte (blue)
can flow. Zooming in on a porous electrode reveals a lot of structure which is
depicted on top. Throughout this manuscript the focus lies on one tiny part of
the electrode (the pore of interest), which we consider our system, and model
it as two parallel charged surfaces separated by a distance $H$ being part of
the same electrode. Hence, the potential on the two surfaces, denoted by
$\Phi_{0}$, is the same and measured w.r.t. the potential in the reservoir,
which we assume to be grounded. Invoking the primitive model reduces the
solvent to a continuous dielectric medium at temperature $T$ and with uniform
dielectric constant $\epsilon$.

We will employ the primitive model (PM) in which the solvent is treated as a
continuous dielectric medium. The properties of the solvent can then be
completely captured by a single parameter: the Bjerrum length
$\lambda_{B}=e^{2}\beta/4\pi\epsilon\epsilon_{0}$, where $e$ is the proton
charge and $\beta=1/k_{B}T$ the inverse temperature. The Bjerrum length of
water at room temperature is $\lambda_{B}=0.73$ nm, which will be kept
constant throughout this work. The ions of species $j$ are modelled as hard
spheres with diameter $d_{j}$ and a centro-symmetric charge $ez_{j}$ with
$z_{j}$ denoting the valency; the pair potential $u_{ij}(r)$ between ions of
species $i$ and $j$ separated by a distance $r$ is described by

$\displaystyle\beta u_{ij}(r)=\begin{cases}\infty,\quad&r<d_{ij};\\\
z_{i}z_{j}\frac{\lambda_{B}}{r},\quad&r\geq d_{ij},\end{cases}$ (1)

where $d_{ij}=(d_{i}+d_{j})/2$. Note that the PM neglects dispersion and
polarizability. We denote the bulk concentration of ion species $j$ in the
macroscopic electrolyte reservoir between the two electrodes by $\rho_{b,j}$,
and define the total ionic packing fraction in the bulk as

$\displaystyle\eta_{b}=\frac{\pi}{6}\sum_{j}d_{j}^{3}\rho_{b,j}.$ (2)

When all particles have the same diameter and valency, the model is referred
to as the Restrictive Primitive Model (RPM). The bulk RPM is completely
characterized by the ionic blk concentration and the dimensionless temperature
$T^{*}=d/\lambda_{B}$. We stress that the medium and the diameters are chosen
such that $T^{*}\gg T_{c}^{*}=0.05$ the critical temperature [24], e.g. in our
manuscript $\lambda_{B}=0.73$ nm is fixed and the ion diameters range from
$d_{j}=0.2$ nm to $d_{j}=1$ nm, for which $T^{*}$ ranges from $T^{*}=0.27$ to
$T^{*}=1.37$, respectively.

The (non-electrostatic part of the) external potential that determines the
pore size reads

$\displaystyle V_{ext}^{j}(x)=\begin{cases}\infty\quad
x\leq\frac{d_{j}}{2};\\\ 0\quad\frac{d_{j}}{2}<x<H-\frac{d_{j}}{2};\\\
\infty\quad x\geq H-\frac{d_{j}}{2},\\\ \end{cases}$ (3)

which describes two planar hard walls positioned at $x=0$ and $x=H$.
Throughout this manuscript, the surface separation is fixed at $H=5$ nm. Note
that the electrostatic part of the external potential, due to the surface
charge on the electrodes, will be taken into account through the Poisson
equation, as discussed in Appendix A.

Within this simplification, we can controllably study the differential
capacitance as a function of the surface potential $\Phi_{0}$, for several
choices of the ion diameter $d_{j}$, valency $z_{j}$, and ionic bulk
concentrations in the reservoir $\rho_{b,j}$.

###  II.2 Thermodynamics in the Pore

Let us now focus on a single pore as depicted in Fig. 1. The remainder of the
electrode and the space between the electrodes we consider to be the ion and
heat reservoir for the pore of interest. Assuming that the space between the
electrodes is large enough, we can treat the pore of interest grand
canonically. This justifies the introduction of a bulk ionic concentration
$\rho_{b,j}$ (corresponding to a fixed ionic chemical potential $\mu_{j}$) in
the reservoir that is independent of the electric potential difference $\Psi$
between the electrodes. Exchange of ions between the system and reservoir
takes place and we assume equilibrium. The surface potential $\Phi_{0}$ on the
surfaces in the pore is fixed, which allows exchange of electronic charge
between the electrode and a charge reservoir; however, no charge exchange
takes place between electrode and ions, i.e. no electrochemistry is
considered. The pore also has a fixed volume $V=AH/2$ (rigid electrode) with
total surface area $A$, and temperature $T$. To summarize, the system is
described by a $\\{\mu\\},V,T,\Phi_{0},A,H$ ensemble, characterized by the
grand potential $\Omega(\\{\mu\\},V,T,\Phi_{0},A,H)$ [10], where the brackets
$\\{\mu\\}$ denote the set of chemical potentials of each of the ion species.
The differential of $\Omega$ for this system reads

$\displaystyle\mathrm{d}\Omega=$
$\displaystyle-p\mathrm{d}V-S\mathrm{d}T-\sum_{j}N_{j}\mathrm{d}\mu_{j}+$ (4)
$\displaystyle-Q\mathrm{d}\Phi_{0}+\gamma\mathrm{d}A-f\mathrm{d}H,$

where $N_{j}$ is the average number of ions of species $j$ in the pore, $p$
the pressure, $Q$ the electronic charge on both surfaces of the pore, $\gamma$
the surface tension, and $f$ the force between the surfaces that are separated
by a distance $H$. It is convenient to separate the volumetric and areal
contribitions in Eq. (4), and write
$\Omega=-p(T,\\{\mu\\})V+\gamma(T,\\{\mu\\},\Phi_{0},H)A$ as two separate
differentials [10]

$\displaystyle\mathrm{d}p$
$\displaystyle=s_{b}\mathrm{d}T+\sum_{j}\rho_{b,j}\mathrm{d}\mu_{j},$ (5)
$\displaystyle\mathrm{d}\gamma$
$\displaystyle=-s_{s}\mathrm{d}T-\sum_{j}\Gamma_{j}\mathrm{d}\mu_{j}-\sigma\mathrm{d}\Phi_{0}-\frac{f}{A}\mathrm{d}H,$
(6)

where $S$ and $N_{j}$ have been separated into a volumetric and areal part
according to $S=Vs_{b}+As_{s}$, and $N_{j}=V\rho_{b,j}+A\Gamma_{j}$ with
$\rho_{b,j}$ the bulk concentration of species $j$ and $\Gamma_{j}$ the
adsorption. The average electronic surface charge density is denoted by
$\sigma=Q/A$ and obeys the charge neutrality condition

$\displaystyle\sigma=-\sum_{j}ez_{j}\Gamma_{j}.$ (7)

Eq. (5) is the Gibbs-Duhem equation and Eq. (6) the Lipmann equation.

Now the question arises: what does this have to do with capacitances? The
short answer to this question is that the (areal) differential capacitance is
defined by

$\displaystyle C_{\mu,T,H}$
$\displaystyle=\left(\frac{\partial\sigma}{\partial\Phi_{0}}\right)_{\mu,T,H}$
(8)
$\displaystyle=-\left(\frac{\partial^{2}\gamma}{\partial\Phi_{0}^{2}}\right)_{\mu,T,H},$
(9)

which opens the door for the longer answer. First of all, it turns out that
there is not one differential capacitance $C$, but a whole set of differential
capacitances, one for each ensemble. This is similar to constant-volume and
constant-pressure heat capacities, and isothermal and isentropic
compressibilities, for instance. Throughout this manuscript we will study the
differential capacitance at fixed temperature $T$ and surface separation $H$,
and simplify the notation: the ensemble in Eq. (4) gives rise to $C_{\mu}$,
while an $\\{N\\},V,T,\Psi,H$ ensemble would give rise to $C_{N}$. Those two
are, in fact, not unrelated and follow similar thermodynamic identities as
those of the heat capacities at constant volume or pressure [25, 10]. We will
focus on $C_{\mu}(\Phi_{0},\\{\mu\\})$.

###  II.3 Differential Capacitance Revisited

Before tackling the primitive-model electrolyte confined between two hard
walls by Density Functional Theory (DFT), let us first discuss the mechanical
equilibrium in a planar geometry. The $xx$ (normal) component of the pressure
tensor within the pore follows from the force balance [26, 27]

$\displaystyle p_{xx}(H)=$
$\displaystyle\sum_{j}\int_{0}^{H}\mathrm{d}x\rho_{j}(x)F_{x}^{j}(x),$ (10)

where $\rho_{j}(x)$ is the local density of species $j$,
$F_{x}^{j}(x)=-\partial V_{j}(x)/\partial x$ the normal force on the surface
due to the full interaction potential $V_{j}(x)$ (including the electrostatic
interactions due to the surface charge). If the integral in Eq. (10) does not
depend on $H$, which is the case when the range of $F_{x}^{j}$ is much smaller
than $H$, then $p_{xx}(H)=p$ is the actual bulk pressure of the electrolyte
and Eq. (10) reduces to [26, 27]

$\displaystyle\beta
p=\sum_{j}\rho_{j}\left(\frac{d_{j}^{+}}{2}\right)-\frac{\beta\sigma^{2}}{2\epsilon\epsilon_{0}},$
(11)

where $d_{j}^{+}=\lim_{x\downarrow d_{j}}x$. This sum rule relates the bulk
pressure $p$ of the ion reservoir to the local contact densities
$\rho_{j}(d_{j}^{+}/2)$ and the surface charge density $\sigma$. In the case
of two charged hard parallel surfaces separated by a finite distance, the
integral in Eq. (10) does depend on $H$, especially at low ionic bulk
concentrations. However, with the parameters chosen throughout this
manuscript, this only plays a minor role and can be neglected for all
practical purposes. Let us note that for soft particle-wall interactions the
integral in Eq. (10) does not reduce to a contact value and needs to be fully
evaluated. Nonetheless, the gradient of the particle-wall interaction is by
construction largest at values close to the contact distance $d_{j}/2$ at
which it gives the largest contribution to the integral.

One interesting feature of Eq. (11) is that the global bulk pressure is
independent of the surface potential $\Phi_{0}$, while the local quantities
are not. Therefore, taking the derivative of Eq. (11) w.r.t. $\Phi_{0}$ leads
to

$\displaystyle
0=\sum_{j}\left(\frac{\partial\rho_{j}(d_{j}^{+}/2)}{\partial\Phi_{0}}\right)_{\mu}-\frac{\beta\sigma}{\epsilon\epsilon_{0}}\left(\frac{\partial\sigma}{\partial\Phi_{0}}\right)_{\mu},$
(12)

where we recognize the differential capacitance given in Eq. (8) in the second
term. Hence, Eq.(12) can be rewritten as

$\displaystyle C_{\mu}$
$\displaystyle=\frac{\epsilon\epsilon_{0}}{\beta\sigma}\sum_{j}\left(\frac{\partial\rho_{j}(d_{j}/2)}{\partial\Phi_{0}}\right)_{\mu},$
(13)

where in the summand one can recognize a fluctuation profile as discussed in
Ref. 28. Within the RPM, Eq. (13) simplifies to

$\displaystyle
C_{\mu}=\frac{\epsilon\epsilon_{0}}{\beta\sigma}\left(\frac{\partial\rho_{N}(d/2)}{\partial\Phi_{0}}\right)_{\mu},$
(14)

where we introduce the total ion density

$\displaystyle\rho_{N}(x)=\sum_{j}\rho_{j}(x).$ (15)

Hence, the differential capacitance is to a large extent dictated by the
response of the number density $\rho_{N}(d/2)$ at contact (first layer of
ions) to the potential $\Phi_{0}$. As stated previously, this key-concept will
stand even for soft particle-wall interactions; the contact value will be
replaced by an integral, whose integrand is peaked near $x=d_{j}/2$.

Interestingly, another observation leads to yet another expression for the
differential capacitance. Due to the ions having a hard core interaction with
the wall and a central charge, the charge density for $x<d/2$ (also for
$x>H-d/2$) vanishes, leading to a linear drop of the electrostatic potential
profile $\Phi(x)$ in that region. Hence, the potential at $x=d/2$ is given by

$\displaystyle\Phi(d/2)=\Phi_{0}+\left.\frac{\partial\Phi(x)}{\partial
x}\right|_{x=0}\frac{d}{2}=\Phi_{0}-\frac{d}{2\epsilon\epsilon_{0}}\sigma,$
(16)

where we used Gauss’ law. Taking again the derivative w.r.t. the surface
potential $\Phi_{0}$ and reordering results in

$\displaystyle
C_{\mu}=\frac{2\epsilon\epsilon_{0}}{d}\left[1-\left(\frac{\partial\Phi(d/2)}{\partial\Phi_{0}}\right)_{\mu}\right].$
(17)

The region between $x=0$ and $x=d/2$ is sometimes referred to as the Stern
layer [29], with the corresponding Stern-layer capacitance [17, 30, 31, 32]
$C_{s}=2\epsilon\epsilon_{0}/d$, although that will not be used in this work
111The Stern layer $C_{EDL}$ separates the diffuse layer from the electrode,
each layer having their own capacitance, which in a serial circuit gives the
total capacitance of the EDL:
$\displaystyle\frac{1}{C_{EDL}}=\frac{1}{C_{s}}+\frac{1}{C_{DL}}.$ (18) With
little effort, one can show using Eqs. (17) and (45) that, by introducing the
Stern-layer capacitance, the capacitance of the diffuse layer reads
$\displaystyle
C_{DL}=C_{s}\left[\left(\frac{\partial\Phi_{0}}{\partial\Phi(d/2)}\right)_{\mu}-1\right],$
(19) which can take negative values in case of overscreening. Hence, splitting
the EDL into a diffuse layer and the stern layer seems not to be desirable for
a generic system, and we will not further distinguish between Stern layers and
other layers, but instead talk about the total capacitance $C_{\mu}$..

The relations in Eqs. (13) and (17) will allow for a new and independent
investigation into the differential capacitance, which helps in understanding
the underlying physics.

###  II.4 Density Functional Theory

We will use classical DFT to study the PM electrolytes. We refer to previous
literature for the details [34, 35, 36]. In short, DFT involves a grand
potential functional, which reads

$\displaystyle\Omega[\\{\rho\\}]=$
$\displaystyle\mathcal{F}_{id}[\\{\rho\\}]+\mathcal{F}_{ex}^{HS}[\\{\rho\\}]+\mathcal{F}_{ex}^{ES}[\\{\rho\\}]-$
$\displaystyle\sum_{j}\int\mathrm{d}\mathbf{r}\rho_{j}(\mathbf{r})\left(\mu_{j}-V_{ext}^{j}(\mathbf{r})\right)-Q\Phi_{0},$
(20)

where $V_{ext}^{j}$ is the (non-electric) external potential of Eq. (3) acting
on particles of species $j$, and $Q=\sigma A$ is the total charge on the
surface in the pore. The first term $\mathcal{F}_{id}[\\{\rho\\}]$ is the
Helmoltz free energy functional for an ideal (non-interacting) system, the
second term $\mathcal{F}_{ex}^{HS}[\\{\rho\\}]$ accounts for the hard-sphere
interactions and is dealt with by Fundamental Measure Theory White-Bear II
[35], and the third term $\mathcal{F}_{ex}^{ES}[\\{\rho\\}]$ describes the
Coulombic interactions for which we invoke the MSAc functional described
extensively in Refs. 37, 38 (see also Appendix A).

The grand potential functional $\Omega[\\{\rho\\}]$ is, in fact, minimized by
the equilibrium density profiles $\\{\rho_{0}\\}$ [34], i.e.

$\displaystyle\left.\frac{\delta\Omega[\\{\rho\\}]}{\delta\rho_{j}(\mathbf{r})}\right|_{\\{\rho\\}=\\{\rho_{0}\\}}=0.$
(21)

Moreover, the grand potential at its minimum is the actual grand potential of
the system, i.e. $\Omega[\\{\rho_{0}\\}]=\Omega(\\{\mu\\},V,T,\Phi_{0},H)$.
Hence, DFT is a powerful theoretical framework to combine thermodynamics and
structure as it gives access to the thermodynamics of the system via the
density profiles. Furthermore, it gives directly both the ionic charge density
profile $q(x)=\sum_{j}z_{j}\rho_{j}(x)$ and via the Poisson equation also the
potential profile $\Phi(x)$ from which the surface charge follows, where we
note that in a planar geometry the profiles only depend on the out-of-plane
coordinate $x$. In principle, we have thus access to $C_{\mu}$ via any of the
routes in Eqs. (8), (9), (13), and (17) laid out in the previous sections;
however, the results that will be presented are calculated from Eq. (8), while
the other definitions are used as a tool to interpret the data. Nevertheless,
we have validated each of the Eqs. (8), (9), (13), and (17) for the sake of
consistency.

###  II.5 Simulations

Although DFT is a very powerful framework, the intrinsic functionals
$\mathcal{F}_{ex}^{HS}$ and $\mathcal{F}_{ex}^{ES}$ are approximate and need
to be tested against simulations or experiments. Although we do not go into
detail here on any specific type of simulation, we wish to present a novel
sumrule that connects quantities that can be measured in simulation but not
with DFT. The quantity of interest throughout this manuscript is the
differential capacitance. From elementary statistical physics one can derive
the relation [39, 40, 25]

$\displaystyle C(\Phi_{0})=\frac{\beta}{A}\left(\langle Q^{2}\rangle-\langle
Q\rangle^{2}\right),$ (22)

where the brackets $\langle\ldots\rangle$ indicate either a time or an
ensemble average.

Inspired by Ref. 28, we also found a new expression to calculate the
differential capacitance, namely

$\displaystyle C$
$\displaystyle=\frac{\epsilon\epsilon_{0}}{\sigma}\sum_{j}\left[\langle\hat{\rho}_{j}(d_{j}/2)Q\rangle-\langle\hat{\rho}_{j}(d_{j}/2)\rangle\langle
Q\rangle\right],$ (23)

where
$\hat{\rho}_{j}(\mathbf{r})=\sum_{i=1}^{N_{j}}\delta(\mathbf{r}-\mathbf{r}_{i})$
denotes the density operator. Note that we use the notation
$\sigma=\langle\sigma\rangle$ throughout the manuscript, but it is important
to keep in mind that the surface charge $\sigma=Q/A$ fluctuates as a
consequence of fixing the surface potential $\Phi_{0}$.

Employing DFT causes the loss of direct information on the fluctuations,
because DFT only returns equilibrium density profiles. Hence, within our DFT
approach one does not have direct access to the quantities such as $\langle
Q^{2}\rangle$ and $\langle\hat{\rho}_{j}(d_{j}/2)Q\rangle$. However, they can
in principle be measured within simulations.

##  III Capacitance Curves

###  III.1 First Glance: Symmetric Electrolyte

Figure 2: The differential capacitance $C_{\mu}$ for the RPM with monovalent
ions with diameter $d=0.5$ nm and for a surface separation $H=5$ nm and
Bjerrum length $\lambda_{B}=0.73$ nm. The curves are shown for concentrations
ranging from $\rho_{b}=1.3$ mM (light-blue) to $\rho_{b}=5.1$ M (green) with
several concentrations in between. The concentrations for the other colored
lines are given in the legend; the concentration for the black lines are
equally spaced to illustrate the trend. The legend also shows the
corresponding bulk packing fractions $\eta_{b}$.

The system that we investigate first is an RPM electrolyte consisting of two
monovalent species of ions (i.e. $z_{+}=-z_{-}=1$), so that
$\rho_{b,+}=\rho_{b,-}\equiv\rho_{b}$. In Fig. 2 we show the differential
capacitance, calculated via Eq. (8), as a function of the applied voltage
$\Phi_{0}$ for ions with $d=0.5$ nm at several ionic bulk packing fractions
ranging from $\eta_{b}=0.0001$ to $\eta_{b}=0.4$, corresponding to bulk ionic
concentrations ranging from $\rho_{b}\approx 1$ mM to $\rho_{b}\approx 5$ M.
As stated previously, the surface separation is fixed at $H=5$ nm and the
solvent is characterized by $\lambda_{B}=0.73$ nm. The foremost noticeable
feature of Fig. 2 is the crossover around a bulk packing fraction
$\eta_{c}\approx 0.12$ (blue line) from the so-called camel-shaped
($\eta_{b}<\eta_{c}$) to the bell-shaped curves ($\eta_{b}>\eta_{c}$). This
crossover has been a subject in many theoretical and experimental studies [17,
41, 42, 43, 44, 45, 46, 47, 16, 48], and is ascribed to excluded-volume
effects. Let us, for convenience, introduce $\Phi^{*}$ as the surface
potential at which the capacitance $C_{\mu}(\Phi_{0})$ takes its maximum value
at a given $\rho_{b}$ or $\eta_{b}$, i.e. $C_{\mu}(\Phi^{*})\geq
C_{\mu}(\Phi_{0})$. Then the camel-shaped capacitance curves correspond to a
finite $\Phi^{*}$, while the bell-shaped capacitance curves are characterized
by a vanishing $\Phi^{*}$.

Figure 3: (a) The number density profiles $\rho_{N}(x)$ for the surface
potentials $\Phi_{0}=0.1$ V (solid), $\Phi_{0}=0.2$ V (dashed), and
$\Phi_{0}=0.5$ V (dotted) at a bulk concentration of $\rho_{b}\approx 0.13$ M
(i.e. $\eta_{b}=0.01$). The inset shows the local packing fraction $\eta(x)$
as defined in Eq. (24). The potential at which the capacitance has its maximum
is $\Phi^{*}\approx 0.2$ V. For $\Phi_{0}<\Phi^{*}$ the profiles only have one
peak at $x=d/2=0.25$ nm, while for $\Phi_{0}>\Phi^{*}$ the profiles have two
peaks; the second located at $x=3d/2=0.75$ nm. Around the maximum at
$\Phi_{0}\approx\Phi^{*}$ the second layer starts forming. (b) The number
density profiles $\rho_{N}(x)$ for the surface potentials $\Phi_{0}=0$ V
(solid) and $\Phi_{0}=0.2$ V (dashed) at a bulk concentration of
$\rho_{b}\approx 3.6$ M (i.e. $\eta_{b}=0.28$). The inset shows the local
packing fraction $\eta(x)$ as defined in Eq. (24). This bulk concentration
$\rho_{b}\approx 3.6$ M corresponds to the same-colored (purple) bell-shaped
capacitance curve in Fig. 2.

To understand this crossover, one has to understand what causes the maximum in
$C_{\mu}(\Phi_{0})$. Consider the orange camel-shaped capacitance curve in
Fig. 2 belonging to $\rho_{b}\approx 0.13$ M (i.e. $\eta_{b}=0.01<\eta_{c}$)
with $\Phi^{*}\approx 0.2$ V. For the same $\rho_{b}$ (and in the same color)
we present in Fig. 3(a) the number density profiles $\rho_{N}(x)$ for
$\Phi_{0}=0.1$ V (solid), $\Phi_{0}=0.2$ V (dashed), and $\Phi_{0}=0.5$ V
(dotted). The inset shows, for the same state points, the profile of the
weighted packing fraction, defined by 222$\eta(x)$ is identical to the
weighted density $n_{3}(x)$ from Fundamental Measure Theory [35].

$\displaystyle\eta(x)=\pi\sum_{j}\int_{x-d_{j}/2}^{x+d_{j}/2}\mathrm{d}x^{\prime}\rho_{j}(x^{\prime})\left[\left(\frac{d_{j}}{2}\right)^{2}-(x-x^{\prime})^{2}\right],$
(24)

which reduces to Eq. (2) in the homogeneous bulk. Interestingly, while the
density profile for $\Phi_{0}<\Phi^{*}$ (solid) only has one (contact) peak at
$x=d/2=0.25$ nm, the one for $\Phi_{0}>\Phi^{*}$ (dotted) has a second peak at
$x=3d/2=0.75$ nm. Hence, the maximum in $C_{\mu}$ indicates a structural
change, analogous to peaks in the heat capacity indicating (smooth) changes of
the thermal occupancy of microstates. When we consider the purple bell-shaped
capacitance curve in Fig. 2, belonging to $\rho_{b}\approx 3.6$ M (i.e.
$\eta_{b}=0.28>\eta_{c}$), and investigate the corresponding $\rho_{N}(x)$ and
$\eta(x)$ in Fig. 3(b) for $\Phi_{0}=0$ V (solid) and $\Phi_{0}=0.2$ V
(dashed), we find the presence of the second peak in $\rho_{N}(x)$ for both
(and in fact for all) surface potentials. Hence, the maximum in
$C_{\mu}(\Phi_{0})$ is clearly related to the onset of layering of counter-
ions near the surface, as is in fact also consistent with Eq. (13). Moreover,
considering the profile of the weighted packing fraction (see insets in Figs.
3(a) and (b)) at $x=d/2=0.25$ nm, one finds that its value is close to or even
larger than the packing fraction at which hard-sphere freezing takes place in
the bulk (i.e. around $\eta_{b}\approx 0.5$). This suggests that the first
layer of ions might in fact even be frozen [50], and invites a study into the
in-plane structure of ions at the charged surface. However, that is beyond the
scope of this manuscript.

To characterize and interpret the crossover from camel- to bell-shaped curves
further, let us consider the crossover potential $\Phi^{*}(\eta_{b})$ as
function of the bulk packing fraction $\eta_{b}$ plotted in Fig. 4, for three
ion diameters. The dashed line at $\eta_{c}\approx 0.12$ represents for
$d=0.5$ nm the capacitance curve crossing over from two maxima at
$\Phi_{0}=\Phi^{*}>0$ (camel) to one maximum at $\Phi_{0}=\Phi^{*}=0$ (bell).
From Fig. 4, it is evident that the crossover is rather gradual and therefore
we refer to a crossover rather than a transition. This is in line with the
maximum in the capacitance curves being finite rather than infinite, the
latter case would imply a thermodynamic phase transition. This graduality is
found in the density profiles as well: where profiles for
$\Phi_{0}\ll\Phi^{*}$ and $\Phi_{0}\gg\Phi^{*}$ can be easily discerned by the
absence and presence of a second peak at $x=1.5d=0.75$ nm (see Fig. 3(a)),
those for potentials close to $\Phi^{*}$ cannot. The same holds for density
profiles at $\Phi_{0}=0$; for low bulk concentrations layering is almost
absent, while at high bulk concentrations layering is clearly present (see
Fig. 3(b)), but one cannot qualitatively discern the two cases close to
$\eta_{b}\approx\eta_{c}$.

Figure 4: The surface potential $\Phi^{*}$ at which $C_{\mu}(\Phi_{0})$ has
its maximum as function of the bulk concentration expressed as the bulk
packing fraction $\eta_{b}$, for the RPM with ion diameter $d=0.2$ nm
(yellow), $d=0.5$ nm (blue), and $d=1$ nm (orange). The vertical dashed line
is located at $\eta_{b}=0.12$ and represents the crossover for ions with
$d=0.5$ nm from camel-shaped capacitance curves to bell-shaped capacitance
curves.

The results so far were all for ions with a common diameter $d=0.5$ nm and
$\lambda_{B}=0.73$ nm. Let us now investigate the effect of ion size, keeping
$\lambda_{B}$ fixed at $\lambda_{B}=0.73$ nm. To this end, we next consider
ions of a common diameter $d=0.2$ nm and $d=1$ nm, for which
$C_{\mu}(\Phi_{0})$ is presented in Fig. 5 for a variety of ionic bulk
concentrations. There is little qualitative difference in $C_{\mu}(\Phi_{0})$
for the different ion sizes. However, quantitatively, one finds that the
smaller the ions, the larger $C_{\mu}(\Phi_{0})$ at a given bulk packing
fraction $\eta_{b}$. It is important to note that ionic bulk concentration
$\rho_{b}$ and packing fraction $\eta_{b}$, although linearly related, are not
the same quantity; small ions have a higher molarity at the same packing
fraction and the balance between packing and electrostatics is quite
different. Experimentally, one often characterizes bulk concentrations in
molarity rather than packing fractions, which renders the ion diameter an
important fit parameter. Fig. 5 shows that the ion size can determine the
shape of $C_{\mu}(\Phi_{0})$ at a given molarity $\rho_{b}$. We note, for
instance, that the camel-bell crossover packing fraction $\eta_{c}$ only
changes slightly from $\eta_{c}\approx 0.1$ for $d=1$ nm to $\eta_{c}\approx
0.15$ for $d=0.2$ nm, while the crossover bulk concentration changes by two
orders of magnitude from $\rho_{c}\approx 32$ M for $d=0.2$ nm to
$\rho_{c}=0.16$ M for $d=1$ nm. Fig. 5 also shows a bell-shaped capacitance
curve for sufficiently concentrated electrolytes, and a camel-to-bell
crossover for dilute electrolytes at a sufficiently high electrode potential,
with a higher $\Phi^{*}$ for lower $\eta_{b}$ and smaller $d$. The generally
lower capacitance for larger ions can be understood qualitatively in terms of
the close proximity of smaller ions to the electrode. This is explicit in Eq.
(17), in which $C_{\mu}$ is inversely proportional to the ion diameter $d$.
Hence, a larger ion diameter results in an overall smaller
$C_{\mu}(\Phi_{0})$.

The results presented in this section evidently show that ion size is a key
parameter that largely determines the magnitude of the differential
capacitance at a given ionic bulk concentration.

Figure 5: The differential capacitance for ions with diameter (a) $d=0.2$ nm
and (b) $d=1$ nm for a surface separation $H=5$ nm and Bjerrum length
$\lambda_{B}=0.73$ nm. The bulk concentration/packing fraction at which the
camel-bell crossover occurs is depicted in dark blue, i.e. $\eta_{c}\approx
0.1$ ($\rho_{c}\approx 32$ M) for $d=0.2$ nm and $\eta_{c}\approx 0.16$
($\rho_{c}\approx 0.16$ M) for $d=1$ nm.

###  III.2 Second Glance: Asymmetric Electrolytes

Let us now consider an asymmetric electrolyte consisting of cations with
diameter $d_{+}=0.4$ nm and anions with diameter $d_{-}=0.6$ nm. The voltage-
dependent differential capacitance for this system is presented in Fig. 6(a),
again for bulk packing fractions ranging from $\eta_{b}=0.0001$
($\rho_{b}\approx 1.1$ mM) to $\eta_{b}=0.4$ ($\rho_{b}\approx 4.5$ M). The
asymmetry in the electrolyte is reflected by the asymmetry of Fig. 6(a) w.r.t.
$\Phi_{0}\rightarrow-\Phi_{0}$; the smaller cations can approach the surface
to closer distances than the anions, causing a larger capacitance for negative
potentials compared to positive potentials. The capacitance
$C_{\mu}(\Phi_{0})$ at its maximum for the camel-shaped curves at $\Phi_{0}<0$
is therefore larger than that at $\Phi_{0}>0$, which motivates the
introduction of the notation $\Phi^{*}_{+}$ and $\Phi_{-}^{*}$ for the maxima
located at positive and negative potentials, respectively. For high bulk
concentrations with only one maximum (bell-shaped curves), we keep using the
notation $\Phi^{*}$, which is not necessarily close to zero now that the
symmetry is broken. In Fig. 6 we see that $\Phi^{*}$ does indeed not vanish in
the asymmetric electrolyte, rather it is located at negative potentials, and
seems for the present choice of parameters to converge to
$\Phi^{*}\approx-0.03$ V at high bulk concentrations. This also implies that
the potential of zero charge takes non-vanishing values.

In Fig. 6(b) the differential capacitance of the asymmetric electrolyte of
Fig. 6(a) is compared to those of the symmetric RPM. We replot the six colored
differential capacitance curves of the asymmetric electrolyte of Fig. 6(a) in
black now, while in blue ($\Phi_{0}>0$) and red ($\Phi_{0}<0$) we plot those
for the RPM with $d=0.6$ nm and $d=0.4$ nm, respectively, at the same bulk
concentration $\rho_{b}$. For the three lowest bulk concentrations
$\rho_{b}=1.1$ mM, $\rho_{b}=11$ mM, and $\rho_{b}=0.11$ M, we see that the
differential capacitance of the asymmetric electrolyte with $d_{+}=0.4$ nm and
$d_{-}=0.6$ nm (black) is indistinguishable from those of the RPM with $d=0.4$
nm (red) at negative potentials and $d=0.6$ nm (blue) at positive potentials.
Hence, $C_{\mu}$ at negative potentials is fully dictated by the cations,
while at positive potentials it is fully dictated by the anions. However, at
higher bulk concentrations, close to the camel-bell transition, one finds
differences between the RPM-approximation and the actual asymmetric
electrolytes, although primarily only for $|\Phi_{0}|<100$ mV. This can be
explained by considering Eq. (13), where packing of asymmetric-sized ions in
the first layer $\sum_{j}\rho(d_{j}/2)$ is evidently different from the RPM.
Nevertheless, at higher surface potentials the differential capacitance of the
asymmetric electrolyte follows again that of the RPM, because the counterions
are fully repelled resulting in a similar packing as in the RPM.

Figure 6: (a) The differential capacitance for an asymmetric electrolyte with
$d_{+}=0.4$ nm and $d_{-}=0.6$ nm for a surface separation $H=5$ nm and
Bjerrum length $\lambda_{B}=0.73$ nm. The thick black curves in (b) correspond
to the colored ones in (a) and are compared with the differential capacitance
of the RPM with $d_{+}=d_{-}=0.4$ nm (red) and $d_{+}=d_{-}=0.6$ nm (blue) at
the same bulk concentration $\rho_{b}$.

Another type of asymmetric electrolyte that we consider is a 1:2 electrolyte
in which $z_{+}=1$ and $z_{-}=-2$, for convenience with equal diameter
$d_{+}=d_{-}=0.5$ nm. The capacitance curves for this system are given in Fig.
7, again for $H=5$ nm and $\lambda_{B}=0.73$ nm. First let us note that the
maximum at positive potentials for bulk concentrations below $\rho_{b,-}=0.85$
M is larger than that of the monovalent RPM at the same packing fraction (see
Fig. 2 or consider the negative potentials). This is a consequence of the
higher valency of the anions; fewer ions are needed to screen the charge on
the surface and therefore layering/packing requires higher surface potentials.
One peculiar feature that was not present in all previous monovalent cases is
the minimum capacitance around $\Phi_{0}\approx 0.4$ V (and therefore also a
second maximum at $\Phi_{0}=\Phi_{++}^{*}\approx 0.6$), which is caused by
overscreening [17].

Figure 7: The differential capacitance for an asymmetric 1:2 electrolyte with
$d_{+}=d_{-}=0.5$ nm. The bulk concentrations for the colored lines are given
in the legend and range from $\rho_{b,-}=\rho_{b,+}/2=0.85$ mM to
$\rho_{b,-}=\rho_{b,+}/2=3.4$ M.

At such relative high surface potentials, the positively charged surface
attracts the divalent anions to such an extent that it creates a layer of
anions with negative charge that exceeds the magnitude of the charge on the
surface. Therefore the ions further away from the electrode perceive the
surface as being negatively charged, rather than positive, which causes a net
attraction of positive ions. This is shown in Fig. 8(b) and (c), which
presents the anion and cation density profiles, respectively, at the surface
potentials indicated by the vertical lines in the capacitance curve of Fig.
8(a), for the bulk concentration of $\rho_{b,-}=\rho_{b,+}/2\approx 85$ mM
($\eta_{b}=0.1$). The inset in Fig. 8(b) shows the local packing fraction as
defined by Eq. (24). For low surface potentials $\Phi_{0}<\Phi^{*}_{+}\approx
0.14$ V (dashed-dotted lines) the behaviour is similar to the RPM, but upon
increasing the surface potential to $\Phi^{*}_{+}<\Phi_{0}=0.3$
V$<\Phi^{*}_{++}$ (dotted line) a clear peak in the cation density profile is
formed, as if the surface is negatively charged. Note that in this case, there
is still only one dense layer of anions near the surface. Increasing the
surface potential to $\Phi_{0}=0.5$ V (dashed lines) causes a dense second
layer of anions, while the cations are repelled. Overscreening is found again
at $\Phi_{0}=0.7$ V $>\Phi^{*}_{++}$ (solid lines). Hence, the structural
changes of the EDL are manifested in the differential capacitance; from a
diffuse EDL to overscreening with one layer of cations and anions (up to the
first maximum in $C_{\mu}$), to two layers of anions and a diffuse cation
layer (between the first maximum and the first minimum in $C_{\mu}$), to
overscreening with two layers of anions (beyond the second maximum in
$C_{\mu}$). This is also visible in the local packing fraction (see inset in
Fig. 8(b)), which shows that a densely-packed first layer is formed for
$\Phi_{0}>\Phi^{*}_{+}$ and a second dense layer for $\Phi_{0}>\Phi^{*}_{++}$.

Figure 8:  In (a) the $C_{\mu}$ curve of the asymmetric 1:2 electrolyte for
which the anion and cations density profiles are portrayed in (b) and (c),
respectively, at the surface potentials $\Phi_{0}=0.05$ V (dashed-dotted),
$\Phi_{0}=0.3$ V (dotted), $\Phi_{0}=0.5$ V (dashed), and $\Phi_{0}=0.7$ V
(solid), as indicated by the vertical lines in (a). The bulk concentration
corresponding to these results is $\rho_{b,+}=\rho_{b,-}/2=85$ mM
($\eta_{b}=0.01$). The inset in (b) shows the local packing fraction. The
$\Phi^{*}_{+}$ and $\Phi^{*}_{++}$ are the surface potentials at which
$C_{\mu}$ has a maximum.

###  III.3 Impurities

It is extremely difficult to experimentally study an aqueous electrolyte in
which only one type of salt is dissolved. Generally, there are always
“impurities", with a much lower bulk concentration than that of the dominant
species. Common cations in water are Calcium (Ca2+), Magnesium (Mg2+), Sodium
(Na+), and Potassium (K+), while the most common anions are Carbonate
(CO${}_{3}^{2-}$), Chloride (Cl-), and Sulfate (SO${}_{4}^{2-}$). Although
there are techniques to obtain purified (dionized, demineralized) water,
extracting all ion types is challenging [51, 52, 53]. One particular example
would be the presence of divalent ions such as Calcium (Ca2+) and Carbonate
(CO${}_{3}^{2-}$), which typically have concentrations ranging from below
$0.01$ mM (purified water) to $0.5$ mM (recommended amount of Calcium in tap
water) to $2$ mM (very hard water) [54, 55].

Figure 9:  (a) The differential capacitance for a mixture of anions, in which
one species is monovalent and the other divalent. The parameter $\nu$ gives
the ratio between the bulk density of those two anion species and ranges from
$\infty$ (one divalent anion species) to 0 (on monovalent anion species),
while the packing fraction for all curves is fixed at $\eta_{b}=0.01$. (b)
Similar as in (a), but there are two monovalent anions with different size;
one being $d_{-}=d_{+}=0.5$ nm, and the other $d^{\mathrm{large}}=2d=1$ nm.
The parameter $\nu$ is the ration between the bulk density of the large anion
w.r.t. small anion, i.e. $\nu=0$ represents the RPM, while $\nu=\infty$
represents the system with one large anion species, while the bulk packing
fraction for all curves is fixed at $\eta_{b}=0.4$.

Let us therefore now consider three-component electrolytes consisting of one
cation species and two anion species, with an asymmetry between the anions.
First, we consider a 1:1:2 electrolyte of equisized ions with diameter $d=0.5$
nm and $\rho_{b,+}=\rho_{b,-}+2\rho_{b,2-}$ where the anion composition is
characterized by $\nu=\rho_{b,2-}/\rho_{b,-}$. The differential capacitance
for this system with a bulk packing fraction $\eta_{b}=0.1$ is presented in
Fig. 9, where $\nu$ takes the values $\infty$, $10^{-1}$, $10^{-2}$,
$10^{-3}$, $10^{-4}$, $10^{-5}$, and $0$, where $\nu=0$ corresponds to a 1:1
electrolyte and $\nu=\infty$ to a 1:2 electrolyte. Interestingly, even for
$\nu=10^{-6}$ (corresponding to $\rho_{b,2-}\approx 1.3\cdot 10^{-6}$ M) we
find the divalent anion species to be predominant for surface potentials
$\Phi_{0}>0.3$ V. Again, this can be explained by using Eq. (13), which for a
1:2 electrolyte can be written in terms of the effective composition
$\tilde{\nu}\equiv\nu\mathrm{e}^{e\beta\Phi(d/2)}$ that indicates whether or
not the divalent ions dominate the surface composition, i.e. for
$\tilde{\nu}\ll 1$ the monovalent anions dominate while for $\tilde{\nu}\gg 1$
the divalent anions dominate. Because $(e\beta)^{-1}\approx 25$ mV, one finds
for $\Phi_{0}=0.3$ V $\approx 12/e\beta$ that Boltzmann factor
$\exp{(e\beta\Phi(d/2))}$ is large, such that $\mathcal{O}(\tilde{\nu})\approx
1$, even if $\nu=10^{-6}$. That is, the two anion species compete with each
other on who may fill the first layer near the surface, and since a divalent
ion has a much stronger interaction with the surface charge it wins over the
monovalent ion for large surface charges, even if their bulk concentration is
much lower. To conclude, even a trace amount of divalent ions, like Ca2+ or
CO${}_{3}^{2-}$, in a predominant monovalent electrolyte can play a dominant
role in the differential capacitance at high surface potentials.

For two species of anions which differ in size rather than valency, the story
is quite different. Let us again consider a three-component system with
monovalent ion species with different ion diameters specified by $d\equiv
d_{+}=0.5$ nm, $d_{-}^{\mathrm{small}}=d$, and $d_{-}^{\mathrm{large}}=2d$.
Fig. 9 shows the differential capacitance for this system for different anion
composition ratios
$\nu=\rho_{b,-}^{\mathrm{large}}/\rho_{b,-}^{\mathrm{small}}$ ranging from
$\nu=0$ (one anion species with $d_{-}=d$) to $\nu=\infty$ (one anion species
with $d_{-}=2d$) at a bulk packing fraction of $\eta_{b}=0.4$. Similar to the
previous three-component case with different valencies there is a competition
between the two anion species; smaller anions can approach the surface to
closer distances than the larger anions, obviously, and therefore the two
anion species again compete on who may fill the first layer. For $\nu<1$, the
smaller anions dominate at all surface potentials, and Fig. 9(b) only shows
one discernible peak in the capacitance profile. When $\nu$ is increased to
$\nu>1$ (the bulk concentration of the larger anions exceeds that of the
smaller anions) the larger anions fill the first layer up to a certain surface
potential beyond which the smaller anions take over. This point is indicated
by a peak in the differential capacitance, after which the smaller anions,
although having a much lower bulk concentration, are the dominating factor for
the differential capacitance. Again, this can be explained from Eq. (13),
where it becomes apparent that the competition is between
$\rho_{-}^{\mathrm{small}}(d/2)$ and $\rho_{-}^{\mathrm{large}}(d)$. For
$\nu>1$ and for surface potentials at which the lager anions make up the first
layer, the quantity $\rho_{-}^{\mathrm{large}}(d)$ is similar to the case of
which there is only one large cation species. However, when the surface
potential is increased such that the smaller ions make up the first layer,
then $\rho_{-}^{\mathrm{large}}(d)$ basically vanishes, due to the presence of
the smaller anions that repel the larger ones. Hence, the differential
capacitance of the impurity system then follows the curve of the differential
capacitance of the system for which there is only one anion species with
$d_{-}=d$. Looking at Fig. 9, it becomes apparent that a camel-shaped curve
can occur at very large packing fractions, although its minimum does not
coincide with the potential of zero charge.

##  IV Conclusion

We have carefully investigated the differential capacitance for primitive
model electrolytes within classical DFT. A new equation (13) to interpret the
differential capacitance was introduced, which greatly helped to understand
its behaviour, because it explicitly showed the importance of the first layer
of ions near the surface. Specifically, it helped to further rationalize the
camel-bell crossover, that has been a topic of interest across many studies,
as a (smooth) structural change in the density profiles near the electrode.
Analogous to peaks in the heat capacity indicating (smooth) changes of the
thermal occupancy of microstates, peaks in the differential capacitance
indicate structural change. This can either be the forming of dense-packed
layers of counterions near the surface, but can also be the reorganisation of
layers of cations and anions as we showed for the 1:2 electrolyte with
divalent anions, in which overscreening causes a rich structural behaviour
that seeped through in the capacitance curves. When the impurities in three-
component electrolytes were studied, we again found that the competition for
the first layer of ions near the surface determined the behaviour of the
capacitance curves. This three-component electrolyte that was studied, with
one species of cations and two species of anions in which the relative
concentration of the anions was changed, gave much richer physics than perhaps
anticipated. Both an asymmetry in the anion valency as well as in the anion
diameter was considered. Surprisingly at first, even when the composition of
the bulk mixture is very asymmetric, both anion species still have their
regime in which they fill up the first layer. For anion asymmetry in the
valency, where the relative concentration of the divalent anions is very
small, we found that the divalent anions can still fill up the first layer
when the surface potential is large enough, due to the much stronger
interaction with the charged surface. For diameter asymmetry of the anions,
with a low relative concentration of the small anions, we found that they can
still fill up the first layer at large surface potential because they can
screen the charge on the surface more efficiently. For both cases it is the
competition between the two anion species that shapes the differential
capacitance. Hence, impurities, as we call them, have a rather strong effect
on the differential capacitance. The differential capacitance can thus be used
as a probe for the electrolyte composition. Knowing that the differential
capacitance is largely determined by the first layer of adsorbed ions, and
that peaks in the differential capacitance indicate structural changes, can
therefore help to distinguish the components in the electrolyte.

Although Eq. (13) for the differential capacitance that we derived strictly
only holds for the primitive model electrolytes in contact with charged hard
walls, the physics that the first layer of ions near the surface dominate its
behaviour still holds when the conditions are loosened up. Important is that
the ion-ion and ion-surface interaction are strongly repulsive at short
separations; our results will be modified if dispersion forces are relevant.

The next step would be to confront the knowledge gained from this analysis to
experimental data, allowing interpretation of the measured differential
capacitance in terms of EDL structure and bulk electrolyte composition. This
could, in turn, lead to a better understanding of the differential
capacitance, and consequently to an improved performance of real EDL
capacitors. An ultimate goal is to provide knowledge to create sustainable
alternatives for Lithium-ion batteries [56, 57], which would contribute
directly to the transition towards renewable energy.

###### Acknowledgements.

We thank Willem Boon for engaging discussions. This work is part of the D-ITP
consortium, a program of the Netherlands Organisation for Scientific Research
(NWO) that is funded by the Dutch Ministry of Education, Culture and Science
(OCW). It also forms part of the NWO programme ‘Datadriven science for smart
and sustainable energy research’, with project number 16DDS014.

##  Appendix A Mean Spherical Approximation

The electrostatic Helmholtz free energy functional is given by

$\displaystyle\mathcal{F}_{ex}^{ES}[\\{\rho\\}]=$
$\displaystyle\mathcal{F}_{ex}^{MF}[\\{\rho\\}]+\mathcal{F}_{ex}^{MSAc}[\\{\rho\\}];$
(25) $\displaystyle\mathcal{F}_{ex}^{MF}[\\{\rho\\}]=$
$\displaystyle\frac{1}{2}\int\mathrm{d}\mathbf{r}q(\mathbf{r})\Phi(\mathbf{r});$
(26) $\displaystyle\mathcal{F}_{ex}^{MSAc}[\\{\rho\\}]=$
$\displaystyle-\frac{k_{B}T}{2}\sum_{ij}\int\mathrm{d}\mathbf{r}\int\mathrm{d}\mathbf{r}^{\prime}\rho_{i}(\mathbf{r})\times$
$\displaystyle\quad\quad\quad\Delta
c^{MSAc}_{ij}(|\mathbf{r}-\mathbf{r}^{\prime}|;\\{\rho_{b}\\})\rho_{j}(\mathbf{r}^{\prime}),$
(27)

where $\Delta c^{MSAc}(\mathbf{r},\\{\rho_{b}\\})$ is a correction on top of
mean-field theory based on the mean spherical approximation (MSA). The
potential profile $\Phi(\mathbf{r})$ is obtained from solving the Poisson
equation

$\displaystyle\nabla^{2}\Phi(\mathbf{r})=-\frac{q_{tot}(\mathbf{r})}{\epsilon\epsilon_{0}},$
(28)

where $q_{tot}(\mathbf{r})=q(\mathbf{r})+q_{f}(\mathbf{r})$ is the total
charge density counting both the ionic charge density
$q(\mathbf{r})=\sum_{j}z_{j}\rho_{j}(\mathbf{r})$ and the fixed charges
$q_{f}(\mathbf{r})$. Solving Eq. (28) for fixed surface potentials in a planar
geometry requires the boundary conditions $\Phi(0)=\Phi(H)=\Phi_{0}$, which
give rise to the surface charge densities via
$\partial_{x}\Phi(0+)=-\partial_{x}\Phi(H^{-})=-4\pi\lambda_{B}\sigma$.
Therefore, the electrostatic part of the external potential arises from the
Poisson equation and is taken care of within
$\mathcal{F}_{ex}^{MF}[\\{\rho\\}]$. In the bulk, $\Delta c^{MSAc}(r)$ is
given by [58]

$\displaystyle\Delta
c_{ij}^{MSAc}(r,\\{\rho_{b}\\})=\begin{cases}c_{ij}^{MSAsh}(r,\\{\rho_{b}\\})+z_{i}z_{j}\frac{\lambda_{B}}{r},\quad&0\leq
r\leq\Delta d_{ij};\\\
c_{ij}^{MSAl}(r,\\{\rho_{b}\\})+z_{i}z_{j}\frac{\lambda_{B}}{r},\quad&\Delta
d_{ij}<r\leq d_{ij};\\\ 0,\quad&r>d_{ij};\end{cases}$ (29)

where $\Delta d_{ij}=|d_{i}-d_{j}|/2$ and $d_{ij}=(d_{i}+d_{j})/2$. For
$d_{i}<d_{j}$ the first term for short (sh) separations reads

$\displaystyle
c_{ij}^{MSAsh}(r;\\{\rho_{b}\\})=2\lambda_{B}\left[z_{i}N_{j}+d_{i}\zeta(X_{i}+\frac{1}{3}d_{i}^{2}\zeta)\right],$
(30)

with

$\displaystyle X_{j}=$ $\displaystyle\frac{z_{j}-d_{j}^{2}\zeta}{1+\Upsilon
d_{j}};$ (31) $\displaystyle N_{j}=$ $\displaystyle\frac{X_{j}-z_{j}}{d_{j}};$
(32) $\displaystyle\Upsilon=$
$\displaystyle\pi\lambda_{B}\sum_{j}\rho_{b,j}X_{j}^{2};$ (33)
$\displaystyle\zeta=$
$\displaystyle\frac{1}{H}\sum_{j}\frac{\rho_{b,j}d_{j}z_{j}}{1+\Upsilon
d_{j}};$ (34) $\displaystyle H=$
$\displaystyle\sum_{j}\frac{\rho_{b,j}d_{j}^{3}}{1+\Upsilon
d_{j}}+\frac{2}{\pi}\left(1-\eta_{b}\right).$ (35)

The second term in Eq. (29) for the longer-ranged part (l) for $\Delta
d_{ij}<r\leq d_{ij}$ is given by

$\displaystyle
c_{ij}^{MSAl}(r;\\{\rho_{b}\\})=\frac{\lambda_{B}}{r}\left[A_{ij}+B_{ij}r+C_{ij}r^{2}+F_{ij}r^{4}\right];$
(36)

with

$\displaystyle A_{ij}=$ $\displaystyle-\Delta
d_{ij}^{2}\left[\zeta(X_{i}+X_{j})+\zeta^{2}d_{ij}^{2}-N_{i}N_{j}\right],$
(37) $\displaystyle B_{ij}=$
$\displaystyle-(X_{i}-X_{j})(N_{i}-N_{j})-(X_{i}^{2}+X_{j}^{2})\Upsilon-$
$\displaystyle 2d_{ij}N_{i}N_{j}+\frac{1}{3}\zeta^{2}(d_{i}^{3}+d_{j}^{3}),$
(38) $\displaystyle C_{ij}=$
$\displaystyle-\zeta(X_{i}+X_{j})+N_{i}N_{j}-\frac{1}{2}\zeta^{2}(d_{i}^{2}+d_{j}^{2}),$
(39) $\displaystyle F_{ij}=$ $\displaystyle\frac{1}{3}\zeta^{2}.$ (40)

These rather involved expressions become more tractable for the monovalent
RPM, where $c_{ij}^{MSAsh}(r;\rho_{b})=0$, $\zeta=0$ and $\Delta d_{ij}=0$,
such that

$\displaystyle c_{ij}^{MSAl}(r;\rho_{b})=$
$\displaystyle\frac{\lambda_{B}}{r}\Big{[}\Big{(}-(X_{i}-X_{j})(N_{i}-N_{j})-2X^{2}\Upsilon$
$\displaystyle-2dN_{i}N_{j}\Big{)}r+N_{i}N_{j}r^{2}\Big{]},$ (41)

which can be simplified by introducing $D=d+1/\Upsilon$, so that with

$\displaystyle X_{j}$ $\displaystyle=\frac{z_{j}}{\Upsilon D},$ (42)
$\displaystyle N_{j}$ $\displaystyle=-\frac{z_{j}}{D},$ (43)

one finds

$\displaystyle
c_{ij}^{MSAl}(r;\rho_{b})=z_{i}z_{j}\frac{\lambda_{B}}{r}\frac{-2Dr+r^{2}}{D^{2}}.$
(44)

For a planar geometry, one integrates out the planar coordinates:
$c(x)=\int\mathrm{d}y\mathrm{d}zc(\sqrt{x^{2}+y^{2}+z^{2}})$.

## References

  * Gielen  _et al._ [2019] D. Gielen, F. Boshell, D. Saygin, M. D. Bazilian, N. Wagner, and R. Gorini, “The role of renewable energy in the global energy transformation,” Energy Strategy Reviews 24, 38–50 (2019).
  * Graabak and Korpås [2016] I. Graabak and M. Korpås, “Variability characteristics of european wind and solar power resources—a review,” Energies 9, 449 (2016).
  * Huber, Dimkova, and Hamacher [2014] M. Huber, D. Dimkova, and T. Hamacher, “Integration of wind and solar power in europe: Assessment of flexibility requirements,” Energy 69, 236 – 246 (2014).
  * Gür [2018] T. M. Gür, “Review of electrical energy storage technologies, materials and systems: challenges and prospects for large-scale grid storage,” Energy Environ. Sci. 11, 2696–2767 (2018).
  * Lai and Locatelli [2021] C. S. Lai and G. Locatelli, “Economic and financial appraisal of novel large-scale energy storage technologies,” Energy 214, 118954 (2021).
  * Matos, Carneiro, and Silva [2019] C. R. Matos, J. F. Carneiro, and P. P. Silva, “Overview of large-scale underground energy storage technologies for integration of renewable energies and criteria for reservoir identification,” Journal of Energy Storage 21, 241 – 258 (2019).
  * Posada  _et al._ [2017] J. O. G. Posada, A. J. Rennie, S. P. Villar, V. L. Martins, J. Marinaccio, A. Barnes, C. F. Glover, D. A. Worsley, and P. J. Hall, “Aqueous batteries as grid scale energy storage solutions,” Renewable and Sustainable Energy Reviews 68, 1174 – 1182 (2017).
  * Bullich-Massagué  _et al._ [2020] E. Bullich-Massagué, F.-J. Cifuentes-García, I. Glenny-Crende, M. Cheah-Mañé, M. Aragüés-Peñalba, F. Díaz-González, and O. Gomis-Bellmunt, “A review of energy storage technologies for large scale photovoltaic power plants,” Applied Energy 274, 115213 (2020).
  * Górka  _et al._ [2010] J. Górka, A. Zawislak, J. Choma, and M. Jaroniec, “Adsorption and structural properties of soft-templated mesoporous carbons obtained by carbonization at different temperatures and koh activation,” Applied Surface Science 256, 5187–5190 (2010), seventh International Symposium Effects of Surface Heterogeneity in Adsorption and Catalysis on Solids - ISSHAC-7.
  * Cats  _et al._ [2021a] P. Cats, R. Sitlapersad, W. den Otter, A. Thornton, and R. van Roij, “Capacitance and structure of electric double layers: Comparing brownian dynamics and classical density functional theory,” J. Solution Chem. (2021a).
  * Bossa  _et al._ [2019] G. V. Bossa, D. L. Caetano, S. J. de Carvalho, and S. May, “Differential capacitance of an electrical double layer with asymmetric ion sizes in the presence of hydration interactions,” Electrochimica Acta 321, 134655 (2019).
  * Kaja  _et al._ [2016] M. Kaja, S. Lamperski, W. Silvestre-Alcantara, L. Bhuiyan, and D. Henderson, “Influence of anisotropic ion shape, asymmetric valency, and electrolyte concentration on structural and thermodynamic properties of an electric double layer,” Condensed Matter Physics 19, 13804 (2016).
  * Henderson  _et al._ [2017] D. Henderson, W. Silvestre-Alcantara, M. Kaja, S. Lamperski, J. Wu, and L. B. Bhuiyan, “Structure and capacitance of an electric double layer of an asymmetric valency dimer electrolyte: A comparison of the density functional theory with monte carlo simulations,” Journal of Molecular Liquids 228, 236–242 (2017), from simple liquids to macromolecular solutions: recent experimental and theoretical developments. In Honor of the 70th birthday of Vojko Vlachy.
  * Biesheuvel and van Soestbergen [2007] P. Biesheuvel and M. van Soestbergen, “Counterion volume effects in mixed electrical double layers,” Journal of Colloid and Interface Science 316, 490–499 (2007).
  * Caetano  _et al._ [2016] D. L. Z. Caetano, G. V. Bossa, V. M. de Oliveira, M. A. Brown, S. J. de Carvalho, and S. May, “Role of ion hydration for the differential capacitance of an electric double layer,” Phys. Chem. Chem. Phys. 18, 27796–27807 (2016).
  * Nakayama and Andelman [2015] Y. Nakayama and D. Andelman, “Differential capacitance of the electric double layer: The interplay between ion finite size and dielectric decrement,” The Journal of Chemical Physics 142, 044706 (2015).
  * Kornyshev [2007] A. A. Kornyshev, “Double-layer in ionic liquids: Paradigm change?” The Journal of Physical Chemistry B 111, 5545–5557 (2007).
  * Frischknecht, Halligan, and Parks [2014] A. L. Frischknecht, D. O. Halligan, and M. L. Parks, “Electrical double layers and differential capacitance in molten salts from density functional theory,” The Journal of Chemical Physics 141, 054708 (2014).
  * Ma, Woodward, and Forsman [2014] K. Ma, C. E. Woodward, and J. Forsman, “Classical density functional study on interfacial structure and differential capacitance of ionic liquids near charged surfaces,” J. Phys. Chem. C 118, 15825–15834 (2014).
  * Vo, Forsman, and Woodward [2020] P. Vo, J. Forsman, and C. E. Woodward, “A semi-gcmc simulation study of electrolytic capacitors with adsorbed titrating peptides,” The Journal of Chemical Physics 153, 174703 (2020).
  * Chen, Paul, and Dai [2017] X. Chen, R. Paul, and L. Dai, “Carbon-based supercapacitors for efficient energy storage,” National Science Review 4, 453–489 (2017).
  * Zhai  _et al._ [2011] Y. Zhai, Y. Dou, D. Zhao, P. F. Fulvio, R. T. Mayes, and S. Dai, “Carbon materials for chemical capacitive energy storage,” Advanced Materials 23, 4828–4850 (2011).
  * Winter and Brodd [2004] M. Winter and R. J. Brodd, “What are batteries, fuel cells, and supercapacitors?” Chemical Reviews 104, 4245–4270 (2004), pMID: 15669155\.
  * Orkoulas and Panagiotopoulos [1994] G. Orkoulas and A. Z. Panagiotopoulos, “Free energy and phase equilibria for the restricted primitive model of ionic fluids from monte carlo simulations,” The Journal of Chemical Physics 101, 1452–1459 (1994).
  * van Roij [2014] R. van Roij, “Statistical thermodynamics of supercapacitors and blue engines,” in  _Electrostatics of Soft and Disordered Matter_ , edited by D. S. Dean, J. Dobnikar, A. Naji, and R. Podgornik (Jenny Stanford Publishing, 2014).
  * Henderson and Blum [1978] D. Henderson and L. Blum, “Some exact results and the application of the mean spherical approximation to charged hard spheres near a charged hard wall,” The Journal of Chemical Physics 69, 5441–5449 (1978).
  * Henderson, Blum, and Lebowitz [1979] D. Henderson, L. Blum, and J. L. Lebowitz, “An exact formula for the contact value of the density profile of a system of charged hard spheres near a charged wall,” Journal of Electroanalytical Chemistry and Interfacial Electrochemistry 102, 315–319 (1979).
  * Eckert  _et al._ [2020] T. Eckert, N. C. X. Stuhlmüller, F. Sammüller, and M. Schmidt, “Fluctuation profiles in inhomogeneous fluids,” Phys. Rev. Lett. 125, 268004 (2020).
  * Stern [1924] O. Stern, “The theory of the electrolytic double layer,” Z. Elektrochem. Angew. Phys. Chem. 30, 508–516 (1924).
  * van Hal, Eijkel, and Bergveld [1996] R. van Hal, J. Eijkel, and P. Bergveld, “A general model to describe the electrostatic potential at electrolyte oxide interfaces,” Advances in Colloid and Interface Science 69, 31–62 (1996).
  * Lockett  _et al._ [2010] V. Lockett, M. Horne, R. Sedev, T. Rodopoulos, and J. Ralston, “Differential capacitance of the double layer at the electrode/ionic liquids interface,” Phys. Chem. Chem. Phys. 12, 12499–12512 (2010).
  * May [2019] S. May, “Differential capacitance of the electric double layer: mean-field modeling approaches,” Current Opinion in Electrochemistry 13, 125–131 (2019), fundamental and Theoretical Electrochemistry Physical and Nanoelectrochemistry.
  * Note [1] The Stern layer $C_{EDL}$ separates the diffuse layer from the electrode, each layer having their own capacitance, which in a serial circuit gives the total capacitance of the EDL:  $\displaystyle\frac{1}{C_{EDL}}=\frac{1}{C_{s}}+\frac{1}{C_{DL}}.$ (45) With little effort, one can show using Eqs. (17) and (45) that, by introducing the Stern-layer capacitance, the capacitance of the diffuse layer reads  $\displaystyle C_{DL}=C_{s}\left[\left(\frac{\partial\Phi_{0}}{\partial\Phi(d/2)}\right)_{\mu}-1\right],$ (46) which can take negative values in case of overscreening. Hence, splitting the EDL into a diffuse layer and the stern layer seems not to be desirable for a generic system, and we will not further distinguish between Stern layers and other layers, but instead talk about the total capacitance $C_{\mu}$.
  * Evans [1979] R. Evans, “The nature of the liquid-vapour interface and other topics in the statistical mechanics of non-uniform, classical fluids,” Advances in Physics 28, 143–200 (1979).
  * Roth [2010] R. Roth, “Fundamental measure theory for hard-sphere mixtures: a review,” Journal of Physics: Condensed Matter 22, 063102 (2010).
  * Härtel [2017] A. Härtel, “Structure of electric double layers in capacitive systems and to what extent (classical) density functional theory describes it,” Journal of Physics: Condensed Matter 29, 423002 (2017).
  * Mier-y-Teran  _et al._ [1990] L. Mier-y-Teran, S. H. Suh, H. S. White, and H. T. Davis, “A nonlocal free-energy density-functional approximation for the electrical double layer,” The Journal of Chemical Physics 92, 5087–5098 (1990).
  * Cats  _et al._ [2021b] P. Cats, R. Evans, A. Härtel, and R. van Roij, “Primitive model electrolytes in the near and far field: Decay lengths from dft and simulations,” The Journal of Chemical Physics 154, 124504 (2021b).
  * Limmer  _et al._ [2013] D. T. Limmer, C. Merlet, M. Salanne, D. Chandler, P. A. Madden, R. van Roij, and B. Rotenberg, “Charge fluctuations in nanoscale capacitors,” Phys. Rev. Lett. 111, 106102 (2013).
  * Scalfi  _et al._ [2020] L. Scalfi, D. T. Limmer, A. Coretti, S. Bonella, P. A. Madden, M. Salanne, and B. Rotenberg, “Charge fluctuations from molecular simulations in the constant-potential ensemble,” Phys. Chem. Chem. Phys. 22, 10480–10489 (2020).
  * Fedorov, Georgi, and Kornyshev [2010] M. Fedorov, N. Georgi, and A. Kornyshev, “Double layer in ionic liquids: The nature of the camel shape of capacitance,” Electrochemistry Communications 12, 296–299 (2010).
  * Fedorov and Kornyshev [2014] M. V. Fedorov and A. A. Kornyshev, “Ionic liquids at electrified interfaces,” Chemical Reviews 114, 2978–3036 (2014), pMID: 24588221.
  * Uematsu, Netz, and Bonthuis [2018] Y. Uematsu, R. R. Netz, and D. J. Bonthuis, “The effects of ion adsorption on the potential of zero charge and the differential capacitance of charged aqueous interfaces,” Journal of Physics: Condensed Matter 30, 064002 (2018).
  * Girotto  _et al._ [2018] M. Girotto, R. M. Malossi, A. P. dos Santos, and Y. Levin, “Lattice model of ionic liquid confined by metal electrodes,” The Journal of Chemical Physics 148, 193829 (2018).
  * Downing, Bossa, and May [2018] R. Downing, G. V. Bossa, and S. May, “The role of ion–ion correlations for the differential capacitance of ionic liquids,” The Journal of Physical Chemistry C 122, 28537–28544 (2018).
  * Islam, Alam, and Ohsaka [2008] M. M. Islam, M. T. Alam, and T. Ohsaka, “Electrical double-layer structure in ionic liquids: A corroboration of the theoretical model by experimental results,” The Journal of Physical Chemistry C 112, 16568–16574 (2008).
  * Jiang, Meng, and Wu [2011] D. Jiang, D. Meng, and J. Wu, “Density functional theory for differential capacitance of planar electric double layers in ionic liquids,” Chemical Physics Letters 504, 153–158 (2011).
  * Jitvisate and Seddon [2018] M. Jitvisate and J. R. T. Seddon, “Direct measurement of the differential capacitance of solvent-free and dilute ionic liquids,” The Journal of Physical Chemistry Letters 9, 126–131 (2018), pMID: 29256620.
  * Note [2] $\eta(x)$ is identical to the weighted density $n_{3}(x)$ from Fundamental Measure Theory [35].
  * Dijkstra [2004] M. Dijkstra, “Capillary freezing or complete wetting of hard spheres in a planar hard slit?” Phys. Rev. Lett. 93, 108303 (2004).
  * Bassyouni  _et al._ [2019] M. Bassyouni, M. Abdel-Aziz, M. S. Zoromba, S. Abdel-Hamid, and E. Drioli, “A review of polymeric nanocomposite membranes for water purification,” Journal of Industrial and Engineering Chemistry 73, 19–46 (2019).
  * Gamaethiralalage  _et al._ [2021] J. G. Gamaethiralalage, K. Singh, S. Sahin, J. Yoon, M. Elimelech, M. E. Suss, P. Liang, P. M. Biesheuvel, R. L. Zornitta, and L. C. P. M. de Smet, “Recent advances in ion selectivity with capacitive deionization,” Energy Environ. Sci. 14, 1095–1120 (2021).
  * Yang  _et al._ [2021] F. Yang, Y. He, L. Rosentsvit, M. E. Suss, X. Zhang, T. Gao, and P. Liang, “Flow-electrode capacitive deionization: A review and new perspectives,” Water Research 200, 117222 (2021).
  * team [2009] W. team,  _Calcium and magnesium in drinking-water_ (World Health Organisation, 2009).
  * Kozisek [2020] F. Kozisek, “Regulations for calcium, magnesium or hardness in drinking water in the european union member states,” Regulatory Toxicology and Pharmacology 112, 104589 (2020).
  * Kaunda [2020] R. B. Kaunda, “Potential environmental impacts of lithium mining,” Journal of Energy & Natural Resources Law 38, 237–244 (2020).
  * Izquierdo  _et al._ [2015] A. Izquierdo, R. Grau, J. Carilla, and E. Casagranda, “Side effects of green technologies: the potential environmental costs of lithium mining on high elevation andean wetlands in the context of climate change,” Newsletter of the Global Land Project 12, 53–56 (2015).
  * Hiroike [1977] K. Hiroike, “Supplement to blum’s theory for asymmetric electrolytes,” Molecular Physics 33, 1195–1198 (1977).

## DATA AVAILABILITY

The data that support the findings of this study are available from the
corresponding authors upon reasonable request.