A semiclassical Thomas-Fermi model to tune the metallicity of electrodes in molecular simulations
Laura Scalfi, Thomas Dufils, Kyle Reeves, Benjamin rotenberg, Mathieu Salanne
AA semiclassical Thomas-Fermi model to tune the metallicity of electrodes inmolecular simulations
Laura Scalfi,
1, 2, 3
Thomas Dufils,
1, 2, 3
Kyle Reeves,
1, 2
Benjamin Rotenberg,
1, 2 and Mathieu Salanne
1, 2 Sorbonne Universit´e, CNRS, Physico-chimie des ´Electrolytes et Nanosyst`emes Interfaciaux, PHENIX,F-75005 Paris R´eseau sur le Stockage Electrochimique de l’Energie (RS2E), FR CNRS 3459, 80039 Amiens Cedex,France These authors contributed equally to this work
Spurred by the increasing needs in electrochemical energy storage devices, the electrode/electrolyte interfacehas received a lot of interest in recent years. Molecular dynamics simulations play a proeminent role inthis field since they provide a microscopic picture of the mechanisms involved. The current state-of-the-art consists in treating the electrode as a perfect conductor, precluding the possibility to analyze the effectof its metallicity on the interfacial properties. Here we show that the Thomas-Fermi model provides a veryconvenient framework to account for the screening of the electric field at the interface and differenciating goodmetals such as gold from imperfect conductors such as graphite. All the interfacial properties are modifiedby screening within the metal: the capacitance decreases significantly and both the structure and dynamicsof the adsorbed electrolyte are affected. The proposed model opens the door for quantitative predictions ofthe capacitive properties of materials for energy storage.
I. INTRODUCTION
The development of constant applied potential meth-ods for simulating electrochemical systems has allowedto solve many outstanding problems in physical electro-chemistry, ranging from the origin of supercapacitance innanoporous electrodes made of carbon or even of MetalOrganic Frameworks to the understanding of the dy-namic aspects of metal surface hydration . These meth-ods are based on the use of an extended Hamiltonian,in which the electrode charges are additional degreesof freedom that obey a constant potential constraint ateach simulation step . They allowed to partly alleviatethe main conceptual difficulty to represent the electrode-electrolyte interface at the molecular scale, which is theneed to account for the electronic structure on the elec-trode side, while the electrolyte is usually better sim-ulated using classical force fields because it requires asampling of the configurational space beyond the reach oftoday’s capabilities with ab initio calculations (see Ref. 6for a recent review of classical molecular simulations ofelectrode-electrolyte interfaces).Despite these successes, the possibility to simulaterealistic systems remains limited by the crudeness ofthe “electronic structure” model, since the electrode istreated as a perfect metal. It is however well knownthat the electronic response of different electrodes ( e.g. graphite vs. gold) to the adsorption of a charge shouldstrongly differ. This was shown in numerous analytical or density functional theory (DFT)-based studies , butalso more recently in an experimental study where strongdifferences in the confinement-induced freezing of ionicliquids were shown depending on the nature of the elec-trode . In the latter study, this effect was interpretedusing analytical developments accounting for the metal-licity of the system, in the framework of the Thomas-Fermi (TF) model . Here we build upon these developments to imple-ment a computational Thomas-Fermi electrode. The TFmodel is based on a local density approximation ofthe free electron gas, limited to its kinetic energy, and itaccounts for the screening of the electrostatic potentialover a characteristic screening length. We consider modelelectrodes with the gold structure and tunable metallic-ity, separated by either vacuum or a simple NaCl aque-ous electrolyte. We show that both the total accumu-lated charge and its distribution within the electrode arestrongly affected. Accounting for screening in the elec-trodes radically changes their response to the adsorptionof the electrolyte, which results in noticeable differencesin the structure of the liquid when a voltage is applied.Screening inside the metal should therefore be accountedfor when simulating electrochemical interfaces, in appli-cations ranging from supercapacitors to Li-ion batteries. II. THE THOMAS-FERMI ELECTRODE MODEL
We consider an electrode composed of N s sites (herethese sites are positioned on the nuclei) with a num-ber density d . Each atom i has Z valence electrons,and we introduce its partial charge q i as a dynamicalvariable accounting for the local excess of electrons. Asshown schematically on Figure 1, in the currently avail-able method the charges fluctuate in time to representperfect metals. The partial charges are calculated at eachsimulation step in order to ensure that the potential isthe same within the whole electrode ; when such anelectrode is put in contact with an electrolyte the screen-ing occurs within a thin layer at the surface only (notethat supercapacitors are often simulated using constantcharge setups, in which the vector { q i } i ∈ [1 ,N s ] containsprescribed (usually identical) values for all the atoms ofeach electrode and does not vary with time, which doesnot correspond to a realistic electrode). Nevertheless, a r X i v : . [ c ond - m a t . m t r l - s c i ] S e p homas-FermiscreeningNon ideal metalPerfectscreening Ideal metal
FIG. 1. Electrode polarization with different simulationmethods. Constant potential simulations (left) correspondto a perfect screening of the charges, hence to the behavior ofan ideal metal, whereas the Thomas-Fermi model introducesa screening length to account for the imperfect screening ofthe charge in a non-ideal metal. many electrode materials have a finite density of statesavailable at the Fermi level. This was sometimes consid-ered in the literature by computing a so-called quantumcapacitance that accounts for the corresponding screen-ing .Here we propose to take these effects into account di-rectly within classical molecular dynamics simulations,by employing the Thomas-Fermi model. It consists in alocal density approximation of the energy of the valenceelectrons. The Thomas-Fermi functional for the kineticenergy reads U T F [ n ( r )] = (cid:90) (cid:126) m e (3 π ) / n ( r ) / d r , (1)where n ( r ) is the local number density of electrons, m e their mass and (cid:126) Planck’s constant. In order to obtaina practical description in molecular simulations, we nowexpress n ( r ) as a sum over discrete atomic sites i , withlocal densities n i = d (cid:104) Z + q i ( − e ) (cid:105) , with e the elementarycharge. If the perturbation in the number of free chargecarriers is small compared to the number of electrons, i.e. | q i | (cid:28) Ze , we can expand the kinetic energy to secondorder in powers of q i as U T F = 35 N s ZE F + E F ( − e ) N s (cid:88) i =1 q i + l T F d (cid:15) N s (cid:88) i =1 q i (2)where E F = (cid:126) k F / m e is the Fermi level of a free-electron gas of density Zd and l T F = (cid:112) (cid:15) (cid:126) π / ( m e e k F )is the Thomas-Fermi length of the material, with the cor-responding Fermi wavevector defined by k F / π = Zd and (cid:15) the vacuum permittivity. The zeroth-order termis the total kinetic energy of an electron gas with N s Z electron (the total number of electrons in the system).The first order corresponds by definition to the chemical potential of the added/removed electrons (depending onthe sign of q i ). The second order term, which is alwayspositive and reaches its minimum when all the partialcharges vanish corresponds to an energy penalty to in-duce non-homogeneous charge distributions.Our system consists of two electrodes, hereafter namedafter their position in the simulation cell: left (L) andright (R). Their atom indices respectively range between[1 , N L ] and [ N L + 1 , N L + N R ], their Thomas-Fermi en-ergies are noted U LT F and U RT F , and they are held at po-tentials Ψ L and Ψ R = Ψ L + ∆Ψ where ∆Ψ is the appliedvoltage. We assume for simplicity that the electrodes aremade of the same material, hence they have the sameFermi level at rest. The total energy of the system reads E tot = K + U C + U vdW + U LT F + U RT F − N L (cid:88) i =1 Ψ L q i − N L + N R (cid:88) i = N L +1 Ψ R q i , (3)where K is the kinetic energy of the electrolyte, U C cor-responds to the Coulombic interactions, U vdW describesthe van der Waals interactions (given by a force field),while the last two terms account for the reversible worknecessary to charge the electrode atoms. U C reads U C = 12 (cid:90) (cid:90) ρ ( r ) ρ ( r (cid:48) )4 π(cid:15) | r − r (cid:48) | d r d r (cid:48) , (4)where the charge distribution ρ ( r ) consists in a col-lection of M point charges for the electrolyte and of N = N L + N R atom-centered Gaussians (with width η − )representing the electrodes: ρ ( r ) = M (cid:88) j =1 q j δ ( r − r j ) + N (cid:88) i =1 q i η π − / e − η | r − r i | , (5)with δ the Dirac function. Note that in Eq. 4 the onlyself-energy to be included is the one due to the Gaussiancharges. By injecting Eq. 2 into Eq. 3 and introducing∆Ψ, the total energy can be rewritten as E tot = K + U C + U vdW + 35 N ZE F + l T F d (cid:15) N (cid:88) i =1 q i − (Ψ L + E F e ) N (cid:88) i =1 q i − ∆Ψ N (cid:88) i = N L +1 q i . (6)As detailed in Ref. 6, in the absence of electrochemi-cal reactions, we impose the electroneutrality constraint (cid:80) Ni =1 q i = 0, so that the electrodes bear opposite chargesand the corresponding term in the reversible work re-duces to the usual Q tot ∆Ψ, with Q tot the total chargeof the positive electrode. As in the constant potentialmethod neglecting the quantum nature of the electrons(corresponding to l T F = 0 . ∂E tot /∂q i = 0 . Compared to this per-fect metal case, the modifications of the algorithm areminimal and virtually don’t add any computational cost.Our approach, which involves fluctuating charges, maybe related to the charge equilibration model , in par-ticular to its extension to electrochemical systems pro-posed by Onofrio et al. . This method is based on twomain chemical quantities, the electronegativity χ and thehardness H of each atomic species. The self-consistentequations to solve are equivalent if we take χ ∼ E F and H ∼ e l T F d/(cid:15) . However, these concepts, whichare related to those of electronic affinity and ionizationenergy , are rooted in the description of the electronicproperties of atoms and molecules, rather than that ofbulk materials, which are more naturally described interms of band structure. The issue of starting fromthe correct reference state for (electro-)chemical potentialequalization methods was already pointed out in Ref. 22,where York and Yang derived a fluctuating charge modelfrom DFT for molecules and underlined the difference be-tween atomic and molecular reference states to determinethe electronegativities and hardnesses. More recently, adetailed discussion on the correspondance between con-stant potential electrode models and the charge equilibra-tion approach was provided in Ref. 23. Another physicalmodel of electrodes was proposed , in which the Hamil-tonian is constructed in the tight-binding approximation. III. EMPTY CAPACITOR
As a first validation of our implementation, we studya model system composed of two planar (100) gold elec-trodes separated by a distance L and held at a constantpotential difference ∆Ψ = 1 V. Each electrode consistsof n atomic planes with an inter-spacing a in the z di-rection. We compare the simulated results against ana-lytical predictions of the corresponding continuum modelwhere the Poisson equation for the one dimensional po-tential Ψ( z ) reads Ψ (cid:48)(cid:48) ( z ) = l − T F Ψ( z ) inside each electrodeand Ψ (cid:48)(cid:48) ( z ) = 0 between them. The total capacitance ofthe system is given by C = Q tot / ∆Ψ.Assuming that the width of the material is large com-pared to the Thomas-Fermi length, the in-plane charge Q k at z = ka ( k ∈ [1 , n ]) can be expressed as Q k Q tot = e − ( k − a/l TF (cid:104) − e − a/l TF (cid:105) , (7)Figure 2a shows a very good agreement between Eq. 7and the simulation for large l T F values. Small deviationsfor large z are due to the finite number of planes. Theabove exponentially decaying charge distribution insidethe metal, due to the screening over the Thomas-Fermilength l T F , results according to the continuous model ina capacitance per unit area1 C EC = 1 C vac + 2 C T F = L vac (cid:15) + 2 l T F (cid:15) , (8) with C vac = (cid:15) /L vac the theoretical capacitance per unitarea for perfect metallic electrodes ( l T F = 0 ˚A) sepa-rated by a vacuum slab of width L vac and C T F = (cid:15) /l T F that for a single Thomas-Fermi electrode. This resultcan be simply understood in terms of the equivalent cir-cuit (hence the subscript C EC ) illustrated in Figure 2b,with three capacitors in series (see Supplementary Sec-tion S1 for a discussion of the continuum descriptions andequivalent circuit models). As shown in Figure 2c, thesimulation results are consistent with the prediction of alinear relation between 1 /C and L/(cid:15) , where L is the dis-tance between the first atomic planes on each electrode,with a constant shift which increases with l T F .However, the width of the vacuum slab between theelectrodes is not exactly the distance between the firstatomic planes. Indeed, each atomic site is surroundedby electrons, and the boundary between the free elec-tron gas inside the electrode and the vacuum (the so-called “Jellium edge” ) is rather shifted half of the inter-plane distance away from the electrode. Since this fea-ture is present on both electrodes, the actual vacuumslab width is more consistent with L vac = L − a . Fig-ure 2d shows that using this prescription, Eq. 8 providesa very good description of the simulated capacitance C over a wide range of distances between the electrodes andThomas-Fermi lengths, which confirms the consistency ofthe present classical model to represent the charge distri-bution within the metal. The decay length of the chargeinside the electrode coincide with l T F within 1 % for allvalues l T F (cid:38) a . The slight deviations from the predic-tions of the continuous theory can be analyzed by intro-ducing an effective length l eff from the measured capac-itance as 1 C = L − a(cid:15) + 2 l eff (cid:15) . (9)The results obtained for various l T F at fixed L , illustratedin Figure 2e, indicate that this effective length deviatesfrom the Thomas-Fermi length only when the latter be-comes comparable to the atomic details of the electrodes(interplane and interatomic distances, width of the Gaus-sian distributions). An additional test was performed byadding a single charge at various distances between theelectrodes and comparing the energy of the system to anapproximate analytical expression . The results, whichare provided in Supplementary Section S2, also show agood agreement over a broad range of l T F values.
IV. IMPACT OF THE THOMAS-FERMI LENGTH ONTHE ELECTROCHEMICAL INTERFACE PROPERTIES
In order to understand the impact of screening insidethe metal on the properties of electrode/electrolyte inter-faces, we study a system consisting of two (100) gold-likeelectrodes in contact with an aqueous solution of NaCl(with concentration 1 mol L − ), illustrated in Figure 3a.The TF length l T F was systematically varied from 0.0 to
IG. 2. Empty Thomas-Fermi capacitor. All results correspond to a (100) gold-like electrode structure with n = 50 atomicplanes and L = 300 ˚A between the electrodes where not stated otherwise. Charges were computed by applying a voltage∆Ψ = 1 V between the electrodes for different Thomas-Fermi lengths l TF ranging from 0.0 to 16.0 ˚A, that are representedboth by different symbols and by different colors indicated in the colorbar. (a) Total charge per plane on the positive electrodeas a function of the position from the surface ( k is the index of the atomic plane), normalized by the total electrode charge Q tot . The symbols are simulated values for different Thomas-Fermi lengths l TF and the lines are the prediction of Eq. 7. (b)Snapshot of the simulated system and its equivalent circuit representation corresponding to the capacitance obtained with thecontinuum theory (see text). (c) Computed reciprocal capacitance as a function of the analytical predictions for perfect metalsusing L vac = L , and (d) for Thomas-Fermi metals using Eq. 8 with L vac = L − a , for varying electrode spacing L (between 10and 200 ˚A). (e) Effective length l eff , defined in Eq. 9, as a function of l TF . ). Simulations were performed for voltages∆Ψ =0, 1 and 2 V between the two electrodes.As a first illustration of the impact of screening on theelectrochemical interface, we compute the Poisson poten-tial across the cell. The results for an applied potentialof 2 V are displayed on Figure 3b. We observe a verydifferent pattern inside the electrode depending on l T F :for the perfect metal the applied potential is reached atpositions corresponding to the first atomic plane, whilefor the TF model we clearly see the desired effect of fieldpenetration with an exponential decay inside the elec-trode.Figure 3c shows that the integral capacitance decreasessignificantly with l T F (note that it remains constant be-tween 1 and 2 V, see Supplementary Figure S2). Theeffect is already non negligible for l T F = 0 . C metal using the equivalent circuit depicted on Figure 3a(see Supplementary Section S1). This approach, used forexample by Gerischer to interpret experimental data ,has been applied in many simulation works where the ad-ditional term due to the screening was computed usingDFT and therefore called “quantum capacitance”, whilethe perfect metal capacitance was computed using eithera mean-field theory or molecular dynamics . Never-theless, it neglects the interplay between the electronicstructure of the electrode and the ionic structure of theadsorbed electrolyte. This coupling is self-consistentlytaken into account in our model, which therefore pro-vides a perfect framework to test this approximation. Ascan be seen on Figure 3c, the equivalent circuit approx- IG. 3. The capacitance decreases significantly with Thomas-Fermi length. (a) Snapshot of the simulated system and itsequivalent circuit representation, where C metal stands for the capacitance computed for the perfect metal simulation. (b)Poisson potential across the simulation cell for a system made of two (100) gold-like electrodes in contact with a NaCl aqueoussolution. The applied voltage is 2 V and different l TF values ranging from 0.0 to 5.0 ˚A are represented by different colorsindicated in the colorbar. The screening of the potential inside the electrodes increases markedly with l TF . (c) Variation of thecapacitance with l TF . The results from the simulations are compared with the equivalent circuit approximation. Error barsare extracted from the standard error of the charge distribution corrected for sample correlations. imation underestimates rather significantly the real ca-pacitance (by 20 to 30 %).At null voltage, the average structure of the liquid doesnot vary significantly with l T F (see Supplementary Fig-ure S3). As shown on Figure 4a-b, it is characterizedby several adsorption layers, mainly consisting of watermolecules. By computing the distribution of the angle θ between the vector normal to the surface and the waterdipole (see dashed black curve on Figure 4c-d) or the O-Hbonds (see Supplementary Figure S4) for molecules in thefirst adsorbed layer, we observe that they mostly lie in aplane parallel to the surface or with one H atom pointingaway from the surface. A small population is orientedtowards the surface, which results in a small shoulder onthe H atoms atomic density profiles.The ions have different adsorption profiles: the Na + density is characterized by a large peak located close tothe one of O atoms, so that they can be considered tobelong to the first layer, while the Cl − ions are locatedfurther away from the electrode surface. Their profiledisplays a small pre-peak in the region where the waterdensity is very low and a peak with a larger intensitylocated in the second hydration layer. Once a potentialis applied, the liquid mainly responds on the two elec-trodes through (i) a stronger orientation of the watermolecules towards/away from the negative/positive elec-trode as shown on Figure 4c-d and Supplementary Fig-ures S4 and S5, and (ii) the appearance of a new adsorp-tion peak for the Na + ions near the negative electrode(Figure 4e) and an increase of the pre-peak intensity inthe Cl − density profiles on the positive electrode side(Figure 4f). In all cases, the modifications in the struc-ture depend strongly on l T F . This shows that depend-ing on the type of material, we can expect all the elec-trochemical double-layer properties to change markedlywith the nature of the chosen electrode.Dynamical properties are particularly important for electrochemical applications. They control for examplethe power delivered by an energy storage device. Theequilibrium fluctuations of the electrode charge at 0 V,which reflect the linear response to a small applied volt-age, are shown on Figure 5 for the various l T F . An in-creased screening yields faster dynamics for the relax-ation of the electrochemical double-layer. Such a differ-ence was somewhat unexpected given that the systemsat null potential have on average the same structural fea-tures, but it can be qualitatively understood as the resultof weaker interactions with the more diffuse charges in-duced within the electrode. This means that the dynam-ics do not only depend on the nature of the electrolyte,but also on the electronic structure of the electrode ma-terial.
V. CONCLUSION
Understanding the electrode/electrolyte interface is aprerequisite not only for the design of more efficient en-ergy storage devices , but also for understanding wet-ting phenomena involved in lubrication or heterogeneouscatalysis . Although in the past decades molecular sim-ulations have provided many insights on the structure ofthe electrochemical double-layer, they still fail at predict-ing quantitatively many experimental quantities, such asthe variation of the differential capacitance with the ap-plied voltage . This is particularly true in the case ofcarbon materials, due to their complex electronic struc-ture properties that deviate largely from the ones of typ-ical metals. Many intriguing experimental observations,such as the capillary freezing of ionic liquids confinedbetween metallic surfaces or the emergence of longer-than-expected electrostatic screening lengths in concen-trated electrolytes , remain to be explained quanti-tatively. The Thomas-Fermi model, by allowing to tune +- + FIG. 4. The structure of the electrochemical interface depends on the Thomas-Fermi length at finite voltages. (a,b) Atomicdensity profiles for the O, H, Na + and Cl − atoms near the electrode at null potential for l TF = 0 . l TF values as shown on Supplementary Figure S3). Note that in the case of H atoms the profile is dividedby two to facilitate the comparison with O atoms. (c,d) Distribution of the adsorbed water molecules orientation with respectto the vector normal to the electrode surfaces for an applied potential of 2 V for the whole range of simulated l TF indicatedby the colorbar; the distribution for 0 V and l TF = 0 . + and Cl − ions for an applied potential of 2 V for the whole range of simulated l TF indicated by thecolorbar. The negative (positive) electrode is located at negative (positive) z . the metallicity of the electrode using a single parame-ter (and without introducing additional computationalcosts) should lead to a more accurate understanding ofthe interfacial properties of such electrodes using molec-ular simulations. The extension of this work to complexmaterials such as nanoporous carbons will require addi-tional efforts, in order to take into account the effect ofthe local environment of each atom on its electronic re-sponse. In that case, it might be relevant to sacrifice someof the simplicity of the TF model by including atom-specific or even bond-specific terms in the energy, fol-lowing e.g. the split charge equilibration approach .In this context, the present work suggests that it couldbe possible to determine the associated parameters from a simplified representation of the underlying electronicdensity. APPENDIX : SIMULATION DETAILS
The TF electrode model was implemented in themolecular dynamics code MetalWalls . All simulationswere run using a matrix inversion method to enforceboth the constant potential and the electroneutrality con-straints on the charges. Electrode atoms have a Gaus-sian charge distribution of width η − = 0 .
56 ˚A centeredon zero and the Thomas-Fermi length l T F ranges from0.0 to 16.0 ˚A for the empty capacitor and from 0.0 to
IG. 5. The relaxation of the electrode charge indicates afaster dynamics of the interfacial electrolyte near screenedmetals. Normalized auto-correlation function of the totalcharge at null potential for varying l TF values ranging from0.0 to 5.0 ˚A indicated by the colorbar. z direction using an accurate 2D Ewaldsummation method to compute electrostatic interactions.A cutoff of 17.0 ˚A was used for both the short range partof the Coulomb interactions and the intermolecular in-teractions. For the latter we used the truncated shiftedLennard-Jones potential. The box length in both the x and y directions was L x = L y = 36 .
630 ˚A with 162atoms per atomic plane. The structure is face-centeredcubic with a lattice parameter of 4 .
07 ˚A and a separa-tion between planes a = 2 .
035 ˚A in the (100) direction(the atomic density d is 0.59 · m − ). The emptycapacitors have 50 planes per electrode whereas the elec-trochemical cells have 10 (leading to a total of 1620 atomsper electrode). In the latter case, the electrolyte is com-posed of 2160 water molecules, modeled using the SPC/Eforce field , and 39 NaCl ion pairs. The Lennard-Jonesparameters for Na + and Cl − were taken from Ref. 38and the ones for the electrode atoms from Ref. 39; theLorentz-Berthelot mixing rules were used. The simu-lation boxes were equilibrated at constant atmosphericpressure for 500 ps by applying a constant pressure forceto the electrodes with l T F = 0 . L = 56.8 ˚A. The simulations were run at298 K with a timestep of 1 fs. Each system was run forat least 8 ns. ACKNOWLEDGEMENTS
The authors thank M. Sprik, P. A. Madden, L. Bocquetand B. Coasne for useful discussions. This project hasreceived funding from the European Research Councilunder the European Union’s Horizon 2020 research and innovation programme (grant agreement No. 771294).This work was supported by the French National Re-search Agency (Labex STORE-EX, Grant ANR-10-LABX-0076, and project NEPTUNE, Grant ANR-17-CE09-0046-02).
SUPPLEMENTARY INFORMATIONS1. CONTINUUM DESCRIPTION AND EQUIVALENTCIRCUIT MODELS
As discussed in the main text, the capacitance of theelectrochemical cell comprising the electrodes, the elec-trolyte and the interfaces between them is often analysedin term of a simple model based on capacitors in series.Considering the symmetry of the problem, the mean-fieldPoisson potential Ψ only depends on the position z in thedirection perpendicular to the electrodes and satisfies oneof the following equations: • Thomas-Fermi: Ψ (cid:48)(cid:48) ( z ) = 1 l T F Ψ( z ) inside the elec-trodes described by the Thomas-Fermi model, with l T F the Thomas-Fermi length of the material. • Debye-H¨uckel: Ψ (cid:48)(cid:48) ( z ) = 1 λ D Ψ( z ), with λ D = (cid:16) e (cid:15) (cid:15) r k B T (cid:80) i c i z i (cid:17) − / the Debye screening lengthfor a dilute electrolyte, with e the elementarycharge, k B T the thermal energy, (cid:15) r the relativepermittivity of the solvent, c i the concentration ofionic species i and z i its valency. This equation isthe linearized version of the more general Poisson-Boltzmann equation and is only valid for small ap-plied voltages. • Poisson: Ψ (cid:48)(cid:48) ( z ) = − ρ q ( z ) (cid:15) in the more general case,with ρ q the charge density which can be obtainedfrom the density profiles of all species (includingthe contributions of the O and H of water in theaqueous systems considered in the present work).These equations need to be solved in the vari-ous regions of the system, with appropriate continu-ity equations at the boundaries between them andoverall boundary conditions lim z →−∞ Ψ( z ) = Ψ L andlim z → + ∞ Ψ( z ) = Ψ R , with Ψ L and Ψ R the potentialsof the left and right electrodes, respectively. The capac-itance of the system, defined as the ratio between thecharge accumulated inside the electrodes and the voltage∆Ψ = Ψ R − Ψ L , can then be expressed as1 C = (cid:88) k C k (S1)with C k the capacitance corresponding to each region k .This rule allows to determine the relative contributionf each region as a function of the composition (elec-trode material, electrolyte) and geometry (distance be-tween the electrodes) of the system.In both the Thomas-Fermi electrode and the Debye-H¨uckel electrolyte, the potential behaves as a sum of twoexponentials. If the width of the corresponding regions(electrode slab or electric double layer) are large com-pared to the associated screening lengths, then only oneexponential contributes and the capacitance is given by C T F = (cid:15) l T F (analogous to the “quantum capacitance”sometimes introduced to capture the contribution of theelectrode) and C DH = (cid:15) λ D , respectively. For an emptycapacitor, the vacuum slab between the electrodes cor-responds to a capacitance C vac = (cid:15) L vac , with L vac thewidth of the slab, i.e. the distance between the surfaceof the electrodes (which slightly differs from the differ-ence between the position of the first atomic planes oneach surface, as explained in the main text).For the pure water case, one should distinguish thecontribution of the bulk dielectric liquid, with capaci- tance C bulk = (cid:15) (cid:15) r L bulk , and that of the structured layers atthe interface, which can be associated with an effectiveinterfacial capacitance C w,int . This interfacial contribu-tion of water is also present in the aqueous electrolyte,and the ionic contribution may be estimated with C DH in the dilute and low potential limits, from the non-linear Poisson-Boltzmann capacitance, or directly fromthe ionic concentration profiles from molecular simula-tions. S2. SINGLE CHARGE BETWEEN TWO ELECTRODES
We consider a model system consisting of a unit pointcharge between two graphite electrodes held at constantpotential difference ∆Ψ = 0 V, and compute the energyof the system as a function of the charge position z . Theenergy of a particle of charge q interacting with an infinitecontinuous Thomas-Fermi metal with a plane surface hasbeen derived in Ref 12, which also provides the followinganalytical ansatz U T F ( z ) = − q π(cid:15) z [1 − (S2)13 . k T F z ) + 37 . k T F z ) + 18 . k T F z ) + 127 . k T F z ) + 73 . k T F z ) + 70 . k T F z ) + 20 . k T F z ) + 1 (cid:21) where k T F = 1 /l T F and z is the distance from the sur-face. Note that for large k T F , we recover the result ob-tained from the image charge method . To adapt thisexpression to our setup, we assumed the lateral dimen-sions of the box to be large enough to neglect the effect ofthe periodic images and the two electrodes to be decor-related in order to simply sum both contributions. Wefirst consider the position of the charge with respect tothe atomic position of the electrodes, leading to the fol-lowing simple expression for the energy of the capacitorcomposed by a single charged particle and the two elec-trodes U c = U T F ( z ) + U T F ( L − z ) (S3)where L is the distance between the electrodes. The re-sults are shown on Figure S1. The agreement between thesimulations and the analytical calculation is good both atlong and short ranges and we recover the energy depen-dence as a function of the Thomas-Fermi length. How-ever, it should be pointed out that such a good compar-ison can only be achieved using a box with large lateraldimensions (around 200 ˚A) compared to the electrodesseparation L = 50 ˚A. This is because, for large distances z , the effect of the periodic images of the point chargealong the lateral directions x and y on the induced chargedistribution within the electrodes cannot be neglected, so that the analytical prediction for an isolated ion is notappropriate. In addition, our analytical expression as-sumes a superposition of the contributions due to eachelectrode treated independently, neglecting in particularthe fact that the total charge induced on each electrode isnot a full elementary charge but depends on the positionof the point charge. Finally, for short distances, thereis a threshold below which the continuum approximationbreaks down and the molecular structure of the electrodeplays a role. S3. CAPACITANCES OF THE GOLD-LIKEELECTRODES AT 1 AND 2 VS4. ADDITIONAL STRUCTURALCHARACTERIZATIONS J. I. Siepmann and M. Sprik, “Influence of Surface-Topology andElectrostatic Potential on Water Electrode Systems,” J. Chem.Phys. , 511–524 (1995). C. Merlet, B. Rotenberg, P. A. Madden, P.-L. Taberna, P. Simon,Y. Gogotsi, and M. Salanne, “On the Molecular Origin of Su-percapacitance in Nanoporous Carbon Electrodes,” Nat. Mater. , 306–310 (2012). S. Bi, H. Banda, M. Chen, L. Niu, M. Chen, T. Wu, J. Wang,R. Wang, J. Feng, T. Chen, M. Dinca, A. A. Kornyshev, and
IG. S1. Electrostatic energy of a point charge between twoelectrodes for a graphitic capacitor with 10 graphene planesper electrode with L x = 191 .
80 ˚A and L y = 196 .
88 ˚A. En-ergies were computed for a range of Thomas-Fermi lengths l TF ranging from 0.0 to 8.0 ˚A, represented both by differentsymbols and by different colors indicated in the colorbar. Thesimulation results (symbols) are compared to the analyticalcalculation (straight lines) derived from Kaiser et al. . Theenergies have been shifted to be equal to zero when the ion isat the center of the box.FIG. S2. Integral capacitances computed at applied volt-ages of 1 and 2 V for the system made of two (100) gold-likeelectrodes in contact with a NaCl aqueous solution. G. Feng, “Molecular understanding of charge storage and charg-ing dynamics in supercapacitors with mof electrodes and ionicliquid electrolytes,” Nat. Mater. , 552–558 (2020). D. T. Limmer, A. P. Willard, P. Madden, and D. Chandler,“Hydration of metal surfaces can be dynamically heterogeneousand hydrophobic,” Proc. Natl. Acad. Sci. U.S.A. , 4200–4205
FIG. S3. Atomic density profiles for the O, H, Na + and Cl − atoms near the electrode for the systems made of two (100)gold-like electrodes in contact with a NaCl aqueous solution atnull potential. The structure is similar for all the l TF values.FIG. S4. Distribution of the water molecules O-H bonds ori-entation with respect to the vector normal to the electrodesurfaces for an applied potential of 2 V for the whole range ofsimulated l TF indicated by the colorbar; the distribution for0 V and l TF = 0 . (2013). L. Scalfi, D. T. Limmer, A. Coretti, S. Bonella, P. A. Mad- + FIG. S5. Atomic density profiles for the H atoms for an ap-plied potential of 2 V for the whole range of simulated l TF indicated by the colorbar. The negative (positive) electrodeis located at negative (positive) z . den, M. Salanne, and B. Rotenberg, “Charge fluctuations frommolecular simulations in the constant-potential ensemble,” Phys.Chem. Chem. Phys. , 10480–10489 (2020). L. Scalfi, M. Salanne, and B. Rotenberg, “Molecular simulationof electrode-solution interfaces,” https://arxiv.org/abs/2008.11967 (2020), arXiv:2008.11967. A. A. Kornyshev and M. A. Vorotyntsev, “Analytic expressionfor the potential energy of atest charge bounded by solid stateplasma,” J. Phys. C: Solid State Phys. , L691–L694 (1978). A. A. Kornyshev, W. Schmickler, and M. A. Vorotyntsev, “Non-local electrostatic approach to the problem of a double layer at ametal-electrolyte interface,” Phys. Rev. B , 5244–5256 (1982). N. B. Luque and W. Schmickler, “The electric double layer ongraphite,” Electrochim. Acta , 82–85 (2012). A. A. Kornyshev, N. B. Luque, and W. Schmickler, “Differentialcapacitance of ionic liquid interface with graphite: the story oftwo double layers,” J. Solid State Electrochem. , 1345–1349(2014). J. Comtet, A. Nigu`es, V. Kaiser, B. Coasne, L. Bocquet, andA. Siria, “Nanoscale capillary freezing of ionic liquids confinedbetween metallic interfaces and the role of electronic screening,”Nat. Mater. , 634–639 (2017). V. Kaiser, J. Comtet, A. Nigu`es, A. Siria, B. Coasne, and L. Boc-quet, “Electrostatic interactions between ions near thomas-fermisubstrates and the surface energy of ionic crystal at imperfectmetals,” Faraday Discuss. , 129–158 (2017). L. Thomas, “The calculation of atomic fields,” Proc. CambridgePhil. Roy. Soc. , 542–548 (1927). E. Fermi, “Un metodo statistico per la determinazione di al-cuna priorieta dell’atome,” Rend. Accad. Naz. Lincei , 602–607(1927). S. K. Reed, O. J. Lanning, and P. A. Madden, “ElectrochemicalInterface Between an Ionic Liquid and a Model Metallic Elec-trode,” J. Chem. Phys. , 084704 (2007). E. Paek, A. J. Pak, and G. S. Hwang, “On the influence ofpolarization effects in predicting the interfacial structure and ca-pacitance of graphene-like electrodes in ionic liquids,” J. Chem.Phys. , 024701 (2015). R. F. Nalewajski, “Electrostatic effects in interactions betweenhard (soft) acids and bases,” J. Am. Chem. Soc. , 944–945(1984). W. J. Mortier, S. K. Ghosh, and S. Shankar, “Electronegativity-equalization method for the calculation of atomic charges inmolecules,” J. Am. Chem. Soc. , 4315–4320 (1986). A. K. Rappe and W. A. Goddard III, “Charge equilibration formolecular dynamics simulations,” J. Phys. Chem. , 3358–3363(1991). N. Onofrio, D. Guzman, and A. Strachan, “Atomic origin ofultrafast resistance switching in nanoscale electrometallizationcells,” Nat. Mater. , 440–446 (2015). M. Buraschi, S. Sansotta, and D. Zahn, “Polarization effects indynamic interfaces of platinum electrodes and ionic liquid phases:A molecular dynamics study,” J. Phys. Chem. C , 2002–2007(2020). D. M. York and W. Yang, “A chemical potential equalizationmethod for molecular simulations,” The Journal of ChemicalPhysics , 159–172 (1996), publisher: American Institute ofPhysics. H. Nakano and H. Sato, “A chemical potential equaliza-tion approach to constant potential polarizable electrodes forelectrochemical-cell simulations,” J. Chem. Phys. , 164123(2019). L. Pastewka, T. T. J¨arvi, L. Mayrhofer, and M. Moseler,“Charge-transfer model for carbonaceous electrodes in polar en-vironments,” Phys. Rev. B , 165418 (2011). N. D. Lang and W. Kohn, “Theory of metal surfaces: inducedsurface charge and image potential,” Phys. Rev. B , 3541–3550(1973). N. V. Smith, C. T. Chen, and M. Weinert, “Distance of theimage plane from metal surfaces,” Phys. Rev. B , 7565–7573(1989). H. Gerischer, “An interpretation of the double layer capacity ofgraphite electrodes in relation to the density of states at theFermi level,” J. Phys. Chem. , 4249–4251 (1985). A. J. Pak, E. Paek, and G. S. Hwang, “Relative contributions ofquantum and double layer capacitance toward the supercapaci-tor performance of carbon nanotubes in an ionic liquid,” Phys.Chem. Chem. Phys. , 19741–19747 (2013). M. Salanne, B. Rotenberg, K. Naoi, K. Kaneko, P.-L. Taberna,C. P. Grey, B. Dunn, and P. Simon, “Efficient Storage Mecha-nisms for Building Better Supercapacitors,” Nat. Energy , 16070(2016). J. Carrasco, A. Hodgson, and A. Michaelides, “A molecularperspective of water at metal interfaces,” Nat. Mater. , 667–674 (2012). M. V. Fedorov and A. A. Kornyshev, “Ionic liquids at electrifiedinterfaces,” Chem. Rev. , 2978—3036 (2014). M. A. Gebbie, M. Valtiner, X. Banquy, E. T. Fox, W. A. Hen-derson, and J. N. Israelachvili, “Ionic liquids behave as diluteelectrolyte solutions,” Proc. Natl. Acad. Sci. U.S.A. , 9674–9679 (2013). A. M. Smith, A. A. Lee, and S. Perkin, “The electrostatic screen-ing length in concentrated electrolytes increases with concentra-tion,” J. Phys. Chem. Lett. , 2157–2163 (2016). R. A. Nistor, J. G. Polihronov, and M. H. M¨user, “A general-ization of the charge equilibration method for nonmetallic mate-rials,” J. Chem. Phys. (2006). R. A. Nistor and M. H. M¨user, “Dielectric properties of solidsin the regular and split-charge equilibration formalisms,” Phys.Rev. B , 104303 (2009). A. Marin-Lafl`eche, M. Haefele, L. Scalfi, A. Coretti, T. Du-fils, G. Jeanmairet, S. Reed, A. Serva, R. Berthin, C. Bacon,S. Bonella, B. Rotenberg, P. A. Madden, and M. Salanne, “Met-alwalls: A classical molecular dynamics software dedicated tothe simulation of electrochemical systems,” ChemRxiv preprint under review (2020). H. J. C. Berendsen, J. R. Grigera, and T. P. Straatsma, “TheMissing Term in Effective Pair Potentials,” J. Phys. Chem. ,6269–6271 (1987). L. X. Dang, J. Am. Chem. Soc. , 6954 (1995). A. Berg, C. Peter, and K. Johnston, J. Chem. Theory Comput.13