Phase diagram of a bulk 1d lattice Coulomb gas

V. Démery, R. Monsarrat, David S. Dean, R. Podgornik

► To cite this version:

V. Démery, R. Monsarrat, David S. Dean, R. Podgornik. Phase diagram of a bulk 1d lattice Coulomb gas. EPL - Europhysics Letters, 2016, 113 (1), pp.18008 (1-6). 10.1209/0295-5075/113/18008. hal-01284082

HAL Id: hal-01284082

https://hal.science/hal-01284082

Submitted on 7 Mar 2016

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution - ShareAlike 4.0 International License
Phase diagram of a bulk 1d lattice Coulomb gas

V. DÉMERY¹, R. MONSARRAT¹, D. S. DEAN² and R. PODGORNIK³,⁴

¹ PCT, UMR Gulliver 7083, ESPCI Paris and CNRS - 10 rue Vauquelin, F-75005 Paris, France

² Université de Bordeaux and CNRS, Laboratoire Ondes et Matière d'Aquitaine (LOMA), UMR 5798 F-33400 Talence, France

³ Department of Physics, Faculty of Mathematics and Physics, University of Ljubljana - SI-1000 Ljubljana, Slovenia

⁴ Department of Theoretical Physics, J. Stefan Institute - SI-1000 Ljubljana, Slovenia

PACS 82.45.Gj - Electrolytes

PACS 05.20.Jj - Statistical mechanics of classical fluids

PACS 05.70.Fh - Phase transitions: general studies

**Abstract** – The exact solution, via transfer matrix, of the simple one-dimensional lattice Coulomb gas (1d LCG) model can reproduce peculiar features of ionic liquid capacitors, such as overscreening, layering, and camel- and bell-shaped capacitance curves. Using the same transfer matrix method, we now compute the bulk properties of the 1d LCG in the constant voltage ensemble. We unveil a phase diagram with rich structure exhibiting low-density disordered and high-density ordered phases, separated by a first-order phase transition at low temperature; the solid state at full packing can be ordered or not, depending on the temperature. This phase diagram, which is strikingly similar to its three-dimensional counterpart, also sheds light on the behaviour of the confined system.

Ionic liquids have intriguing properties when used in capacitors. These include self-organisation in alternating positively and negatively charged layers close to a charged surface [1–4], overscreening of the surface charge [5] and camel- and bell-shaped capacitance curves depending on temperature [6]. The first attempt to explain these features was based on a model with electrostatic interactions and steric repulsion treated at the mean-field level. This treatment successfully predicted the camel- and bell-shaped capacitance curves [7,8], however additional terms are required to explain the layering [9]. Later, a 1d lattice model was solved exactly, showing that electrostatic interaction and steric repulsion alone can qualitatively reproduce the observed features [10,11]. A continuum version of the one-dimensional model, modeling the ions as charged rods, showed a better qualitative agreement with experimental electrode separation-pressure curves [12]. However, in this case, an exact solution was not found, and a numerical computation is only possible for small systems. All these studies of one-dimensional systems focused on the capacitor behaviour while their bulk properties remained unexamined. We now extend this analysis to the bulk system with the rationale that the understanding of bulk properties may potentially provide valuable insight into the behaviour of small systems.

The basic ingredients of ionic liquids, electrostatic interaction and steric repulsion, are incorporated in the restricted primitive model (charged hard spheres), the phase diagram of which has been extensively studied and is now well established [13–20]. It features a phase coexistence between a low-density gas and a high-density ordered liquid at low temperature, a Néel line of second-order transitions at larger temperatures, and a solid phase at high density. The lattice version of this model (LRPM) implements the steric exclusion by allowing at most one ion per site. Its phase diagram has the same structure as its continuous counterpart, but also exhibits significant differences stemming from the fact that the lattice facilitates organisation [20–29]. The bulk properties of RPM and LRPM are thus well understood, but the link with the peculiar behaviour of an ionic liquid capacitor remains to be made.

Here, we use the exact transfer matrix solution of the one-dimensional LRPM (or one-dimensional lattice Coulomb gas) developed in refs. [10,11] in the context of a fluid capacitor to compute its bulk phase diagram. We show that, remarkably, this simple model produces a rich phase diagram showing striking similarities to its 3d counterpart. Specifically, we find that the charge density oscillations observed close to a charged surface correspond
<figure><img src="image_2.png" /><figcaption>Fig. 1: (Colour online) Notations for the one-dimensional lattice Coulomb gas: sites can either be empty, host a cation (red) or an anion (blue). Electrodes bear charges $-Q$ and $Q$.</figcaption></figure>

to an ordered bulk phase, and the abrupt change of the capacitance at the point of zero charge [11] of an ionic liquid capacitor corresponds to a first-order bulk phase transition.

**Model.** – We begin by defining the one-dimensional lattice Coulomb gas model, shown in fig. 1. The ions are confined between an anode at $i_a = \lceil-(N+1)/2\rceil$ and a cathode at $i_c = \lceil(N+1)/2\rceil$; $N$ sites are thus available. Each site $i \in (i_a, i_c)$ can be empty ($s_i = 0$), host a cation ($s_i = 1$) or an anion ($s_i = -1$). In the constant charge ensemble, the cathode and anode bear a charge $Q \in \mathbb{R}_+$ and $-Q$, respectively. The grand partition function can be written as

$$ \Xi = \sum_{\langle s_i \rangle} \delta_{0, \sum_i s_i} \mu^{\sum_i |s_i|} \exp \left( \frac{1}{4T} \sum_{i,j} |i-j| s_i s_j \right), \quad (1) $$

where the index $i$ runs over $(i_a, i_c)$, except in the argument of the exponential where the sum includes the electrodes, with $s_{i_a} = -Q$ and $s_{i_c} = Q$. In what follows, we use a unit charge, unit lattice spacing, and unit Boltzman constant: the only free parameters are thus the fugacity $\mu$ and the temperature $T$. Electroneutrality is enforced via the Kronecker $\delta$. If the electrodes can exchange charges to maintain a potential difference $\Delta V$ between them, the relevant ensemble is the constant voltage ensemble. In this ensemble, the grand partition function is given by [11]

$$ \hat{\Xi}_{\Delta V} = \int_{-\infty}^{\infty} \exp(Q\Delta V) \Xi_Q dQ. \quad (2) $$

**Transfer matrix solution.** – The model defined above is exactly solvable using a Hubbard-Stratonovitch transformation, which allows the rewriting of the partition function via the transfer matrix formalism [10,11]. The constant charge partition function thus reads

$$ \Xi = \langle \psi_b | K^N | \psi_b \rangle, \quad (3) $$

where the vector $\psi_b$ is a function and the transfer matrix $K$ is a linear integral operator defined on functions in $\mathbb{R}$:

$$ \psi_{b,k} = \exp\left(-\frac{k^2}{4T}\right) \delta(k - Q) \quad (4) $$

$$ K_{kl} = \exp\left(-\frac{k^2+l^2}{4T}\right) \times (\delta(k-l) + \mu[\delta(k-l+1) + \delta(k-l-1)]). \quad (5) $$

The moments of the density can be computed using the cation (+) and anion (-) density operators defined as

$$ R_{\pm,kl} = \mu \exp \left( -\frac{k^2 + l^2}{4T} \right) \delta(k - l \pm 1). \quad (6) $$

The charge and occupancy operators are given by $R = R_+ - R_-$ and $R_o = R_+ + R_-$, respectively. In the following, we will be interested in the average density and the charge-charge correlation function, which are given by

$$ \langle |s_i| \rangle = \Xi^{-1} \langle \psi_b | K^{i-i_c-1} R_o K^{i_a-i-1} | \psi_b \rangle, \quad (7) $$

$$ \langle s_i s_j \rangle = \Xi^{-1} \langle \psi_b | K^{i-i_c-1} R K^{j-i-1} R K^{i_a-j-1} | \psi_b \rangle. \quad (8) $$

In the thermodynamic limit, $N \to \infty$, the maximal eigenvalue $\lambda^*$ of the transfer matrix and its associated eigenvector $\psi^*$ dominate, we can thus write $K^N \sim \lambda^{*N}|\psi^*\rangle\langle\psi^*|$ and it follows that:

$$ \Xi \sim \lambda^{*N} |\langle \psi_b | \psi^* \rangle|^2, \quad (9) $$

$$ \langle |s_i| \rangle \to \lambda^{*-1} \langle \psi^* | R_o | \psi^* \rangle = \rho, \quad (10) $$

$$ \langle s_i s_j \rangle \to \lambda^{*- (j-i+1)} \langle \psi^* | R K^{j-i-1} R | \psi^* \rangle = \rho c_{j-i}, \quad (11) $$

where we have defined the average occupancy $\rho$ and charge correlation function $c_i$.

The transfer matrix $K$ and density operators $R_\pm$ (eqs. (5), (6)) map the subspaces $\mathcal{E}_\theta = \{\psi_k, k \in \mathbb{Z} + \theta\}$, where $\theta \in [0, 1)$, onto themselves. Hence, for each $\theta$, a discrete (but still infinite) matrix can be associated with the transfer matrix restricted to the subspace $\mathcal{E}_\theta$, $K_\theta = K|_{\mathcal{E}_\theta}$. For constant charges $\pm Q$ on the electrodes, the boundary vector is given by eq. (4), and only the subspace $\mathcal{E}_{Q-\lfloor Q \rfloor}$ is explored. This means that the eigenvalue (and corresponding eigenvector) entering eqs. (9)–(11) is actually the maximal eigenvalue of $K_{Q-\lfloor Q \rfloor}$, which we denote $\lambda^*_{Q-\lfloor Q \rfloor}$ (it is positive from the Perron-Frobenius theorem). This dependence of the bulk behaviour of the system on the boundary conditions of a finite system is related to the $\theta$-vacua encountered in quantum field theory [30].

On the other hand, in the constant voltage ensemble, the partition function (eq. (2)) obviously involves all the possible boundary charges. As a consequence, the maximal eigenvalue $\lambda^*$ which enters eqs. (9)–(11) is the maximal eigenvalue of the full transfer matrix $K$. The charges $Q \in \mathbb{Z} + \theta$ such that $\lambda^*_{\theta} = \lambda^*$ are said to be “selected” (we also say that a value of $\theta$ is selected).

For each $\theta$, the transfer matrix $K_\theta$ has infinite dimension, so that a cut-off $k_c$ has to be introduced on its indices in order to diagonalise it numerically. The spectrum of $K_{\theta=0}$ for $T = 10$ and $\mu = 1$ is shown in fig. 2; increasing the number of modes beyond $k_c \approx 5$ increases the density of eigenvalues around 0, but does not change the maximal (and minimal) eigenvalues, neither the spectral gap between the first and second eigenvalues. This rapid convergence with the cutoff is due to the rapid decay of the transfer matrix coefficients with the indices (see eq. (5)). The accuracy of the evaluation is better than the standard
<figure><img src="image_3.png" /><figcaption>Fig. 2: (Colour online) Spectrum of the transfer matrix $K_{\theta=0}$ as a function of the numerical cutoff $k_c$ for $T = 10$, $\mu = 1$. The $2k_c + 1$ eigenvalues corresponding to each value of $k_c$ are plotted along the corresponding y-axis and are colour coded, blue denoting the smallest and red the largest. One clearly sees the relative insensitivity of $\lambda^*$ to $k_c$ for $k_c > 5$.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 3: (Colour online) Dependence of the maximal eigenvalue $\lambda_{\theta}^*$ on $\theta$ at $T = 0.25$ (top) and $T = 0.2$ (bottom) around a value of the fugacity $\mu$ where the selected $\theta$ goes from 0 to $1/2$ (top: $\mu \in \{1.122, 1.123, 1.12404664, 1.125, 1.126\}$; bottom: $\mu \in \{1.324, 1.325, 1.32559787, 1.326, 1.327\}$). For $T = 0.25$ the selected value of $\theta$ goes continuously from 0 to $1/2$ while the transition is sharp for $T = 0.2$.</figcaption></figure>

precision of Python, $10^{-15}$, for the cut-off $k_c = 25$ used throughout the manuscript.

The relations $\lambda_{\theta}^* = \lambda_{-\theta}^*$ and $\lambda_{\theta}^* = \lambda_{1+\theta}^*$ mean that $\theta = 0$ and $1/2$ are stationary points of $\lambda_{\theta}^*$ with respect to $\theta$. Numerically, we find that integer ($\theta = 0$) or half-integer ($\theta = 1/2$) boundary charges correspond to the maximal

<figure><img src="image_5.png" /><figcaption>Fig. 4: (Colour online) Phase diagram in the density-temperature plane. Blue dotted lines delineate regions of integer ($\theta = 0$) and half-integer ($\theta = 1/2$) selected charges. Solid blue lines delineate the coexistence regions. Red triangles denote the onset of charge oscillation, i.e., ordering. The solid phase at $\rho = 1$ is ordered for $\theta = 1/2$ (gray lines, os) and disordered for $\theta = 0$ (ds).</figcaption></figure>

value of $\lambda_{\theta}^*$ almost everywhere. Curves of $\lambda_{\theta}^*$ are shown for values of the fugacity close to a transition in fig. 3. Depending on the shape of these curves, the selected $\theta$ can either change continuously (e.g., for $T = 0.25$) or discontinuously (e.g., for $T = 0.2$); the latter corresponds, as we shall see, to a first-order phase transition. On the contrary, no such transition exists in the constant charge ensemble. We are interested in the bulk phase transitions that may take place in the system, so that we restrict our study to the constant voltage ensemble; in the thermodynamic limit, bulk averages do not depend on the imposed voltage (see eqs. (10), (11)). Note that the occurrence of the free energy minima at $\theta = 0$ for low-density systems and $\theta = 1/2$ for high-density systems can be partially explained by mapping the model here onto a discrete interface model [11].

**Density.** — When the selected charges switch discontinuously between integers and half-integers, as the fugacity increases, the density changes discontinuously (from $\rho_a$ to $\rho_b$). At imposed density, this discontinuity becomes a coexistence region for $\rho \in [\rho_a, \rho_b]$. Coexistence should be understood in the statistical sense: the statistical state of the system is a combination of low- and high-density states. In this region, the maximal eigenvalue $\lambda^*$ of the transfer matrix and the pressure $p = T \log(\lambda^*)$ do not depend on the density. The density-temperature phase diagram, shown in fig. 4, displays three coexistence regions. The main one, at low temperature, is reminiscent of the one present in the RPM [17]. The one at higher temperature along the same transition line and the one around $T \approx 1$, $\rho \approx 0.8$ are more exotic.

**Spatial organisation.** — Spatial organisation of the charge is encoded in the charge correlation function $c_i$. The correlation function is shown in fig. 5 for several densities and temperatures; it may decay monotonically or
<figure><img src="image_6.png" /><figcaption>Fig. 5: (Colour online) Charge correlation function as a function of the distance. Top: $T = 1$ and $\rho \in \{0.1, 0.7, 0.9\}$; inset: focus on $\rho \in \{0.7, 0.9\}$ at large distance. Bottom: $T = 0.1$ and $\rho \in \{0.1, 0.9, 0.99\}$.</figcaption></figure>

with oscillations. We focus on the long-distance decay to understand these two behaviours.

The long-distance behaviour of the correlation function can be investigated by keeping only the two maximal eigenvalues of $K$ (in absolute value), $\lambda^*$ and $\lambda_2$, corresponding to the eigenvectors $\psi^*$ and $\psi_2$, respectively:

$$c_i \sim \frac{\langle \psi^* | R (\lambda^{*i-1} |\psi^*) \langle \psi^* | + \lambda_2^{i-1} |\psi_2\rangle \langle \psi_2 | \rangle |\psi^* \rangle}{\lambda^{*i+1} \rho} \quad (12)$$

The system is neutral, $0 = \langle s_0 \rangle = \langle \psi^* | R | \psi^* \rangle / \lambda^*$, and $R$ is antisymmetric, so that

$$c_i \sim \frac{|\langle \psi_2 | R | \psi^* \rangle|^2}{\lambda^* \rho} \left( \frac{\lambda_2}{\lambda^*} \right)^{i-1} \quad (13)$$

Note that this relation breaks down if $\langle \psi_2 | R | \psi^* \rangle = 0$; in this case the largest eigenvalue such that the corresponding eigenvector $\psi_{\text{eig}}$ satisfies $\langle \psi^* | R | \psi_{\text{eig}} \rangle \neq 0$ should be considered. If $\lambda_2 < 0$ (i.e., $\lambda_2$ is the minimal eigenvalue of $K$) and $\langle \psi_2 | R | \psi^* \rangle \neq 0$, the correlation function oscillates. We find numerically that this last condition is never satisfied in the phase $\theta = 0$, but is satisfied in the phase $\theta \neq 1/2$ at large enough density. The charge oscillations are seen in fig. 5 for $T = 1$ and $\rho = 0.7$ ($\theta = 1/2$ for this density, while $\theta = 0$ for $\rho \in \{0.1, 0.9\}$) and for $T = 0.1$ and $\rho \in \{0.9, 0.99\}$. We denote the onset of charge oscillation, and thus organisation, by red triangles

<figure><img src="image_7.png" /><figcaption>Fig. 6: (Colour online) Long-range charge order as a function of temperature for $\rho = 1$. The inset is a zoom on the $T > 1$ region, where a small long-range order exists.</figcaption></figure>

in the phase diagram, fig. 4. The correlation length is given by $\xi = \log(|\lambda^*/\lambda_2|)^{-1}$, which appears to diverge as the density approaches 1 (see fig. 5 for $T = 0.1$). This is investigated later using an expansion around $\rho = 1$.

**Dense packing limit.** — The dense packing limit $\rho \to 1$ corresponds to $\mu \to \infty$. In this limit, some analytical results may be obtained from the properties of $\bar{K} = K/\mu$ and $\bar{R}_\pm = R_\pm/\mu$; the eigenvalues now refer to eigenvalues of $\bar{K}$. First, we focus on the case of dense packing, i.e., $\rho = 1$, where the system is in a solid state. From the eigenvector $\psi^*$, we construct the vector $\psi_-$ with components $(\psi_-)_k = (-1)^{\lfloor k \rfloor} \psi_k^*$. From $\bar{K}_{kl}(-1)^{\lfloor l \rfloor} = -(\bar{1})^{\lfloor k \rfloor} \bar{K}_{kl}$, we deduce that $\bar{K} |\psi_-\rangle = -\lambda^* |\psi_-\rangle$: $-\lambda^*$ is also an eigenvalue of $\bar{K}$. Hence $\lambda_2 = -\lambda^*$ in eq. (13), and we get

$$c_i \sim \frac{|\langle \psi_- | \bar{R} | \psi^* \rangle|^2}{\lambda^* \rho} (-1)^i \quad (14)$$

The parameter $\omega(T) = \lambda^{-2} |\langle \psi_- | \bar{R} | \psi^* \rangle|^2$ characterises the long-range charge order. It depends only on temperature. If $\omega(T) > 0$, the system is an ordered solid. The long-range charge order is plotted in fig. 6: it is zero for $\theta = 0$ and non-zero for $\theta = 1/2$, although it decays rapidly with increasing temperature; the ordered and disordered solid regions are shown in the phase diagram, fig. 4. While ordered and disordered solid regions exist in the RPM [15,17,19], the peculiar dependence on the temperature is specific to the one-dimensional system considered here. Finally, we note that the limits $N \to \infty$ and $\rho \to 1$ do not commute. Taking $\rho \to 1$ first, $-\lambda^*$ should be included in eqs. (9)–(11). Here, we take the thermodynamic limit at $\rho < 1$ and then take the limit $\rho \to 1$.

Next, we show that the long-range charge order is necessarily zero for $\theta = 0$. For $\theta = 0$ and $\theta = 1/2$, the transfer matrix $\bar{K}_{kl}$ is even (i.e., $\bar{K}_{-k,-l} = \bar{K}_{k,l}$), so that the maximal eigenvector $\psi_k^*$ is also even. On the other hand, the function $k \mapsto (-1)^{\lfloor k \rfloor}$ is even for $\theta = 0$ and odd for $\theta = 1/2$, and so is $(\psi_-)_k$. Finally, the operator $\bar{R}$ is
odd: $\bar{R}_{-k,-l} = -\bar{R}_{k,l}$. Putting these relations together im-
plies that $\langle \psi_- | \bar{R} | \psi^* \rangle = 0$ for $\theta = 0$: there is no long-range
order and the system is in the disordered solid phase. On
the other hand, long-range order can occur for $\theta = 1/2$,
and more generally for $\theta \neq 0$.

We have shown that the correlation length is infinite
in the regions where $\theta = 1/2$ for $\rho = 1$. Here, we use an
expansion around $\rho = 1$ to determine its dependence on $\rho$.
We expand the operators, eigenvalues and eigenvectors as
a function of the small parameter $\epsilon = 1/\mu$. The occupancy
operator is $\bar{R}_o(\epsilon) = \bar{K}(\epsilon = 0)$, and the density is given by

$$
\rho(\epsilon) = \frac{\langle \psi^*(\epsilon) | \bar{K}(0) | \psi^*(\epsilon) \rangle}{\lambda^*(\epsilon)} = 1 - \epsilon \frac{\lambda^{*\prime}(0)}{\lambda^*(0)} + \mathcal{O}(\epsilon^2). \quad (15)
$$

On the other hand, the correlation length satisfies the
relation

$$
\xi(\epsilon)^{-1} = \log \left( - \frac{\lambda^*(\epsilon)}{\lambda_2^*(\epsilon)} \right) = \epsilon \frac{\lambda^{*\prime}(0) + \lambda_2^*(0)}{\lambda^*(0)} + \mathcal{O}(\epsilon^2), \quad (16)
$$

where we have used $\lambda_2(0) = -\lambda^*(0)$. The variation of an
eigenvalue, e.g. $\lambda_2$, is given by $\lambda'_2 = \langle \psi_2 | \bar{K}' | \psi_2 \rangle$. Using
that $[\psi_2(\epsilon=0)]_k = (-1)^{\lfloor k \rfloor} [\psi^*(\epsilon=0)]_k$ and $\bar{K}'(0) =
\exp(-k^2/[2T])\delta(k-l)$, we get that $\lambda'_2(0) = \lambda^{*\prime}(0)$. Insert-
ing this equality in eq. (16) and comparing with eq. (15)
leads to

$$
\xi \underset{\rho \to 1}{\sim} \frac{1}{2} (1 - \rho)^{-1}, \qquad (17)
$$

the correlation length is therefore half of the average dis-
tance between vacancies.

**Discussion.** – In this letter, we have computed ex-
actly the bulk phase diagram of the 1d LCG (or LRPM)
in the constant voltage ensemble (fig. 4). It features a
low-density disordered phase and a high-density ordered
phase, separated by a first-order phase transition at low
temperatures. The solid state at full occupancy may be
charge ordered or charge disordered, depending on the
temperature.

We note that the presence of a phase transition in
our one-dimensional system does not violate the Landau-
Peierls argument since the electrostatic interaction is long
range (or, more precisely, does not decay faster than
$1/x^2$) [31,32]). In the transfer matrix formalism, this ar-
gument is rephrased in the van Hove theorem [33]. For
the theorem to apply, the transfer matrix should be irre-
ducible, which is not the case here because the transfer
matrix can be decomposed on the subspaces $\mathcal{E}_\theta$.

Some of the bulk properties observed here may be re-
lated to the capacitor behaviour of the 1d LCG. Firstly,
the discontinuous jump of the capacitance at the point of
zero charge (PZC), noticed in ref. [11], which corresponds
to the transition between camel- and bell-shaped capaci-
tance curves [5,8], is clearly related to the first-order phase
transition that we observe. In the low-density phase, inte-
ger charges are selected, so that the PZC lies in a plateau
of the voltage-charge curve and the differential capacitance
is small. In the high-density phase, half-integer charges

are selected and the PZC lies at the transition between
two plateaus, hence the differential capacitance is large.
Secondly, charge density oscillations were observed in the
constant charge ensemble for non-integer surface charges.
Spatial organisation of the charge at an interface thus cor-
responds directly to charge correlation in the bulk; our
results suggest that charge density oscillations seen at a
charged interface are also present in the constant voltage
ensemble in the regions where half-integer charges are se-
lected in the bulk.

Whether such a close relationship between the bulk and
capacitor properties is also present in more realistic models
remains an open question. Up till now, numerical simu-
lations of ionic liquids at interfaces have not investigated
this question [6]. However, we expect the relationship to
be more subtle in three-dimensional systems than in our
one-dimensional model. For example, the observed layer-
ing where charge oscillates only in the direction orthogonal
to the bounding surface is not directly compatible with the
ordered state of the bulk liquid in three dimensions.

We thank A. C. MAGGS for useful discussions about the
properties of the transfer matrix and the van Hove theo-
rem. RP acknowledges support from the ARRS through
the program P1-0055 and from ESPCI ParisTech during
his stay as a visiting *Joliot Chair* Professor.

REFERENCES

[1] MEZGER M., SCHRÖDER H., REICHERT H., SCHRAMM S., OKASINSKI J. S., SCHRÖDER S., HONKIMÄKI V., DEUTSCH M., OCKO B. M., RALSTON J., ROHWERDER M., STRATMANN M. and DOSCH H., *Science*, **322** (2008) 424.

[2] HAYES R., BORISENKO N., TAM M. K., HOWLETT P. C., ENDRES F. and ATKIN R., *J. Phys. Chem. C*, **115** (2011) 6855.

[3] PERKIN S., CROWHURST L., NIEDERMEYER H., WELTON T., SMITH A. M. and GOSVAMI N. N., *Chem. Commun.*, **47** (2011) 6572.

[4] PERKIN S., *Phys. Chem. Chem. Phys.*, **14** (2012) 5052.

[5] FEDOROV M. V. and KORNYSHEV A. A., *Chem. Rev.*, **114** (2014) 2978.

[6] FEDOROV M. V. and KORNYSHEV A. A., *J. Phys. Chem. B*, **112** (2008) 11868.

[7] BORUKHOV I., ANDELMAN D. and ORLAND H., *Electrochem. Acta*, **46** (2000) 221.

[8] KORNYSHEV A. A., *J. Phys. Chem. B*, **111** (2007) 5545.

[9] BAZANT M. Z., STOREY B. D. and KORNYSHEV A. A., *Phys. Rev. Lett.*, **106** (2011) 046102.

[10] DÉMERY V., DEAN D. S., HAMMANT T. C., HORGAN R. R. and PODGORNIK R., *EPL*, **97** (2012) 28004.

[11] DÉMERY V., DEAN D. S., HAMMANT T. C., HORGAN R. R. and PODGORNIK R., *J. Chem. Phys.*, **137** (2012) 064901.

[12] LEE A. A., VELLA D., PERKIN S. and GORIELY A., *J. Chem. Phys.*, **141** (2014) 094904.
[13] STELL G., *Phys. Rev. A*, **45** (1992) 7628.

[14] ORKOULAS G. and PANAGIOTOPOULOS A. Z., *J. Chem. Phys.*, **101** (1994) 1452.

[15] VEGA C., BRESME F. and ABASCAL J. L. F., *Phys. Rev. E*, **54** (1996) 2746.

[16] STELL G., *New Results on Some Ionic-Fluid Problems*, in *New Approaches to Problems in Liquid State Theory*, edited by CACCAMO C., HANSEN J.-P. and STELL G., *NATO Sci. Ser.*, Vol. **529** (Springer, Netherlands) 1999, pp. 71-89.

[17] BRESME F., VEGA C. and ABASCAL J. L. F., *Phys. Rev. Lett.*, **85** (2000) 3217.

[18] LUIJTEN E., FISHER M. E. and PANAGIOTOPOULOS A. Z., *Phys. Rev. Lett.*, **88** (2002) 185701.

[19] VEGA C., ABASCAL J. L. F., MCBRIDE C. and BRESME F., *J. Chem. Phys.*, **119** (2003) 964.

[20] PANAGIOTOPOULOS A. Z., *J. Phys.: Condens. Matter*, **17** (2005) S3205.

[21] WALKER A. B. and GILLAN M. J., *J. Phys. C: Solid State Phys.*, **16** (1983) 3025.

[22] PANAGIOTOPOULOS A. Z. and KUMAR S. K., *Phys. Rev. Lett.*, **83** (1999) 2981.

[23] DICKMAN R. and STELL G., *AIP Conf. Proc.*, **492** (1999) 225.

[24] BROGNARA A., PAROLA A. and REATTO L., *Phys. Rev. E*, **65** (2002) 066113.

[25] KOBELEV V., KOLOMEISKY A. B. and FISHER M. E., *J. Chem. Phys.*, **116** (2002) 7589.

[26] CIACH A. and STELL G., *Phys. Rev. Lett.*, **91** (2003) 060601.

[27] DIEHL A. and PANAGIOTOPOULOS A. Z., *Phys. Rev. E*, **71** (2005) 046118.

[28] REN R., O'KEEFFE C. J. and ORKOULAS G., *J. Chem. Phys.*, **125** (2006) 124504.

[29] BARTSCH H., DANNENMANN O. and BIER M., *Phys. Rev. E*, **91** (2015) 042146.

[30] AIZENMAN M. and FRÖHLICH J., *J. Stat. Phys.*, **26** (1981) 347.

[31] LANDAU L. D. and LIFSHITZ E. M., *Statistical Physics, Course of Theoretical Physics*, Vol. **5** (Butterworth-Heinemann) 1980.

[32] THOULESS D. J., *Phys. Rev.*, **187** (1969) 732.

[33] CUESTA J. and SÁNCHEZ A., *J. Stat. Phys.*, **115** (2004) 869.