# Brownian dynamics simulations of electric double-layer capacitors with tunable metallicity

Paul Desmarchelier, Alexandre P dos Santos, Yan Levin, Benjamin Rotenberg

► To cite this version:

Paul Desmarchelier, Alexandre P dos Santos, Yan Levin, Benjamin Rotenberg. Brownian dynamics simulations of electric double-layer capacitors with tunable metallicity. The Journal of Chemical Physics, 2026, 164 (7), pp.074106. 10.1063/5.0314749. hal-05517000

HAL Id: hal-05517000

https://hal.sorbonne-universite.fr/hal-05517000v1

Submitted on 18 Feb 2026

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 CC BY 4.0 - Attribution - International License
# Brownian dynamics simulations of electric double-layer capacitors with tunable metallicity

Paul Desmarchelier,¹ Alexandre P. dos Santos,² Yan Levin,² and Benjamin Rotenberg¹, ³, a)

¹) Sorbonne Université, CNRS, Laboratoire PHENIX (Physicochimie des Electrolytes et Nanosystèmes Interfaciaux), 4 place Jussieu, 75005 Paris, France

²) Instituto de Física, Universidade Federal do Rio Grande do Sul, Caixa Postal 15051, CEP 91501-970 Porto Alegre, RS, Brazil

³) Réseau sur le Stockage Electrochimique de l'Energie (RS2E), FR CNRS 3459, 80039 Amiens Cedex, France

(Dated: 19 January 2026)

We introduce an efficient description of electrodes, characterized by their Thomas-Fermi screening length $l_{TF}$ inside the metal, for Brownian dynamics (BD) simulations of capacitors. Within a Born-Oppenheimer approximation for the electron charge density inside the electrodes, we derive the effective many-body potential for ions in an implicit solvent between Thomas-Fermi electrodes, taking into account the constraints of applied voltage and of global electro-neutrality of the system, as well as the 2D periodic boundary conditions along the electrode surfaces. We derive the average charge and the fluctuation-dissipation relation for the differential capacitance, highlighting the contribution of the fluctuations of the net ionic dipole moment, as well as those from the solvent polarization and of the electron density, whose fluctuations are suppressed within the Born-Oppenheimer description. We demonstrate the relevance of this model by validating its predictions against known results for the force on ions as a function of the ion-surface distance in simple geometries. The equilibrium ionic density profiles from BD simulations are in excellent agreement with those from an explicit electrode model for perfect metals, and are obtained at a significantly lower computational cost. Finally, we discuss with the present model the effect of the Thomas-Fermi screening length on the equilibrium ionic density profiles and the capacitance. While limited to parallel plate capacitors, the present simulation method allows to consider larger systems, lower concentrations, and longer time scales concentrations than molecular simulations in order to predict the electrochemical properties of Thomas-Fermi capacitors and correlate them with the ion dynamics.

## I. INTRODUCTION

The interface between metals and liquid electrolytes plays a key role in electrochemical energy storage devices (batteries, electric double layer capacitors), electrocatalysis, or electro-chemical sensing. In these contexts, the interfacial properties result from the mutual influence of the electron density inside the metal and of the charge distribution due to the solvent molecules and the ions inside the liquid. This interplay controls in particular the equilibrium properties such as the ionic densities and electrostatic potential profiles, or the capacitance that links the charge accumulated on both sides of the interface to the potential drop across it. From the theoretical point of view, the predictions of equilibrium properties generally rely on the mean-field Poisson-Boltzmann theory, its linearized limit (Debye-Hückel theory), or its extensions¹⁻⁵. More advanced liquid state theories such as integral equations, field theory or classical Density Functional Theory (cDFT) have been developed to include in particular electrostatic correlations⁶,⁷, the finite size of the ions (e.g. within the primitive model of electrolytes consisting of charged hard spheres)⁸,⁹, the explicit polarization of the solvent in addition to the ionic densities¹⁰⁻¹³ or even molecular features of the solvent around the ions¹⁴,¹⁵. For dynamical properties, most results have been obtained at the same level of description as Poisson-Boltzmann theory (point ions in an implicit solvent interacting via the mean-field electrostatic potential),

within Poisson-Nernst-Planck theory. This model allowed in particular to investigate the charging dynamics in capacitors¹⁶⁻¹⁹ or their frequency-dependent impedance²⁰,²¹. Here again, progress e.g. with time-dependent cDFT allowed to analyze the effect of some physical ingredients such as the finite size of the ions or electrokinetic couplings on the dynamical response of capacitors²²⁻²⁵.

In the last decades, these theoretical predictions were complemented by simulations at various levels of description. Ab initio or hybrid QM/MM molecular dynamics (MD) simulations provide the most accurate description but entail a large computational cost that limits the system size and trajectory length that can be simulated²⁶⁻³³. Classical all-atom and coarse-grained MD simulations offer a good trade-off between accuracy of the description and computational cost, and the implementation in generic or dedicated open source simulation packages³⁴⁻³⁸ of methods to sample systems in a statistical ensemble where the voltage between the electrodes is imposed³⁹,⁴⁰ have enabled their use for a variety of electrode/electrolyte interfaces. We refer the reader e.g. to Refs. 41 and 42 for recent reviews on these topics. Such simulations allow in particular to compute the capacitance or the frequency-dependent impedance of electrochemical cells and to correlate these observables with the microscopic properties of the confined liquid⁴³⁻⁴⁶.

In order to further decrease the computational cost of particle based simulations, implicit solvent models, where the effect of the solvent on ion-ion interactions is captured by the solvent permittivity and the ions evolve according to Langevin or Brownian dynamics (BD), instead of Newton's equation of

<footnote>a) benjamin.rotenberg@sorbonne-universite.fr</footnote>
<figure><img src="image_2.png" /></figure>

motion. The effects of the confining dielectric walls or perfect conductors and of voltage are then introduced as external potentials or enter in the ion-ion interactions, in order to model driven electrolytes confined in slit pores or capacitors under an applied voltage³⁵,⁴⁷⁻⁵³. For perfect conductors, efficient simulations at this level of description, taking into account periodic boundary conditions inherent to particle-based simulations in condensed matter systems, have been proposed within the Green's function formalism to compute effective ion-ion interactions⁵⁴⁻⁵⁷. An alternative approach for BD simulations, combining the description of explicit constant-potential electrodes using in MD simulations and ions with an implicit solvent described by its permittivity has also been proposed by Cats et al.⁵⁸. Dynamical properties, such as the electric double layer relaxation after a charge transfer event⁵⁹, the charge fluctuations of an electrode due to the interfacial ion dynamics and redox reaction⁶⁰,⁶¹, or the frequency- and field-dependent response of confined electrolytes⁶² have also been investigated using BD simulations.

Even though most theoretical and simulation studies to date have been limited to perfect metals, the importance of charge screening inside the metal has long been recognized. At the fundamental level, this phenomenon arises from the kinetic energy of the electrons, which has to be described using quantum mechanics. Thomas-Fermi (TF) theory<sup>63,64</sup> introduces this within a simple kinetic energy functional. Within a linearization approximation, the electrostatic potential then satisfies an equation identical to the Debye-Hückel equation for electrolytes, with the screening length replaced by the so-called Thomas-Fermi screening length, *l*<sub>TF</sub>, that depends on the density of states of the metal. This behavior can also be understood in terms of non-local electrostatics, via a wave vector dependent permittivity that depends on *l*<sub>TF</sub>. TF theory has been used to investigate the consequence of such screening inside the metal on ion-wall, ion-ion interactions or interfacial capacitance<sup>65-70</sup> (note that the electrostatic problem is the same as the one for an ion in a slit pore surrounded by an electrolyte treated at the Debye-Hückel level, that was solved analytically in Ref. 71). The consequences on the collective interfacial behavior has been demonstrated experimentally, for example on the capillary freezing of ionic liquids confined depending on the metallicity of the substrate<sup>72</sup>, which could be rationalized within TF theory<sup>73</sup>, or the kinetics of electron transfer at twisted graphene bilayers<sup>74</sup>. Several methods have also been proposed to introduce TF theory in atomistic simulations, to examine the consequences of screening inside the metal on the interfacial properties<sup>44,75-79</sup>. They suffer however from the above-mentioned limitations in terms of system size, salt concentration and trajectory length due to the computational cost of molecular simulations.

In the present work, we introduce an efficient description of electrodes, characterized by their Thomas-Fermi screening length *l*<sub>TF</sub> inside the metal, for Brownian dynamics (BD) simulations of capacitors. In Section II, we derive within a Born-Oppenheimer approximation for the electron charge density inside the electrodes, the effective many-body potential for ions in an implicit solvent between Thomas-Fermi electrodes, taking into account the constraints of applied voltage and of

global electro-neutrality of the system, as well as the 2D periodic boundary conditions along the electrode surfaces. We derive the average charge and the fluctuation-dissipation relation for the differential capacitance, highlighting the contribution of the fluctuation of the net ionic dipole moment, as well as those from the solvent polarization and of the electron density, whose fluctuations are suppressed within the Born-Oppenheimer description. In section III, we demonstrate the relevance of this model by validating its predictions against known results for the force on ions as a function of the ion-surface distance in simple geometries. The equilibrium ionic density profiles from BD simulations are in excellent agreement with those from an explicit electrode model for perfect metals, and are obtained at a significantly lower computational cost. Finally, we discuss with the present model the effect of the Thomas-Fermi screening length on the equilibrium ionic density profiles and the capacitance. While limited to parallel plate capacitors, the present simulation method allows to consider larger systems, lower concentrations, and longer time scales concentrations than molecular simulations in order to predict the electrochemical properties of Thomas-Fermi capacitors and correlate them with the ion dynamics.

## II. THEORY

### A. System

<figure><img src="image_3.png" /><figcaption>FIG. 1. Capacitor consisting of two Thomas-Fermi electrodes separated over a distance *L* by an electrolyte solution, under an applied voltage *ΔΨ*. The electrolyte consists of explicit ions in an implicit solvent with permittivity *ε*₀*ε*<sub>s</sub>, while the electrodes are characterized by a permittivity *ε*₀ and a Thomas-Fermi screening-length *l*<sub>TF</sub>. Periodic boundary conditions in the *x* and *y* directions along the electrode-electrolyte interfaces, with box dimensions *L*<sub>x</sub> and *L*<sub>y</sub>, are indicated by dashed lines. The dotted lines indicate the connexion with an electric circuit for *z* → ±∞ that imposes the external voltage.</figcaption></figure>

We consider a capacitor consisting of two electrodes labeled left, *l*, and right, *r*, separated by an electrolyte solution
(s), as illustrated in Fig. 1. These regions are characterized by their position $z < -L/2$, $z > L/2$ and $z \in [-L/2, L/2]$, respectively, while the system is periodic in the x and y directions. Note that the reference to left and right is only for convenience and is not related to a specific orientation of the capacitor. The electrodes are modeled by a continuous background of positive charges, with charge density $\rho_{l,r}^{\text{back}}(\mathbf{r}) = +en_{l,r}^0$, with $n_{l,r}^0$ the average electron density in each electrode, and continuous electron distributions, with charge density $\rho_{l,r}^e(\mathbf{r}) = -e[n_{l,r}^0 + \delta n_{l,r}(\mathbf{r})]$, while the electrolyte is modeled by point charges $q_i$ at positions $\mathbf{r}_i$ in a solvent with permittivity $\epsilon_0\epsilon_s$. The corresponding charge densities are

$$ \rho_{l,r}(\mathbf{r}) = \rho_{l,r}^{\text{back}}(\mathbf{r}) + \rho_{l,r}^{e}(\mathbf{r}) = -e\delta n_{l,r}(\mathbf{r}) \quad (1) $$

$$ \rho_{\text{ions}}(\mathbf{r}) = \sum_{i=1}^{N_{\text{ions}}} q_i \delta(\mathbf{r} - \mathbf{r}_i) \quad (2) $$

The total charge density can be written as $\rho^{\text{tot}}(\mathbf{r}) = \rho_l(\mathbf{r}) + \rho_{\text{ions}}(\mathbf{r}) + \rho_r(\mathbf{r})$, even though only one of the three terms contributes in each region. It is the source of an electrostatic potential $\phi(\mathbf{r})$ satisfying the Poisson equation $\Delta\phi = -\rho^{\text{tot}}/\epsilon_0\epsilon$ in each region of space ($\epsilon = 1$ inside the electrodes and $\epsilon_s$ in the electrolyte), with $\Delta$ the Laplacian operator, as well as continuity of $\phi$ and $\epsilon\phi'$ at the interfaces between regions. The values of the potentials for $z \to \pm\infty$ will be discussed in the following. Since both ions and electrons are part of the system, the total electrostatic energy is

$$ \begin{aligned} U_{\text{Coul}}[n_l, n_r, \rho_{\text{ions}}] &= \frac{1}{2} \int \rho(\mathbf{r})\phi(\mathbf{r}) \, d\mathbf{r} \\ &= \frac{1}{2} \int \rho_l \phi \, d\mathbf{r} + \frac{1}{2} \int \rho_{\text{ions}} \phi \, d\mathbf{r} + \frac{1}{2} \int \rho_r \phi \, d\mathbf{r}, \end{aligned} \quad (3) $$

where in the second line we have separated the integral over the three relevant regions of space. The electrostatic potential can be written as a $\phi(\mathbf{r}) = \int \rho^{\text{tot}}(\mathbf{r}')G(\mathbf{r}, \mathbf{r}')d\mathbf{r}'$, where the Green's function $G(\mathbf{r}, \mathbf{r}')$ describes the electrostatic interaction between unit point charges at $\mathbf{r}$ and $\mathbf{r}'$ and the integral runs over the whole space occupied by the electrodes and the electrolyte. This highlights the quadratic form of the electrostatic energy and the fact that a factor of $1/2$ is needed to avoid double-counting. Note that $G(\mathbf{r}, \mathbf{r}')$ diverges when $\mathbf{r}' \to \mathbf{r}$, for example as $1/4\pi\epsilon_0\epsilon_s|\mathbf{r}' - \mathbf{r}|$ in the electrolyte solution (see Section II B 2 and Appendix A for more details). This diverging self-interaction term should be removed when computing the energy.

In addition to the Coulomb energy, we consider the kinetic energy of the electrons described by

$$ T[n_l, n_r] = \int_l \tau(n_l(\mathbf{r})) d\mathbf{r} + \int_r \tau(n_r(\mathbf{r})) d\mathbf{r} \quad (4) $$

where

$$ \tau(n) = An^{5/3} = \frac{3h^2}{40m_e} \left(\frac{3}{\pi}\right)^{2/3} n^{5/3}, \quad (5) $$

with $h$ Planck's constant and $m_e$ the mass of the electron, the Thomas-Fermi (TF) kinetic energy functional. This term,

which is of quantum-mechanical origin, results together with the Coulomb energy in screened electrostatic interactions that can be described within the classical Thomas-Fermi model.

While the system is canonical for the electrolyte (fixed number of ions $N_{\text{ions}}$) and for the background positive charges in the electrodes, it is grand-canonical for the electrons: the electron densities in both electrodes fluctuate by exchanging electrons with reservoirs setting their electrochemical potentials $\mu_{l,r}$ (one value per electrode). In addition, we impose the global electroneutrality of the system:

$$ \int_l (-e) \delta n_l d\mathbf{r} + \int_s \rho_{\text{ions}} d\mathbf{r} + \int_r (-e) \delta n_r d\mathbf{r} = 0. \quad (6) $$

In practice, this amounts to imposing that electrons are exchanged between the electrodes, rather than with two independent reservoir, hence (see below) to imposing and a finite electrochemical potential difference $\Delta\mu = \mu_r - \mu_l$ or voltage $\Delta\Psi$ between the electrodes instead of independent electrochemical potentials. Assuming that the electrolyte is also overall neutral, this implies that the total charge of the two electrodes are opposite to each other. The corresponding thermodynamic potential for the electrons is:

$$ \begin{aligned} \Omega[n_l, n_r | \rho_{\text{ions}}, \mu_l, \mu_r] = & U_{\text{Coul}}[n_l, n_r, \rho_{\text{ions}}] + T[n_l, n_r] \\ & - \mu_l \int_l (n_l^0 + \delta n_l) d\mathbf{r} - \mu_r \int_r (n_r^0 + \delta n_r) d\mathbf{r} \\ & - \lambda \left[ \int_l (-e) \delta n_l d\mathbf{r} + \int_s \rho_{\text{ions}} d\mathbf{r} + \int_r (-e) \delta n_r d\mathbf{r} \right] \end{aligned} \quad (7) $$

where $\lambda$ is a Lagrange multiplier enforcing global electroneutrality, and the total Coulomb energy and the kinetic energy of the electrons are given by Eqs. 3 and 4, respectively. Note that we have made explicit here that we consider a functional of the electron densities $n_{l,r}(\mathbf{r})$ that depends parametrically on the ionic distribution $\rho_{\text{ions}}(\mathbf{r})$ and the electrochemical potentials $\mu_{l,r}$.

## B. Born-Oppenheimer description

For a given configuration of the ions and corresponding charge distribution, $\rho_{\text{ions}}(\mathbf{r})$, we look for the electron density that minimizes the total energy under the constraints of fixed electrochemical potentials $\mu_{l,r}$ and global electroneutrality. This is achieved by minimizing the thermodynamic potential Eq. 7 with respect to the charge densities $n_{l,r}$ in the electrodes, i.e. solving the Euler-Lagrange equations:

$$ \left. \frac{\delta \Omega[n_l, n_r | \rho_{\text{ions}}, \mu_l, \mu_r]}{\delta n_{l,r}(\mathbf{r})} \right|_{n_{l,r}^{\text{BO}}} = 0, \quad (8) $$

where the BO superscript refers to the Born-Oppenheimer distribution of the electrons satisfying these constraints. In the following, we will not always include the BO superscript explicitly but one should remember that the equations below hold only for $n_{l,r}^{\text{BO}}$. Once these distributions are determined,
one can define an effective (many-body) potential of the system consisting only of ions as:

$$V^{\text{eff}}[\{\mathbf{r}_i\}|\mu_l, \mu_r] = \Omega[n_1^{\text{BO}}, n_r^{\text{BO}}] |\rho_{\text{ions}}, \mu_l, \mu_r| \quad (9)$$

where we have again highlighted the parametric dependence on the electrochemical potentials $\mu_{l,r}$. The corresponding force acting on the ions can then be computed as the gradient with respect to the ionic positions

$$\mathbf{F}_i^{\text{eff}} = -\nabla_{\mathbf{r}_i} V^{\text{eff}}[\{\mathbf{r}_i\}|\mu_l, \mu_r] \quad (10)$$

which accounts for the reorganization of the charge distribution within the electrodes as the ions move. This is analogous to what is done in ab initio molecular dynamics, via the Hellmann-Feynman theorem, or to the case of explicit colloids with small ions treated at the DFT level, see Ref. 80.

## 1. Euler-Lagrange equations

In order to proceed further, we assume that $|\delta n_{l,r}| \ll n_{l,r}^0$ and approximate the kinetic energy Eq. 4, up to quadratic order in the electrode charge densities for consistency with the Coulomb term, using

$$\tau(n_{l,r}) = \tau(n_{l,r}^0) + \mu_{l,r}^T \delta n_{l,r} + \frac{1}{2} \alpha_{l,r} \delta n_{l,r}^2 + \mathcal{O}(\delta n_{l,r}^3) \quad (11)$$

where (see Eq. 5)

$$\mu_{l,r}^T = \left. \frac{\partial \tau}{\partial n} \right|_{n_{l,r}^0} = \frac{5A}{3} (n_{l,r}^0)^{2/3} \quad (12)$$

$$\alpha_{l,r} = \left. \frac{\partial^2 \tau}{\partial n^2} \right|_{n_{l,r}^0} = \frac{10A}{9} (n_{l,r}^0)^{-1/3} \quad (13)$$

The corresponding local chemical potentials are

$$\frac{\delta T[n_l, n_r]}{\delta n_{l,r}(\mathbf{r})} = \mu_{l,r}^T + \alpha_{l,r} \delta n_{l,r}(\mathbf{r}) + \mathcal{O}(\delta n_{l,r}^2) \quad (14)$$

Minimizing the thermodynamic potential with respect to the electronic densities in the left and right electrodes yields (using Eq. 11 and the fact that the Coulomb energy is quadratic in the densities):

$$\frac{\delta \Omega[n_l, n_r | \rho_{\text{ions}}, \mu_l, \mu_r]}{\delta n_{l,r}(\mathbf{r})} = 0 \\ = \mu_{l,r}^T + \alpha_{l,r} \delta n_{l,r}(\mathbf{r}) - e\phi(\mathbf{r}) - \mu_{l,r} + \lambda e. \quad (15)$$

In particular, in the bulk of the electrodes $\delta n_{l,r}(\mathbf{r}) = 0$ so that

$$\mu_{l,r} = \mu_{l,r}^T - e(\Psi_{l,r} - \lambda) \quad (16)$$

where $\Psi_l = \lim_{z \to -\infty} \phi$ and $\Psi_r = \lim_{z \to +\infty} \phi$ are such that $\Delta\Psi = \Psi_r - \Psi_l$ is the voltage imposed between the two electrodes, and $\lambda$ amounts to a constant potential shift of both potentials, which does not change the voltage, as expected. As anticipated in Section II A, Eq. 16 also highlights the fact

that in practice, one imposes the difference $\Delta\mu = \mu_r - \mu_l$ rather than the two values independently. The Euler-Lagrange equations can be rewritten as

$$\alpha_{l,r} \delta n_{l,r}(\mathbf{r}) = e[\phi(\mathbf{r}) - \Psi_{l,r}] \quad (17)$$

Since inside the electrode, the charge densities $\delta n_{l,r}$ are also related to the Poisson equation $\Delta\phi = -(e)\delta n/\epsilon_0$, the potential satisfies the Thomas-Fermi equation in each electrode:

$$\Delta\phi = k_{TF}^2 [\phi(\mathbf{r}) - \Psi_{l,r}] \quad (18)$$

where we have introduced the inverse of the Thomas-Fermi screening length (possibly different in both electrodes, even though here we consider identical electrodes)

$$k_{TF}^2 = \frac{1}{l_{TF}^2} = \frac{e^2}{\epsilon_0 \alpha_{l,r}} \quad (19)$$

where $\alpha_{l,r}$ is given in Eq. 13. Compared to the familiar form of the TF equation, usually written for a single electrode, Eq. 18 differs in the constant potential shifts $\Psi_{l,r}$ in each electrode that enforce the applied voltage between electrodes. As shown below (see Section II B 4), the constant potential shift $\lambda$ does not enter in the effective many-body potential between the ions, and the electrode potentials only appear as the voltage $\Delta\Psi$.

## 2. Electrostatic potential

The electrostatic potential in the electrodes satisfies the above Thomas-Fermi equations Eq. 18, together with the boundary conditions $\Psi_l = \lim_{z \to -\infty} \phi$ and $\Psi_r = \lim_{z \to +\infty} \phi$ such that $\Delta\Psi = \Psi_r - \Psi_l$. In the electrolyte, it satisfies the Poisson equation $\Delta\phi = -\rho_{\text{ions}}^{\text{tot}}/\epsilon_0\epsilon_s$ with $\epsilon_s$ the solvent permittivity. At the electrode/electrolyte interfaces, $\phi$ is continuous and its derivative is such that $\epsilon\phi'$ is continuous.

The full electrostatic problem is conveniently expressed as the sum of two contributions $\phi = \phi_{\Delta\Psi} + \phi_{\text{ions}}$, where $\phi_{\Delta\Psi}$ is the solution in the presence of the voltage $\Delta\Psi$ but in the absence of ions, while $\phi_{\text{ions}}$ is the solution in the presence of ions and in the absence of voltage. The former only depends on the z coordinate:

$$\phi_{\Delta\Psi}(z) = \begin{cases} \Psi_l + \Delta\Psi \frac{\epsilon_s l_{TF}}{L_{\text{eff}}} e^{\frac{1}{2}k_{TF}(L+2z)} & -\infty < z < -\frac{L}{2} \\ \frac{\Psi_l + \Psi_r}{2} + \Delta\Psi \frac{z}{L_{\text{eff}}} & -\frac{L}{2} < z < \frac{L}{2} \\ \Psi_r - \Delta\Psi \frac{\epsilon_s l_{TF}}{L_{\text{eff}}} e^{\frac{1}{2}k_{TF}(L-2z)} & \frac{L}{2} < z < +\infty \end{cases} \quad (20)$$

where we have introduced the effective length (further discussed in Section II B 3):

$$L_{\text{eff}} = L + 2\epsilon_s l_{\text{TF}}. \quad (21)$$

The potential due to the ions can further be separated into their individual contributions as $\phi_{\text{ions}} = \sum_i^{N_{\text{ions}}} \phi_i$, where each
term is the potential satisfying the problem where the ion distribution is limited to a single ion with charge $q_i$ at position $\mathbf{r}_i$ and its periodic images, in the absence of voltage:

$$ \rho_i(\mathbf{r}) = \sum_{m_x=-\infty}^{\infty} \sum_{m_y=-\infty}^{\infty} q_i \delta(\mathbf{r}-\mathbf{r}_i + m_x L_x \mathbf{e}x + m_y L_y \mathbf{e}_y). \quad (22) $$

Following Ref. 54 for the case of a perfect metal ($l_{TF} = 0$), we solve the problem in reciprocal space in the x and y directions by looking for a Green's function of the form

$$ \phi_i(\mathbf{r}) = \frac{q_i}{L_x L_y} \sum_{\mathbf{k}} g_{\mathbf{k}}(z_i, z) \cos[\mathbf{k} \cdot (\mathbf{r} - \mathbf{r}_i)], \quad (23) $$

with $\mathbf{k} = (\frac{2\pi m_x}{L_x}, \frac{2\pi m_y}{L_y})$. This leads to solutions of the form $g_k = Ae^{-kz} + Be^{kz}$ with $k = |\mathbf{k}|$ between the electrodes and $g_k = C_{\pm}e^{\pm\chi_{TF}z}$ where

$$ \chi_{TF} = \sqrt{k^2 + k_{TF}^2} \quad (24) $$

inside the electrodes, the constants being determined by the boundary conditions. In addition to the one mentioned above-mentioned for $z \to \pm\infty$ and at $z = \pm L/2$, one also needs to account for continuity of $g_k(z_i, z)$ at $z = z_i$ and the discontinuity of its derivative:

$$ \epsilon_s \frac{\partial g_k(z_i, z)}{\partial z}\bigg|_{z=z_i^+} = \epsilon_s \frac{\partial g_k(z_i, z)}{\partial z}\bigg|_{z=z_i^-} - \frac{q_i}{\epsilon_0} \quad (25) $$

The complete solutions are provided in Appendix A. It can be checked in particular that $\phi_i(\mathbf{r})$ diverges as $1/4\pi\epsilon_0\epsilon_s|\mathbf{r} - \mathbf{r}_i|$ for $\mathbf{r} \to \mathbf{r}_i$, as expected. As already mentioned in Section II A, this diverging self-interaction term should be removed when computing the energy.

### 3. Electrode charge

The electrostatic potential also allows to express the local charge induced inside each electrode via the TF Eq. 18 and the Poisson equation. In turn, this provides the total charge of the electrodes $Q_{l,r} = (-e) \int_{l,r} \delta n_{l,r}(\mathbf{r}) d\mathbf{r}$ induced by the applied voltage and the presence of ions. When the electrolyte is neutral, the electrodes are oppositely charged and we can introduce the charge of the right electrode $Q = Q_r = -Q_l$. Using Eqs. 20 and A1 the electrode charge can be expressed as:

$$ \begin{aligned} Q[\{\mathbf{r}_i\}|\Delta\Psi] &= L_x L_y \frac{\epsilon_0 \epsilon_s k_{TF} \Delta\Psi}{(2\epsilon_s + k_{TF} L)} - \frac{k_{TF}}{(2\epsilon_s + k_{TF} L)} \sum_{i}^{N_{\text{ions}}} q_i z_i \\ &= C_0 \Delta\Psi - \frac{M_{\text{ions}}}{L_{\text{eff}}} \end{aligned} \quad (26) $$

where we have introduced

$$ M_{\text{ions}} = \sum_{i} q_i z_i \quad (27) $$

the dipole (along the z direction) of the ion distribution and the capacitance of the ion-free capacitor such that

$$ \frac{L_x L_y}{C_0} = \frac{L}{\epsilon_0 \epsilon_s} + \frac{2l_{TF}}{\epsilon_0} = \frac{L_{\text{eff}}}{\epsilon_0 \epsilon_s} \quad (28) $$

and we also used the electroneutrality of the electrolyte to introduce $\phi_{\Delta\Psi}$ in the second equality. Note that Eq. 28 can be interpreted as an equivalent circuit, with two capacitors each with capacitance per unit area $\epsilon_0/l_{TF}$ accounting for the effect of the TF metals in series with one with capacitance per unit area $\epsilon_0\epsilon_s/L$ accounting for the effect of the dielectric slab (see also Ref. 81). For $l_{TF} = 0$, Eq. 26 reduces to the result derived for an ideal metal in Ref. 55.

The effective length $L_{\text{eff}}$ defined in Eq. 21 can thus be understood as the width of an equivalent dielectric slab, and $E_{\text{eff}} = -\Delta\Psi/L_{\text{eff}}$ is the corresponding uniform electric field in the slab induced by the applied voltage. This effective length is also consistent with the one introduced to describe the lateral decay of the charge induced by an ion in a single TF electrode, proposed recently in Ref. 77 as an extension of earlier work by Vorotyntsev and Kornyshev$^{66}$, even though here the distance is the one between the two electrodes (hence the factor of two for their contribution) and not an ion-surface distance. Finally, Eq. 26 can be rearranged as

$$ M_{\text{elec}} + M_{\text{ions}} + M_{\text{solv}} = 0, \quad (29) $$

where we have introduced the dipole of the electrode charge distribution $M_{\text{elec}} = QL_{\text{eff}}$ and the solvent dipole $M_{\text{solv}} = \epsilon_0\epsilon_s L_x L_y LE$, with $E = -\Delta\Psi/L$. This underlines the various contributions to the screening of the electric field throughout the system, as discussed in the case of $l_{TF} = 0$ for explicit electrodes and solvent (without ions)$^{45}$.

### 4. Effective many-body potential between ions

As explained in Eq. 9, we evaluate the thermodynamic potential $\Omega[n_l^{\text{BO}}, n_r^{\text{BO}}, \rho_{\text{ions}}]$ for the charge density inside the electrodes satisfying the above Euler-Lagrange equations. As a result, it only depends on the ion distribution and, parametrically, on the electron chemical potentials. Using Eqs. 7 and 11, as well as Eq. 16, Eq. 17 and the electroneutrality condition, we obtain:

$$ \begin{aligned} V^{\text{eff}}[\{\mathbf{r}_i\}|\mu_l, \mu_r] &\equiv \Omega[n_l^{\text{BO}}, n_r^{\text{BO}}|\rho_{\text{ions}}, \mu_l, \mu_r] \\ &= \int_{l,r} [\tau(n_{l,r}^0) - \mu_{l,r} n_{l,r}^0] d\mathbf{r} + \frac{1}{2} \int_s \rho_{\text{ions}} \phi d\mathbf{r} \\ &\quad + \int_{l,r} \left[ \mu_{l,r}^T - \mu_{l,r} + \frac{1}{2} (\alpha_{l,r} \delta n_{l,r} - e\phi) \right] \delta n_{l,r} d\mathbf{r} \\ &= \Omega_0 + \frac{1}{2} \int_s \rho_{\text{ions}} \phi d\mathbf{r} + \int_{l,r} \frac{e}{2} (\Psi_{l,r} - \lambda) \delta n_{l,r} d\mathbf{r} \\ &= \Omega_0 + \frac{1}{2} \int_s \rho_{\text{ions}} \phi d\mathbf{r} - \frac{1}{2} [(\Psi_{r} - \lambda)Q_{r} + (\Psi_{l} - \lambda)Q_{l}]. \end{aligned} \quad (30) $$

where $\Omega_0 = \int_{l,r} [\tau(n_{l,r}^0) - \mu_{l,r} n_{l,r}^0] d\mathbf{r}$ is the electronic grand-potential in the absence of ions and voltage, which does not
depend on the ionic configuration and therefore doesn't contribute to the force on the ions. We remind the reader that Eq. 30 corresponds to an expansion of the kinetic energy to second order in $\delta n_{l,r}$, whereas the Coulomb contribution is exactly quadratic in these quantity. The last term can be rewritten, when the electrolyte is neutral ($\int_s \rho_{ions} d\mathbf{r} = \sum_i q_i = 0$), by introducing the charge of the right electrode $Q = Q_r = -Q_l$, whose expression is provided in Eq. 26, the voltage $\Delta\Psi$, and $\phi = \phi_{\Delta\Psi} + \phi_{ions}$, as:

$$ 
\begin{aligned}
V^{\text{eff}}[\{\mathbf{r}_i\}|\Delta\Psi] &= \Omega_0 + \frac{1}{2} \int_s \rho_{\text{ions}}(\mathbf{r})\phi(\mathbf{r}) d\mathbf{r} - \frac{1}{2} Q \Delta\Psi \\
&= \Omega_0 - \frac{1}{2} C_0 \Delta\Psi^2 + \frac{1}{2} \sum_i q_i \phi_{\text{ions}}(\mathbf{r}_i) + \sum_i q_i \phi_{\Delta\Psi}(\mathbf{r}_i) \\
&= \Omega_0 - \frac{1}{2} C_0 \Delta\Psi^2 + U_{\text{Coul}}^{\text{ions}} + \frac{\Delta\Psi}{L_{\text{eff}}} M_{\text{ions}},
\end{aligned}
\quad (31) $$

with

$$ U_{\text{Coul}}^{\text{ions}}[\{\mathbf{r}_i\}] = \frac{1}{2} \sum_i q_i \phi_{\text{ions}}(\mathbf{r}_i) \quad (32) $$

where, as mentioned previously, the diverging self-energy $q_i^2/4\pi\epsilon_0\epsilon_s|\mathbf{r} - \mathbf{r}_i|$ is implicitly removed when evaluating $\phi_{ions}(\mathbf{r}_i)$.

Note that a constant shift of all electrostatic potentials does not change the r.h.s. of Eq. 31 since it would cancel in $\Delta\Psi$ and would only add a contribution proportional to $\sum_i q_i = 0$ due to the electroneutrality of the electrolyte. In fact, this remark also holds more generally in Eq. 30 due to the global electroneutrality ($Q_l + Q_r + \sum_i q_i = 0$). In addition, the first two terms in Eq. 31 do not depend on the ionic positions and therefore do not contribute to the forces acting on them (see Eq. 10). The effect of ion-ion interactions in the presence of the implicit solvent and Thomas-Fermi electrodes will be discussed in more detail in Section II C 2. Finally, the force acting on ion $i$ due to the applied voltage (negative gradient with respect to $\mathbf{r}_i$ of the last term) is simply $q_i E_{eff} e_z$ with the effective field $E_{eff} = -\Delta\Psi/L_{eff}$ introduced in Section II B 3.

### 5. Partition function, average charge and capacitance

Brownian dynamics simulations of ions using the effective many-body potential Eq. 31 sample the statistical ensemble corresponding to fixed number of ions $N_{ions}$, volume $V$, temperature $T$ and voltage $\Delta\Psi$. The corresponding partition function is

$$ \mathcal{L}[\Delta\Psi] = \int d\mathbf{r}^{N_{ions}} e^{-\beta V^{\text{eff}}[\{\mathbf{r}_i\}|\Delta\Psi]} \quad (33) $$

where $\beta = 1/k_B T$ with $k_B$ Boltzmann's constant, and the corresponding thermodynamic potential is the free energy $\mathcal{F}[\Delta\Psi] = -k_B T \ln \mathcal{L}[\Delta\Psi]$ (see Ref. 44). The average charge, for a given voltage, follows straightforwardly from Eq. 26:

$$ 
\begin{aligned}
\langle Q \rangle &= \frac{1}{\mathcal{L}} \int d\mathbf{r}^{N_{\text{ions}}} Q[\{\mathbf{r}_i\} | \Delta\Psi] e^{-\beta V^{\text{eff}}[\{\mathbf{r}_i\} | \Delta\Psi]} \\
&= C_0 \Delta\Psi - \frac{\langle M_{\text{ions}} \rangle}{L_{\text{eff}}},
\end{aligned}
\quad (34) $$

where $\langle M_{ions} \rangle$ is the ensemble-averaged dipole of the ionic distribution. Furthermore, the charge for fixed configuration of the ions and voltage, $Q[\{\mathbf{r}_i\}|\Delta\Psi]$ is equal to $-\partial V^{\text{eff}}[\{\mathbf{r}_i\}|\Delta\Psi]/\partial\Delta\Psi$ (see Eqs. 31 and 26), so that:

$$ 
\begin{aligned}
\langle Q \rangle &= \frac{1}{\mathcal{L}} \left( -\frac{1}{\beta} \frac{\partial}{\partial\Delta\Psi} \int d\mathbf{r}^{N_{\text{ions}}} e^{-\beta V^{\text{eff}}[\{\mathbf{r}_i\}|\Delta\Psi]} \right) \\
&= \frac{1}{\beta} \frac{\partial \ln \mathcal{L}}{\partial\Delta\Psi} = -\frac{\partial \mathcal{F}}{\partial\Delta\Psi},
\end{aligned}
\quad (35) $$

similarly to what was found for simulations with explicit electrodes in the constant-potential ensemble$^{44}$. This shows that the free energy change during the charge of the capacitor from zero voltage to finite voltage $\Delta\Psi_{max}$ is

$$ \Delta \mathcal{F} = - \int_{0}^{\Delta \Psi_{\text{max}}} \langle Q \rangle_{\Delta \Psi} d\Delta \Psi, \quad (36) $$

i.e. the reversible electrical work exchanged with the charge reservoir. The subscript highlights here the fact that the ensemble average of the charge is made at fixed voltage $\Delta\Psi$.

The differential capacitance is the derivative of the average electrode charge with respect to voltage:

$$ 
\begin{aligned}
C_{\text{diff}} = \frac{\partial \langle Q \rangle}{\partial \Delta \Psi} &= \beta [\langle Q^2 \rangle - \langle Q \rangle^2] \\
&= C_0 + \frac{\beta}{L_{\text{eff}}^2} [\langle M_{\text{ions}}^2 \rangle - \langle M_{\text{ions}} \rangle^2].
\end{aligned}
\quad (37) $$

For a purely capacitive system, this quantity is independent of voltage and equal to the integral capacitance $C_{int} = \langle Q \rangle / \Delta\Psi$. However this is not necessarily the case in general. Eq. 37 is a fluctuation-dissipation relation relating the response of the average electrode charge to a change in voltage to the equilibrium fluctuations of the electrode charge at fixed voltage. These fluctuations arise from two contributions: the thermal fluctuations of the ions, as obvious from the second term in the last line of Eq. 37, and the thermal fluctuations of the implicit electrons and the solvent, which are suppressed in the Born-Oppenheimer description and embedded in $C_0$, via the Thomas-Fermi screening length $l_{TF}$ and the relative permittivity $\epsilon_s$, respectively (see Eq. 28). This result is analogous to the one obtained in Ref. 44 with explicit electrode atoms and electrolyte (including the solvent). In the latter case, the explicit expression for the contribution of the electrolyte-free capacitor has the same physical origin as $C_0$ but is not as simple as Eq. 28, as it depends on the microscopic details of the electrodes.

### C. Brownian Dynamics simulations

The ions are described as Brownian particles with a diffusion coefficient $D_i$. Their positions $\mathbf{r}_i$ evolve according to the overdamped Langevin equation:

$$ \dot{\mathbf{r}}_i = \beta D_i \mathbf{F}_i + \sqrt{2D_i} \xi_i \quad (38) $$

where $\xi_i$ is a Gaussian white noise, and the force $\mathbf{F}_i$ acting on the ions includes the interactions between them and with
the walls. They include short-range forces, described in Section II C 1, as well as the one deriving from the effective many-body potential (see Eq. 10), described further in Section II C 2.

## 1. Short-range interactions

Short-range repulsion interactions between ions i and j are described by Weeks-Chandler-Andersen (WCA) potentials:

$$v_{ij}^{\text{WCA}}(r) = \begin{cases} v_{ij}^{\text{LJ}}(r) - v_{ij}^{\text{LJ}}(r^*) & , r \le r^*, \\ 0 & , r > r^*, \end{cases} \quad (39)$$

with the Lennard-Jones (LJ) potential

$$v_{ij}^{\text{LJ}}(r) = 4\epsilon_{ij} \left[ \left( \frac{\sigma_{ij}}{r} \right)^{12} - \left( \frac{\sigma_{ij}}{r} \right)^{6} \right], \quad (40)$$

and $r^* = 2^{1/6}\sigma_{ij}$ the position of the minimum of $v_{ij}^{\text{LJ}}$. The LJ energy and diameter $\epsilon_{ij}$ and $\sigma_{ij}$ are computed from the corresponding parameters for ions i and j using the Lorentz-Berthelot mixing rules.

Two descriptions are considered for the short-range interactions between ions and the electrodes. For the comparison with explicit electrodes (see Section II C 3), we use the same short-range interactions, i.e. pairwise WCA interactions with explicit atoms, so that only the many-body effective potential describing Coulomb interactions and the effect of the kinetic energy of the electrons differ from the explicit case. For the other cases, the short-range interactions with the wall are also treated implicitly via a potential of the form

$$U(z) = V_w(z+L/2) + V_w(z-L/2), \quad (41)$$

which depends only on z, where $V_w$ is the so-called Steele potential, obtained by integrating LJ interactions over atomic planes. The resulting potential includes short-range repulsion and an attractive well leading to adsorption$^{82–85}$.

$$V_w^{\text{Steele}}(z) = 2\pi \rho_{\text{surf}} \epsilon_w \sigma_w \Delta \left[ \frac{2}{5} \left( \frac{\sigma_w}{z} \right)^{10} - \left( \frac{\sigma_w}{z} \right)^4 - \frac{\sigma_w^4}{3\Delta(z+0.61\Delta)^3} \right], \quad (42)$$

with $\rho_{\text{surf}}$ the surface number density of electrode atoms, $\Delta$ the distance between atomic planes within the electrode, and where $\epsilon_w$ and $\sigma_w$ tune the strength and range of the ion-wall interaction. In order to model purely repulsive walls, we use

the same truncation and shifting as for the short-range ion-ion interactions, i.e.

$$V_w(z) = \begin{cases} V_w^{\text{Steele}}(z) - V_w^{\text{Steele}}(z^*) & , z \le z^*, \\ 0 & , z > z^*, \end{cases} \quad (43)$$

with $z^* \approx 0.986\sigma_w$. More details on interactions are provided in Section II C 3.

## 2. Effective many-body potential

The force arising from electrostatic interactions and the kinetic energy of the electrons is computed as the negative gradient with respect to $\mathbf{r}_i$ Eq. 10 of the effective many-body potential. As already mentioned in Section II B 4, the first two terms in Eq. 31 do not contribute to the force acting on ion i. The last term provides the force due to the applied voltage and is simply $q_i E_{\text{eff}} \mathbf{e}_z$ with the effective field $E_{\text{eff}} = -\Delta\Psi/L_{\text{eff}}$ introduced in Section II B 3. The force due to the effect of ion-ion interactions in the presence of the implicit solvent and Thomas-Fermi electrodes (third term in Eq. 31) is the negative gradient of $U_{\text{Coul}}^{\text{ions}}[\{\mathbf{r}_i\}]$ defined in Eq. 32 where the diverging self-energy $q_i^2/4\pi\epsilon_0\epsilon_s|\mathbf{r}-\mathbf{r}_i|$ is implicitly removed when evaluating $\phi_{\text{ions}}(\mathbf{r}_i)$.

The electrostatic energy $U_{\text{Coul}}^{\text{ions}}[\{\mathbf{r}_i\}]$ can be split into two contributions: one due to direct Coulomb interactions between the ions without electrodes (corresponding to the limit of infinite interelectrode distance $L \to \infty$), which can be efficiently calculated under 2D periodic boundary conditions using Ewald summation techniques$^{86,87}$, and another due to the polarization of the electron charge distribution inside the metal. Specifically, we write

$$U_{\text{Coul}}^{\text{ions}} = U_{\text{Coul}}^{\infty} + \delta U_{\text{Coul}}^{\text{ions}} \quad (44)$$

with (see also Appendix A):

$$\begin{aligned} U_{\text{Coul}}^{\infty}[\{\mathbf{r}_i\}] &\equiv \lim_{L \to \infty} \frac{1}{2} \sum_i q_i \phi_{\text{ions}}(\mathbf{r}_i) \\ &= \frac{1}{2} \sum_i \sum_{j \neq i} \sum_{\mathbf{n}} \frac{q_i q_j}{4\pi \epsilon_0 \epsilon_s |\mathbf{r}_j - \mathbf{r}_i + \mathbf{n}|}, \end{aligned} \quad (45)$$

where $\mathbf{n} = m_x L_x \mathbf{e}_x + m_y L_y \mathbf{e}_y$ accounts for the sum over periodic images and the prime now explicitly indicates the removal of the term $j=i$ when $\mathbf{n}=\mathbf{0}$. Using Eq. A1 for the potential $\phi_{\text{ions}}(\mathbf{r})$ with Eq. A3 for the Green's function in the electrolyte, as well as the electroneutrality of the electrolyte ($\Sigma_i q_i = 0$), we obtain the following expression for the difference:

$$\delta U_{\text{Coul}}^{\text{ions}} = U_{\text{Coul}}^{\text{ions}} - U_{\text{Coul}}^{\infty} = \frac{1}{2L_x L_y} \left\{ \sum_i \sum_j \sum_{\mathbf{k} \neq \mathbf{0}} q_i q_j \delta g_{\mathbf{k}}^s (z_i, z_j) \cos[\mathbf{k} \cdot (\mathbf{r}_j - \mathbf{r}_i)] + \frac{M_{\text{ions}}^2}{\epsilon_0 \epsilon_s L_{\text{eff}}} \right\} \quad (46)$$
with (see Eqs. A3 and A5)

$$ \delta g_k^s(z_i, z_j) = \frac{\left[ e^{-k(z_j+z_i+L)} + e^{k(z_j+z_i-L)} \right] (\epsilon_s^2 k^2 - \chi_{TF}^2) + \left[ e^{k(z_i-z_j-2L)} + e^{k(z_j-z_i-2L)} \right] (\epsilon_s k - \chi_{TF})^2}{2\epsilon_0\epsilon_s k \left[ (1-e^{-2kL})(\epsilon_s^2 k^2 + \chi_{TF}^2) + 2(1+e^{-2kL})\epsilon_s k \chi_{TF} \right]} \quad (47) $$

where $\chi_{TF}$ is defined in Eq. 24 and the last term in Eq. 46 corresponds to $\mathbf{k} = 0$. Since $z_i$ and $z_j$ are in the interval $]-L/2, L/2]$, the sum over modes $\mathbf{k}$ converges exponentially fast, regardless of the positions of the two ions. Note that there is no diverging self-energy term to be subtracted in Eq. 46.

### 3. Simulation details

The coupled equations Eq. 38 are solved numerically using the MetalWalls simulation package$^{36,37}$, in which we implemented the Euler-Maruyama integrator, a rescaling of the 2D Ewald summation by the permittivity $\epsilon_s$ to compute $U_{Coul}^{ions}$ (see Eq. 32) and the new term $\delta U_{Coul}^{ions}$ (see Eq. 46) to account for the effect of the implicit Thomas-Fermi electrodes, as well as the corresponding forces.

The parameters for the short-range interactions (see Section II C 1) are identical for cations, anions and wall atoms, namely $\sigma_{ij} = \sigma_w = 5$ Å and $\epsilon_{ij} = \epsilon_w = 2.477$ kJ/mol. The cut-off for ion-ion interactions is thus $r^* = 5.61$ Å, while that for the Steele potential is $z^* = 4.92$ Å. For the definition of the latter, the graphite structure corresponds to $\rho_{surf} = 0.38$ Å$^{-2}$ and $\Delta = 3.354$ Å. The repulsive walls (explicit atom sites or zero of the Steele potential) are placed at the same position as the dielectric interface, i.e. $z = \pm L/2$.

For the Coulomb interactions, we use a real-space cut-off of $r_c = 15.9$ Å in the Ewald summation and choose the number of wave vectors to compute $U_{Coul}^{\infty}$ and $\delta U_{Coul}^{ions}$ to achieve a prescribed tolerance $\delta U_{tol} = 3 \times 10^{-5}$ eV. For $\delta U_{Coul}^{ions}$, which is computed only in reciprocal space, this leads to $k_{max} > \log\left[\frac{q_{max}^2}{8\pi d \epsilon_0 \delta U_{tol}}\right]/2d$, with $d$ the distance of closest approach of ions from the interface. For isolated ions (Section III 1) and ion pairs (Section III 2) we can choose this according to the range of considered ion-distance surfaces. For the Brownian dynamics simulations of Sections III 3 and III 4, we choose $d = \sigma/2$, with $\sigma$ the common diameter describing short-range interactions between ions and with the wall. Except for the results in vacuum for validation on isolated ions and ion pairs, the simulated systems correspond to a solvent with the relative permittivity of water, $\epsilon_s = 78$.

For BD simulations, we consider ions with a diffusion coefficient $D_{\pm} = 1.12 \times 10^{-9}$ m$^2$s$^{-1}$ and a temperature $T = 298$ K. The overdamped Langevin equations 38 are integrated with a time step of 5 fs. The systems are first equilibrated for 100 ns, before production runs of 400 ns, with electrode charges (computed from Eq. 26) sampled every 0.25 ps and ion positions every 0.5 ps. The reported uncertainty on ionic concentration profiles is the standard error over 10 consecutive blocks of the trajectory. The differential capacitance is computed from the variance of the charge distribution (first

line in Eq. 37). The corresponding uncertainty is obtained as in Ref. 44 following Zwanzig and Ailawadi$^{88}$: The standard error estimated as $\langle \delta Q^2 \rangle \times \sqrt{4\tau_c/\tau_s}$ where $\tau_c$ is the sampling time and $\tau_c$ is the correlation time computed as $\tau_c = \int_0^\infty \langle \delta Q(t) \delta Q(0) \rangle^2 / \langle \delta Q^2 \rangle^2 dt$. The uncertainty over the integral capacity is computed via the same method with $\langle Q \rangle \times \sqrt{4\tau_c^2/\tau_s}$ with $\tau_c = \int_0^\infty \langle \delta Q(t) \delta Q(0) \rangle / \langle \delta Q^2 \rangle dt$. We refer the reader to Ref. 46 for more details on the computation of these correlation times.

### III. RESULTS

We now demonstrate the relevance of the proposed implicit electrode description, by first validating its predictions for the force on the ions as a function of the ion-surface distance in simple geometries: we compare our results with (semi-)analytical for isolated ions near electrodes for various $l_{TF}$ in Section III 1) and with numerical results under periodic boundary conditions for ion pairs next to perfect metals ($l_{TF} = 0$) in Section III 2. We then turn in Section III 3 to the equilibrium ionic density profiles and compare our results with Brownian Dynamics simulations with an explicit electrode model for perfect metals proposed in Ref. 58, before using our model that also applies to finite $l_{TF}$ to discuss the effect of the Thomas-Fermi screening length on the ionic density profiles in Section III 4. Section III 5 finally discusses the capacitance as a function of $l_{TF}$ and compares, when possible, several levels of description.

#### 1. Isolated ions

As a first test of the present method, we compare our numerical predictions for a periodic system with two electrodes and several ions to the semi-analytical results from the literature for a single charge in a dielectric next to a Thomas-Fermi metal. We use the result for this non-periodic case with a single solid/liquid interface in a superposition approximation to compare with our case with two electrodes and several ions, without considering the effect of periodic boundary conditions. Specifically, we compare $U_{Coul}^{ions}[\{\mathbf{r}_i\}]$ for a simple geometry with the reference:

$$ U_{Coul}^{\text{ref}}[\{\mathbf{r}_i\}] = \frac{1}{2} \sum_i q_i \phi_i^{\text{ref}}(\mathbf{r}_i) = \sum_i U_i^{\text{ref}}(\mathbf{r}_i) \quad (48) $$
with

$$U_i^{\text{ref}}(\mathbf{r}_i) = \frac{q_i^2}{8\pi\varepsilon_0\varepsilon_s} \int_0^\infty dk \frac{\varepsilon_s k - \sqrt{k^2 + k_{TF}^2}}{\varepsilon_s k + \sqrt{k^2 + k_{TF}^2}} e^{-2kd_i} \quad (49)$$

where $d_i$ is the distance of ion i from the nearest electrode. Eq. 49 was first derived by Kornyshev et al. (see Eqs. 6 and following of Ref. 65) and used as reference in several subsequent works (see in e.g. Eq. S19 of Ref. 75 or more recently Eq. 69 of Ref. 70). This superposition approximation assumes in particular that $d_i \ll L$ for all ions and that the lateral dimensions are sufficiently large $d_i \ll L_x, L_y$. Eq. 49 reduces to the well-known limits:

$$\lim_{l_{TF} \to 0} U_i^{\text{ref}}(\mathbf{r}_i) = - \frac{q_i^2}{16\pi\varepsilon_0\varepsilon_s d_i} \quad (50)$$

and

$$\lim_{l_{TF} \to \infty} U_i^{\text{ref}}(\mathbf{r}_i) = \left( \frac{\varepsilon_s - 1}{\varepsilon_s + 1} \right) \frac{q_i^2}{16\pi\varepsilon_0\varepsilon_s d_i} \quad (51)$$

for perfect conductors and insulating walls, respectively.

Figure 2 compares the energy $U_{\text{Coul}}^{\text{ions}}$ with $U_{\text{Coul}}^{\text{ref}}$ for a system with an interelectrode distance $L = 159$ Å and lateral dimensions $L_x = L_y = L$, and two ions with opposite charges $q = \pm e$ placed symmetrically at varying distances $d$ from each electrode, i.e. positions $x_i = y_i = 0$ and $z_i = \pm(\frac{L}{2} - d)$, in vacuum (panel 2a) or in a solvent with relative permittivity $\varepsilon_s = 78$ (panel 2b). In both cases, the numerical results for all considered TF screening length $l_{TF}$ ranging from 0 (perfect metal) to $5a_0$, with $a_0 \approx 0.53$ Å the Bohr radius, are in excellent agreement with the semi-analytical prediction. This simultaneously shows that the box sizes and considered $d$ range are consistent with the superposition approximation and that the proposed scheme provides the expected results. We observe in particular that, while in vacuum the force is always attractive, decays with the ion-surface distance $d$ and with increasing $l_{TF}$, the dielectric contrast with the solvent results in a cross-over from purely attractive for a perfect metal ($l_{TF} = 0$) to repulsive for large $l_{TF}$, as expected for insulating walls ($l_{TF} \to \infty$). For intermediate $l_{TF}$, the force profile may be non-monotonic, with both repulsive and attractive forces depending on $d$, as observed here for $l_{TF} = 0.1a_0$.

## 2. Ion pair next to a perfect metal

We then compare the present method with the existing one for perfect metals ($l_{TF} = 0$) described in Ref. 57 for a simple system where an ion pair is approached from the same electrode. Specifically, we consider a pair with charges $+e$ at $(0, 0, -\frac{L}{2} + d)$ and $-e$ at $(x, 0, -\frac{L}{2} + d)$ with $x = 0.714$ Å, in vacuum. The box dimensions are $L_x = 67.69$ Å, $L_y = 36.64$ Å, $L = 39.72$ Å. Fig. 3 reports the components of the electrostatic force on the positive charge in the directions perpendicular and parallel to the electrodes as a function of the distance $d$ from them (the force on the negative charge is simply equal

<figure><img src="image_4.png" /><figcaption>FIG. 2. Electrostatic force in the direction perpendicular to the implicit electrodes, as a function of the distance $d$ from them, for a system consisting of two “isolated” ions with charges $+e$ at $(0, 0, -\frac{L}{2} + d)$ and $-e$ at $(0, 0, +\frac{L}{2} - d)$ in vacuum (a) or in a solvent with relative permittivity $\varepsilon_s = 78$ (b), in the absence of voltage. The reported force is that on the positive charge; the force on the negative charge is simply the opposite. The numerical results with the present method (lines) for various Thomas-Fermi screening lengths $l_{TF}$ (in atomic units, with $a_0 \approx 0.53$ Å the Bohr radius) indicated by colors are compared with the force corresponding to Eqs. 48 and 49, for independent ions next to an infinite interface between a dielectric medium and a Thomas-Fermi metal (symbols). The dotted line in panel b indicates the insulating limit $l_{TF} \to \infty$ (see Eq. 51).</figcaption></figure>

for the $z$ component and opposite for the $x$ one). The results for both components as a function of distance from the interface are in excellent agreement with the ones obtained by the method of Ref. 57. This proves that the present approach, which is implemented as a correction to the case of ions in the absence of walls under 2D periodic boundary conditions (see Eq. 44) correctly describes the effect of the interface in the case of a perfect metal ($l_{TF} = 0$).

## 3. Comparison with explicit electrodes and implicit solvent

After the above validation for specific configurations of the ions by comparison with previous analytical and numerical results, we now turn to the study of equilibrium properties of capacitors obtained by Brownian dynamics simulations. We compare the results obtained with the present model describing both the polar solvent and the electrodes implicitly via the Green's functions to that obtained in the case of perfect metals ($l_{TF} = 0$) by Cats et al. using explicit electrodes 58. In this work, the authors computed the forces on ions to be used in Brownian dynamics (in their case underdamped Langevin dynamics) by adapting the method introduced for molecular dynamics simulations with explicit electrode atoms maintained
<figure><img src="image_5.png" /><figcaption>FIG. 3. Electrostatic force in the directions perpendicular (a) and parallel (b) to the implicit electrodes as a function of the distance $d$ from them, for a pair of ions with charges $+e$ at $(0,0,-\frac{L}{2}+d)$ and $-e$ at $(x,0,-\frac{L}{2}+d)$ with $x = 0.714$ Å, in vacuum, for perfect metals ($l_{TF} = 0$). The reported force is that on the positive charge; the force on the negative charge is simply equal for the $z$ component and opposite for the $x$ one. The numerical results with the present method (lines) are compared with that obtained with the method of Ref. 57 (symbols), which only applies to perfect metals.</figcaption></figure>

at a constant potential (with a potential difference between the two electrodes) and an explicit electrolyte (ions and solvent molecules). They suggested that for perfect metals, an implicit solvent can be modeled by simultaneously rescaling the ionic charges by a factor $1/\sqrt{\epsilon_s}$ and the voltage by a factor $\sqrt{\epsilon_s}$. While this rescaling in fact amounts to considering that the electrode atoms are also embedded in a medium with permittivity $\epsilon_s$, this approach resulted in an excellent agreement for the considered systems (monovalent ions in a very polar solvent) and voltages with classical Density Functional Theory calculations.

Here, we consider exactly the same system as in Ref. 58, namely a capacitor consisting of two graphite electrodes separated by $L = 39.72$ Å. Each electrode is modeled by a single atomic layer of graphite with unit cell parameter $a = 2.46$ Å and total dimensions $L_x \times L_y = 67.69 \times 36.64$ Å$^2$, with a total of 960 atoms. We compare the results of overdamped Langevin dynamics simulations (Eq. 38) with forces computed as in Ref. 58 or with the present method, which does not rely on explicit electrode atoms to compute the electrostatic forces on the ions. For the short-range interactions, we consider either explicit atomic sites (interacting with the ions via pairwise WCA potentials, see Eq. 39) to isolate the effect of the electrostatic interactions, or implicit walls (see Eqs. 41-43), to isolate the effect of the short-range ones.

The results obtained with the three methods, shown in Fig. 4, are in excellent agreement with each other. As ex-

<figure><img src="image_6.png" /><figcaption>FIG. 4. Equilibrium anion (red solid lines) and cation (blue dashed lines) density profiles from Brownian dynamics simulations of ions in an implicit solvent between perfect metals ($l_{TF} = 0$), with (a) the explicit electrode model from Ref. 58 (see text), (b) the present model of implicit electrodes and short-range interactions computed with the same atomic sites as in the previous panel and (c) the present fully implicit model where short-range interactions are computed via the Steele potential (see Eqs. 41-43). In all panels, the density profiles are shown for $\Delta\Psi = 0, 0.1$ and $0.2$ V from light to dark lines, respectively.</figcaption></figure>

pected for an electrolyte consisting of ions with opposite charges and identical short range interactions, the profiles are always symmetric with respect to $z = 0$. They are identical in the absence of voltage, and display a small maximum near the electrodes, which results not only from the weak attractive interactions of each ion with the metal (see the orange line in Fig. 2b for $l_{TF} = 0$) but also from the pressure inside the relatively dense liquid due to short-range repulsion (see also Section III 4). As voltage is applied, cations (resp. anions) accumulate at the electrode with the lower (resp. higher) potential and are depleted from the other one, as expected. The accumulation/depletion of ions near the corresponding electrodes increases with applied voltage.
The agreement between the various methods show that the same prediction for the ionic density profiles can be obtained with a significantly simpler model and, accordingly, a reduced computational cost. Indeed, for the considered system with a single layer of electrode atoms in each electrode (as in Ref. 58), with a total of 1920 electrode atoms in the system, the computational cost of electrostatic interactions is reduced by a factor of ~ 6 (all comparisons reported here are obtained from simulations with a single CPU on the same computer). For a system with three atomic planes in each electrodes, the speed-up is by a factor of ~ 60. For the same systems, the speed-up achieved using an implicit wall for short-range interactions instead of explicit sites is of approximately 25 and 80, respectively. For these systems, the computation of electrostatic interactions represents ~ 95% of the cost per step in the fully implicit model. Finally, while we compared the present implicit electrode model to one with explicit electrodes, in both cases the solvent is treated implicitly.

The speed-up with respect to a fully atomistic simulation is even more dramatic. As an example, for a system considered in Ref. 46 with a 0.5 M NaCl aqueous solution between model gold electrodes separated by a distance L ≈ 5 nm (with 40 ions, 2160 water molecules and 3240 electrode atoms), the present BD simulations (with only 40 ions and implicit solvent and electrodes) is ≈ 6 times faster on a single CPU than the MD simulations on a GPU. The present model only approximately describes the effect of the solvent, neglecting in particular subtle effects that may occur in the first molecular layers at the interface, and a detailed comparison remains out of the scope of the present work (see some discussion in Section IV). Nevertheless, it offers a computationally efficient approach to investigate the effect of the TF screening length on the properties of capacitors with dimensions much larger than the ones currently within reach of MD simulations. This includes the possibility to consider lower salt concentrations (typically above 0.1 mol/L in MD simulations), by extending the lateral dimensions to ensure a sufficient number of ions in the system, or larger inter-electrode distances (typically limited to 10-20 nm in MD simulations). Importantly, the reduced computational cost per step and the larger time steps used in Brownian dynamics offer the possibility to also investigate the dynamics over much longer times scales (typically limited to a few tens of ns in MD simulations).

### 4. Effect of the Thomas-Fermi screening length

We now use the present implicit model to investigate the effect of the Thomas-Fermi screening length on the ionic distributions and the capacitance. We recall that this was not possible with previous methods based on Green's functions, which only applied so far to the perfect metal case⁵⁴,⁵⁵,⁵⁷. In the following, we consider the same system as in Section III 3, except for the variable $l_{TF}$.

Fig. 5 shows the effect of Thomas-Fermi screening lengths $l_{TF}$ on the ion density profile ($c_{ions} = c_{+} + c_{-}$), in the absence of voltage ($\Delta\Psi = 0$). Specifically, it reports the difference between the profile for a given $l_{TF}$ with respect to the profile for

perfect metals ($l_{TF} = 0$). This reference profile can be seen in Fig. 4, where the cation and anion density profiles overlap for $\Delta\Psi = 0$ V. The negative values of the difference next to the wall results reflect the cross-over from attractive to repulsive ion-surface interactions, illustrated for isolated ions in panel Fig. 2b, upon increasing $l_{TF}$. Note that the decrease in the local ion concentration at the surface remains small (-0.3 mol/L at most for the largest considered $l_{TF}$, comparable to an insulating wall) compared to the maximum of the concentration ($c_{max} \sim 2.9 \pm 0.1$ mol/L, as can be deduced from the ionic density profiles of Fig. 4): the balance of interactions, including the short-range repulsion within the electrolyte, still results in density maxima near the walls in that case (see also Fig. 6 below). Far from the surface, there is no significant effect of $l_{TF}$ on the salt concentration, as expected.

Fig. 6a then shows the anion density profiles for a finite voltage $\Delta\Psi = 0.1$ V, for various $l_{TF}$. The cation density profiles, not shown for clarity, are simply symmetric with respect to $z = 0$. In all cases, this voltage corresponding to 4 times the thermal potential $k_B T/e \sim 25$ mV induces an asymmetry in the anion density profile, with a decrease (resp. increase) with respect to $\Delta\Psi = 0$ V next to the left (resp. right) electrode at the lower (resp. larger) potential. This asymmetry is less pronounced as $l_{TF}$ increases. The effect of voltage is further emphasized in panel 6b, which reports the difference with respect to the cation and anion density profiles for $\Delta\Psi = 0$ V, across the left half of the capacitor ($z < 0$). The difference highlights the fact that the decrease in cation concentration and increase in anion concentration induced by voltage are symmetric when $l_{TF}$ is large: in that case the effective electric field inside the capacitor, corresponding to the effective length $L_{eff}$ (see Eq. 21) is weaker due to screening within the electrodes; this results in small changes in the ionic density profiles, symmetric and decaying approximately exponentially consistently with Debye-Hückel (linearized Poisson-

<figure><img src="image_7.png" /><figcaption>FIG. 5. Difference between the ion concentration profiles ($c_{ions} = c_{+} + c_{-}$) for various Thomas-Fermi screening lengths $l_{TF}$ (in units of the Bohr radius $a_0$) and the concentration profiles for $l_{TF} = 0$, in the absence of voltage ($\Delta\Psi = 0$). The system is identical to that of Fig. 4, except for the change in $l_{TF}$.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>FIG. 6. (a) Anion density profile under an applied voltage $\Delta\Psi = 0.1$ V for various Thomas-Fermi screening lengths $l_{TF}$. (b) Difference between the anion (solid lines) and cation (dashed lines) concentration profiles at $\Delta\Psi = 0.1$ V with respect to the profiles at $\Delta\Psi = 0$ V, for various $l_{TF}$. The system is identical to that of Fig. 4, except for the change in $l_{TF}$.</figcaption></figure>

Boltzmann) theory, even though one should not expect it to be accurate at such a high concentration. For smaller $l_{TF}$ values, the lack of screening within the metal results in large and asymmetric depletion/enrichment of the corresponding ions.

## 5. Capacitance

Finally, Fig. 7 reports the integral and differential capacitance of capacitors measured at $\Delta\Psi = 0.1$ V, for capacitors with electrodes characterized by various Thomas-Fermi screening $l_{TF}$, using the implicit electrode model for electrostatic interactions and either the explicit or implicit description of short-range interactions. For $l_{TF} = 0$, results for the explicit electrode model are also reported. In all cases, the integral and differential capacitances are in good agreement, as expected for a purely capacitive system. A similarly good agreement is observed between the capacitances predicted

<figure><img src="image_9.png" /><figcaption>FIG. 7. Integral ($C_{int}$) and differential ($C_{diff}$) capacitance of capacitors measured at $\Delta\Psi = 0.1$ V, for capacitors with electrodes characterized by various Thomas-Fermi screening $l_{TF}$ (in units of the Bohr radius $a_0$), using the implicit electrode model for electrostatic interactions and either explicit or implicit (Steele) description of short-range interactions. For $l_{TF} = 0$, results for the explicit electrode model are also reported. In each case, the contribution $C_0$ of the ion-free capacitor to the differential capacitance (see Eqs. 28 and 37) is indicated as hatched areas.</figcaption></figure>

with the various descriptions of the same system. This highlights in particular the relevance of the simplest (fully implicit electrodes) introduced in the present work, which provides the same results at a significantly reduced computational cost.

The capacitance decreases upon increasing $l_{TF}$, consistently with the cross-over from perfect to insulating walls. Fig. 7 also indicates the contribution $C_0$ of the ion-free capacitor to the differential capacitance (see Eq. 37). Both $C_0$, corresponding to the fluctuations of the solvent polarization in the absence of ions (see Eq. 28) that are taken into account implicitly via the permittivity $\epsilon_s$, and the remaining contribution corresponding to the ionic contribution to the fluctuations of the dipole of the liquid slab, decrease with increasing $l_{TF}$. However, the ionic contribution, which dominates for small $l_{TF}$, becomes negligible for larger values.

## IV. CONCLUSION AND PERSPECTIVES

We have introduced an efficient description of electrodes, characterized by the Thomas-Fermi screening length $l_{TF}$ inside the metal, for Brownian dynamics simulations of capacitors. Within a Born-Oppenheimer approximation for the electron charge density inside the electrodes, we derived the effective many-body potential for ions in an implicit solvent between Thomas-Fermi electrodes, taking into account the constraints of applied voltage between them and of global electro-neutrality of the system, as well as the 2D periodic boundary conditions along the electrode surfaces. Following previous work in the case of perfect metals ($l_{TF} = 0$), the problem is solved using Green's functions. The final result is con-
veniently expressed as a correction to the Ewald summation used to compute ion-ion interactions in standard simulation packages and fully computed in reciprocal space. By considering the statistical ensemble corresponding to fixed number of ions, volume, temperature and voltage, we derived using the effective many-body potential the expressions of the average charge and the fluctuation-dissipation relation for the differential capacitance, highlighting the contributions of ion distribution inside the capacitor and is fluctuations, as well as those from the solvent polarization and of the electron density, whose fluctuations are suppressed within the Born-Oppenheimer description.

We implemented this implicit description of the electrostatic interactions between ions in the presence of the solvent and the electrodes in the Metal Walls simulation package, considering both explicit and implicit descriptions of the short-range interactions between ions and the walls. We demonstrated the relevance of this model by first validating its predictions for the force on the ions as a function of the ion-surface distance in simple geometries, by comparing with results from the literature in limit cases (isolated ions or perfect metals). We then compared the equilibrium ionic density profiles from Brownian Dynamics simulations with those obtained using an explicit electrode model for perfect metals, finding excellent agreement and demonstrating the benefit of using the present implicit model in terms of computational cost. Finally, we used the present model to discuss the effect of the Thomas-Fermi screening length on the equilibrium ionic density profiles and the capacitance.

The present model is limited to parallel plate capacitors, although the various steps of Section II could in principle be extended to undulating surfaces, following Ref. 53 where Green's functions were derived for undulating dielectric membranes. But it can be used to investigate the effect of the TF screening length on the properties of capacitors with dimensions much larger than the ones currently within reach of molecular dynamics simulations. This would offer the possibility to consider lower salt concentrations (typically above 0.1 mol/L in MD simulations), by extending the lateral dimensions to ensure a sufficient number of ions in the system, or larger inter-electrode distances (typically limited to 10-20 nm in MD simulations). Importantly, the reduced computational cost per step and the larger time steps used in Brownian dynamics offer the possibility to also investigate the dynamics over much longer times scales (typically limited to a few tens of ns in MD simulations). It would in particular be interesting to derive the fluctuation-dissipation relation for the frequency-dependent admittance (inverse of the impedance), already successfully applied in MD simulations to link the electrochemical response with the dynamics within the liquid slab<sup>45,46</sup>, in the case of Brownian dynamics, for example following the work of Ref. 62 for the frequency-dependent conductivity of confined electrolytes. Beyond the study of these effects, the present BD approach can provide reference data to assess the predictions of simpler models such as PNP theory and its extensions, e.g. the ones mentioned in the introduction.

In the present work, we only considered systems without redox reactions at the electrode surfaces. In this purely capac-

itive case, one can assume that the electrolyte is overall neutral (hence that the two electrode charges are opposite) to obtain simple expressions for the electrode charge (Eq. 26) and the many-body effective potential (Eq. 31). However many steps of the derivation in Section II do not rely on this assumption. In particular, Eq. 30 remains valid and one could obtain an expression when the electrolyte is not neutral and its charge is compensated by the net charge of the electrodes. The result could then be used in hybrid simulations where the Brownian dynamics of ions is combined with a Monte Carlo scheme for tentative instantaneous charge transfers between ions and electrodes, as done in Ref. 89 for perfect metals. Such an extension would allow to investigate the effect of the Thomas-Fermi length e.g. the relaxation of the electric double layer after an electron transfer<sup>59</sup> or on the correlation between electrodes in electrochemical cells<sup>90</sup>.

Even though the TF model was successfully used to predict the effect of the screening length on a variety of interfacial properties, the accuracy of its predictions for a specific material can be limited, especially for properties that crucially depend on short-range features (expectedly so, since it is the low-wavevector and low-frequency limit of Lindhard theory of screening). Subtle effects, such as the electron spillover at the surface of the metal, which can enhance the capacitance (sometimes described as negative contributions to the inverse capacitance, e.g. in jellium models<sup>91,92</sup>), could be approximately captured at the present level by a careful definition of an effective sharp interface between the jellium-like metal and its neighborhood (vacuum or liquid), see e.g. Refs. 44,77. In order to go beyond, one could follow the steps described in Section II with a better kinetic energy functional. Adding the von Weizsäcker correction, that involves the gradient of the electronic density<sup>93,94</sup>, to the TF term Eq. 5 would be a natural step in this direction.

Future possible improvements also pertain to the description of the solvent. The effective interaction between an ion embedded in a dielectric and a Thomas-Fermi metal (see Fig. 2b) only partly reflects the actual potential of mean force (PMF), which could in principle be computed from all atoms MD simulations with constant-potential TF electrodes<sup>81</sup>. While a comparison of the present model with an explicit solvent is beyond the scope of the present work, we might anticipate that features such as the molecular layering of the solvent can lead to differences with the predictions of the present implicit model. Some of them can be captured in an effective way by tuning the position of the interface between the metal and the implicit solvent, using the concept of Dielectric Dividing Surface introduced by Netz an co-workers<sup>95</sup> (see also Ref. 96), as illustrated for the frequency-dependent permittivity of nanocapacitors<sup>45,46</sup>.

However others, in particular related to the solvation of interfacial ions, might require more elaborate descriptions such as molecular Density Functional Theory<sup>14,15,78</sup>. Such molecular descriptions would however not allow to simply express the many-body effective potential between ions. This could however be achieved by introducing the solvent polarization as an additional field in the free energy functional and treating this polarization at the Born-Oppenheimer level, as the elec-
tron charge density. In fact, the present model can be recovered in such a perspective using the simplest functional corresponding to a linear dielectric medium characterized by its permittivity. More elaborate functionals for the solvent (see e.g. Refs. 70,97,98) could then be introduced naturally at this level of description. In practice, a simpler approach would be to introduce molecular features through a one-body ion-wall PMF, derived or inspired from MD simulations⁹⁹. Such terms (e.g. with damped oscillations to capture the solvent layering), have already been used in extensions of Poisson-Boltzmann theory¹⁰⁰ or in simulations of electrokinetic phenomena¹⁰¹ and could straightforwardly be included in the present BD simulations via additional contributions to Eq. 41.

## ACKNOWLEDGMENTS

This article is dedicated to Christoph Dellago, whose work on developing simulation methods based on advanced Statistical Mechanics has been an inspiration over the years, on the occasion of his 60th birthday. The authors acknowledge discussions with Lydéric Bocquet, Benoit Coasne, David Limmer, Leonardo Coello, Mathieu Salanne and Andrea Grisafi. This project received funding from the European Research Council under the European Union's Horizon 2020 research and innovation program (grant agreement no. 863473). APdS acknowledges financial support from CNPq under grant 303310/2025-1.

## AUTHOR DECLARATIONS

### Conflict of interest

There is no conflict of interest to declare.

### Author contributions

**Paul Desmarchelier:** Conceptualization (equal); Formal analysis (equal); Investigation (lead); Methodology (supporting); Software (lead); Writing/Original Draft Preparation (supporting); Writing – review & editing (equal). **Alexandre Dos Santos:** Conceptualization (equal); Formal analysis (equal); Investigation (supporting); Methodology (equal); Writing – review & editing (equal). **Yan Levin:** Conceptualization (equal); Formal analysis (equal); Investigation (supporting); Methodology (equal); Writing – review & editing (equal). **Benjamin Rotenberg:** Conceptualization (equal); Formal analysis (equal); Funding Acquisition (lead); Investigation (supporting); Methodology (equal); Supervision (lead); Writing/Original Draft Preparation (lead); Writing – review & editing (equal).

## DATA AVAILABILITY

The original data presented in this study are openly available in Zenodo at [DOI/URL] inserted after final acceptance.

## Appendix A: Electrostatic potential

The electrostatic potential at position **r** due to the ion distribution $\sum_i \rho_i(\mathbf{r})$ where $\rho_i(\mathbf{r})$ given in Eq. 22 corresponds to ion **i** with charge $q_i$ at position $\mathbf{r}_i$ and its periodic images is:

$$ \phi_{\text{ions}}(\mathbf{r}) = \left\{ \begin{array}{ll} \frac{1}{L_x L_y} \sum_{i}^{N_{\text{ions}}} q_i \left\{ \sum_{\mathbf{k} \neq 0} g_{\mathbf{k}}^{l}(z_i, z) \cos[\mathbf{k} \cdot (\mathbf{r} - \mathbf{r}_i)] + \frac{e^{\frac{1}{2}k_{\text{TF}}(L+2z)}(2\varepsilon_s + k_{\text{TF}}(L-2z_i))}{2\varepsilon_0 k_{\text{TF}}(2\varepsilon_s + k_{\text{TF}}L)} \right\} & -\infty < z < -\frac{L}{2} \\ \\ \frac{1}{L_x L_y} \sum_{i}^{N_{\text{ions}}} q_i \left\{ \sum_{\mathbf{k} \neq 0} g_{\mathbf{k}}^{s}(z_i, z) \cos[\mathbf{k} \cdot (\mathbf{r} - \mathbf{r}_i)] + \frac{[2\varepsilon_s + k_{\text{TF}}(L+2z)][2\varepsilon_s + k_{\text{TF}}(L-2z_i)]}{4\varepsilon_0 k_{\text{TF}}(2\varepsilon_s + k_{\text{TF}}L)} \right\} & -\frac{L}{2} < z < z_i \\ \\ \frac{1}{L_x L_y} \sum_{i}^{N_{\text{ions}}} q_i \left\{ \sum_{\mathbf{k} \neq 0} g_{\mathbf{k}}^{s}(z_i, z) \cos[\mathbf{k} \cdot (\mathbf{r} - \mathbf{r}_i)] + \frac{[2\varepsilon_s + k_{\text{TF}}(L-2z)][2\varepsilon_s + k_{\text{TF}}(L+2z_i)]}{4\varepsilon_0 k_{\text{TF}}(2\varepsilon_s + k_{\text{TF}}L)} \right\} & z_i < z < \frac{L}{2} \\ \\ \frac{1}{L_x L_y} \sum_{i}^{N_{\text{ions}}} q_i \left\{ \sum_{\mathbf{k} \neq 0} g_{\mathbf{k}}^{r}(z_i, z) \cos[\mathbf{k} \cdot (\mathbf{r} - \mathbf{r}_i)] + \frac{e^{\frac{1}{2}k_{\text{TF}}(L-2z)}(2\varepsilon_s + k_{\text{TF}}(L+2z_i))}{2\varepsilon_0 k_{\text{TF}}(2\varepsilon_s + k_{\text{TF}}L)} \right\} & \frac{L}{2} < z < +\infty \end{array} \right. \quad (\text{A1}) $$

where the sums run over non-zero wave vectors corresponding to the periodic boundary conditions in the x and y directions along the surface, $\mathbf{k} = (\frac{2\pi m_x}{L_x}, \frac{2\pi m_y}{L_y})$, and for the left electrode,

inter-electrode and right electrode regions:
$$g_k^l(z_i, z) = \frac{e^{\frac{1}{2}(\chi_{TF}(L+2z) - k(L+2z_i))} (e^{k(2z_i - L)} (\epsilon_s k - \chi_{TF}) + (\epsilon_s k + \chi_{TF}))}{\epsilon_0 [(1 - e^{-2kL})(\epsilon_s^2 k^2 + \chi_{TF}^2) + 2(1 + e^{-2kL})\epsilon_s k \chi_{TF}]} \quad (A2)$$

$$g_k^s(z_i, z) = \frac{\left(e^{-k(z+z_i+L)} + e^{k(z+z_i-L)}\right) (\epsilon_s^2 k^2 - \chi_{TF}^2) + e^{k(|z-z_i|-2L)} (\epsilon_s k - \chi_{TF})^2 + e^{-k|z_j-z_i|} (\epsilon_s k + \chi_{TF})^2}{2\epsilon_0 \epsilon_s k [(1 - e^{-2kL})(\epsilon_s^2 k^2 + \chi_{TF}^2) + 2(1 + e^{-2kL})\epsilon_s k \chi_{TF}]} \quad (A3)$$

$$g_k^s(z_i, z) = \frac{e^{\frac{1}{2}(\chi_{TF}(L-2z) - k(L-2z_i))} (e^{-k(2z_i + L)} (\epsilon_s k - \chi_{TF}) + (\epsilon_s k + \chi_{TF}))}{\epsilon_0 [(1 - e^{-2kL})(\epsilon_s^2 k^2 + \chi_{TF}^2) + 2(1 + e^{-2kL})\epsilon_s k \chi_{TF}]} \quad (A4)$$

In practice, in order to compute the effective many-body potential one only requires the potential experienced by the ions, i.e. $g_k^s(z_i, z)$ given in Eq. A3. The terms outside the sums over modes in Eq. A1 correspond to $k = 0$ and can be obtained as the $k \to 0$ limits of the coefficients in Eq. A1, so that one could in fact write the sum over all modes including $k = 0$ for compactness. However, it is convenient to write this term separately to consider the limit $L \to \infty$, which is used to compute the energy $U_{\text{Coul}}^{\text{ions}}$ (see Eq. 44). In this limit, for $k \neq 0$ one has

$$g_k^\infty(z_i, z) \equiv \lim_{L \to \infty} g_k^s(z_i, z) = \frac{e^{-k|z-z_i|}}{2\epsilon_0 \epsilon_s k}, \quad (A5)$$

which corresponds to the 2D Fourier expansion of the Green's function for bulk ions in an implicit solvent, $1/4\epsilon_0\epsilon_s|\mathbf{r}-\mathbf{r}_i|$. In the same limit, the $k=0$ term in $\phi_{\text{ions}}(\mathbf{r})$ diverges but the corresponding $\mathcal{O}(L^2)$ and $\mathcal{O}(L)$ terms in the energy $\sum_i q_i \phi_{\text{ions}}(\mathbf{r}_i)$ vanish for any finite $L$ due to the electroneutrality of the electrolyte, so that the limit $L \to \infty$ of the energy is well defined (see Eq. 46).

1. R. Parsons, Chemical Reviews 90, 813 (1990).

2. A. A. Kornyshev, The Journal of Physical Chemistry B 111, 5545 (2007).

3. Z. A. H. Goodwin, G. Feng, and A. A. Kornyshev, Electrochimica Acta 225, 190 (2017).

4. A. Baskin and D. Prendergast, Journal of The Electrochemical Society 164, E3438 (2017).

5. D. A. Mazur, P. E. Brandyshev, S. V. Doronin, and Y. A. Budkov, ChemPhysChem 25, e202400650 (2024).

6. R. Kjellander and S. Marcélja, Chemical Physics Letters 112, 49 (1984).

7. T. Zhou, D. Bruch, and Z.-G. Wang, Phys. Rev. E 110, 044129 (2024).

8. P. Cats, R. Evans, A. Härtel, and R. van Roij, The Journal of Chemical Physics 154, 124504 (2021).

9. P. Cats and A. Härtel, The Journal of Chemical Physics 159, 184707 (2023).

10. A. Levy, D. Andelman, and H. Orland, Physical Review Letters 108, 227801 (2012).

11. A. Levy, D. Andelman, and H. Orland, The Journal of Chemical Physics 139, 164909 (2013).

12. N. Bruch, T. Binninger, J. Huang, and M. Eikerling, The Journal of Physical Chemistry Letters 15, 2015 (2024).

13. N. Bruch, T. Binninger, J. Huang, and M. Eikerling, The Journal of Chemical Physics 162, 224703 (2025).

14. G. Jeanmairet, M. Levesque, R. Vuilleumier, and D. Borgis, The Journal of Physical Chemistry Letters 4, 619 (2013).

15. G. Jeanmairet, B. Rotenberg, D. Borgis, and M. Salanne, The Journal of Chemical Physics 151, 124111 (2019).

16. M. Z. Bazant, K. Thornton, and A. Ajdari, Physical Review E 70, 021506 (2004).

17. M. Janssen and M. Bier, Physical Review E 97, 052616 (2018).

18. I. Palaia, A. J. Asta, M. Dutta, P. B. Warren, B. Rotenberg, and E. Trizac, Physical Review Letters 135, 148002 (2025).

19. I. Palaia, A. J. Asta, M. Dutta, P. B. Warren, B. Rotenberg, and E. Trizac, Physical Review E 112, 035417 (2025).

20. G. Barbero, F. Batalioto, and A. M. Figueiredo Neto, Applied Physics Letters 92, 172908 (2008).

21. A. Antonova, G. Barbero, L. R. Evangelista, and P. Tilli, Journal of Statistical Mechanics: Theory and Experiment 2020, 043202 (2020).

22. M. S. Kilic, M. Z. Bazant, and A. Ajdari, Phys. Rev. E 75, 021503 (2007).

23. J. Jiang, D. Cao, D. E. Jiang, and J. Wu, Journal of Physics Condensed Matter 26, 284102 (2014).

24. A. J. Asta, I. Palaia, E. Trizac, M. Levesque, and B. Rotenberg, The Journal of Chemical Physics 151, 114104 (2019).

25. K. Ma, M. Janssen, C. Lian, and R. van Roij, The Journal of Chemical Physics 156, 084101 (2022).

26. N. Bonnet, T. Morishita, O. Sugino, and M. Otani, Phys. Rev. Lett. 109, 266101 (2012).

27. J. D. Elliott, A. Troisi, and P. Carbone, Journal of Chemical Theory and Computation 16, 5253 (2020).

28. J.-B. Le, A. Chen, L. Li, J.-F. Xiong, J. Lan, Y.-P. Liu, M. Iannuzzi, and J. Cheng, JACS Au 1, 569 (2021).

29. K. Takahashi, H. Nakano, and H. Sato, The Journal of Chemical Physics 157, 234107 (2022).

30. A. Groß and S. Sakong, Chemical Reviews, acs.chemrev.1c00679 (2022).

31. S. Sakong, J. Huang, M. Eikerling, and A. Groß, Current Opinion in Electrochemistry 33, 100953 (2022).

32. A. Groß, Current Opinion in Electrochemistry 40, 101345 (2023).

33. A. Grisafi, A. Bussy, M. Salanne, and R. Vuilleumier, Phys. Rev. Mater. 7, 125403 (2023).

34. F. Weik, R. Weeber, K. Szuttor, K. Breitsprecher, J. de Graaf, M. Kuron, J. Landsgesell, H. Menke, D. Sean, and C. Holm, The European Physical Journal Special Topics 227, 1789 (2019).

35. T. D. Nguyen, H. Li, D. Bagchi, F. J. Solis, and M. Olvera De La Cruz, Computer Physics Communications 241, 80 (2019).

36. A. Marin-Laflèche, M. Haefele, L. Scalfi, A. Coretti, T. Dufils, G. Jeanmairet, S. Reed, A. Serva, R. Berthin, C. Bacon, S. Bonella, B. Rotenberg, P. Madden, and M. Salanne, Journal of Open Source Software 5, 2373 (2020).

37. A. Coretti, C. Bacon, R. Berthin, A. Serva, L. Scalfi, I. Chubak, K. Goloviznina, M. Haefele, A. Marin-Laflèche, B. Rotenberg, S. Bonella, and M. Salanne, The Journal of Chemical Physics 157, 184801 (2022).

38. L. J. V. Ahrens-Iwers, M. Janssen, S. R. Tee, and R. H. Meißner, The Journal of Chemical Physics 157, 084801 (2022).

39. J. I. Siepmann and M. Sprik, The Journal of Chemical Physics 102, 511 (1995).

40. S. K. Reed, O. J. Lanning, and P. A. Madden, The Journal of Chemical Physics 126, 084704 (2007).

41. L. Scalfi and B. Rotenberg, Proceedings of the National Academy of Sciences 118, e2108769118 (2021).

42. G. Jeanmairet, B. Rotenberg, and M. Salanne, Chemical Reviews 122, 10860 (2022).

43. D. T. Limmer, C. Merlet, M. Salanne, D. Chandler, P. A. Madden, R. Van Roij, and B. Rotenberg, Physical Review Letters 111, 106102 (2013).

44. L. Scalfi, D. T. Limmer, A. Coretti, S. Bonella, P. A. Madden, M. Salanne, and B. Rotenberg, Physical Chemistry Chemical Physics 22, 10480 (2019).
(2020).

⁴⁵G. Pireddu and B. Rotenberg, Physical Review Letters **130**, 098001 (2023).

⁴⁶G. Pireddu, C. J. Fairchild, S. P. Niblett, S. J. Cox, and B. Rotenberg, Proceedings of the National Academy of Sciences **121**, e2318157121 (2024).

⁴⁷A. Arnold, J. De Joannis, and C. Holm, The Journal of Chemical Physics **117**, 2496 (2002).

⁴⁸S. Tyagi, M. Süzen, M. Sega, M. Barbosa, S. S. Kantorovich, and C. Holm, The Journal of Chemical Physics **132**, 154112 (2010).

⁴⁹K. Breitsprecher, K. Szuttor, and C. Holm, The Journal of Physical Chemistry C **119**, 22445 (2015).

⁵⁰O. Maxian, R. P. Peláez, L. Greengard, and A. Donev, The Journal of Chemical Physics **154**, 204107 (2021).

⁵¹F. Jiménez-Ángeles, A. Ehlen, and M. Olvera De La Cruz, Faraday Discussions **246**, 576 (2023).

⁵²A. P. Dos Santos, F. Jiménez-Ángeles, A. Ehlen, and M. Olvera De La Cruz, Physical Review Research **5**, 043174 (2023).

⁵³N. Pogharian, A. P. Dos Santos, A. Ehlen, and M. Olvera De La Cruz, The Journal of Chemical Physics **160**, 094704 (2024).

⁵⁴A. P. Dos Santos, M. Girotto, and Y. Levin, The Journal of Chemical Physics **147**, 184105 (2017).

⁵⁵M. Girotto, A. P. Dos Santos, and Y. Levin, The Journal of Chemical Physics **147**, 074109 (2017).

⁵⁶R. M. Malossi, M. Girotto, A. P. Dos Santos, and Y. Levin, The Journal of Chemical Physics **153**, 044121 (2020).

⁵⁷I. M. Telles, A. P. Dos Santos, and Y. Levin, The Journal of Chemical Physics **161**, 214102 (2024).

⁵⁸P. Cats, R. S. Sitlapersad, W. K. Den Otter, A. R. Thornton, and R. Van Roij, Journal of Solution Chemistry **51**, 296 (2022).

⁵⁹F. Grün, M. Jardat, P. Turq, and C. Amatore, The Journal of Chemical Physics **120**, 9648 (2004).

⁶⁰E. Kätelhön and B. Wolfrum, Physica Status Solidi (a) **209**, 881 (2012).

⁶¹K. J. Krause, K. Mathwig, B. Wolfrum, and S. G. Lemay, The European Physical Journal Special Topics **223**, 3165 (2014).

⁶²T. Hoang Ngoc Minh, G. Stoltz, and B. Rotenberg, The Journal of Chemical Physics **158**, 104103 (2023).

⁶³L. Thomas, Proc. Cambridge Phil. Roy. Soc. **23**, 542 (1927).

⁶⁴E. Fermi, Rend. Accad. Naz. Lincei **6**, 602 (1927).

⁶⁵A. A. Kornyshev, A. I. Rubinshtein, and M. A. Vorotyntsev, Physica Status Solidi (b) **84**, 125 (1977).

⁶⁶M. A. Vorotyntsev and A. A. Kornyshev, Soviet Journal of Experimental and Theoretical Physics **51**, 509 (1980).

⁶⁷C. C. Rochester, A. A. Lee, G. Pruessner, and A. A. Kornyshev, ChemPhysChem **14**, 4121 (2013).

⁶⁸A. A. Lee, D. Vella, and A. Goriely, EPL (Europhysics Letters) **113**, 38005 (2016).

⁶⁹A. A. Lee and S. Perkin, The Journal of Physical Chemistry Letters **7**, 2753 (2016).

⁷⁰J. G. Hedley, K. K. Bhatt, H. Berthoumieux, and A. A. Kornyshev, The Journal of Chemical Physics **162**, 114703 (2025).

⁷¹Y. Levin, Europhysics Letters (EPL) **76**, 163 (2006).

⁷²J. Comtet, A. Niguès, V. Kaiser, B. Coasne, L. Bocquet, and A. Siria, Nature Materials **16**, 634 (2017).

⁷³V. Kaiser, J. Comtet, A. Niguès, A. Siria, B. Coasne, and L. Bocquet, Faraday Discussions **199**, 129 (2017).

⁷⁴Y. Yu, K. Zhang, H. Parks, M. Babar, S. Carr, I. M. Craig, M. Van Winkle, A. Lyssenko, T. Taniguchi, K. Watanabe, V. Viswanathan, and D. K. Bediako, Nature Chemistry **14**, 267 (2022).

⁷⁵A. Schlaich, D. Jin, L. Bocquet, and B. Coasne, Nature Materials **21**, 237 (2022).

⁷⁶K. Goloviznina, J. Fleischhaker, T. Binninger, B. Rotenberg, H. Ers, V. Ivanistsev, R. Meissner, A. Serva, and M. Salanne, Advanced Materials **24**, 05230 (2024).

⁷⁷S. Nair, G. Pireddu, and B. Rotenberg, Molecular Physics **123**, e2365990 (2025).

⁷⁸S. Nair, G. Jeanmairet, and B. Rotenberg, The Journal of Chemical Physics **163**, 014107 (2025).

⁷⁹L. Coello Escalante and D. T. Limmer, Nano Letters (2024), 10.1021/acs.nanolett.4c04690, publisher: American Chemical Society.

⁸⁰H. Löwen, P. A. Madden, and J.-P. Hansen, Physical Review Letters **68**, 1081 (1992).

⁸¹L. Scalfi, T. Dufils, K. G. Reeves, B. Rotenberg, and M. Salanne, The Journal of Chemical Physics **153**, 174704 (2020).

⁸²W. A. Steele, Surface Science **36**, 317 (1973).

⁸³W. A. Steele, The Journal of Physical Chemistry **82**, 817 (1978).

⁸⁴J. J. Magda, M. Tirrell, and H. T. Davis, The Journal of Chemical Physics **83**, 1888 (1985).

⁸⁵J. J. Magda, M. Tirrell, and H. T. Davis, The Journal of Chemical Physics **84**, 2901 (1986).

⁸⁶M. Mazars, Molecular Physics **105**, 1909 (2007).

⁸⁷M. Mazars, Molecular Physics **105**, 1927 (2007).

⁸⁸R. Zwanzig and N. K. Ailawadi, Phys. Rev. **182**, 280 (1969).

⁸⁹K. A. Dwelle and A. P. Willard, The Journal of Physical Chemistry C **123**, 24095 (2019), https://doi.org/10.1021/acs.jpcc.9b06635.

⁹⁰J. Huang, Y. Chen, and M. Eikerling, Proceedings of the National Academy of Sciences **120**, e2307307120 (2023), https://www.pnas.org/doi/pdf/10.1073/pnas.2307307120.

⁹¹W. Schmickler and D. Henderson, The Journal of Chemical Physics **80**, 3381 (1984).

⁹²E. Leiva and W. Schmickler, Journal of Electroanalytical Chemistry and Interfacial Electrochemistry **229**, 39 (1987).

⁹³C. F. v. Weizsäcker, Zeitschrift für Physik **96**, 431 (1935).

⁹⁴R. Benguria, H. Brezis, and E. H. Lieb, Commun. Math. Phys. **79**, 167 (1981).

⁹⁵A. Schlaich, E. W. Knapp, and R. R. Netz, Physical Review Letters **117**, 048001 (2016).

⁹⁶S. J. Cox and P. L. Geissler, Chemical Science **13**, 9102 (2022).

⁹⁷H. Berthoumieux and F. Paillusson, The Journal of Chemical Physics **150**, 094507 (2019).

⁹⁸M. R. Becker, R. R. Netz, P. Loche, D. J. Bonthuis, D. Mouhanna, and H. Berthoumieux, Physical Review Letters **134**, 158001 (2025).

⁹⁹P. Loche, L. Scalfi, M. Ali Amu, O. Schullian, D. J. Bonthuis, B. Rotenberg, and R. R. Netz, The Journal of Chemical Physics **157**, 094707 (2022).

¹⁰⁰D. Horinek and R. R. Netz, Physical review letters **99**, 226104 (2007).

¹⁰¹L. Joly, C. Ybert, E. Trizac, and L. Bocquet, The Journal of Chemical Physics **125**, 204716 (2006).