supporting_info

# The Role of Energy Scales for the Structure of Ionic Liquids at Electrified
Interfaces - A Theory-Based Approach

Max Schammer  German Aerospace Center, Pfaffenwaldring 38-40, 70569 Stuttgart,
Germany  Helmholtz Institute Ulm, Helmholtzstraße 11, 89081 Ulm, Germany
Arnulf Latz<EMAIL_ADDRESS>German Aerospace Center, Pfaffenwaldring
38-40, 70569 Stuttgart, Germany  Helmholtz Institute Ulm, Helmholtzstraße 11,
89081 Ulm, Germany  Universität Ulm, Albert-Einstein-Allee 47, 89081 Ulm,
Germany  Birger Horstmann<EMAIL_ADDRESS>German Aerospace Center,
Pfaffenwaldring 38-40, 70569 Stuttgart, Germany  Helmholtz Institute Ulm,
Helmholtzstraße 11, 89081 Ulm, Germany  Universität Ulm, Albert-Einstein-Allee
47, 89081 Ulm, Germany

###### Abstract

Ionic liquids offer unique bulk and interfacial characteristics as battery
electrolytes. Our continuum approach naturally describes the electrolyte on a
macroscale. An integral formulation for the molecular repulsion, which can be
quantitatively determined by both experimental and theoretical methods, models
the electrolyte on the nanoscale. In this article, we perform a systematic
series expansion of this integral formulation, derive a description of
chemical potentials in terms of higher-order concentration gradients, and
rationalize the appearance of fourth-order derivative-operators in modified
Poisson equations, recently proposed in this context. In this way, we
formulate a rigorous multi-scale methodology from atomistic quantum chemistry
calculations to phenomenological continuum models. We apply our generalized
framework to ionic liquids near electrified interfaces and perform analytic
asymptotic analysis. Three energy scales describing electrostatic forces
between ions, molecular repulsion, and thermal motion determine the shape and
width of the long-ranging charged double layer. We classify the charge
screening mechanisms dependent on the system parameters dielectricity, ion
size, interaction strength, and temperature. We find that the charge density
of electrochemical double layers in ionic liquids either decays exponentially,
for negligible molecular repulsion, or oscillates continuously. Charge
ordering across several ion-diameters occurs if the repulsion between
molecules is comparable with thermal energy and Coulomb interaction.
Eventually, phase separation of the bulk electrolyte into ionic layers emerges
once the molecular repulsion becomes dominant. Our framework predicts the
exact phase boundaries among these three phases as a function of temperature,
dielectricity and ion-sizes.

##  I Introduction

Strong electrostatic correlations in crowded environments play an important
role in biology, chemistry and physics.1, 2, 3 For example, in molecular
biology, they account for DNA packing,4 which is crucial for the
compactification of genetic materials in viruses,5 impact the cytoskeleton
organization,6 and influence transport in ion channels.7 Furthermore, such
correlations explain the thermodynamic stability of plasmas,8, 9 and charged
colloidal suspensions.10, 11

Surprisingly, the complexity of these phenomena can be understood to a large
degree by models derived initially for electrolyte solutions.12 Starting from
the fundamental Debye-Hückel theory for dilute solutions,13 increasingly
accurate models for concentrated electrolytes were developed,14 taking more
complex Coulomb correlations into account.

Because ionic liquids (ILs) consist only of positive and negative ions without
neutral solvent, they constitute the extreme limit for the examination of
electrostatic correlations in electrolytic solutions. Indeed, ILs possess
characteristic properties in the bulk-regime,15, 16 but also near electrified
interfaces.17 This makes them highly attractive from both fundamental and
applied perspectives.18, 19, 20, 21, 22, 23, 24 The study of interfacial
electrochemistry is of wide-ranging interest. For example, the behavior of ILs
near electrified interfaces has paramount importance for their performance as
battery electrolytes.25, 26

Theoretical studies of ILs near electrified interfaces discuss the structure
of charged electrochemical double layers (EDL) on atomistic/molecular scales.
These include classical density functional theory (cDFT) simulations, and
molecular-dynamics (MD) simulations. cDFT gives detailed insights into the
arrangement of molecules in the EDL.27, 28, 29 MD resolves the molecular
motion and can elucidate the EDL structure.30, 31, 32

However, cDFT/MD simulations are limited by their computational costs.
Simulations at length-scales above the nanometer scale are hardly accessible
to the atomistic / molecular approach. Thus, continuum theories, and mean-
field-theories (MFT) provide a complementary methodology for the simulation of
larger systems, where the microscopic details can be neglected, or are used as
averaged parameters (e.g. constant dielectric parameters).

Usually, MFTs for electrolytes are based on lattice gas models of ions, first
proposed by Bikermann.33 Recently, MFTs have attracted great interest fot the
study for ILs. As proposed by Santangelo for aqueous systems,34 the extension
of MFTs by higher-order electrostatic correlations is useful for the
description of long-ranged structures emerging in electrolytes. Bazant, Storey
and Kornyshev (BSK) applied this approach to ILs near electrified
interfaces.35 By using a phenomenological model, which is based on a
generalized Ginzburg-Landau functional, BSK describe charge oscillations known
as overscreening and charge saturation known as crowding. Yochelis et al.
rationalized this approach, and extended it to bulk properties.36, 37, 38, 38,
39, 40 However, MFT models are usually restricted to equilibrium effects of
binary ILs with structureless bulk, although rare MFT models, complemented by
continuum methods 41, 39 and extended to the ternary case,42 exist.

This highlights the advantage of continuum frameworks, which describe
dynamical transport processes. In addition, continuum models based on rigorous
physical assumptions43 identify coupled phenomena arising from the interplay
of mechanics, thermodynamics, and electromagnetic theory.44 Furthermore, this
approach allows to develop a unified, thermodynamically consistent framework
that provides the common theoretical basis for the description of different
electrochemical systems.45, 46, 47, 48, 49 Continuum models are not restricted
to binary or ternary systems, because they can be formulated for arbitrary
many species, charged and uncharged. Thus, they apply to more realistic
electrolytes.

Recently, we proposed such a novel continuum transport theory.45 In this
theory, we take account for steric effects via the mean-volume, which is due
to finite molar volumes of the ion species. For this purpose, we impose a
volume constraint on the electrolyte. This mechanism stabilizes the bulk
structure against Coulomb collapse50 and leads to charge-saturation near
electrified interfaces. Thus our theory resolves the deficiencies of the
classical Poisson-Boltzmann (PB) theory, which predicts unrealistically high
interface-concentrations.14 Furthermore, the existence of finite molar volumes
of the electrolyte species leads to a pressure dependence of the chemical
potentials (in accordance with thermodynamic arguments).51

However, this bulk-framework cannot describe the emergence of long-range
structures in ILs near electrified interfaces. Therefore, in a joint
experimental / theoretical work, we extended our framework with non-local
interactions and validated it with results obtained from atomic-force-
microscopy.52 Thus, we extended the mean-volume-effect of the bulk theory with
molecular volume exclusion due to hardcore repulsion. Our holistic framework
allows us to couple dynamic transport processes occurring in the bulk-
electrolyte with interfacial electrochemical processes. Thus, we provide a
continuum model which bridges the length-scales from nano-meters, e.g., EDL,
to millimeters, e.g., battery cells. Moreover, our framework allows us to
connect the continuum description with correlation functions generated by MD.

However, the dependence of EDL structures on molecular repulsion, molecular
size, temperature, and dielectricity is still unknown. In this paper, we
derive such an understanding with asymptotic analysis. To this aim, we present
our thermodynamically consistent transport theory with an integral formulation
of non-local interactions in section II.1. These correlations represent
atomistic volume exclusion and lead to modified constitutive equations
(LABEL:eq:SI_const_eq_entro, LABEL:eq:SI_const_eq_elf,
LABEL:eq:SI_const_eq_magf and LABEL:eq:SI_const_eq_chempot). Moreover, the
interactions impose contributions to the stress-tensor, and thus modify the
mechanical coupling to the transport equations
(LABEL:eq:SI_const_eq_stress_tensor). In section II.2, we approximate the
interaction functional with a gradient expansion, which facilitates the
analytic asymptotic analysis of the EDL structure. In section II.3, we apply
our extended framework to study the EDL structure for neat ILs. When we non-
dimensionalize our dynamical description in section II.4, three competing
energy scales describing electrostatic forces among ions, molecular repulsion,
and thermal motion appear in the theory. Because our focus lies on the
formation of equilibrium structures, we discuss the stationary state in
section II.5. In section II.6, we discuss limiting cases of our stationary
theory.

We perform numerical simulations and analytic asymptotic analysis to study the
interplay and the effect of the competing energy scales on EDL structures.
First, in section III, we discuss the EDL structure for the mean volume
constraint. Second, in section IV, we incorporate molecular repulsion into our
analysis and classify the EDL structure dependence on the relation between
competing energy scales.

##  II Theory

###  II.1 Generalized Transport Theory

Recently, we have proposed a free energy functional
$F^{\ce{b}}=\int\ce{d}V\,\uprho\upvarphi_{\textnormal{H}}$ for the dynamical
description of ionic liquids in the bulk-phase.45 In this bulk-model, the
Helmholtz free energy density
$\upvarphi_{\textnormal{H}}(\varUpsilon)=\upvarphi_{\textnormal{H}}(T,c_{1},\ldots,c_{\textnormal{N}},\boldsymbol{D},\boldsymbol{B},\boldsymbol{\upkappa})$
is a function of the variables temperature $T$, concentrations $c_{\alpha}$,
dielectric displacement $\boldsymbol{D}$, magnetic field $\boldsymbol{B}$, and
strain-rate tensor $\boldsymbol{\upkappa}$. This variable-set $\varUpsilon$
constitutes material-specific properties of multi-component, viscous and
polarizable media in the liquid state.

In contrast, models describing non-local interactions rely on functionals
$F^{\ce{int}}[\varUpsilon]$, such that the free energy takes the form

$F[\varUpsilon]=F^{\ce{int}}[\varUpsilon]+F^{\ce{b}}(\varUpsilon)=F^{\ce{int}}[\varUpsilon]+\int\ce{dV}\,\uprho\upvarphi_{\textnormal{H}}.$
(1)

This functional approach constitutes a more general description for
electrolyte materials and allows the incorporation of non-local correlations
between field quantities. Such correlations typically arise from microscopic
effects occurring on the nano-meter scales, e.g., in the vicinity of
electrified interfaces. Despite the conceptual difference between the
functional approach eq. 1 and the canonical bulk approach, the derivation of
the resulting transport theory is rather similar to the rationale outlined in
great detail in Ref. 45 for the free energy
$F^{\ce{b}}=\int\ce{dV}\uprho\upvarphi_{\textnormal{H}}$. We present a
detailed derivation for the functional approach in the Supporting Information
(see LABEL:sec:SItheory).

The extension of the free energy according to eq. 1 leads to modified
constitutive equations for entropy density $s$, electric field strength
$\boldsymbol{\mathscr{E}}$, magnetic field $\boldsymbol{\mathscr{H}}$,
chemical potentials $\upmu_{\alpha}$, and the stress tensor
$\boldsymbol{\upsigma}$ in the form of functional-derivatives (see
LABEL:eq:SI_const_eq_entro, LABEL:eq:SI_const_eq_elf,
LABEL:eq:SI_const_eq_magf, LABEL:eq:SI_const_eq_chempot and
LABEL:eq:SI_const_eq_stress_tensor in LABEL:sec:SItransport-theory). We
evaluate this framework for the bulk-energy of a linear dielectric medium
discussed in Ref. 45, see LABEL:eq:SI_model_bulk_energy. The resulting forces
are supplemented by contributions stemming from the non-local correlations
(LABEL:eq:SI_chemical_force).

For the remaining part of this work, we neglect thermal driving forces by
setting the temperature equal to constant values, and assume the electrostatic
limit, $\boldsymbol{B}=0$ and $\boldsymbol{\mathscr{H}}=0$. This determines
the electric field $\boldsymbol{E}=\boldsymbol{\mathscr{E}}$ by the
electrostatic potential, $\boldsymbol{E}=-\nabla\varPhi$.

###  II.2 Gradient Expansion of Molecular Interactions

We proposed a model for hardcore-interactions based on a convolution-
functional for the interaction free energy in Ref. 52,

$F^{\ce{int}}=\frac{1}{2}\sum_{\alpha,\,\beta}^{\textnormal{N}}\iint\ce{d}x^{3}\ce{d}y^{3}\
\mathscr{F}_{\alpha\beta}(|\mathbf{x}-\mathbf{y}|)c_{\alpha}(\mathbf{x})c_{\beta}(\mathbf{y}),$
(2)

leading to transport contributions in the form of (see LABEL:sec:SIfunct-
deriv)

$\frac{\updelta F^{\ce{int}}}{\updelta
c_{\alpha}}(\mathbf{x})=\sum_{\beta}^{\textnormal{N}}\int\ce{d}y^{3}\
\mathscr{F}_{\alpha\beta}(\lvert\mathbf{x}-\mathbf{y}\rvert)c_{\beta}(\mathbf{y}).$
(3)

The symmetric potentials $\mathscr{F}_{\alpha\beta}$ determine the correlation
length $\ell_{\ce{int}}$, and the magnitude of the interaction. The number of
additional parameters describing this interaction depends upon the model for
$\mathscr{F}_{\alpha\beta}$. In a previous publication we used a Lennard-
Jones-type force-field for $\mathscr{F}_{\alpha\beta}$.52 Such potentials are
often used in the literature.53, 54, 55, 56, 57 Furthermore, because

$\frac{\updelta^{2}F^{\ce{int}}}{\updelta c_{\beta}(\mathbf{z})\updelta
c_{\alpha}(\mathbf{x})}=\mathscr{F}_{\alpha\beta}(\lvert\mathbf{x}-\mathbf{z}\rvert),$
(4)

the potentials $\mathscr{F}_{\alpha\beta}$ determine the direct pair
correlation functions used in liquid state theory.50

Experimental results suggest that such interactions typically decay after some
ionic diameters.52 Thus, we focus on potentials $\mathscr{F}_{\alpha\beta}$
ranging over the size of one molecule. Their extend $\ell_{\ce{int}}$ is large
compared to the exponential decay of the electric field, i.e., the Debye-
length,58, 59 yet small compared to the battery cell.

In the SI (LABEL:sec:SIgradient-expansion), we show that such a convolution
functional $F^{\ce{int}}$ can be approximated in power-series of concentration
gradients when $\mathscr{F}_{\alpha\beta}$ are short-ranged,

$F^{\ce{int}}[c_{\gamma}]=\frac{1}{2}\sum_{\alpha,\,\beta}^{\textnormal{N}}\sum_{n=0}^{\infty}\varGamma^{2n}_{\alpha\beta}\int\ce{d}^{3}y\
c_{\alpha}(y)\cdot\nabla^{2n}c_{\beta}(y),$ (5)

where

$\varGamma^{n}_{\alpha\beta}(\mathscr{F}_{\alpha\beta})=\frac{1}{n!}\cdot\int\ce{d}^{3}x\,\mathscr{F}_{\alpha\beta}(|x|)\cdot
x^{n}.$ (6)

Here, $\varGamma^{n}_{\alpha\beta}$ are symmetric perturbation coefficients of
dimension
$[\varGamma^{2n}_{\alpha\beta}]=$\mathrm{J}\text{\,}\mathrm{m}${}^{3+2n}${\mathrm{mol}}^{-2}$$.
We state the complete free energy functional for IL electrolytes in the SI
(LABEL:eq:SI_const_eq_stress_tensor).

The excess chemical potentials are determined via functional derivatives
(LABEL:eq:SI_const_eq_chempot). In LABEL:sec:SIfunct-deriv, we show that this
Ansatz leads to force contributions

$\frac{\updelta F^{\ce{int}}}{\updelta
c_{\alpha}(z)}=\sum_{\beta=1}^{\textnormal{N}}\sum_{n=0}^{\infty}\varGamma^{2n}_{\alpha\beta}\cdot\nabla^{2n}c_{\beta}(z).$
(7)

The corresponding electrochemical potentials for ionic liquid electrolytes are

$\nabla\upmu^{\ce{el}}_{\alpha}=\left(Fz_{\alpha}-\upnu_{\alpha}\varrho\right)\nabla\varPhi+RT\cdot(\nabla
c_{\alpha})/c_{\alpha}-RT\upnu_{\alpha}\nabla c\\\
+\sum_{n=0}^{\infty}\sum_{\beta,\,\gamma}^{\textnormal{N}}\left(\updelta^{\alpha}_{\beta}-\upnu_{\alpha}c_{\beta}\right)\varGamma^{2n}_{\beta\gamma}\cdot\nabla^{2n+1}c_{\gamma}.$
(8)

We specify our electrolyte model and assume a one-dimensional Gaussian
interaction potential for symmetric ions,

$\mathscr{F}_{\alpha\beta}=\left(2\sqrt{2\pi}\right)^{3}\cdot\mathscr{V}^{0}\cdot(N_{A})^{2}\cdot\exp\left[-2\left(\frac{x\pi}{a}\right)^{2}\right].$
(9)

Here, $\mathscr{V}^{0}$ denotes the characteristic interaction energy, and $a$
is the extension of ion-pairs. This material parameter determines the
correlation length $\ell_{\ce{int}}=a/\pi\sqrt{2}$ of the interaction. We
assume that $a$ emerges naturally from the common molar volume
$\upnu=\upnu_{+}+\upnu_{-}$ via $\upnu=N_{A}a^{3}$, which is justified below
(see section IV.1.2 and eq. 48). Thus, only $\mathscr{V}^{0}$ is introduced as
a novel independent material parameter. In contrast, potentials of Lennard
Jones type need at least one more parameter for the well depth. In the case of
a binary system, the lowest order expansion coefficients (eq. 6) of the inter-
species correlations for the Gaussian model in eq. 9 are (LABEL:sec:SIgauss-
model-inter)

$\displaystyle\varGamma^{0}_{12}=\mathscr{V}^{0}N_{A}\upnu=\mathscr{V}^{0}\left(N_{A}\right)^{2}a^{3},$
(10)
$\displaystyle\varGamma^{2}_{12}=\frac{\ell_{\ce{int}}^{2}}{2}\cdot\varGamma^{0}_{12}=\frac{a^{2}}{4\pi^{2}}\cdot\varGamma^{0}_{12}.$
(11)

###  II.3 Binary IL

In this section, we apply our formalism to binary ILs at electrified
interfaces. Thus, we use the extended electrochemical forces eq. 8 in our
multicomponent framework derived in Ref. 45.

Figure 1: Scheme of the as-modelled set-up. The binary IL-electrolyte
PYR[1,4]TFSI is subject to the negatively charged interface at the left, which
causes the formation of an electrochemical double layer (EDL). Charge-ordering
diminishes with increasing distance from the interface (towards the right
side), and the electrolyte is electroneutral in the bulk.

As discussed in section II.1, binary electrolytes are described with the
variables electric potential $\varPhi$, charge density $\varrho$, and center-
of-mass convection velocity $\mathbf{v}$. Furthermore, the electric
conductivity $\kappa$ is the only independent transport parameter in this
case.45 The dynamical transport equations are

$\displaystyle\partial_{t}\varrho=-\nabla\left(\varrho\mathbf{v}\right)-\nabla\boldsymbol{\mathcal{J}},$
(12)
$\displaystyle\varrho=-\varepsilon_{0}\nabla\varepsilon_{\ce{R}}\nabla\varPhi,$
(13)
$\displaystyle\nabla\mathbf{v}=(M_{+}\upnu_{-}-M_{-}\upnu_{+})/Fz_{+}M_{\ce{IL}}\cdot\nabla\boldsymbol{\mathcal{J}}.$
(14)

Here, $M_{\pm}$ are the molar masses of the ionic species which sum to
$M_{\ce{IL}}$, and $\boldsymbol{\mathcal{J}}$ is the electric current relative
to the center-of-mass motion,
$\boldsymbol{\mathcal{J}}=\kappa(M_{+}\cdot\nabla\upmu^{\ce{el}}_{-}-M_{-}\nabla\upmu^{\ce{el}}_{+})/M_{\ce{IL}}Fz_{+}$.

Solutions to eqs. 12, 13 and 14 determine the ionic concentrations via
$\varrho=Fz_{+}(c_{+}-c_{-})$ (charge-conservation), and via the Euler
equation for the volume LABEL:eq:SI_elyte_eq_of_state.

We restrict our setup to one spatial dimension and assume that the inert
electrified interface is located at $x=$0$$. The electroneutral boundary
condition $\varrho(x{\to}\infty)=$0$$ implies that the bulk concentration
$c^{\ce{b}}=c_{\pm}(x{\to}\infty)$ is completely determined by the total
partial molar volume $\upnu=\upnu_{+}+\upnu_{-}$ via
$c^{\ce{b}}\cdot\upnu=$1$$. Because binary ILs are electrically neutral,
$z_{-}=-z_{+}$, and we choose $z_{+}>0$.

We neglect viscous forces in our discussion of the EDL
($\nabla\boldsymbol{\tauup}=$0$$). Therefore, the Gibbs-Duhem relation (see
LABEL:eq:SI_trivial_gibbs_duhem) becomes
$c_{+}\nabla\upmu^{\ce{el}}_{+}+c_{-}\nabla\upmu_{-}^{\ce{el}}=$0$$, and the
expression for the electric flux simplifies to
$\boldsymbol{\mathcal{J}}=-\kappa\upnu\uprho/Fz_{+}M_{\ce{IL}}\cdot\nabla\upmu_{\ce{IL}}$,
where we use the chemical potential of the anion-species to determine the IL-
electrolyte ($\nabla\upmu_{\ce{IL}}=\nabla\upmu^{\ce{el}}_{-}$),

$\nabla\upmu_{\ce{IL}}=\nabla\left(Fz_{+}\varPhi-\gammaup_{+}\frac{\updelta
F^{\ce{int}}}{\updelta c_{-}}+\gammaup_{-}\frac{\updelta
F^{\ce{int}}}{\updelta c_{+}}-\right.\\\
\left.-RT\gammaup_{+}\ln\left[\frac{c_{-}}{c^{\ce{b}}}\right]+RT\gammaup_{-}\ln\left[\frac{c_{+}}{c^{\ce{b}}}\right]\right).$
(15)

Here, we introduced the relative magnitude of the molar volumes
$\gammaup_{\pm}=\upnu_{\pm}/\upnu$. Thus, the forces given by eq. 15 depend
upon the model for $F^{\ce{int}}$. Furthermore, by using the Gauss-model (see
eq. 9), we can either close the forces via the “complete” integral equation
(eq. 3) or we can use the gradient expansion (eq. 7).

We want to apply the half-cell potential $\varDelta\phi$. Because the electric
potential $\varPhi$ is continuous across the electrode-electrolyte interface,
$\varPhi(0)$ in the electrolyte is subject to the boundary condition

$\varPhi(0)-\varPhi(x\to\infty)=\varDelta\phi.$ (16)

Without loss of generality, we set the electrolyte potential in the bulk to
zero, $\lim_{x\to\infty}\varPhi=$0$$. Hence, $\varDelta\phi=\varPhi($0$)$ is
the potential applied to the electrode.

We perform one-dimensional numerical simulations of this system of equations
12, and 13 in the completely dissociated state, subject to an inert
electrified interface (for more details, see LABEL:sec:comp-deta). This
electrolyte is part of the IL-family composed of TFSI anions and PYR cations.
Because of their excellent electrochemical properties, these ILs are widely
studied and used for applications in lithium-ion batteries.60 We state the
electrolyte parameters in the SI (see LABEL:sec:SI_binary-ionic-liquid).

###  II.4 Energy Scales and Dimensions

In this section, we clarify the notation and state the non-dimensional form of
principal quantities appearing in our theory. For a complete discussion, we
refer to the SI (see LABEL:sec:SInonimensionalization), where we share the
motivation for our choices.

We introduce dimensionless variables for electric potential
$\tilde{\varPhi}=\varPhi\cdot Fz_{+}/RT$, charge density
$\tilde{\varrho}=\varrho\cdot\upnu\tilde{c}^{\ce{b}}/Fz_{+}$, and
concentration $\tilde{c}_{\alpha}=c_{\alpha}\cdot\upnu\tilde{c}^{\ce{b}}$. As
a consequence, the Euler equation for the volume becomes
$\tilde{c}^{\ce{b}}=\gammaup_{+}\tilde{c}_{+}+\gammaup_{-}\tilde{c}_{-}$, with
the dimensionless molar volumes $\gammaup_{\pm}=\upnu_{\pm}/\upnu$ and the
dimensionless bulk concentration $\tilde{c}^{\ce{b}}$. The Poisson equation
suggests defining the generalized Debye-length

$L_{\ce{D}}=\sqrt{\frac{k_{\ce{B}}T\varepsilon_{0}\varepsilon_{\ce{R}}a^{3}\tilde{c}^{\ce{b}}}{(ez_{+})^{2}}}.$
(17)

Here, we used $F=eN_{\ce{A}}$ for the Faraday constant, where $e$ is the
elementary charge, and $N_{\ce{A}}$ is Avogadro’s constant, and the model
$\upnu=N_{\ce{A}}a^{3}$ for the partial molar volumes introduced above. This
Debye-length differs from the canonical definition by the asymmetry-factors
$\gammaup_{\pm}$.61 However, it reproduces the textbook definition for
symmetric ions ($\gammaup_{\pm}=$0.5$$). $L_{\ce{D}}$ becomes minimal for
$\gammaup_{\pm}=$0.5$$ because the mixing entropy of a binary electrolyte is
extremal for equal ion-size. Thus, asymmetry increases the Debye-screening
length.

With this length scale, we non-dimensionalize our grid, viz.
$\tilde{x}=x/L_{\ce{D}}$ and $\tilde{\nabla}=L_{\ce{D}}\cdot\nabla$, and
obtain the dimensionless Poisson-equation,

$\tilde{\varrho}=-\tilde{\nabla}^{2}\tilde{\varPhi}.$ (18)

In the SI (see LABEL:sec:SInonimensionalization), we non-dimensionalize the
transport equations (see eqs. 12 and 14) for binary symmetric ILs. Because we
neglect convective effects in our EDL discussion, the complete set of
equations consists of the Poisson equation and one transport equation for the
charges. By substituting eq. 3 into eq. 15, we find for the integral
description

$\partial_{\tilde{t}}\tilde{\varrho}=\tilde{\nabla}\Biggl{[}\left(1+\chiup\tilde{\varrho}\right)\cdot\tilde{\nabla}\Biggl{(}\tilde{\varPhi}-\gammaup_{+}\ln\left[\frac{\tilde{c}_{-}}{\tilde{c}^{\ce{b}}}\right]+\gammaup_{-}\ln\left[\frac{\tilde{c}_{+}}{\tilde{c}^{\ce{b}}}\right]-\\\
-\int\ce{d}\tilde{x}^{3}\tilde{\mathscr{F}}_{\alpha\beta}(|\tilde{\mathbf{x}}-\tilde{\mathbf{y}}|)\tilde{\varrho}(\tilde{\mathbf{y}})\Biggr{)}\Biggr{]},$
(19)

where
$\chiup=(\gammaup_{-}M_{+}/M_{\ce{IL}}-\gammaup_{+}M_{-}/M_{\ce{IL}})/\tilde{c}^{\ce{b}}$
measures the “asymmetry” of the ion-species, and
$\partial_{\tilde{t}}=\varepsilon_{0}\varepsilon_{\ce{R}}/\kappa\cdot\partial_{t}$.

The interaction potential is non-dimensionalized (see
LABEL:eq:SI_dimless_interaction_pot) by two energy-scales for thermal energy
$E_{\ce{th}}$ and electrostatic energy $E_{\ce{el}}$,

$\displaystyle
E_{\ce{th}}=k_{\ce{B}}T\cdot\frac{\tilde{c}^{\ce{b}}}{2\gammaup_{+}\gammaup_{-}},$
(20) $\displaystyle
E_{\ce{el}}=\frac{(ez_{+})^{2}}{$4$\pi\varepsilon_{0}\varepsilon_{\ce{R}}}\cdot\frac{$1$}{$4$\gammaup_{+}\gammaup_{-}}\cdot\frac{$1$}{a},$
(21)

such that
$\tilde{\mathscr{F}}_{\alpha\beta}=\mathscr{F}_{\alpha\beta}/(N_{\ce{A}})^{2}E_{\ce{th}}\cdot(L_{\ce{D}}/a)^{3}$
In the case of symmetric ions $\gammaup_{\pm}=$0.5$$, these energy scales take
the textbook form for thermal energy and Coulomb energy of charges at distance
$a$. Apparently, both energy scales are coupled by the generalized Debye-
length $L_{\ce{D}}$,

$E_{\ce{th}}/E_{\ce{el}}=8\pi\left(L_{\ce{D}}/a\right)^{2}.$ (22)

The integral form eq. 19 for the transport equation allows to relate our
continuum framework to MD-simulations as discussed in section V.1.

In this article we restrict the gradient-expansion of the interaction to the
trivial and first non-trivial modes
($\tilde{\varGamma}^{0}_{+-}=\mathscr{V}^{0}/E_{\ce{th}}$ and
$\tilde{\varGamma}^{2}_{+-}=2/\pi\cdot\mathscr{V}^{0}/E_{\ce{th}}\cdot
E_{\ce{el}}/E_{\ce{th}}$, see LABEL:eq:SI_nondim_pertmodes) and obtain

$\partial_{\tilde{t}}\tilde{\varrho}=\tilde{\nabla}\Biggl{[}\left(1+\chiup\tilde{\varrho}\right)\cdot\tilde{\nabla}\Biggl{(}\tilde{\varPhi}-\gammaup_{+}\ln\left[\frac{\tilde{c}_{-}}{\tilde{c}^{\ce{b}}}\right]+\gammaup_{-}\ln\left[\frac{\tilde{c}_{+}}{\tilde{c}^{\ce{b}}}\right]\\\
-\frac{\mathscr{V}^{0}}{E_{\ce{th}}}\left(1+\frac{2}{\pi}\frac{E_{\ce{el}}}{E_{\ce{th}}}\cdot\tilde{\nabla}^{2}\right)\tilde{\varrho}\Biggr{)}\Biggr{]}.$
(23)

###  II.5 Stationary State

Because our focus is the formation of equilibrium structures, we discuss the
system of equations in the stationary limit. This allows us to integrate the
differential equations using electroneutral boundary conditions, which results
in a simplified description susceptible to analytic techniques.

Stationarity ($\partial_{t}\varrho{=}$0$$) implies that all fluxes are
constant. Here, we have no flux conditions
$\tilde{\boldsymbol{\mathcal{J}}}=\tilde{\mathbf{v}}=$0$$, which implies that
both species are in equilibrium,
$\nabla\upmu^{\ce{el}}_{+}=\nabla\upmu^{\ce{el}}_{-}$, i.e.,

$\tilde{\nabla}\tilde{\upmu}_{\ce{IL}}=$0$.$ (24)

Thus, the stationary state for the binary electrolyte is described by the
Poisson equation and eq. 24. Here, we evaluate the equilibrium condition using
the gradient description (eq. 15) in the non-dimensionalized form (see
LABEL:eq:SI_nondim_pertmodes and eq. 23), and integrate the result using
electroneutral boundary conditions in the bulk,

$$0$=\tilde{\varPhi}-\frac{\mathscr{V}^{0}}{E_{\ce{th}}}\left($1$+\frac{$2$}{\pi}\frac{E_{\ce{el}}}{E_{\ce{th}}}\tilde{\nabla}^{2}\right)\tilde{\varrho}-\left(\gammaup_{+}\ln\left[\frac{\tilde{c}_{-}}{\tilde{c}^{\ce{b}}}\right]-\right.\\\
\left.-\gammaup_{-}\ln\left[\frac{\tilde{c}_{+}}{\tilde{c}^{\ce{b}}}\right]\,\right).$
(25)

Apparently, in contrast to the dynamical case where electrolyte momentum is
important, the molar masses appearing as parameters in the fluxes $\mathbf{v}$
and $\boldsymbol{\mathcal{J}}$ become irrelevant in the stationary limit.
Instead, the relative magnitude of the molar volumes $\gammaup_{\pm}$ enters
the system of equations. This highlights the principal role of molar volumes
as parameters in the stationary state and is a consequence of the Euler
equation for the volume, LABEL:eq:SI_elyte_eq_of_state.

For completeness, we state the integral transport equation 19 in the
stationary limit (see LABEL:sec:inter-contr-chem)

$$0$=\tilde{\varPhi}-\int\ce{d}\tilde{y}\tilde{\mathscr{F}}_{+-}\left(\tilde{x},\tilde{y}\right)\tilde{\varrho}(\tilde{y})-\left(\gammaup_{+}\ln\left[\frac{\tilde{c}_{-}}{\tilde{c}^{\ce{b}}}\right]-\right.\\\
\left.-\gammaup_{-}\ln\left[\frac{\tilde{c}_{+}}{\tilde{c}^{\ce{b}}}\right]\right).$
(26)

###  II.6 Small and Large Potentials

Equation 18 and eq. 25 (or eq. 26) constitute the complete set of equations,
necessary to describe a binary IL-electrolyte in a stationary state. In
section III.1 and section IV.3, we solve these equations using numerical
methods. Our goal is to supplement these numerical methods by an analytic
examination of the gradient description. However, the analytic solution of the
gradient description is hindered by the higher order gradients appearing in
eq. 25 and by the different prefactors of the logarithmic terms (in general,
$\gammaup_{+}\neq\gammaup_{-}$). Therefore, we distinguish different limiting
cases in our analysis in sections III.2 and IV.1. In the SI, we describe the
special case of symmetric ion-species (see LABEL:sec:integr-transp-equat and
LABEL:sec:SImean-volume-effect).

In section III.2 we show that the limiting case of small charge-densities,
$\tilde{\varrho}\ll$1$$, is useful. In this case, we can expand the
logarithmic term in eq. 25 and eq. 26 around the electroneutral state,

$\gammaup_{+}\ln\left[\tilde{c}_{-}/\tilde{c}^{\ce{b}}\right]-\gammaup_{-}\ln\left[\tilde{c}_{+}/\tilde{c}^{\ce{b}}\right]\approx-\tilde{\varrho},$
(27)

such that eqs. 18 and 25 become

$\displaystyle\tilde{\varrho}=-\tilde{\nabla}^{2}\tilde{\varPhi},$ (28)
$\displaystyle$0$=\tilde{\varPhi}+\hat{\varepsilon}_{\ce{R}}\tilde{\varrho},$
(29)

where $\hat{\varepsilon}_{\ce{R}}$ is defined as the dielectric operator

$\hat{\varepsilon}_{\ce{R}}=1-\sum_{n=0}^{1}\tilde{\varGamma}^{2n}_{+-}\cdot\tilde{\nabla}^{2n}=$1$-\frac{\mathscr{V}^{0}}{E_{\ce{th}}}-\frac{2}{\pi}\frac{\mathscr{V}^{0}}{E_{\ce{th}}}\frac{E_{\ce{el}}}{E_{\ce{th}}}\tilde{\nabla}^{2}.$
(30)

In the absence of molecular repulsion, $\mathscr{V}^{0}=$0$$, the dielectric
operator reduces to the canonical, scalar-valued dielectric parameter
$\hat{\varepsilon}_{\ce{R}}\to$1$$.

Furthermore, quantities similar to $\hat{\varepsilon}_{\ce{R}}$ also arise in
the liquid state theory of classical statistical mechanics. This expansion
corresponds to a small wave vector expansion of the dielectric function
expressed as a correlation function of the molecular dipole densities (see,
e.g., Ref.50).

In the following sections, we show that the gradient expansion, eqs. 29 and
28, allows significant insights into the competing effects of the interactions
$\mathscr{V}^{0},E_{\ce{el}},E_{\ce{th}}$, and predicts the EDL structure as a
function of the temperature, dielectricity, ion-size, and interaction
strength.

##  III Mean Steric Effect: Charge Saturation

Figure 2: Simulation results of the EDL structure perpendicular to the
electrode-electrolyte interface for a binary IL for $\mathscr{V}^{0}=$0$$ (see
eqs. 18 and 23 ). If not mentioned otherwise, $T=$300\text{\,}\mathrm{K}$$,
$\varepsilon_{\ce{R}}=$15$$, and $\varDelta\phi=$-0.1\text{\,}\mathrm{V}$$. a)
Profiles of electric potential and charge density (inset) for different
electrode potentials $\varDelta\phi$. b) Concentration profiles of the anions
and cations for different electrode potentials $\varDelta\phi$. c)
Concentration profiles for different volume ratios $\gammaup_{+}$. The inset
shows the corresponding electric potential. d) Concentration profiles for
varying dielectric constants (dashed lines) and temperatures (solid lines).

In this section, we neglect non-local interactions, $\mathscr{V}^{0}=0$, and
discuss the EDL structure of the electrolyte due to bulk effects alone based
on $F^{\ce{b}}$ (see eq. 1). In this way, we reveal the competition between
Coulombic ordering and entropic disordering, i.e., diffusion.

Toward this aim, we consider the system of equations constituted by the
Poisson eq. 18 and eq. 25 subject to $\mathscr{V}^{0}=0$,

$0=\tilde{\varPhi}-\gammaup_{+}\ln\left[\tilde{c}_{-}/\tilde{c}^{\ce{b}}\right]\\\
+\gammaup_{-}\ln\left[\tilde{c}_{+}/\tilde{c}^{\ce{b}}\right].$ (31)

First, in section III.1, we solve this system of equations numericaly. We
supplement this investigation by an analytic analysis, and focus on the two
limiting regimes of large and small electric potentials. In section III.2.1 we
discuss the case $\tilde{\varPhi}\ll$1$$, and in section III.2.2, we discuss
the case $\tilde{\varPhi}\gg$1$$. For the special case of symmetric ion-
species ($\gammaup_{\pm}=0.5$), we derive analytic solutions for the electric
field $\tilde{\boldsymbol{E}}(\tilde{\varPhi})$, and for the charge density
$\tilde{\varrho}(\tilde{\varPhi})$ as functions of the electrolyte electric
potential in the SI (see LABEL:sec:SImean-volume-effect).

###  III.1 Simulations

Figure 2 shows numeric results for the system of eqs. 18 and 31. Figure 2a and
fig. 2b illustrate screening profiles of the electric potential, the charge
density, and the ion concentrations for varying electrode potentials
$\varDelta\phi$.

Apparently, the application of a negative electrode potential
($\varDelta\phi<$0$$) polarizes the electrolyte. The electric potential (see
fig. 2a) is continuous across the interface and decays smoothly towards the
electroneutral bulk region. The inset of fig. 2a shows that, for low electrode
potentials, the charge density decays exponentially. A similar behavior can be
observed in fig. 2b for the concentrations. The concentration of positive
counter-ions increases toward the interface, whereas negative ions get
depleted. Apparently, the electrolyte screens the electrode potential by
accumulation of counter-ions. However, above
$\varDelta\phi\approx$-0.05\text{\,}\mathrm{V}$$, the counter-ion
concentration saturates near the interface. A further increase of
$\varDelta\phi$ broadens the EDL.

This behavior can be explained by the mean volume effect. The application of a
negative potential $\varDelta\phi$ implies that positive ions accumulate near
the interface, and negative ions deplete. However, the Euler equation for the
volume, LABEL:eq:SI_elyte_eq_of_state, implies the saturation concentration
$c_{+}^{\ce{sat}}=$1$/\upnu_{+}$. Once the accumulated species reaches this
saturation, the screening mechanism transitions from increasing the
concentration at the interface to broadening the width of the EDL. The
simulated EDL approaches a thickness of $0.6\text{\,}\mathrm{nm}$ at
$\varDelta\phi\approx$-0.05\text{\,}\mathrm{V}$$, which is significantly wider
than predicted by the canonical Debye-Hückel-theory with the Debye length
$L_{\ce{D}}\approx$0.7\text{\,}\mathrm{\SIUnitSymbolAngstrom}$$ (see eq. 17).
This phenomenon is typically denoted “crowding”.35

Because the saturation concentration depends upon the molar volume,
$c_{\alpha}^{\ce{sat}}=$1$/\upnu_{\alpha}$, the partial molar volumes directly
affect the screening behavior. Figure 2c shows numerical results for the ionic
concentrations at different volume ratios $\gammaup_{\pm}=\upnu_{\pm}/\upnu$
(in which $\upnu$ is kept fixed). The EDL is thinner for smaller
$\gammaup_{+}$, as this allows for tighter packing of cations.

The effect of temperature $T$ and dielectricity $\varepsilon_{\ce{R}}$ on the
EDL structure is illustrated in fig. 2d. The screening is more effective for
smaller values of $\varepsilon_{\ce{R}}$ and the EDL width increases with
increasing magnitude of $\varepsilon_{\ce{R}}$. This is in qualitative
agreement with the screening behavior for dilute solutions as predicted by the
Debye-Hückel theory. Likewise, the EDL becomes more diffuse with increasing
temperature because of the disordering effect of thermal motion. The observed
effects of $T$ and $\varepsilon_{\ce{R}}$ highlight the competing interplay
between the electrostatic effect of charge ordering and the disordering effect
of entropy.

To summarize, the simulations show two distinct regimes of EDL structure.
First, for large electrode potentials, $\Delta\tilde{\phi}\gg$1$$, the charge
is saturated near the interface. Second, near the electroneutral bulk region,
where $\Delta\tilde{\phi}\ll$1$$, the charge density decays exponentially.
These two distinct EDL structures, charge saturation and exponential decrease,
correspond to two disjointed electrolyte regimes,
$\lvert\tilde{\varPhi}\rvert\gg$1$$ and $\lvert\tilde{\varPhi}\rvert\ll$1$$,
respectively.

###  III.2 Asymptotic Analysis

Figure 3: Comparison of the asymptotic analysis of the EDL structure with
numerical results obtained from eqs. 18 and 23 ($T=$300\text{\,}\mathrm{K}$$,
$\mathscr{V}^{0}=$0$$ and $\varepsilon_{\ce{R}}=$15$$). We consider two
different values of the interface potential, corresponding to the regimes of
small and large potentials at $\varDelta\phi=$-0.01\text{\,}\mathrm{V}$$
($\Delta\tilde{\phi}=$-0.4$$), and $\varDelta\phi=$-0.35\text{\,}\mathrm{V}$$
($\Delta\tilde{\phi}=$-13.5$$). a) Profiles of the electric potentials as
predicted by the analytical approximation (see eqs. 33 and 36). The inset
highlights the region close to the interface for the case
$\varDelta\phi=$-0.01\text{\,}\mathrm{V}$$. b) Profiles of the charge density
as predicted by the asymptotic analysis in eqs. 34 and 38. The inset compares
numerical and analytical values for the total charge in the EDL (see
LABEL:sec:SI_asympt-analys-large in the SI).

The simulation results in section III.1 motivate our procedure for analysing
the EDL structure. First, we study the EDL far away from the electrode close
to the electroneutral bulk (large $\tilde{x}$). For this purpose, we expand
the stationary equations around the electroneutral bulk for small charge
densities, $\tilde{\varrho}\ll$1$$. According to eq. 31, this regime
corresponds to small dimensionless potentials
$\lvert\tilde{\varPhi}\rvert\ll$1$$. Note that this coincides with the high-
temperatures regime because $\tilde{\varPhi}=\varPhi\cdot Fz_{+}/RT$ (see
section III.2.1).

Second, we analyze the behavior close to the electrode (small $\tilde{x}$) at
large electrode potentials $\lvert\Delta\tilde{\phi}\rvert\gg$1$$, where the
electrolyte potential satisfies $\lvert\tilde{\Phi}\rvert\gg$1$$. This
corresponds to low temperatures (see section III.2.2).

####  III.2.1 Asymptotic Analysis: Small Potentials
$\lvert\tilde{\varPhi}\rvert\ll$1$$

As outlined above, we begin our analytic examination of the EDL in the limit
of small dimensionless potentials $\lvert\tilde{\varPhi}\rvert\ll$1$$. Our
idea is to approach the EDL from the electroneutral bulk region along the
direction of decreasing $\tilde{x}$. Thus, we use the expansion of ionic
concentrations around the bulk electrolyte in eq. 27,
$\tilde{c}_{\pm}=\tilde{c}^{\ce{b}}\pm\tilde{\varrho}\gammaup_{\mp}\approx\tilde{c}^{\ce{b}}$.

For this aim, we insert $\tilde{\varrho}$ from eq. 28 into eq. 29 for
$\hat{\varepsilon}_{\ce{R}}=$1$$, yielding

$\tilde{\nabla}^{2}\tilde{\varPhi}=\tilde{\varPhi}.$ (32)

With the boundary conditions discussed above,
$\tilde{\varPhi}($0$)=\varDelta\tilde{\phi}=Fz_{+}\varDelta\phi/RT$ and
$\tilde{\varPhi}(\tilde{x}\to\infty)=$0$$, we obtain the solution

$\displaystyle\tilde{\varPhi}=\varDelta\tilde{\phi}\cdot\exp(-\tilde{x}),$
(33)
$\displaystyle\tilde{\varrho}=-\varDelta\tilde{\phi}\cdot\exp(-\tilde{x}).$
(34)

The corresponding dimensionalized electrolyte potential,

$\varPhi=\varDelta\phi\cdot\exp(-x/L_{\ce{D}}),$ (35)

decays exponentially with damping parameter $1/L_{\ce{D}}$. Thus, the decay
length in the limit $\lvert\tilde{\varPhi}\rvert\ll$1$$ is the Debye length
$L_{\ce{D}}$ (see eq. 17).

In fig. 3, we compare the analytic predictions for this limit (dashed green
lines) with the numeric results (blue lines) for different electrode
potentials. Apparently, the analytical and numerical results for the electric
potential and the charge density are in excellent agreement for small
electrode potentials, $\Delta\tilde{\phi}\approx$-0.4$$, when the condition
$\lvert\tilde{\varPhi}\rvert\ll$1$$ is fulfilled everywhere.

In the SI, we derive the expressions for the total EDL-surface-charge-density
and the differential capacitance (see LABEL:sec:SI_asympt-analys-small).

####  III.2.2 Asymptotic Analysis: Large Potentials
$\lvert\tilde{\varPhi}\rvert\gg$1$$

Next, we discuss the EDL in the limit of large potentials
$\lvert\tilde{\phi}\rvert\gg$1$$. This regime can be found for large electrode
potentials $\Delta\tilde{\varPhi}\gg$1$$ close to the electrode/electrolyte
interface. As $\tilde{\varPhi}=Fz_{+}\varDelta\phi/RT$, this analysis is exact
in the limit of zero temperature, $T=$0$$.

In this case, the logarithmic terms in eq. 31 must compensate the diverging
potential term $\tilde{\varPhi}$. Because of the mean volume constraint (see
LABEL:eq:SI_elyte_eq_of_state), one of the logarithmic terms is diverging if
one species is depleted and the other species saturates,
$\varrho(x=$0$)=Fz_{\alpha^{\ast}}/\upnu_{\alpha^{\ast}}$. Here, we denote the
saturating species by the index $\alpha^{\ast}$. Because electric potentials
are continuous across interfaces, the saturation species $\alpha^{\ast}$ is
uniquely determined by the sign of the electrode potential,
$\ce{sign}(z_{\alpha^{\ast}})=-\ce{sign}\,(\varDelta\tilde{\phi})$.

Therefore,
$\tilde{\varrho}^{\ce{sat}}=-\ce{sign}(\varDelta\tilde{\phi})\cdot\tilde{c}^{\ce{sat}}$
solves eq. 31, where
$\tilde{c}^{\ce{sat}}=\tilde{c}^{\ce{b}}/\gammaup_{\alpha^{\ast}}$. Upon
integration of the Poisson eq. 28, we find

$\tilde{\varPhi}(\tilde{x})=\varDelta\tilde{\phi}\cdot\left($1$-\tilde{x}/\tilde{L}_{\ce{EDL}}\right)^{2},$
(36)

where the width of the EDL depends on the electrode potential,
$\tilde{L}_{\ce{EDL}}=\sqrt{$2$\gammaup_{\alpha^{\ast}}\lvert\varDelta\tilde{\phi}\rvert/\tilde{c}^{\ce{b}}}$.
Thus, the dimensionalized EDL length is

$L_{\ce{EDL}}=L_{\ce{D}}\tilde{L}_{\ce{EDL}}=\sqrt{$2$a^{3}\lvert\varDelta\tilde{\phi}\rvert\gammaup_{\alpha^{\ast}}k_{\ce{B}}T\varepsilon_{0}\varepsilon_{\ce{R}}/(ez_{+})^{2}}.$
(37)

Apparently, the decay length increases with ion size because of the mean-
volume effect. Also, it scales with the asymmetry
$\sqrt{\gammaup_{\alpha^{\ast}}}$; i.e., it is small for small screening
species. Comparison with the Debye screening-length shows that
$L_{\ce{EDL}}>L_{\ce{D}}$ in the limit of small temperatures $T$ or large
potentials $\lvert\varDelta\tilde{\phi}\rvert\gg$1$$.

In the limit of vanishing temperature, $T=$0$$, the charge profile is box
shaped and is determined by the screening-length $L_{\ce{EDL}}$,

$\tilde{\varrho}=-\theta(\tilde{L}_{\ce{EDL}}-\tilde{x})\cdot\ce{sign}(\varDelta\tilde{\phi})\cdot\tilde{c}^{\ce{sat}},$
(38)

where $\theta$ is the Heaviside function. In the SI, we calculate the analytic
expressions for the total charge in the EDL and for the differential
capacitance in this limit (see LABEL:sec:SI_asympt-analys-large).

In fig. 3, we compare the analytical predictions for this limit (dashed yellow
lines) with the numerical results (solid red lines) for
$\Delta\tilde{\phi}=$-13.5$$. The box-shaped charge-profile is in good
qualitative agreement with the numerical results because it almost predicts
the correct width of the EDL. However, the transition from saturation to the
bulk-state is more diffuse in the numerical profile. This is due to the
entropic, thermal influence, which “washes out” the box. Nevertheless, the
inset shows that the charge in the EDL, as predicted by the analytical
approximation, is quantitatively in good agreement with the numerical results.
We note that this profile and its temperature dependence are reminiscent of
the Fermi distribution.

##  IV Non-Local Interactions: Charge Oscillations

In this section, we discuss the influence of non-local interactions
($\mathscr{V}^{0}\neq 0$) on the structure of the EDL. As in section III, we
discuss the two limiting cases of small and large potentials separately.

###  IV.1 Static Asymptotic Analysis

####  IV.1.1 Asymptotic Analysis: Large Potentials
$\lvert\tilde{\varPhi}\rvert\gg$1$$

Let us first discuss the regime of diverging electrolyte potentials
$\lvert\tilde{\varPhi}\rvert\to\infty$. In this limit, the interaction
contribution cannot compensate the diverging electrolyte potential in eq. 25.
The logarithmic terms are diverging if one species is depleted, and we recover
the same results as described in section III.2.2 for the case of vanishing
interaction contributions $\mathscr{V}^{0}=0$.

####  IV.1.2 Asymptotic Analysis: Small Potentials
$\lvert\tilde{\varPhi}\rvert\ll$1$$

Figure 4: Phase spaces for EDL structure as functions of $T$ and
$\varepsilon_{\ce{R}}$ for equally sized ions ($\gammaup_{\pm}=$0.5$$), see
eq. 45. a) Critical interaction energies $\mathscr{V}^{0}_{\pm}$ as function
of temperature (here, $\varepsilon_{\ce{R}}=$15$$,
$a=$1.3\text{\,}\mathrm{nm}$$). b) Critical interaction energies
$\mathscr{V}^{0}_{\pm}$ as function of dielectricity (here,
$T=$300\text{\,}\mathrm{K}$$ and $a=$1.3\text{\,}\mathrm{nm}$$). Three phases
are present: exponentially damped charge density (shaded red), decaying
oscillatory charge density (shaded gray), quasi-crystalline (shaded blue).

In this section, we consider the full theory with molecular repulsion in the
regime of small potentials $\lvert\tilde{\varPhi}\rvert\ll$1$$. As outlined
above, we expand the interaction free energy in gradients of the charge
density, and restrict our analysis to the first two perturbation modes $n=$0$$
and $n=$1$$. With the assumption of small charge densities $\tilde{\varrho}\ll
1$, we derived linear equations above (eqs. 28 and 29), which we rephrase in
matrix form,

$\displaystyle\left(\begin{matrix}\tilde{\nabla}^{2}\tilde{\varrho}\\\
\tilde{\nabla}^{2}\tilde{\varPhi}\end{matrix}\right)=\tilde{\boldsymbol{\mathscr{A}}\
}\cdot\left(\begin{matrix}\tilde{\varrho}\\\
\tilde{\varPhi}\end{matrix}\right),$ (39) where
$\displaystyle\tilde{\boldsymbol{\mathscr{A}}\
}=\left(\begin{matrix}\displaystyle\frac{\pi
E_{\ce{th}}^{2}}{$2$E_{\ce{el}}\mathscr{V}^{0}}\left($1$-\frac{\mathscr{V}^{0}}{E_{\ce{th}}}\right)&\quad\quad\quad&\displaystyle\frac{\pi
E_{\ce{th}}^{2}}{$2$E_{\ce{el}}\mathscr{V}^{0}}\\\ \\\
-$1$&&0\end{matrix}\right).$ (40)

We solve eq. 39 via the eigenvalue decomposition with the eigenvalues
$\tilde{\upalpha}_{1,2}$ of the matrix in eq. 40. These are determined by the
relative magnitudes of the three competing energy scales, $E_{\ce{th}}$,
$E_{\ce{el}}$, and $\mathscr{V}^{0}$,

$\tilde{\upalpha}_{1,2}=-\frac{\pi}{$4$}\cdot\frac{E_{\ce{th}}}{E_{\ce{el}}}\cdot\left[$1$-\frac{E_{\ce{th}}}{\mathscr{V}^{0}}\mp\sqrt{\left(1-\frac{E_{\ce{th}}}{\mathscr{V}^{0}}\right)^{2}-\frac{$8$}{\pi}\cdot\frac{E_{\ce{el}}}{\mathscr{V}^{0}}}\right].$
(41)

Each eigenvalue $\tilde{\upalpha}_{i}$ gives rise to a dimensionless wave-
vector,

$\tilde{k}_{1,2}=\sqrt{\tilde{\upalpha}_{1,2}}.$ (42)

These determine the general solution of eq. 39 together with the eigenvectors
$\tilde{\boldsymbol{a}}^{\tilde{\upalpha}_{i}}=(\tilde{a}^{\tilde{\upalpha}_{i}}_{1},\tilde{a}^{\tilde{\upalpha}_{i}}_{2})^{T}=(-\tilde{\upalpha}_{i},1)^{T}$,

$\left(\begin{matrix}\tilde{\varrho}\\\
\tilde{\varPhi}\end{matrix}\right)=\left(\begin{matrix}\boldsymbol{a}^{\tilde{\upalpha}_{1}}&\boldsymbol{a}^{\tilde{\upalpha}_{2}}\end{matrix}\vphantom{\left(\begin{matrix}A_{1}e^{\tilde{k}_{1}\tilde{x}}+A_{2}e^{-\tilde{k}_{1}\tilde{x}}\\\
A_{3}e^{\tilde{k}_{2}\tilde{x}}+A_{4}e^{-\tilde{k}_{2}\tilde{x}}\end{matrix}\right)}\right)\cdot\left(\begin{matrix}A_{1}e^{\tilde{k}_{1}\tilde{x}}+A_{2}e^{-\tilde{k}_{1}\tilde{x}}\\\
A_{3}e^{\tilde{k}_{2}\tilde{x}}+A_{4}e^{-\tilde{k}_{2}\tilde{x}}\end{matrix}\right).$
(43)

The expansion coefficients $A_{i}$ are determined by boundary conditions and
physical arguments. Apparently, the corresponding wave-vectors are functions
$\tilde{k}_{1,2}(T,z_{\alpha},\varepsilon_{\ce{R}},\upnu_{\alpha},\mathscr{F}_{\alpha\beta})$,
which determine the structure of the EDL,

$\tilde{k}_{1,2}\in\begin{cases}\mathbb{R},&\textnormal{exponential
damping},\\\ \mathbb{R}+\mathrm{i}\cdot\mathbb{R},&\textnormal{damped
oscillations},\\\
\mathrm{i}\cdot\mathbb{R},&\textnormal{oscillations}.\end{cases}$

Thus, the EDL structure depends upon the relative magnitudes of the energies
$E_{\ce{th}}$, $E_{\ce{el}}$, and $\mathscr{V}^{0}$ via eq. 41. In particular,
the classification of $\tilde{k}_{\pm}$ depends upon the sign of the root

$\mathscr{W}=\left($1$-E_{\ce{th}}/\mathscr{V}^{0}\right)^{2}-$8$E_{\ce{el}}/\mathscr{V}^{0}\pi$
(44)

appearing in eq. 41. Thus, the critical values $\mathscr{V}^{0}_{\pm}$,
defined by the condition $\mathscr{W}(\mathscr{V}^{0}_{\pm})=0$, determine the
thresholds for the transition between the phases of EDL structure,

$\mathscr{V}^{0}_{\pm}=E_{\ce{th}}+$4$E_{\ce{el}}/\pi\pm$2$\sqrt{$2$E_{\ce{el}}($2$E_{\ce{el}}/\pi+E_{\ce{th}})/\pi}.$
(45)

Thus, with eq. 45 we can draw the phase diagram for the EDL structures. Since
$$0$<\mathscr{V}^{0}_{-}<\mathscr{V}^{0}_{+}$, there are three phases. In the
SI, we discuss each case in great detail (see LABEL:sec:SIphase-analysis).
Next, we give a short description of each phase.

Phase 1: $$0$\leq\mathscr{V}^{0}\leq\mathscr{V}^{0}_{-}$. In this regime,
$\tilde{\upalpha}_{1,2}\geq$0$$, which implies a real-valued wave-vector.
Thus, this phase corresponds to exponentially damped profiles,
$\tilde{\varrho}\propto\exp(-\tilde{k}_{\mathbb{R}}\tilde{x})$. A harmonic
analysis of the root $\sqrt{\mathscr{W}}$ reveals that (see
LABEL:sec:SIlimiting-cases)

$\lim_{\mathscr{V}^{0}\to$0$}\tilde{k}_{1}=\infty\quad\text{and}\quad\lim_{\mathscr{V}^{0}\to
0}\tilde{k}_{2}=$1$.$ (46)

Thus, solutions with the damping parameter $\tilde{k}_{1}$ vanish quickly and
are rendered as unphysical, whereas the limit of vanishing interactions for
$\tilde{k}_{2}$ reproduces the “bulk”-expansion for $\tilde{\varPhi}\ll 1$
from section III.2.1, see eq. 33.

Phase 2: $\mathscr{V}^{0}_{-}<\mathscr{V}^{0}<\mathscr{V}^{0}_{+}$. In this
regime, $\mathscr{W}<0$, and, thus, the root eq. 44 becomes complex.
Therefore, the wave-vector has non-vanishing real and imaginary parts,
$\tilde{k}_{1,2}\in\mathbb{R}\times\mathrm{i}\cdot\mathbb{R}$. This
corresponds to charge-profiles of exponentially damped oscillations,

$\tilde{\varrho}\propto\exp(-\tilde{k}_{\mathbb{R}}\tilde{x})\cdot\cos(\tilde{k}_{\mathbb{C}}\tilde{x}).$
(47)

Figure 5: Real and imaginary parts of the non-dimensionalized wave-vector
$\tilde{k}=k\cdot L_{\ce{D}}$ as function of the relative magnitude of the
energy scales $\mathscr{V}^{0}$ and $E_{\ce{th}}$ (see eqs. 41 and 42). Here,
$T=$300\text{\,}\mathrm{K}$$, $\varepsilon_{\ce{R}}=$15$$, and
$a=$1.3\text{\,}\mathrm{nm}$$.

Phase 3: $\mathscr{V}^{0}_{+}\leq\mathscr{V}^{0}$. In this regime, both
eigenvalues are real but negative, $\tilde{\upalpha}_{1,2}\leq 0$. Therefore,
$k_{\pm}\in\mathrm{i}\cdot\mathbb{R}$, which corresponds to undamped
oscillatory profiles. The limiting case for indefinitely strong interactions
yields

$\lim_{\mathscr{V}^{0}\to\infty}k_{1,2}=\lim_{\mathscr{V}^{0}\to\infty}\tilde{k}_{1,2}/L_{\ce{D}}=\begin{cases}0,\\\
\pm\mathrm{i}\cdot 2\pi/a.\end{cases}$ (48)

Thus, the result for $\lim_{\mathscr{V}^{0}\to\infty}k_{2}$ reproduces the
experimental findings obtained by AFM-measurements that the wave-lengths
$\uplambda\approx a/$2$\pi$ of the observed oscillations scale with the size
of molecules $a$.52, 62 Apparently, the incompressibility of ions in our model
prevents a further decrease of the wavelength.

Thus, the critical values $\mathscr{V}^{0}_{\pm}$ constitute exactly the
boundaries between the different EDL phases.

In fig. 4, we illustrate the phase space of EDL structures as functions of
temperature and dielectricity (in LABEL:sec:phase-space-scre we also show the
phase space as function of ion-size, ion-asymmetry and valency). Apparently,
three distinct phases of EDL structures are present. The exponentially damped
EDL phase corresponds to the regions below $\mathscr{V}^{0}_{-}$ (red line),
whereas the damped-oscillatory EDL phase corresponds to the regions between
the blue and red lines. Finally, the undamped oscillatory EDL phase
corresponds to the regions above $\mathscr{V}^{0}_{+}$ (blue line). Figure 4a
illustrates that temperature $T$ and hardcore-interactions $\mathscr{V}^{0}$
are in competition and that the critical interaction strengths
$\mathscr{V}^{0}_{\pm}$ increase with increasing temperature $T$, i.e.,
thermal energy $E_{\text{th}}$.

Figure 4b reveals the influence of dielectricity, i.e., electro-static forces,
on the EDL phases. Apparently, the damped oscillatory phase becomes narrower
for ILs with larger dielectricity $\varepsilon_{\ce{R}}$, i.e., smaller
electric energy $E_{\text{el}}$.

We note that the exponentially damped regime for small interaction strengths
(small compared to electrostatic and thermal energy) corresponds to the EDL
structure found in section III.2.1 in the absence of hardcore-interactions.
However, as can be inferred from fig. 4a, this phase is hardly present for
reasonable temperatures.

Figure 5 shows the non-dimensionalized wave-vector as function of the relative
energy scale $\mathscr{V}^{0}/E_{\ce{th}}$ (where
$T{=}$300\text{\,}\mathrm{K}$$, $\varepsilon_{\ce{R}}{=}$15$$, and
$a{=}$1.3\text{\,}\mathrm{nm}$$). For small interaction energies
$\mathscr{V}^{0}<\mathscr{V}^{0}_{-}$ the wave-vector is real,
$\tilde{k}=\tilde{k}_{\mathbb{R}}$, which corresponds to exponentially damped
profiles. In particular, the static profile at $\mathscr{V}^{0}\to 0$
reproduces the case described in section III.2 where the exponential profile
is determined by the scales $E_{\ce{th}}$ and $E_{\ce{el}}$ alone, i.e.
$\tilde{k}=1$. Apparently, $\tilde{k}_{\mathbb{R}}$ increases with
$\mathscr{V}^{0}$ up to the threshold $\mathscr{V}^{0}_{-}$, beyond which
which it starts to decrease. Thus, the EDL has minimal extension at
$\mathscr{V}^{0}=\mathscr{V}^{0}_{-}$. This suggests that the increasing
strength of the repulsive ion-correlations compresses the screening-layer.
Once the hardcore potential exceeds $\mathscr{V}^{0}_{-}$, the system
overscreens, i.e., the ion-layers begin to oscillate. The damping parameter
$\text{Re}(\tilde{k})$ vanishes exactly when
$\mathscr{V}^{0}=\mathscr{V}^{0}_{+}$, i.e., when the system transitions into
nano-segregation of the ion species. Interestingly, the frequency of the
oscillations $\text{Im}(\tilde{k})$ exhibits a local maximum and minimum in
the regime of damped oscillations. Furthermore, $\text{Im}(\tilde{k})$ attains
its maximal value in the limit of prevailing interaction strength
$\mathscr{V}^{0}\to\infty$.

In the SI (see LABEL:sec:SItriv-order-expans and LABEL:sec:SIline-order-
expans) we investigate the influence of the individual perturbation modes
$\varGamma^{0}_{12}$ and $\varGamma^{2}_{12}$ on the phase diagram. As it
turns out, neglecting all but the zeroth-order correction $\varGamma^{0}_{12}$
results in a binary phase diagram comprising only exponentially damped
profiles and undamped, oscillatory profiles. In contrast, neglecting the
zeroth-order correction, and taking only the first non-trivial order
$\varGamma^{2}_{12}$ into account, results in a binary phase-diagram
comprising only exponentially damped profiles and damped oscillatory profiles.
This is the case for MFTs based on the BSK-framework. Thus, for the “complete”
set of the three different phases, both perturbation modes
$\varGamma^{0}_{12}$ and $\varGamma^{2}_{12}$ are neccessary.

Interestingly, for the pathologic case of negative interaction-strengths
$\mathscr{V}^{0}<0$, the phase space reduces to the two screening phases of
exponentially damped profiles and undamped oscillatory profiles. This follows
straightforwardly from eq. 41 (see also the discussions in LABEL:sec:SItriv-
order-expans and LABEL:sec:SIline-order-expans).

###  IV.2 Dynamic Asymptotic Analysis: Linear Stability Analysis

In this section, we complement the static analysis of section IV.1 by an
analytic analysis of the dynamic transport equation in the gradient
description (eq. 23).

For this purpose, we perform a linear stability analysis and consider the
limit of small potentials, $\lvert\tilde{\varPhi}\rvert\ll 1$. Thus, the
logarithmic terms can approximated as in eq. 27, and eq. 23 becomes

$\partial_{\tilde{t}}\tilde{\varrho}=-\tilde{\nabla}^{2}\left(\tilde{\varPhi}+\hat{\varepsilon}_{\ce{R}}\tilde{\varrho}\right).$
(49)

We expand the electric potential around an uniform bulk-state
$\tilde{\varPhi}^{\ce{b}}$,

$\tilde{\varPhi}=\tilde{\varPhi}^{\ce{b}}+\sum_{i=1}^{\infty}\epsilon^{i}\cdot\tilde{\varPhi}_{i}.$
(50)

Here, the equilibrium state is determined by the electroneutral bulk-condition
$\tilde{\varPhi}^{\ce{b}}=0$ and $\tilde{\varrho}^{\ce{b}}=0$. Thus, the first
order perturbation takes the form

$\tilde{\varPhi}_{1}=\exp\left[\tilde{s}\tilde{t}\right]\cdot\cos\left[\tilde{k}\tilde{x}\right].$
(51)

Here, the wave-number $\tilde{k}$ determines the spatial distribution of the
dimensionless perturbation $\tilde{\epsilon}^{1}\ll 1$, and the parameter
$\tilde{s}$ measures the temporal growth rate of this perturbation.

We restrict our analysis to probing the linear stability and substitute eq. 50
and the Poisson equation into eq. 49. Next, we collect terms up to the first
order in the perturbation mode $\tilde{\epsilon}^{1}$, which yields a
dispersion relation for the growth rate of the perturbation,

$\tilde{s}(\tilde{k})=-1-\left(1-\frac{\mathscr{V}^{0}}{E_{\ce{th}}}+\frac{2}{\pi}\frac{\mathscr{V}^{0}}{E_{\ce{th}}}\frac{E_{\ce{el}}}{E_{\ce{th}}}(\tilde{k})^{2}\right)\cdot(i\tilde{k})^{2}.$
(52)

The uniform state is stable under perturbations if and only if $\tilde{s}<0$.
This defines an instability onset $\tilde{k}^{\ce{c}}$ for the wave-numbers

$\tilde{k}^{\ce{c}}_{1,2}=\pm\frac{1}{2}\sqrt{\pi(\mathscr{V}^{0}E_{\ce{th}}-E_{\ce{th}}^{2})/E_{\ce{el}}\mathscr{V}^{0}}.$
(53)

The corresponding stability criterion
$\tilde{s}(\tilde{k}^{\ce{c}}_{1,2}){<}0$ determines the phase boundary at
which the bulk of the IL-electrolyte becomes unstable. This stability
threshold exactly equals the phase boundary between the damped oscillatory
phase and the nano-segregated phase (see eq. 45),

$\mathscr{V}^{0}_{\ce{+}}=E_{\ce{th}}+4E_{\ce{el}}/\pi+2\sqrt{2E_{\ce{el}}(2E_{\ce{el}}/\pi+E_{\ce{th}})/\pi}.$
(54)

Thus, for interaction energies $\mathscr{V}^{0}{>}\mathscr{V}^{0}_{+}$ the
bulk state of the system becomes unstable and phase separation emerges. The
initial cause for the structure-formation can be driven by external agents, or
boundary conditions, e.g., by the application of an electric potential to an
IL/electrode interface.

This stability analysis complements the static analysis, and rationalizes the
emergence of phase separation into ionic layers occurring at interaction
energies above $\mathscr{V}^{0}_{+}$.

###  IV.3 Validation With Simulation

Our goal in this section is to compare the results of our asymptotic analysis
(see eqs. 41, 42 and 43 in section IV.1.2), with numeric simulations of the
completely-coupled system subject to the two theoretical descriptions
(integral description eqs. 18 and 19, and gradient description eqs. 18 and
23).

We start our numerical investigations with an overview of the screening
profiles for the charge density at different values $\mathscr{V}^{0}$, as
obtained from the integral description eqs. 18 and 19 (see fig. 6). Next, we
compare our different EDL-descriptions in detail for two different energies
$\mathscr{V}^{0}$, see figs. 7 and 8. Finally, we generalize these exemplary
findings via a systematic study over the complete phase-space of interaction
energies, see fig. 9. This provides a clear illustration of the complete set
of phase-transitions which the system undergoes, and highlights the
consistency between the three descriptions.

All simulations were performed for a symmetric cell set-up, where the IL-
electrolyte is located within two oppositely charged, blocking interfaces
separated by a distance of $L_{\ce{cell}}=$60\text{\,}\mathrm{nm}$$. The
electrode on the left side is negatively charged with
$\varDelta\phi=$-100\text{\,}\mathrm{mV}$$, whereas, on the right side, the
electrode is positively charged with
$\varDelta\phi=$100\text{\,}\mathrm{mV}$$. Since charge saturation begins
roughly at $\varDelta\phi=$70\text{\,}\mathrm{mV}$$, see figs. 2 and 3, the
charge distribution can safely be assumed saturated adjacent to the interfaces
(i.e. $\lvert\tilde{\varrho}\rvert=1$). The electrolyte is considered to
consist of symmetric ions ($\gammaup_{\pm}=$0.5$$) of size
$a=$1.3\text{\,}\mathrm{nm}$$. Hence, the cell-geometry allows a “maximal”
number of roughly $90$ ions. In addition, we assume room temperature,
$T=$300\text{\,}\mathrm{K}$$, and $\varepsilon_{\ce{R}}=$15$$. The phase
boundaries corresponding to these parameters as predicted by our analytic
description are $\mathscr{V}^{0}_{-}=$3\text{\,}\mathrm{meV}$$ and
$\mathscr{V}^{0}_{+}=$253\text{\,}\mathrm{meV}$$ (see eq. 45).

Figure 6: Screening profiles of the charge density $\tilde{\varrho}$, as
obtained from numerical simulations of the integral description (see eqs. 18
and 19) for different values $\mathscr{V}^{0}$. The $y$-axis is scaled from
$-1$ to $1$, where $\lvert\tilde{\varrho}\rvert=1$ corresponds to charge
saturation.

Figure 6 shows the numerical results of the charge density for the integral
description (eqs. 18 and 19), where $\mathscr{V}^{0}$ takes values across two
orders of magnitude. First, at $\mathscr{V}^{0}=$1\text{\,}\mathrm{meV}$$ the
profile shows charge saturation near the two electrified electrodes,
$\tilde{\varrho}(x{=}0)=1$ and $\tilde{\varrho}(x{=}L_{\ce{cell}})=-1$. Near
both electrodes, the profile decays exponentially towards the electroneutral
bulk ($\tilde{\varrho}=0$). This corresponds to the profiles which we
discussed in great detail in section III.1. Because
$\mathscr{V}^{0}<\mathscr{V}^{0}_{-}$, this is in accordance with the analytic
prediction. The next two profiles show results for interaction energies within
the intermediate phase,
$\mathscr{V}^{0}_{-}<\mathscr{V}^{0}<\mathscr{V}^{0}_{+}$. Both simulations
show damped oscillatory profiles, where the long-ranged oscillatory profiles
span many nanometers. Apparently, the oscillations in the profile for
$\mathscr{V}^{0}=$196\text{\,}\mathrm{meV}$$ extend almost across the entire
cell. A slight increase by $2\text{\,}\mathrm{meV}$ to
$\mathscr{V}^{0}=$198\text{\,}\mathrm{meV}$$ causes the profile to transition
into a crystalline phase with undamped oscillatory shape. Note that the
amplitudes between the electrodes are smaller than unity, i.e. the bulk region
consists of mixed ion layers with one dominant ion species. An increase to
$\mathscr{V}^{0}=$204\text{\,}\mathrm{meV}$$ enhances the amplitudes of the
oscillations further, i.e. enhances segregation of ion-species. The last plot
shows the corresponding profile for a significantly enhanced interaction
energy ($\mathscr{V}^{0}=$500\text{\,}\mathrm{meV}$$). Here, the amplitudes of
the oscillations have reached saturation ($\lvert\tilde{\varrho}|=1$), and the
electrolyte has transitioned into a crystalline phase consisting of
alternating pure ion layers. In LABEL:fig:detail_crystalline_phase (see
LABEL:sec:SI-crystalline-limit), we highlight that the ionic layers coincide
exactly with the ion size $a$. Thus, with increasing energies
$\mathscr{V}^{0}$, the interfacial structure increases into the bulk
electrolyte, until the bulk itself gets nanostrucured by the layering of the
ion species. This phase transition occurs rapidly within a few
$\text{\,}\mathrm{meV}$.

Apparently, the numerical results for the integral description confirm the
existence of three different screening phases. However, quantitative
deviations between our descriptions are present. As we show in the SI (see
LABEL:fig:SI-peak-numbers) the phase transitions from exponential decay to
damped oscillations occurs roughly at
$\mathscr{V}^{0}=$2\text{\,}\mathrm{meV}$$. In addition, as can be inferred
from fig. 6, the transition from damped oscillations to undamped oscillations
appears at $\mathscr{V}^{0}=$200\text{\,}\mathrm{meV}$$. Hence, both phase
boundaries are slightly shifted to smaller values, compared to the analytical
predictions $\mathscr{V}^{0}_{-}=$3\text{\,}\mathrm{meV}$$ and
$\mathscr{V}^{0}_{+}=$253\text{\,}\mathrm{meV}$$ (see eq. 45). Thus, the
analytical prediction, which is based on the gradient description, slightly
underestimates the influence of $\mathscr{V}^{0}$, when compared with
$E_{\ce{th}}$ and $E_{\ce{el}}$. This can be attributed to the fact that the
gradient description is an approximation based on only the first two
perturbation modes, whereas the integral description comprises all modes.

Figure 7: Results for electric potential and charge density as obtained from
numerical simulations of the integral description (eqs. 18 and 19), and of the
gradient description (eqs. 18 and 23), and as predicted by the analytic
predictions (see eqs. 41, 42 and 43) at
$\mathscr{V}^{0}=$180\text{\,}\mathrm{meV}$$.

Next, we give a quantitative comparison between the numerical results of the
two theoretical descriptions, and the analytical predictions. Here, we
restrict our discussion to the interaction energy
$\mathscr{V}^{0}=$180\text{\,}\mathrm{meV}$$, i.e. the intermediate phase of
damped oscillations. Figure 7 shows the profiles for the charge density and
electrolyte electric potential as obtained from the numeric simulations, and
as predicted by the analytic description for the first few nanometers of the
left half-cell. Figure 7a illustrates the charge distribution adjacent to the
negatively charged electrode. The dashed blue line shows the screening profile
obtained from the gradient description, which exhibits a damped oscillatory
shape. This confirms the analytic prediction for this interaction energy. The
dashed yellow line shows the resulting analytic profile. Note that the
analytic prediction in section IV does not capture charge-saturation, but only
determines the damping-parameter and the oscillation-frequency of the
screening-profile. However, in section III.2.2, we derived an analytic
prediction for the saturation-width $L_{\ce{EDL}}$, valid close to the
interface (see eq. 37). Hence, in order to reconstruct the "complete" profile,
we supplement the contribution emerging from the bulk eq. 43, valid far away
from the interface, by constant charge-saturation $\tilde{\varrho}=1$ spanning
over the width $L_{\ce{EDL}}$. Apparently, the analytic and numeric results of
the gradient description are quantitatively in very good agreement. Finally,
the solid red line in fig. 7a shows the numerical results for the integral
description. In accordance with the results shown in fig. 6, these results
reproduce the analytically predicted screening phase, but the oscillations are
more pronounced. Hence, the influence of the interaction energy
$\mathscr{V}^{0}$ is more dominant in the integral description than in the
gradient description. Next, in fig. 7b) we show the profiles for the
normalized electrolyte electric potential. The dashed blue line shows the
profile due to the gradient description. It is in accordance with the charge-
profile shown in fig. 7a, see also eqs. 28 and 27. Again, we reconstruct the
analytic profile by supplementing the profile eq. 36, valid close to the
interface, by the profile eq. 43, valid towards the electroneutral bulk.
Apparently, the analytic results are quantitatively in very good agreement
with the results stemming from the gradient description. The red line shows
the profile as obtained from the integral description. Like in fig. 7a for the
charge density, the oscillations are slightly enhanced when compared with the
gradient description. In addition, the brown solid lines show the analytic
envelopes for the screening. Apparently, it captures both numerical results
qualitatively very well. Interestingly, the differences in electrolyte
potential between the three descriptions depend on electrode potential near
the electrode as shown in LABEL:fig:figure_SI_phidependence_overscreening (see
LABEL:sec:SI_phi_dependence_overscreening). However, the qualitative agreement
is independent from the boundary conditions.

Figure 8: Screening profile for the charge density $\tilde{\varrho}$ obtained
from numerical simulations with respect to the gradient description (eqs. 18
and 23), and according to the analytic description (see eqs. 41, 42 and 43),
at $\mathscr{V}^{0}=$250\text{\,}\mathrm{meV}$$.

In fig. 8, we show results for the charge distribution at enhanced interaction
energy $\mathscr{V}^{0}=$250\text{\,}\mathrm{meV}$$, i.e., close to the phase
boundary $\mathscr{V}^{0}_{+}=$253\text{\,}\mathrm{meV}$$. As can be inferred
from fig. 6, the integral description has already transitioned to the phase of
undamped oscillations for this interaction energy. Hence, we only show the
screening profile as obtained from numerical simulation of the gradient
description (solid blue line), and compare it with the analytic prediction
(yellow line). Apparently, in accordance with the analytic prediction, the
numerical profile has a damped oscillatory shape, where the oscillations
extend over roughly $40$ ion sizes. This highlights the influence of the
enhanced interaction energy, see also fig. 7. Overall, the analytic profile
shown here is in nice agreement with the numeric results.

Figure 9: Meta analysis of the interfacial profiles for some thousand
simulations. The dashed vertical yellow lines show the phase boundaries
$\mathscr{V}^{0}_{\pm}$ (see eq. 45). The dashed and solid red / blue lines
show the peak-variance of the complete set of simulations as defined by eq.
55, with respect to the integral description (see eqs. 18 and 19), and with
respect to the gradient-description (see eqs. 18 and 23). The left inset shows
the onset of the oscillations at small interaction energies. The right inset
shows the variance in a non-logarithmic setting, which highlights the
occurrence of phase-transitions.

Finally, we conduct a quantitative comparison between the two EDL-discriptions
and the analytic description across multiple orders of magnitude of
$\mathscr{V}^{0}$. To address this goal, we examine simulation results of
roughly 4000 EDL-simulations across the parameter range from
$0.1\text{\,}\mathrm{meV}$ up to $500\text{\,}\mathrm{meV}$. As above, we
apply $\Delta\phi=$\pm 100\text{\,}\mathrm{meV}$$ at the electrodes such that
we can safely assume charge saturation near the interfaces.

We evaluate the simulation results by extracting two characteristic
properties. First, we count the number of peaks
$\textnormal{N}_{\ce{peaks}}(\mathscr{V}^{0})$ appearing in each screening
profile. Due to charge saturation, a minimal number of two peaks always
occurs. At most, roughly 90 ion-layers fit into the cell geometry of length
$L_{\ce{cell}}=$60\text{\,}\mathrm{nm}$$. We present the number of peaks
occurring in the full cell as function of the interaction strengths in the SI
(see LABEL:fig:SI-peak-numbers in LABEL:sec:simulation-results).

However, beyond the number of peaks, we want to evaluate the peak amplitudes,
too. For this purpose, we investigate the peak variance
$\sigma(\mathscr{V}^{0})$ of the left half-cell, defined by

$\sigma^{2}=\frac{\sum_{i=1}^{\textnormal{N}_{\ce{peaks}}}\tilde{\varrho}_{i}\cdot\left(x_{i}\right)^{2}}{\sum_{j}^{\textnormal{N}_{\ce{peaks}}}\tilde{\varrho}_{j}}.$
(55)

Here, $x_{i}$ is the discrete location of the $i$-th peak
$\tilde{\varrho}_{i}=\lvert\tilde{\varrho}(\tilde{x}_{i})\rvert$ (such that
$0\leq\tilde{\varrho}_{i}\leq 1$, where $\tilde{\varrho}_{i}=1$ corresponds to
a saturated peak, i.e. a pure ion-layer) appearing in the profile of the
charge density. In the SI (see LABEL:sec:SI_meta_anlss_peak_nmbrs), we show
analytically that $\sigma$ converges to $L_{\ce{cell}}/\sqrt{3}$ if the set of
simulation-energies comprises energies
$\mathscr{V}^{0}\gg\mathscr{V}^{0}_{+}$. For such interaction energies, the
bulk electrolyte has transitioned into a crystalline phase composed of nano-
segregated ion-layers (see fig. 6).

Figure 9 shows the results for variance $\sigma$ normalized to its maximum
$L_{\ce{cell}}/\sqrt{3}$ in logarithmic scale. In this figure, the vertical
dashed yellow lines indicate the phase-boundaries $\mathscr{V}^{0}_{\pm}$, as
predicted by the analytic description (see eq. 45). The left inset shows the
simulation-results for small values $\mathscr{V}^{0}$, and the right inset
comprises the overall results in a non-logarithmic representation highlighting
the transition. The blue dashed line shows the results for $\sigma$ according
to the gradient description (see eqs. 18 and 23). At small interaction
energies $\mathscr{V}^{0}<\mathscr{V}^{0}_{-}$, the variance is zero. This
corresponds to an exponentially damped, non-oscillatory screening profiles
(note that the only peak, due to charge saturation, is located at
$\tilde{x}_{i}=0$). The variance starts increasing exactly at
$\mathscr{V}^{0}_{-}$ (see also the left inset). This corresponds to an
increasing number of damped oscillations, where the amplitudes of the peaks
also increase with $\mathscr{V}^{0}$. Finally, at $\mathscr{V}^{0}_{+}$, the
variance converges to it’s constant limiting value $L_{\ce{cell}}/\sqrt{3}$
(see also the right inset). In this energy-regime, the bulk electrolyte
consists completely of pure ion-layers. Altogether, these results reproduce
exactly the phases as predicted by the analytic description.

The red curve shows the results for $\sigma$ according to the integral
description (see eqs. 18 and 19). In contrast to the gradient description, the
variance starts increasing from zero at roughly
$\mathscr{V}^{0}=$1\text{\,}\mathrm{meV}$$, i.e. before the analytically
predicted phase boundary $\mathscr{V}^{0}_{-}$ (see also the left inset).
Hence, the phase transition from exponentially damped screening-profiles to
damped oscillatory screening-profiles is slightly shifted to smaller energies.
Next, the variance increases exponentially up to roughly
$\mathscr{V}^{0}=$200\text{\,}\mathrm{meV}$$, above which it transitions into
the constant limiting value $L_{\ce{cell}}/3$. Altogether, the phase
boundaries of the integral description still exhibit a qualitatively good
agreement with the analytic predictions, although being slightly shifted to
smaller values. Apparently, this behaviour is due to the cumulative effect of
the integral-term in eq. 19, which comprises all interaction modes. In
contrast, we consider only the first two modes (n=0 and n=1) of the gradient
expansion in eq. 23.

##  V Multi-Scale Methodology

In this section, we highlight the relation of our model to theories on smaller
and larger length scales. We discuss in section V.1, based on basic concepts
from liquid state theory, how atomistic simulations can directly parametrize
our theory. Next, in section V.2 we sketch the phenomenologic BSK continuum
approach for the description of ILs near electrified interfaces and illustrate
its relation to our work. In addition, we state the relation of our framework
to AFM-experiments in the LABEL:sec:relat-exper-meth ( see also Ref.52).

###  V.1 From Molecular Dynamics to Non-Equilibrium Thermodynamics

Here, we explain how the parameters of our continuum theory can be rigorously
calculated with quantum chemistry, i.e., DFT and MD.

Ab-initio DFT calculations predict the forces between ions and molecules by
calculating their electronic structure. The DFT-generated force-fields are the
focal quantity for MD simulations, 63 which calculate the classical
trajectories of ions and molecules. Results from MD simulations are often
interpreted via profiles of the radial distribution function $g(r)$.

Liquid state theory,50 connects this atomistic description to thermodynamic
concepts and scattering experiments.64 On the one hand, the radial
distribution function allow a straightforward comparison with the structure
factor $S$ from scattering experiments.65, 66 On the other hand, the density
distribution function $g(r)$ can be used to calculate different correlation
functions. By subtracting its asymptotic value follows the so-called total
correlation function used in integral equation theories (IETs),
$h(r)=g(r)-1$.3 In IETs, the pairwise total correlation function $h$ relates
to the the direct correlation function $c^{(2)}$, used in classical density
functional theory (cDFT), via the Ornstein-Zernike relation,67

$h(r)=c^{(2)}+\uprho_{\ce{b}}\int\ce{d}r^{\prime}c^{(2)}(\lvert
r-r^{\prime}\rvert)\cdot h(r^{\prime}).$ (56)

In cDFT, the direct pair correlation functions $c^{(2)}_{\alpha\beta}$ account
for pairwise interactions between two ions of species $\alpha$ and $\beta$,
i.e., the excess free energy due to pairwise ion-interactions.68 Thus, they
can be obtained via the twofold functional derivative of $F^{\ce{int}}$,64
i.e., via our interaction potential $\mathscr{F}$ (see eq. 4),

$c^{(2)}_{\alpha\beta}(\lvert
r-r^{\prime}\rvert)=-\frac{1}{k_{\ce{B}}T(N_{\ce{A}})^{2}}\cdot\mathscr{F}_{\alpha\beta}(\lvert
r-r^{\prime}\rvert).$ (57)

To summarize, DFT determines force fields for MD, MD determines $g(r)$ for
liquid state theory, $g(r)$ determines $c^{(2)}$ via the Ornstein-Zernike
relation, $c^{(2)}$ determines $F^{\ce{int}}$ and generates our non-
equilibrium thermodynamic theory.

The dynamic properties of our theory can be determined from atomistic
simulations, too. These dynamic properties are encoded in the Onsager
coefficients,45 which can be measured experimentally.69 The Onsager
coefficients can be determined by MD simulations (“Green Kubo relations”).70,
71, 72

###  V.2 From Non-Equilibrium Thermodynamics to Phenomenologic BSK Theory

Now, we compare our thermodynamically consistent continuum approach with the
phenomenologic theory proposed by Bazant, Storey, and Kornyshev (BSK), a
seminal MFT-approach for ILs near electrified interfaces.35 In their continuum
model of the EDL, BSK incorporate ion-correlations using a modified linear
dielectric relation
$\bar{\boldsymbol{D}}=\hat{\bar{\varepsilon}}\boldsymbol{E}$ between the
electrostatic fields $\bar{\boldsymbol{D}}$ and $\boldsymbol{E}$, where
$\hat{\bar{\varepsilon}}=\varepsilon_{\ce{R}}\varepsilon_{0}(1-\ell_{\ce{c}}\nabla^{2})$
is their dielectric operator. The second order gradient term in
$\bar{\varrho}$ accounts for non-local ion-interactions, being effectively
short-ranged with correlation length $\ell_{\ce{c}}$. This Ansatz yields a
modified Poisson equation,
$\hat{\bar{\varepsilon}}\nabla^{2}{\varPhi}=-\bar{\varrho}$. The chemical
potential connects electric potential and charge density. Finally,

$\hat{\bar{\varepsilon}}\nabla^{2}{\varPhi}=\varPhi.$ (58)

holds in the limit of small potentials $\tilde{\varPhi}$.

Our model conceptually differs from BSK theory. Since we incorporate
electrostatic correlations in the free energy, non-local ion-interactions
enter the set of equations via the chemical potentials. This implies that the
MFT-quantities appearing in the BSK description, $\bar{\boldsymbol{D}}$ and
$\bar{\varrho}=\nabla\bar{\boldsymbol{D}}$, differ from the corresponding
quantities $\varrho$ and $\boldsymbol{D}$ appearing in our formalism. In
contrast to the "mean field charge density" $\bar{\varrho}$, the charge
density $\varrho$ relates to the "bulk"-quantity $\boldsymbol{D}$, which does
not incorporate ion correlations.

Despite these differences, the resulting model equations are very similar.
This can be seen as follows. The complete set of equations eqs. 28 and 29 for
the limit of small potentials can be cast into one equation for the electric
potential alone,

$\hat{\varepsilon}_{\ce{R}}\tilde{\nabla}^{2}\tilde{\varPhi}=\tilde{\varPhi},$
(59)

where the dielectric operator $\hat{\varepsilon}_{\ce{R}}$ is defined in eq.
30. Noting the conceptual similarity between the dielectric operators
$\hat{\varepsilon}_{\ce{R}}$ and $\hat{\bar{\varepsilon}}$, the similarity
between our model and BSK theory becomes apparent. In this way, we give
physical meaning to the correlation length $\ell_{\ce{c}}$ in BSK theory and
outline its calculation.

Finally, we emphasize that the higher-order gradient-terms, which are
phenomenologically incorporated in the BSK approach, emerge naturally within
our rigorous continuum model. In particular, they merely constitutes the
limiting case for small potentials of the more fundamental integral
formulation eq. 26. Furthermore, in contrast to the phenomenological BSK
model, our order-expansion comprises also a zero-order correction in the
dielectric operator, see eq. 30. This mode is mandatory to realize the
“complete” phase-space of interfacial profiles (see LABEL:sec:SItriv-order-
expans and LABEL:sec:SIline-order-expans in the SI).

###  V.3 Outlook

In this section, we discuss how our framework can be extended to account for
additional microscopic IL effects.

In this work, we have supplemented our bulk description for ILs and highly
correlated electrolytes, recently presented in Ref. 45, by non-local
interactions. Furthermore, we applied the resulting framework for the case of
short-ranged hardcore interactions. However, the generality of our framework
based on the modelling of the free energy, offers the possibility to
incorporate a wide range of non-local effects into our framework.

This includes properties like ion asymmetry, ion geometry, polarization, and
charge delocalization, which have a significant influence on the structure of
ILs near electrified interfaces.73, 62, 74, 75, 76 These effects result partly
from the relative orientation between the ions, which makes a one-dimensional
approach challenging. Nevertheless, assuming a highly symmetric set-up, the
one-dimensional description might still capture some basic consequences of
these effects.

Similar to the force fields used in atomistic simulations, the short-ranged
repulsive interaction can be supplemented by a longer-ranged attractive tail,
taking account for higher-order electrostatic effects of van-der-Waals type,
or for larger ions of complex geometry, _i.e._ , long alkyl chains.53, 54, 63
Also, by refining the short-ranged repulsive interaction potential, more
detailed models for the ion geometry and ion asymmetry can be incorporated
into our model. However, the strong influence of these microscopic properties
on the EDL structure may lead to some novel features within our framework. For
example, the three energy scales which determine the screening profile might
transition into field quantities which exhibit spatial variation. Also, the
phase space of screening profiles might become higher-dimensional, which can
lead to a more complex set of phase boundaries allowing for "mixed" screening
types.

Non-trivial polarization effects could be incorporated into our linear
constitutive model for the coupling between the electric field and the
dielectric displacement. This would result in a spatially varying dielectric
function $\varepsilon_{\ce{R}}(x)$ appearing in the electrostatic energy scale
$E_{\ce{el}}$ (see eq. 21) and a direct coupling of the chemical potentials
with the ion polarization. For small dielectric perturbations, we hypothesis
that the electrostatic energy scale gets more "diffuse", which has the effect
that the phase boundaries between the screening profiles wash out. Only in the
case of large dielectric variations, we expect the phase space of screening
profiles to be altered significantly.

Our dynamical theory offers the possibility to investigate transport processes
occurring in electrochemical devices, _e.g._ , the influence of EDL-charging
on the electrolyte performance, or the influence of the EDL structure on the
electrode transfer kinetics. However, electrochemical devices have some
characteristic properties which must be carefully taken account for when they
are modelled. For example, overlapping double layers in nano-porous electrodes
could be taken account for.77

In our description we assumed an ideally flat surface, which can be a bad
approximation for many electrochemical systems.78. The influence of interface
roughness on the EDL structure can be modelled by modifying the entropic
contributions in the free energy functional.79 In our analysis, this would
alter the thermal energy scale $E_{\ce{th}}$, see eq. 20. Depending upon the
surface morphology, this would enhance the disordering effect of the thermal
energy on the EDL structure. As result, the formation of crystalline phases
might become suppressed at rougher surfaces, similar to increasing the
temperature.

##  VI Conclusion

In this work, we complement our thermodynamically consistent continuum
framework for IL electrolytes by non-local molecular repulsion. Our integral
formulation can be determined by ab-initio MD simulations. Assuming short-
ranged interactions, we expand the interaction free energy in gradients of the
concentrations and adjust the dynamic equations for transport. The resulting
equations connect to the phenomenologic approach of BSK theory. We validate
our approach by simulations and find a remarkable agreement between the
different variants of our theory.

In this way, we develop a predictive multi-scale approach to the theory of ILs
at electrified interfaces. Atomistic density functional theory calculations
parameterize MD simulations, MD simulations yield an integral formulation for
molecular repulsion in our thermodynamic consistent transport theory, our
theory can be expanded to give the phenomenological BSK theory.

The expanded continuum approach allows to perform analytic asymptotic analysis
which creates deeper insights into parameter dependence of EDL structure as we
demonstrate for the example of binary ILs. First, we have neglected molecular
repulsion. We can analytically describe both limits, the dilute Debye limit,
where charge density is exponentially decaying, and the concentrated crowding
limit, where charge is saturated due to steric effects. Second, we have taken
into account molecular repulsion. We discuss the structure of the EDL
dependent on energy scales for thermal motion, molecular repulsion, and
electric Coulomb forces and find three different phases. For small
interactions, we recover the dilute Debye limit. For intermediate
interactions, a multi-layer structure of ions emerges which is washed out over
several atom layers. For very large interactions, the analysis predicts a
long-ranged, non-decaying crystalline order of the EDL. In simulations of our
full theory, we eventually observe charge ordering of quasi-crystalline multi-
layers in this case.

In summary, we have proposed a thermodynamic consistent description of ILs at
electrified interfaces that closes a gap in their multi-scale understanding.
This makes possible a predictive theoretical approach for tailoring ILs. We
proof that the inter-molecular forces determine the EDL structure of binary
ILs. Future works should extend the work to ternary mixtures of ILs and
incorporate the shape of molecules into the theory.

## Supporting Information

Transport Theory of Interacting Electrolytes (including Gradient Expansion,
Functional Derivative); Non-Dimensional Formulation; Binary Ionic Liquid
(including Chemical Forces, Symmetric Ion Species, Charge Saturation, Charge
Oscillations); Computational Details; Parameters Binary Ionic Liquid; Relation
to Experimental Methods; Simulation Results

##  VII Acknowledgment

This work was supported by the European Union’s Horizon 2020 research and
innovation program via the “Si-DRIVE” project (grant agreement No 814464).

The authors acknowledge support by the German Research Foundation (DFG)
through grant no INST 40/575-1 FUGG (JUSTUS 2 cluster) and the state of Baden-
Württemberg through bwHPC.

## References

  * Eisenberg [2013] B. Eisenberg, Biophysical journal 104, 1849 (2013).
  * Grosberg  _et al._ [2002] A. Y. Grosberg, T. Nguyen, and B. Shklovskii, Reviews of modern physics 74, 329 (2002).
  * Henderson [2009] D. Henderson, Interdisciplinary sciences: computational life sciences 1, 1 (2009).
  * Chen and Weeks [2006] Y.-G. Chen and J. D. Weeks, Proceedings of the National Academy of Sciences 103, 7560 (2006).
  * Bloomfield [1991] V. A. Bloomfield, Biopolymers: Original Research on Biomolecules 31, 1471 (1991).
  * Wong and Pollack [2010] G. C. Wong and L. Pollack, Annual review of physical chemistry 61, 171 (2010).
  * Eisenberg  _et al._ [2010] B. Eisenberg, Y. Hyon, and C. Liu, The Journal of Chemical Physics 133, 104104 (2010).
  * Fisher and Levin [1993] M. E. Fisher and Y. Levin, Physical review letters 71, 3826 (1993).
  * Levin and Fisher [1996] Y. Levin and M. E. Fisher, Physica A: Statistical Mechanics and its Applications 225, 164 (1996).
  * van Roij and Hansen [1997] R. van Roij and J.-P. Hansen, Physical review letters 79, 3082 (1997).
  * Hansen and Löwen [2000] J.-P. Hansen and H. Löwen, Annual Review of Physical Chemistry 51, 209 (2000).
  * Levin [2002] Y. Levin, Reports on progress in physics 65, 1577 (2002).
  * Hückel [1924] E. Hückel, in  _Ergebnisse der exakten naturwissenschaften_ (Springer, 1924) pp. 199–276.
  * Kilic  _et al._ [2007] M. S. Kilic, M. Z. Bazant, and A. Ajdari, Phys. Rev. E 75, 021502 (2007).
  * Hayes  _et al._ [2015] R. Hayes, G. G. Warr, and R. Atkin, Chemical Reviews 115, 6357 (2015).
  * Bruce  _et al._ [2017] D. W. Bruce, C. P. Cabry, J. N. Canongia Lopes, M. L. Costen, L. D’Andrea, I. Grillo, B. C. Marshall, K. G. McKendrick, T. K. Minton, S. M. Purcell, S. Rogers, J. M. Slattery, K. Shimizu, E. Smoll, and M. A. Tesa-Serrate, The Journal of Physical Chemistry B 121, 6002 (2017).
  * Fedorov and Kornyshev [2014] M. V. Fedorov and A. A. Kornyshev, Chemical reviews 114, 2978 (2014).
  * Plechkova and Seddon [2008] N. V. Plechkova and K. R. Seddon, Chemical Society Reviews 37, 123 (2008).
  * Wishart [2009] J. F. Wishart, Energy Environ. Sci. 2, 956 (2009).
  * Torimoto  _et al._ [2010] T. Torimoto, T. Tsuda, K.-i. Okazaki, and S. Kuwabata, Advanced Materials 22, 1196 (2010).
  * Werner  _et al._ [2010] S. Werner, M. Haumann, and P. Wasserscheid, Annual review of chemical and biomolecular engineering 1, 203 (2010).
  * Welton [1999] T. Welton, Chemical reviews 99, 2071 (1999).
  * Endres [2004] F. Endres, Zeitschrift für Physikalische Chemie 218, 255 (2004).
  * Armand  _et al._ [2011] M. Armand, F. Endres, D. R. MacFarlane, H. Ohno, and B. Scrosati, in  _Materials For Sustainable Energy_ (World Scientific, 2011) pp. 129–137.
  * Endres  _et al._ [2010] F. Endres, O. Höfft, N. Borisenko, L. H. Gasparotto, A. Prowald, R. Al-Salman, T. Carstens, R. Atkin, A. Bund, and S. Z. El Abedin, Physical Chemistry Chemical Physics 12, 1724 (2010).
  * MacFarlane  _et al._ [2010] D. R. MacFarlane, J. M. Pringle, P. C. Howlett, and M. Forsyth, Physical Chemistry Chemical Physics 12, 1659 (2010).
  * Wu  _et al._ [2011] J. Wu, T. Jiang, D.-e. Jiang, Z. Jin, and D. Henderson, Soft Matter 7, 11222 (2011).
  * Henderson and Wu [2012] D. Henderson and J. Wu, The journal of physical chemistry B 116, 2520 (2012).
  * Jiang  _et al._ [2014] J. Jiang, D. Cao, D.-e. Jiang, and J. Wu, Journal of Physics: Condensed Matter 26, 284102 (2014).
  * Sharma  _et al._ [2016] S. Sharma, A. Gupta, and H. K. Kashyap, The Journal of Physical Chemistry B 120, 3206 (2016).
  * Sharma and Kashyap [2015] S. Sharma and H. K. Kashyap, The Journal of Physical Chemistry C 119, 23955 (2015).
  * Hu  _et al._ [2013] Z. Hu, J. Vatamanu, O. Borodin, and D. Bedrov, Physical Chemistry Chemical Physics 15, 14234 (2013).
  * Bikerman [1942] J. Bikerman, The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science 33, 384 (1942).
  * Santangelo [2006] C. D. Santangelo, Physical Review E 73, 041512 (2006).
  * Bazant  _et al._ [2011] M. Z. Bazant, B. D. Storey, and A. A. Kornyshev, Phys. Rev. Lett. 106, 046102 (2011).
  * Yochelis [2014a] A. Yochelis, Physical Chemistry Chemical Physics 16, 2836 (2014a).
  * Yochelis [2014b] A. Yochelis, The Journal of Physical Chemistry C 118, 5716 (2014b).
  * Yochelis  _et al._ [2015] A. Yochelis, M. B. Singh, and I. Visoly-Fisher, Chemistry of Materials 27, 4169 (2015).
  * Gavish and Yochelis [2016] N. Gavish and A. Yochelis, The Journal of Physical Chemistry Letters 7, 1121 (2016).
  * Bier  _et al._ [2017] S. Bier, N. Gavish, H. Uecker, and A. Yochelis, Phys. Rev. E 95, 060201 (2017).
  * Lee  _et al._ [2015] A. A. Lee, S. Kondrat, D. Vella, and A. Goriely, Phys. Rev. Lett. 115, 106101 (2015).
  * Gavish  _et al._ [2018] N. Gavish, D. Elad, and A. Yochelis, The Journal of Physical Chemistry Letters 9, 36 (2018).
  * Bothe and Dreyer [2015] D. Bothe and W. Dreyer, Acta Mechanica 226, 1757 (2015).
  * Kovetz [2000] A. Kovetz,  _Electromagnetic theory_, Vol. 975 (Oxford University Press Oxford, 2000).
  * Schammer  _et al._ [2021] M. Schammer, B. Horstmann, and A. Latz, Journal of the Electrochemical Society 168, 026511 (2021).
  * Latz and Zausch [2015] A. Latz and J. Zausch, Beilstein journal of nanotechnology 6, 987 (2015).
  * Schmitt  _et al._ [2019] T. Schmitt, T. Arlt, I. Manke, A. Latz, and B. Horstmann, Journal of Power Sources 432, 119 (2019).
  * Becker-Steinberger  _et al._ [2020] K. Becker-Steinberger, S. Schardt, B. Horstmann, and A. Latz, arXiv , 2101.10294 (2020), arXiv:2101.10294 .
  * Single  _et al._ [2019] F. Single, B. Horstmann, and A. Latz, The Journal of Physical Chemistry C 123, 27327 (2019).
  * Hansen and McDonald [2006] J. Hansen and I. McDonald,  _Theory of Simple Liquids_ (Elsevier Science, 2006).
  * Landstorfer  _et al._ [2016] M. Landstorfer, C. Guhlke, and W. Dreyer, Electrochimica Acta 201, 187 (2016).
  * Hoffmann  _et al._ [2018] V. Hoffmann, G. Pulletikurthi, T. Carstens, A. Lahiri, A. Borodin, M. Schammer, B. Horstmann, A. Latz, and F. Endres, Phys. Chem. Chem. Phys. 20, 4760 (2018).
  * Bedrov  _et al._ [2019] D. Bedrov, J.-P. Piquemal, O. Borodin, A. D. MacKerell, B. Roux, and C. Schröder, Chemical Reviews 119, 7940 (2019).
  * Heinz  _et al._ [2013] H. Heinz, T.-J. Lin, R. Kishore Mishra, and F. S. Emami, Langmuir 29, 1754 (2013).
  * Eisenberg  _et al._ [2011] B. Eisenberg, Y. Hyon, and C. Liu, Communications in Mathematical Sciences 9, 459 (2011).
  * Lin and Eisenberg [2014] T.-C. Lin and B. Eisenberg, Communications in Mathematical Sciences 12, 149 (2014).
  * Lin and Eisenberg [2015] T.-C. Lin and B. Eisenberg, Nonlinearity 28, 2053 (2015).
  * Lee  _et al._ [2017] A. A. Lee, C. Perez-Martinez, A. M. Smith, and S. Perkin, Faraday Discussions 199, 239 (2017), 1701.08151 .
  * Rotenberg  _et al._ [2018] B. Rotenberg, O. Bernard, and J.-P. Hansen, Journal of Physics: Condensed Matter 30, 054005 (2018).
  * Karuppasamy  _et al._ [2020] K. Karuppasamy, J. Theerthagiri, D. Vikraman, C.-J. Yim, S. Hussain, R. Sharma, T. Maiyalagan, J. Qin, and H.-S. Kim, Polymers 12, 918 (2020).
  * Schmickler and Santos [2010] W. Schmickler and E. Santos,  _Interfacial electrochemistry_ (Springer Science & Business Media, 2010).
  * Li  _et al._ [2013] H. Li, F. Endres, and R. Atkin, Physical Chemistry Chemical Physics 15, 14624 (2013).
  * González, M.A. [2011] González, M.A., JDN 12, 169 (2011).
  * Beijeren and Felderhof [1979] H. V. Beijeren and B. Felderhof, Molecular Physics 38, 1179 (1979).
  * Billinge [2019] S. J. L. Billinge, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 377, 20180413 (2019).
  * Kats and Muratov [2018] E. I. Kats and A. R. Muratov, Phys. Rev. E 97, 012610 (2018).
  * Archer  _et al._ [2017] A. J. Archer, B. Chacko, and R. Evans, The Journal of Chemical Physics 147, 034501 (2017).
  * Cao  _et al._ [2019] S. Cao, K. A. Konovalov, I. C. Unarta, and X. Huang, Advanced Theory and Simulations 2, 1900049 (2019).
  * Fong  _et al._ [2021] K. D. Fong, J. Self, B. D. McCloskey, and K. A. Persson, Macromolecules 54, 2575 (2021).
  * Fong  _et al._ [2020] K. D. Fong, H. K. Bergstrom, B. D. McCloskey, and K. K. Mandadapu, AIChE Journal 66, e17091 (2020).
  * Götze and Latz [1989] W. Götze and A. Latz, Journal of Physics: Condensed Matter 1, 4169 (1989).
  * Latz [1991] A. Latz, Verallgemeinerte konstituierende gleichungen und formfaktoren für einfache glasbildner. ph.d. dissertation, tu münchen (1991).
  * de Souza  _et al._ [2022] J. P. de Souza, K. Pivnic, M. Z. Bazant, M. Urbakh, and A. A. Kornyshev, The Journal of Physical Chemistry B (2022).
  * Wei  _et al._ [2021] C. Wei, K. Jiang, T. Fang, and X. Liu, Green Chemical Engineering 2, 402 (2021).
  * Izgorodina  _et al._ [2009] E. I. Izgorodina, M. Forsyth, and D. R. MacFarlane, Physical chemistry chemical physics 11, 2452 (2009).
  * Paek  _et al._ [2015] E. Paek, A. J. Pak, and G. S. Hwang, The Journal of Chemical Physics 142, 024701 (2015).
  * Lian  _et al._ [2021] Z. Lian, H. Chao, and Z.-G. Wang, ACS nano 15, 11724 (2021).
  * Jänsch  _et al._ [2015] T. Jänsch, J. Wallauer, and B. Roling, The Journal of Physical Chemistry C 119, 4620 (2015).
  * Aslyamov  _et al._ [2021] T. Aslyamov, K. Sinkov, and I. Akhatov, Physical Review E 103, L060102 (2021).