Capacitance and Structure of Electric Double Layers: Comparing Brownian Dynamics and Classical Density Functional Theory

We present a study of the structure and differential capacitance of electric double layers of aqueous electrolytes. We consider Electric Double Layer Capacitors (EDLC) composed of spherical cations and anions in a dielectric continuum confined between a planar cathode and anode. The model system includes steric as well as Coulombic ion-ion and ion-electrode interactions. We compare results of computationally expensive, but"exact", Brownian Dynamics (BD) simulations with approximate, but cheap, calculations based on classical Density Functional Theory (DFT). Excellent overall agreement is found for a large set of system parameters $-$ including variations in concentrations, ionic size- and valency-asymmetries, applied voltages, and electrode separation $-$ provided the differences between the canonical ensemble of the BD simulations and the grand-canonical ensemble of DFT are properly taken into account. In particular a careful distinction is made between the differential capacitance $C_N$ at fixed number of ions and $C_\mu$ at fixed ionic chemical potential. Furthermore, we derive and exploit their thermodynamic relations. In the future these relations are also useful for comparing and contrasting.

experimental data with theories for supercapactitors and other systems. The quantitative agreement between simulation and theory indicates that the presented DFT is capable of accounting accurately for coupled Coulombic and packing effects. Hence it is a promising candidate to cheaply study room temperature ionic liquids at much lower dielectric constants than that of water.
Keywords Electrolytes · Electric Double Layer · Density Functional Theory · Brownian Dynamics · Differential Capacitance · Capacitors

Introduction
Electric double layer capacitors (EDLCs) are promising energy storage devices, in which electric energy is stored in the net ionic charge that is present in the vicinity of an electrode-electrolyte interface. In EDLCs the cathode attracts cations and repels anions and vice versa for the anode; more so the higher the applied voltage between the cathode and the anode [1]. This energy storage mechanism leads to much higher power densities than those of batteries; the discharge of the so-called Electric Double Layer (EDL) of an EDLC can be much faster than the redox reactions in batteries [2,3]. However, the energy densities of EDLCs are much lower than those of batteries [2].
One of the factors that contributes to the low energy density in EDLCs is the limited potential window in which conventional electrolytes are stable with respect to detrimental chemical reactions. Conventional electrolytes in EDLCs consist of a salt (e.g. tetraethylammonium tetrafluoroborate) dissolved in a solvent (e.g. propylene carbonate or acetonitrile) [4]. In order to maximise the energy density of EDLCs, one could use alternative electrolytes with a larger potential window [5]. Room temperature Ionic liquids (ILs) are potential alternatives to conventional electrolytes in EDLCs, since they can have a potential window of up to 6 V [6]; conventional electrolytes in EDLCs have a potential window of only 2.5 to 2.8 V [5]. Other advantages of ILs over conventional electrolytes are their stability at high temperatures and a low vapour pressure (low volatility and non-flammability), which makes IL-based EDLCs much safer [5,7]. Modelling and theoretically understanding concentrated ILs is difficult, because of the steric repulsions at short ionic separations and the Coulombic interactions at longer ranges are simultaneously at play. While in experiments dispersion forces, polarisation, and orientation degrees of freedom often also play a role, we restrict attention to the combined effects of packing and electrostatics. In order to prepare for the challenges posed by ILs, we here focus on the parameter regime of aqueous systems.
Several methods have been applied to investigate the electric double layer of EDLCs. On the one hand there are continuum methods, such as the mean field Gouy-Chapman-Stern (GCS) theory for point ions and classical Density Functional Theory (DFT) that can include steric effects. On the other hand there are methods in which each ionic constituent is treated explicitly, such as in Molecular Dynamics (MD) and Brownian Dynamics (BD) simulations.
All have their advantages and disadvantages. GCS theory can be solved analytically, but is rather inaccurate for larger ionic concentrations and surface charges, DFT is computationally fast but is an approximate theory for a given model. MD simulations might be considered as 'exact' and provide dynamics at the molecular scale but are computationally expensive. DFT calculations [8][9][10][11][12][13][14] and MD simulations are extensively used to study double layers in ILs and aqueous electrolytes [15][16][17][18][19][20][21][22][23][24][25][26][27][28][29][30][31][32][33]. As an alternative to the MD simulation method, BD is less accurate but computationally cheaper. In BD the explicit solvent of MD is eliminated by including solvent effects -like friction, Brownian noise and the dielectric constant of the medium -in an approximate way in the equations of motion of the ions [34].
The main aim of the paper is to explore whether DFT can be used to efficiently model electrolytes; therefore we are interested in the ability of DFT to quantitatively match BD simulations of a primitive model electrolyte. To this end, we study the EDLs of an aqueous electrolyte confined between a planar anode and cathode using BD simulations and classical DFT. A given potential difference is applied between the two electrodes, which are modeled as graphene-like electrodes in the BD simulations and as impenetrable planar walls in the DFT calculations. We compare both the ion concentration profiles between the electrodes (see also Refs. [35,36]) and the differential capacitance C. The differential capacitance characterises the (additional) charge per (additional) applied voltage, a quantity that is experimentally measurable and often used as a characteristic of the energy storage qualities of a capacitor.
Also it is shown that a careful distinction is required between the capacitance C N at constant number of ions and C µ at constant ion chemical potential, where the former follows naturally from the canonical BD simulations and the latter from the grand-canonical DFT calculations. However, they are related and we show the connecting expressions. We start by a detailed study of our reference system: a 1 mol·L −1 1:1 electrolyte of equal-sized ions of diameter d = 0.5nm in the 4nm gap between two electrodes at a potential difference of 0.2V. Then we vary the salt concentration, the ion valencies, the diameter ratio, the electrode-electrode distance, and the applied voltage. Throughout, the whole explored 5-dimensional parameter space /design of experiments we find excellent agreement between our BD and DFT results.

Model
We consider an aqueous electrolyte confined by two planar electrodes at fixed surface potential Φ L and Φ R , separated by a distance H. The electrolyte contains spherical cations (+) and anions (-) with a diameter d ± and valency Z ± dissolved in a structureless medium with dielectric constant ε = 78 at room temperature T = 298 K (see Fig. 1). The medium is fully characterized by its Bjerrum length λ B = βe 2 /4πεε 0 = 0.72 nm, where β = 1/k B T and ε 0 the dielectric permittivity of free space. The pair potential u ij (r) between a pair of ions of species i and j separated by a distance r is composed of a steric Fig. 1 Illustration of the electrolyte with cations (+) and anions (-) dissolved in a dielectric medium characterized by the Bjerrum length λ B = βe 2 /4πǫ 0 ǫ confined between two planar electrodes separated by a distance H at which a potential Φ L and Φ R is applied on the left and right electrode, respectively. The gradient background color indicates the charge density.
repulsions, characterized by the diameter d, and the Coulombic interaction, i.e.
In the DFT calculations we describe the steric repulsions with a hard-sphere potential whereas in the BD simulations we employ the Weeks-Chandler-Andersen (WCA) pair potential Here ǫ is the interaction parameter that we set to ǫ = k B T , and d ij = (d i + d j )/2. Note that the WCA potential, which is just the repulsive part of the Lennard-Jones potential, is only slighty softer for r < d than the hardsphere potential used in the DFT calculations. However, we will show remarkable agreement between the DFT and BD results in the parameter regime of study, indicating the limited sensitivity of the functional form of the repulsive interaction. By confining the system, we also introduce an external potential, i.e. the interaction of the ions with the fixed particles in the electrode. This ion-electrode interaction is described by the WCA potential and Coulombic interactions, in both DFT and BD. In the DFT description we integrate out the in-plane dimensions, finding for the WCA part of the external potential where ρ w is the surface density of wall particles in number of particles per unit area; ǫ w the interaction strength between the wall particles and ions, for which we take ǫ w = ǫ = k B T ; the contact distance is d wj = (d w + d j )/2 with d w and d j the diameter of the wall particles and ions of species j, respectively. For z > d wj the external potential vanishes. Note that the second electrode is described by the same interaction potential with z replaced by H − z. Our primitive model captures the key features of an electrolyte, within the practical conditions posed by DFT and BD respectively, thereby enabling a quantitative comparison between both methods. For simplicity, we use the same diameter for all ions and electrode particles in the reference system, namely d = 0.5 nm.

Brownian Dynamics
We focus here on the description of the reference system used in the BD simulations, as illustrated in Fig. 2; detailed information on the variations to this system are provided in Section 5. The simulations were performed in LAMMPS [37]. The salt concentration in the slit between the two electrodes is set to what it would have been if the slit was in equilibrium with a reservoir, with a salt concentration of 1 mol·L −1 at 0V. The salt concentration in the slit was determined by DFT calculations, which corresponds to an equal amount of anions and cations, N + = N − = 51. As described in Section 2, the pairwise excluded volume interactions are described by a WCA potential. The ions have valencies of Z + = +1 and Z − = −1, while the variable charges of the wall particles are determined by the Constant Potential Method (see below). The long-ranged Coulombic interactions are evaluated using the particle-particleparticle-mesh method (PPPM) [38], with a cut-off distance of 12 nm and a relative accuracy of 10 −6 in the forces. A correction term allows application of this 3D Ewald summation technique to the current slab geometry [39,40].
In the BD simulation the solvent is implicit and accounted for in the equation of motion. That is, the second order Langevin equation of motion of the ions reads as where x i is the position of the i th particle, m = 50 atomic mass unit (a.m.u.) denotes the mass of the ions (equal for all ions), ξ the friction constant, U the where δ ij is the Kronecker delta, δ(t − t ′ ) the Dirac delta and the angular brackets denote an average. This equation of motion is integrated in LAMMPS by combining the velocity-Verlet scheme [34,41], using a time step of 5 fs, with the Langevin option. Simulations were initiated by placing the ions, stacked into a simple cubic crystal lattice, in the slit. A simulation typically lasted for 200 ns, requiring about 2-3 days on 32 cores in parallel, with the first 10% serving as equilibration phase and the remainder as production run. The two parallel electrodes each consist of one graphene-like layer of 960 particles, covering an area of A = 100.6 d 2 , see Fig. 2. The inter-particle bond lengths in these hexagonal layers are taken as the usual carbon-carbon distance, 0.142 nm = 0.284d. The distance between the electrodes, as measured between the centres of the constituent particles, is H = 8d. All wall particles are frozen, i.e. are excluded from the equation of motion, because the elimination of their rapid vibrations permits the use of a larger time step. Their charges are calculated using the Constant Potential Method (CPM) by imposing a constant voltage difference of Ψ = 0.2 V between the two walls [26,42]. We use the implementation provided by Yang et al. [24] which can be used as a plug-in for LAMMPS. In brief, the CPM determines the charges of all wall particles at every simulation step by solving the linear set of equations that determines the potential at every wall particle, given the positions of all wall particles and ions. As the implementation of CPM for LAMMPS [24] does not take the relative permittivity of the medium into account, the charges of all ions were multiplied by 1/ √ ǫ and the potentials on the walls were multiplied with √ ǫ to reach the same effect -all charges reported below are corrected to refer to an aqueous medium with ǫ = 78. A side effect to this pragmatic inclusion of the solvent's dielectric constant is that not only ion-ion and ionelectrode Coulomb forces are scaled, but so are the Coulomb forces between electrode atoms; the latter is of no consequence, however, since these atoms are immobilized in the simulation. A brief comparison with the fixed charge method, in which all charges are permanently fixed, is provided in Appendix D for the system at zero voltage. The system is periodically repeated in space, with the box lengths in the two directions parallel to the walls dictated by the geometry of the lattice and the height perpendicular to the walls taken as three times the width of the slit. Adding additional layers to the electrodes, creating thin slabs of graphite, does not significantly affect the ion density profiles nor the average total charge of the electrodes; the differences are within the accuracy of the calculation.

Density Functional Theory
The starting point of classical DFT is the grand potential functional Ω of the density profiles ρ j (r) [43], which reads in our case where F id is the intrinsic Helmholtz free energy functional of the ideal gas, F HS ex the excess (over-ideal) Helmholtz free energy functional that deals with the hard-sphere interactions, F ES ex the excess Helmholtz free energy functional that deals with the electrostatic interactions, µ j the chemical potential of species j, and ρ j (r) the local density of species j. Here, is an intrinsic property of the system which depends on the temperature and the interparticle interactions, but not on µ j − V j ext (r). This grand potential functional has the property that it is minimized for a given µ j − V j ext (r) by the equilibrium density profile ρ j,0 (r), i.e. δΩ/δρ j | ρj,0 = 0, at which it equals the actual grand potential Ω introduced later in section 4. Minimizing the grand potential functional w.r.t. the density profiles results in the Euler-Lagrange equations Therefore, once an explicit form of F is constructed, one can find the equilibrium density profiles {ρ 0 } by solving Eq. (8). Although the ideal Helmholtz free energy functional is known exactly and given by with Λ j the thermal wavelength, the excess functional hinges on approximations. One excellent approximation for the hard-sphere functional F HS ex has been developed and goes by the name Fundamental Measure Theory (FMT) [44][45][46], of which we apply the White-Bear II version [47]. The functional that deals with the electrostatics F ES ex is in general more difficult due to the long-range nature of the interactions [8]. However, we use the functional based upon the Mean-Spherical-Approximation (MSA), which for the restrictive primitive model (RPM) reads [48] where c MSA is given by [49] c MSA (r) = and the charge density is defined by q(r) = j Z j ρ j (r). The electrostatic potential makes its entrance upon writing Eq. (10) as the sum of the meanfield electrostatic free energy and MSA corrections, details of which can be found in Refs. [8,48,50]. The electrostatic potential is consequently determined by the Poisson equation, where the constant surface potential is enforced as a boundary condition.
The parameter D = d+ 1/Γ is a length scale that results from MSA, where Γ = ( √ 1 + 2dκ − 1)/2d with κ = √ 8πλ B ρ r the inverse Debye length and 2ρ r the total ion concentration in the reservoir at a given chemical potential µ. For the expression of c MSA beyond the RPM, we refer to [51,52]. Note that c MSA depends on the reservoir concentration ρ r through the parameter D and therefore it also depends on the chemical potential µ. Due to the approximation for the Helmholtz free energy functional F ES ex , it now also depends on µ, which it formally should not depend on. As a result, the Maxwell relation introduced in the next section does not hold exactly and we have two 'routes' to calculate the adsorption Γ (see Appendix C). When we need to calculate the adsorption Γ we use the route Eq. (15) given in the next section for the remainder of the manuscript.
As a final note we mention that our DFT calculations assume planar symmetry in which the in-plane coordinates can again be integrated out, leaving the normal coordinate z perpendicular to the electrodes as the only spatial variable in the numerical calculations.

Thermodynamics
We will treat the above-mentioned model theoretically via DFT, and by BD simulations. However, firstly we discuss the thermodynamics of both methods.
Let us start by considering the ensemble of the BD simulations, which is a closed system with a fixed number of N + cations and N − anions in a volume V at temperature T confined between two planar electrodes of equal area A separated by a distance H and held at a surface potential difference Ψ (see Fig. 3). The system, i.e. electrolyte and the electrodes, together with the charge reservoirs is charge neutral j=± eZ j N j + Q L + Q R = 0. Applying a potential difference Ψ creates an electric field across the system. Since the ions in the system are mobile they will respond to this electric field, and because the electrolyte is confined they will create a charge density near both electrodes. These electrodes are connected to charge reservoirs with which they can exchange charges, such that the charge on the electrode is balanced with the charge density in the electrolyte. In other words, the charge in the electrolyte together with the average charge on the electrodes must vanish, i.e. j=± eZ j N j + Q L + Q R = 0, where Q L/R denotes the average charge. The corresponding thermodynamic potential for this system is the free energy F (N + , N − , V, T, Ψ, A, H), for which the differential form reads with S the entropy, p the pressure, µ ± the chemical potential of the cations and anions, respectively, A the surface area, and γ is the total surface tension (which has contributions from the electrode-electrolyte and in the case of EDLoverlap also from electrode-electrode interactions). We also introduced the force f between the two planar electrodes, were f /A is also referred to as the disjoining pressure. We will call F the Helmholtz free energy, even though it is only a Helmholtz free energy for the ionic species while it is actually a grand-canonical potential for the charge carriers in the electrodes. The DFT calculations, on the other hand, are performed at constant chemical potential µ ± instead of constant number of ions N ± . This implicitly means that the system can freely exchange ions with an ion reservoir. Also, both the electrode potentials are defined w.r.t. a grounded reservoir (see Fig. 3), i.e. the two electrodes are connected to separate charge reservoirs that are independently held at a constant potential. Hence, the potential Φ L and Φ R System System Fig. 3 Illustration of the thermodynamic ensembles applicable to the BD simulations (left) and DFT (right). The simulations are performed with a fixed number of particles N in a fixed volume V at a fixed temperature T , while the DFT calculations employ a fixed chemical potential µ or reservoir concentration ρr. In the simulations the potential difference between the electrodes is fixed at Ψ, while in the DFT calculations the potentials of both electrodes relative to the reservoir are fixed at Φ L/R for the left and right electrode, respectively. In both cases, the electrodes can exchange charge with charge reservoirs to maintain the imposed potential difference(s).
on the left and right electrode, respectively, generate an independent electric field, not only between the electrodes but also between the electrodes and the reservoir. The ions both within the system and in the reservoir respond to this electric field. The role of the charge reservoirs is the same in both ensembles. Since global charge neutrality of the system plus reservoirs still holds, one finds in equilibrium the system charge neutrality condition j=± eZ j N j + Q L + Q R = 0. The corresponding thermodynamic potential is the grand potential Ω(µ ± , V, T, Φ L , Φ R , A, H) with differential The distinction between F (N ± , V, T, Ψ, A, H) and Ω(µ ± , V, T, Φ L , Φ R , A, H) is crucial when comparing results form DFT (constant chemical potential) with BD simulations (constant number of ions).
For macroscopically large systems we can use volumetric and areal extensivity arguments to write Ω = −p(µ + , µ − )V + γ(µ + , µ − , Ψ L , Ψ R , H)A, where we drop the T dependence for convenience as we keep the temperature fixed throughout. Combining the resulting differential dΩ = −pdV − V dp + γdA + Adγ with Eq.(13) gives the Gibbs-Duhem equation for the volumetric terms and, for dH = 0, the Lipmann equation where σ L/R = Q L/R /A denotes the surface charge density and Γ j the adsorption of ions of species j onto both electrodes defined by Here, z denotes the coordinate describing the distance perpendicular to the parallel electrodes, ρ j (z) the local number density of species j at position z, and ρ j,r the reservoir concentration of species j which is dictated by the chemical potentials µ ± . See Appendix A for a more detailed derivation.
The main observable that we focus on in this manuscript is the differential capacitance (per unit area), which can either be obtained at constant number of particles C N or at constant chemical potential C µ , i.e.
where, for simplicity, we consider the RPM where N + = N − ≡ N and dµ + = dµ − ≡ dµ/2, which allows us to write Γ = (Γ + + Γ − )/2. On top of that, we also apply the same (but opposite) potential on both electrodes within the Ω ensemble, i.e. Φ R = −Φ L ≡ Ψ/2, which leads to the same (but opposite) surface charge σ L = −σ R ≡ σ on both electrodes. Interestingly, C µ and C N are related via expressions that are very similar to those between the constant-volume and constant-pressure heat capacities, namely (see Appendix A) where α can be taken equal either to α µ or α Ψ defined by In an exact theory the Maxwell relation α µ = α Ψ is satisfied identically, however our excess free-energy functional is approximate and does not identically satisfy the Maxwell relation. As a consequence our conversion between C N and C µ depends on the choice for α, and thus leads to an inconsistency. However, the numerical differences are limited, as we will see in section 5.7. We also defined which resemble (osmotic) compressibilities of the ions at constant σ and Ψ , respectively. The relations in this section allow us to compare and convert the capacitances at constant N (natural to BD simulation) and at constant µ (natural to DFT).
In the BD simulations, the differential capacitances were determined by three routes. Running a set of equilibrium simulations at a range of potential differences Ψ between the electrodes yielded the mean total charge difference between the electrodes, Q N,Ψ = (1/2) Q R − Q L N,Ψ . The angular brackets denote the canonical ensemble average at the indicated potential difference, which is evaluated in simulations as a time-average [34,41]. The differential capacitance at constant numbers of ions, see Eq. 16, is obtained by numerically differentiating the surface charge with respect to the potential difference, using the central difference formula. Alternatively, the same differential capacitance is extracted from the thermal fluctuations of the wall charge around the average, δQ = Q − Q N,Ψ , over a single equilibrium simulation [53,54], The intrinsic capacitance C 0 , which is a constant independent of Ψ and N , accounts for the thermal fluctuations of the atomic charges around the idealized constant-potential charges calculated by CPM. The numerical value of C 0 is obtained by the fitting procedure discussed in Appendix B.
The differential capacitance at constant chemical potential, see Eq. 16, is obtained from the BD simulations as where each simulation was preceded by a DFT calculation to establish the required numbers of ions N (µ, Ψ ) under the prevailing conditions. The same finite difference relations were used to calculate the capacitance by DFT. The calculations of C N via Eq.(18) will be referred to as C αµ N and C αΨ N , when using α µ and α Ψ , respectively.

Results
In this section, the density profiles and the differential capacitances obtained with BD and DFT are presented and compared. The reference system will be considered first, followed by an exploration of the impacts of its various parameters by varying them one by one. Various authors have previously used DFT and MD to study systems similar to the systems discussed below [8, 31-33, 48, 50, 55, 56]. Although those systems are not identical to the systems we consider, they do show very similar density profiles.

Varying ion concentrations
The first parameter to be varied is the concentration of the electrolyte in the reservoir, from ρ r = 0.1 mol·L −1 to ρ r = 5 mol·L −1 , corresponding to 9 and 252 ion pairs in the BD simulation, respectively. The latter concentration is similar to that of the ion concentration in ILs, while the dielectric constant and the shapes of the ions are markedly different in ILs. The results, see Fig. 5(a) and Fig. 5(b), reveal a slight difference at the peaks for both cases. The lower density shows only one peak, so we can consider this to be in the diluted limit. The higher density shows strong oscillations, due to the layering commonly found at high packing fractions, like the current η = 0.394. The wavelength and amplitude of the oscillations are well captured with DFT, as indicated by the good agreement between the two methods. Nevertheless, there are small discrepancies in the heights of the peaks and valleys. These are mainly due to the treatment of the hard-core potential by the FMT part of the functional [57] and the difference in the representation of the repulsive interactions as WCA in BD and as hard spheres in DFT. The difference between the two plots is in the ensemble being used to control the ion number densities in the slit. Because the anions are smaller than the cations, they both come closer to the electrode -thereby lowering their electrostatic energy -and pack at a higher density. Consequently, as illustrated in Fig. 5(c), when using DFT to impose reservoir concentration of ρ r = 1 mol·L −1 for both ions, the heights of the density peaks adjacent to both electrodes become unequal at the reference potentials of Φ L = −0.1 V and Φ R = 0.1 V, hence Ψ = 0.2 V. From these density profiles, the numbers of ions were calculated as N + = 51 and N − = 55, to the nearest integer, and these numbers were used in the BD simulations shown in the same plot. Figure 5(d) presents results for equal numbers of cations and anions in the slit, taken as the average of the two previous values: N ± = 53. To obtain the desired numbers of ions of each type in DFT, the potential of one electrode was fitted, at constant potential difference Ψ to arrive at Φ L = −0.1116 V and Φ R = 0.0884 V. The density peaks at both electrodes now resemble each other. Note the asymmetry in both plots in the density of the cations at the positive electrode versus that of the anions at the negative electrode, both in their distance from the wall relative to the other ion at the same wall and the distance to the electrode before reaching the constant density plateau in the center of the slit. For both equal and unequal numbers of ions, BD and DFT show good agreement in the density profiles. Figure 5 (e) and (f) present the comparison for smaller Ψ = 0.02 V and larger Ψ = 2 V potential differences with respect to the reference system. Because convergence of the density profile at the lower potential required a very long BD production run, the simulation was performed instead using Newtonian mechanics in combination with a Nosé-Hoover thermostat [34,41]. The thermostat works by rescaling the velocities of all ions at every time step, in such a ways as to recover the correct mean kinetic energy and kinetic energy fluctuations at the desired temperature, and therefore samples the Boltzmann equilibrium distribution also obtained by the BD simulation. The thermostatted method samples configuration space more efficiently by ignoring the slow Brownian motion, which affects the dynamical properties of the system but not the thermodynamic properties studied in this work.

Changing the Surface Potential
Applying a small potential on the electrodes causes a smaller charge density, i.e. the density profiles of the cations and anions are more similar. At high potentials, a large portion of the ions are adsorbed onto the electrodes, causing strong layering. The agreement between the simulations and DFT is excellent for the smaller potentials, while small deviations are observed for the larger potential near the second peak. The latter are a result of strong packing, where, as mentioned before, the WCA potential differs from the FMT approximation to the hard sphere potential in DFT.

Different Valencies
In Fig. 5(g) and (h), the valency of the cation is doubled to twice that of the anion, Z + = −2Z − = 2. Charge neutrality of the reservoir, j Z j ρ r,j = 0, implies that the anion concentration in the reservoir must be double the cation concentration. Here we choose the cation reservoir concentration to be the same as in the reference system, i.e. 2ρ r,+ = ρ r,− = 2 mol·L −1 . Because the symmetry between cations and anions is broken, ions exchange with the reservoir resulted in distinct number of ions in the DFT calculations of Fig. 5(g). The most interesting difference between the reference system and this case is the little hump in the anion density profiles around z = 1 nm. This effect is referred to as overscreening [58], where one finds a negatively charged layer of anions next to the positively charged first layer of cations adjacent to the cathode. Again, excellent agreement is observed between DFT and BD.

Slit width
Lastly, the distance between the electrodes H was varied from H = 1.5 nm in Fig. 5(i) to H = 12 nm in Fig. 5(h). Because a very long production run was required in the BD simulations, these were performed using the Nosé-Hoover thermostat rather than by Langevin Dynamics, as explained in subsection 5.4. At the lower slit width the EDLs overlap substantially, while for the large width the electrolyte acquires the flat distribution of a bulk fluid in the middle of the system. The DFT and BD results are again in excellent agreement. Fig. 6 (a) The average surface charge on the electrodes (σ L −σ R )/2 and (b) the corresponding differential capacitance, as function of the potential difference between the electrodes, for a system in thermal equilibrium with a reservoir at a salt concentration of 1 mol·L −1 , by DFT and BD calculations. The number of ion pairs in the BD simulations was determined by DFT, and varies with the potential.

Capacitance
The capacitance was calculated as a function of the surface potential difference Ψ , at both a constant reservoir concentration and at constant number of ion pairs in the slit. In the former case, which comes naturally to DFT, the DFT calculations at a concentration in the reservoir of ρ r = 1 mol·L −1 were used to determine the numbers of ions in the BD simulations. In the latter case, which comes naturally to BD, the number of ion pairs was fixed at N = 156 and the concentration in the reservoir was varied to reach the desired number of ions in DFT. The charge on the electrode surface in the former case is shown in Fig. 6(a), where the line represents the DFT calculations and the markers the BD simulations. Note that each simulation, although performed at constant numbers of ions, belongs to the same chemical potential. As expected, the charge on the wall and the number of ions in the slit increase with the potential difference between the electrodes. The two methods are in good agreement. The corresponding capacitance is presented in Fig. 6(b), where several calculation methods have been used. The blue solid line and the blue circles represent C ∆ µ Eq. (22) using the data in Fig. 6(a). Because the number of time-consuming BD simulations is necessarily low, the numerical derivative is limited in its accuracy, especially for the last data point at Ψ = 2 V. Nevertheless, the agreement is satisfactory and both methods yield similar camel-shaped curves [59]. Also shown in Fig. 6(b) are calculations of C N , where it should be emphasized that N is not constant across the plot but varies with Ψ . The capacitance C δ N (orange squares) is based on the charge fluctuations in the BD simulations, given in Eq. (21) where also C 0 appears. The value for C 0 , which depends neither on the number of particles nor on the potential difference, is found to be C 0 = 17.6 µF cm −2 (see Appendix B). The DFT calculations of C αµ N (green dotted line) and C αΨ N (green dashed line) are based on the relations in Eq. (17). The approximation made in Sec. 3.2, namely the assumed dependence of the direct correlation function c MSA in Eq. (11) on the chemical potential, resurfaces at this point. The adsorption Γ Fig. 7 The differential capacitance as a function of potential difference between the electrodes, using BD and DFT, for a system containing 156 ion pairs. Because the number of ion pairs is fixed, their chemical potential varies with the electrostatic potential difference.
can be obtained either from Eq. (15) or from the derivative of Eq. (14) with respect to the chemical potential, as derived in Appendix C. In the latter case, the derivative of F ES ex [{ρ}] does not vanish, though in principle it should have. The capacitance C N can therefore be calculated from C µ using either α Ψ or α µ , where the expression for Γ in Eq. (15) was used to calculate α Ψ . Note, however, that the calculation of the surface charge density σ is consistent by construction, since charge neutrality is imposed.
As expected from Eq. (16), in both cases C N is smaller than C µ . Both DFT calculations are in reasonable agreement with the BD results; notably, all three show a bell-shaped curve. The plot shows a substantial difference between C N and C µ , and although the various calculations do not exactly match quantitatively, they agree reasonably well and support the qualitative difference. The reason for the rather large difference between C µ and C N is due to the small electrode-electrode separation, where the region of the EDLs contribute substantially to the total number of ions in the system. In the limit where the separation between the electrodes is infinite, the difference between C µ and C N disappears.
Lastly, we consider the system with a fixed number of ion pairs, N = 156, and vary the surface potential difference Ψ . Shown in Fig. 7 are the capacitances using the same colour and line coding as in Fig. 6, with the addition of an orange curve for C ∆ N using DFT. For Ψ = 0 V to Ψ = 0.3 V, the simulations were run for 800 ns, treating the first 200 ns as equilibration phase, since for lower potential differences the simulations required a long production run for the capacitances to converge. Both differential capacitances are bell-shaped. The agreement between simulations and DFT is remarkably good, and compared to the results in Fig. 6(b), there is little to no qualitative difference between C N and C µ . The small qualitative difference, especially at small potential differences, is mainly due to the relative large number of ions in the system and therefore a corresponding large reservoir concentration. For comparison, the number of ions at Ψ = 0.2 V in Fig. 6 is 51, while only at Ψ = 2 V it is 156. The camel shaped curve in C µ is only existent for small reservoir concentrations ρ r < 1.5 mol·L −1 and bell shaped otherwise. Hence, no camel-shaped capacitance curve is observed within these parameters.

Discussions, Conclusions and Outlook
We presented ionic density profiles for a broad range of parameters applicable to aqueous electrolytes confined between a planar cathode and anode, and found very good agreement between results from DFT calculations and BD simulations. Both methods were also used to calculate differential capacitances, either C µ at constant ionic chemical potential µ or C N at constant number N of ions, via several routes. For a fixed chemical potential of mono-valent ions, at which the ionic reservoir concentration equals ρ r = 1 mol·L −1 , the capacitance curves obtained from DFT and BD are overall in good agreement. The DFT prediction for the capacitance at fixed N , however, gave two somewhat different results due to the approximation for the electrostatic part of the employed functional. Nevertheless, the DFT predictions bracket those of the BD simulations, except at potential differences between cathode and anode below 0.3 V where the simulations were extremely slow. Interestingly, the qualitative difference between C µ and C N is substantial, where C µ is camel shaped and C N is bell shaped. This has to do with the nonlinear relation between µ and N. Furthermore, C µ is found from a linear cut through the landscape in the three-dimensional space spanned by {σ, Ψ, µ}, whereas C N is the result of a non-trivial path through this landscape. We also considered capacitance curves at constant numbers of ions, N = 156, and found excellent agreement between DFT and BD simulations. In this case there is no qualitative difference between C µ and C N . Let us stress the time it takes to obtain the results from DFT and BD simulations. A typical BD simulation of a state point took 2 to 3 days on 32 cores, whereas the DFT calculations took not even 2 seconds on a regular laptop, which amounts to a difference of about 7 orders of magnitude. The accuracy that is lost by applying DFT on these systems is very small, as we have shown throughout this manuscript.
We conclude that with DFT one can obtain the same accuracy in structural and thermodynamic quantities as in BD simulations, at least for aqueous systems. This allows one to explore parameter space much more effectively and to study the properties of these systems thoroughly. A drawback of DFT is that it gives only an equilibrium description of the system, whereas BD simulations also provide the dynamics. We furthermore conclude that one needs to be careful and specify the differential capacitance that is being studied, e.g. C µ or C N , because they can differ both qualitatively and quantitatively. Also the natural choice changes depending on the method employed i.e. C µ for DFT and C N for BD, meaning the direct comparison of results from different methods is not straightforward.
A natural next step will be to divert from aqueous systems, to study systems with a lower dielectric constant. An interesting direction will be to study room temperature ionic liquid (ILs). Although the concentration of ρ r = 5 mol·L −1 in Fig. 5(b) is comparable to that of ionic liquids, the dielectric constant here is considerably higher due to the solvent. It is not sufficient to simply reduce the dielectric constant and redo the calculations, since the electrostatic correlations become much stronger at the low dielectric constants of ILs: a cation-anion pair of sub-nm diameter will bind at contact by Coulombic attractions of several tens of k B T . It is therefore not evident whether DFT or BD simulation will work in this regime. Moreover, the ionic shape in ILs is often non-spherical and needs to be accounted for in DFT. Interestingly, there have been developments in DFT to account for chain-like ions and molecules [60], and these have been applied to some extent to study ILs [9][10][11][12][13][14]. Besides chain-like ions, another approach to implement shape and polarizability is via molecular DFT [61][62][63][64]. Although molecular DFT has been mostly applied to model water, it might prove worthwhile to use this approach for ionic liquids in the future. Not only DFT is challenging at lower dielectric constant, but also simulations become much more challenging due to clustering that occurs at low dielectric constants as a result of the stronger electrostatic interactions. This leads to longer simulation times, which were already nonnegligible in the aqueous systems. Hence, a proper functional can provide the means to study ILs and electrolytes at low dielectric media effectively.
A different line of investigation would be the study of the differences of the differential capacitance C N and C µ . Until now the distinction has not often been made explicitly and further studies are needed to map out the properties and relations between both.

A Thermodynamics Derivation
For details see Ref. [65]. The differential for surface grand potential γ (or also called the surface tension) can be derived from taking the differential of Ω = −pV + γA and equating it with Eq. (13), i.e. dΩ = −pdV − V dp + γdA + Adγ One now needs to separate the volumetric bulk terms from the surface terms. For convenience, let us define N j = V ρ r,j + AΓ j and S = V sr + Ass, where ρ r,j and sr denote the particle density and the entropy density in the reservoir, respectively, while Γ j denotes the adsorption of species j and ss the areal excess entropy. This separation into volumetric and surface terms allows us to properly gather the volume terms on the left and the surface terms on the right of the equation, i.e.
Given that the surface has no influence on the volume term, both sides vanish and we obtain both the Gibbs-Duhem as well as the full Lipmann equation The systems that we considered when calculating the capacitances were at constant temperature T and constant electrode separation H. Moreover, in those calculations we considered the RPM such that Γ + = Γ − ≡ Γ/2 and dµ + = dµ − ≡ dµ and we symmetrized the surface potentials so that dΦ L = −dΦ R ≡ dΨ/2 and σ L = σ R ≡ σ. The Lipmann equation in this situation simplifies to Within the Ω ensemble σ(µ, Ψ) is a function of µ and Ψ, hence Because N (µ) is a function of µ, we therefore find that and the Maxwell relation Fig. 8 The difference in BD between the capacitance C N as obtained by numerical differentiation of the surface charge with potential, see Eq. (16), and as obtained from the thermal charge fluctuations only, see Eq. (21), confirms that the constant C 0 in the latter expression is indeed independent of the potential difference.
one can rewrite Eq. (29) as These relations allows us to relate the differential capacitance obtained from simulations C N to the differential capacitance from DFT Cµ. Within DFT one can calculate Cµ, χ Ψ and α Ψ , which through Eq. (17) gives access to C N . Notice the similarity with the relations between the heat capacity at constant pressure cp and the heat capacity at constant volume c V (see e.g. the thermodynamics book [66]). Although this derivation was done for the symmetric RPM, one can generalize these equations for any system. This is necessary when considering unequal ion sizes/valencies, but also when the potential on both electrodes differ. Hence, in general one needs to consider both electrodes separately.

B Calculation of the Capacitance
The calculation of a differential capacitance C N in BD using the fluctuation expression of Eq. (21) requires the evaluation of the constant C 0 accounting for the neglected thermal charge fluctuations around the idealized charges calculated by CPM. Because C 0 is a property of the electrodes that depends neither on the number of ions nor on surface potential, its value was determined as the difference between the C N obtained from the conventional charge-potential relation, see Eq. 16, and that obtained from the fluctuating contribution in Eq. (21) for C 0 = 0. This difference, plotted in Fig. 8 for a range of potentials, appears indeed to be a constant, with a value of C 0 = 17.6 µF cm −2 . The noise in the data is higher at low potentials, because the simulations converge more slowly at low potentials as well as due to the increased signal-to-noise ratio at the smaller step sizes in numerical differentiation.
(a) (b) Fig. 9 The inconsistency in the adsorption at (a) a fixed reservoir concentration of ρr = 1 mol·L −1 and at (b) a fixed number of ion pairs N = 156 (b).

C Adsorption Inconsistency
The adsorption can be calculated via the two routes Γ = − ∂γ ∂µ T,Ψ,H , and which are both plotted in Fig. 9 at a reservoir concentration of ρr = 1 mol·L −1 . As a short note, any functional that is based upon a bulk expansion like the one we use for the electrostatics suffers from this inconsistency.

D Constant potential versus fixed charge
The constant potential method (CPM) and fixed charge method (FCM) were compared by running simulations at zero voltage and zero charge, respectively. Based on DFT calculations, the equilibrium with a 1 mol·L −1 reservoir results in 47 ion pairs in the simulated slit. In the BD simulations the total charges on the left and right electrodes fluctuate around averages of ±7 nC·cm −2 , which is less then 1% of their standard deviations of 0.8 µC·cm −2 . Hence, the mean total charge on the electrodes is essentially zero. It should be noted that the slab option in LAMMPS can handle non-neutral systems [39,40]. The ionic number densities in both CPM and FCM simulations are similar, see Fig. 10. The slightly higher density peaks near the electrodes for CPM are probably caused by the ions inducing a mirror charge in the electrode and being attracted by that mirror charge; this effect evidently does not occur at fixed wall charges.