Laura Scalfi

Mathieu Salanne

Benjamin Rotenberg
email: benjamin.rotenberg@sorbonne-universite.fr

Molecular Simulation of Electrode-Solution Interfaces

Keywords: numbers: Contents I, INTRODUCTION Interfacial structure C, Interfacial dynamics D, Towards electrochemistry

Many key industrial processes, from electricity production, conversion and storage to electrocatalysis or electrochemistry in general, rely on physical mechanisms occurring at the interface between a metallic electrode and an electrolyte solution, summarized by the concept of electric double layer, with the accumulation/depletion of electrons on the metal side and of ions on the liquid side. While electrostatic interactions play an essential role on the structure, thermodynamics, dynamics and reactivity of electrode-electrolyte interfaces, these properties also crucially depend on the nature of the ions and solvent, as well as that of the metal itself. Such interfaces pose many challenges for modeling, because they are a place where Quantum Chemistry meets Statistical Physics. In the present review, we explore the recent advances on the description and understanding of electrodeelectrolyte interfaces with classical molecular simulations, with a focus on planar interfaces and solvent-based liquids, from pure solvent to water-in-salt-electrolytes.

I. INTRODUCTION

Many key industrial processes, from electricity production, conversion and storage [1], to electrocatalysis or electrochemistry in general [2], rely on physical mechanisms occurring at the interface between a metallic solid (electrode), allowing the transport of electrons, and an electrolyte solution, in which electric currents may arise from the transport of ionic species or the orientation of polar molecules. The most interesting features of such electrode-solution interfaces emerge from the coupling between the charges accumulated on both sides, summarized by the concept of "Electric double layer" (EDL) [3]. From a classical perspective, in perfect metals the electronic charge is localized at the surface and the electric potential is uniform inside the solid, so that the interfacial properties are essentially governed by the ionic densities and electrostatic potential profiles near the interface, which result from the balance between energetic considerations (attraction of counterions to the surface charge and repulsion of co-ions, favoring the build-up of charge) and entropic ones (di↵usion leading to uniform concentrations). Of particular interest are of course the charge accumulated on the electrode and the capacitance of the interface, i.e. the derivative of the charge with respect to the potential drop across the interface, which can be measured in electrochemical experiments. Voltage also provides a handle on wetting properties, by changing surface free energies and the resulting contact angle.

While electrostatic interactions play an essential role on the structure, thermodynamics, dynamics and reactivity of electrode-electrolyte interfaces, these properties also crucially depend on the nature of the ions and solvent, as well as that of the metal itself. The finite size of the ions and of solvent molecules results in a layered structure near solid walls, and their ability to form hydrogen bonds (e.g. for water) constrains their orientation. Such features, which can be investigated using spectroscopic techniques, depend on the atomic lattice of the metal and a given fluid behaves di↵erently on di↵erent faces of the same crystal of a given metal. Finally, many real materials cannot be considered as perfect metals, and the charge and potential distribution within the electrode and their coupling with the liquid must also be taken into account. This e↵ect of the metallic character of the electrode on the properties of the interfacial liquid is a good illustration of the challenge that such interfaces pose for modeling, because by bridging electrons in a solid and ions in a solvent, they are a place where Quantum Chemistry meets Statistical Physics.

On the theoretical side, much progress has been made since the pioneering works of Gouy, Chapman and Stern. At the same continuous level of description, extensions of the mean-field Poisson-Boltzmann theory have been proposed to capture the e↵ects of electrostatic correlations and excluded volume or solvent polarization [4,[START_REF] Mceldrew | Theory of the Double Layer in Water-in-Salt Electrolytes[END_REF] on the structure and capacitance of the EDL, with a low computational cost compatible with routine use in engineering applications. Even the charging dynamics can be investigated at this level [START_REF] Bazant | Di↵use-charge dynamics in electrochemical systems[END_REF], even though the e↵ects of ionic correlations or of the coupling with the solvent dynamics are more accurately described by mesoscopic simulations with explicit or implicit ions [START_REF] Netz | Electrofriction and Dynamic Stern Layers at Planar Charged Surfaces[END_REF][START_REF] Grün | Relaxation of the electrical double layer after an electron transfer approached by Brownian dynamics simulation[END_REF][START_REF] Lobaskin | Di↵usive-convective transition in the non-equilibrium charging of an electric double layer[END_REF][START_REF] Asta | Lattice Boltzmann electrokinetics simulation of nanocapacitors[END_REF]. At the other extreme, quantum calculations, usually based on electronic Density Functional Theory (even though Quantum Monte Carlo can now provide even more accurate results e.g. on water-carbon interactions [START_REF] Striolo | The carbon-water interface: Modeling challenges and opportunities for the waterenergy nexus[END_REF]), allow to capture the density of states of the metal as well as a detailed description of a few interfacial molecules [START_REF] Taylor | First principles reaction modeling of the electrochemical interface: Consideration and calculation of a tunable surface potential from atomic and electronic structure[END_REF][START_REF] Lautar | Modeling interfacial electrochemistry: concepts and tools[END_REF]. However, their computational cost prevents a fully molecular description of the EDL, and resort to a simplified description of the solvent (polarizable continuum) is the rule rather than the exception. As a result, classical molecular simulations have emerged as a powerful compromise between an atomic description and a computational cost allowing a su cient sampling of relevant electrolyte configurations.

In the present review, we explore the recent advances on the description and understanding of electrode-electrolyte interfaces with classical molecular simulations. While many applications involve porous electrodes with disordered structures and complex electrolytes such as room temperature ionic liquids, we restrict ourselves to the simpler yet practically relevant and physically rich case of planar interfaces and solvent-based liquids, from pure solvent to waterin-salt-electrolytes (we refer the readers to e.g. Refs. [START_REF] Merlet | On the Molecular Origin of Supercapacitance in Nanoporous Carbon Electrodes[END_REF][START_REF] Simoncelli | Blue energy and desalination with nanoporous carbon electrodes: Capacitance from molecular simulations to continuous models[END_REF][START_REF] Merlet | The electric double layer has a life of its own[END_REF][START_REF] Fedorov | Ionic liquids at electrified interfaces[END_REF][START_REF] Burt | A review of molecular modelling of electric double layer capacitors[END_REF] on the porous and/or ionic liquid cases). In Section II, we discuss the description of electrode-solution interfaces, emphasizing the choice of the models to represent the metallic character of the electrode and its interactions with the electrolyte solution. Section III then presents the various strategies to simulate electrochemical systems, with a potential di↵erence between two electrodes, and fundamental issues related to such simulations. Finally, Section IV illustrates a selection of properties which can be investigated with such classical molecular simulations, with examples on the capacitance, the interfacial structure and dynamics, as well as steps towards electrochemistry.

The question of how to model metallic electrodes in classical molecular simulations is related to the more general question of the electrostatic response of a medium to an electric charge (ion, or partial charges from molecules), encountered not only in electrochemistry, but in all systems involving interfaces, including e.g. biological macromolecules. From a quantum mechanical perspective, electric conduction in a medium is related to the position of the Fermi level relative to the system's energy levels: In a metal, the conduction band is thermally accessible to the electrons, which can be delocalized over the whole material, whereas in an insulator a band gap hinders the conduction of electrons, which remain localized on atomic sites. From a classical continuum perspective, an essential feature of these interfaces is a contrast in the polarization response of the various media, quantified by their dielectric constant ✏ r , ranging from 1 for vacuum to ⇡ 80 for liquid water and 1 for a perfect metal. Such a dielectric contrast has strong consequences on the behavior of a charge distribution close to the interface. This is usually expressed in terms of image charges. For a sharp flat interface between two media 1 (polar solvent) and 2 (solid wall) with dielectric constants ✏ 1 and ✏ 2 , the electrostatic potential arising from a set of charges (ions) q ext = {q ext 1 , . . . , q ext N } embedded in medium 1 is identical, within this medium, to that arising from a fictitious system in which medium 2 is assigned a dielectric constant ✏ 1 and a set of image charges are placed symmetrically with respect to the boundary (see Figure 1a below), with magnitudes:

q im = ✏ 1 ✏ 2 ✏ 1 + ✏ 2 q ext (1)
For an insulating interface such as water-vacuum (✏ 1 ✏ 2 ) the image charges are similar to the source charges q im ⇡ q ext , whereas for a perfect metallic interface (✏ 2 ! 1) the images have opposite charges q im = q ext . It follows a radically di↵erent electrostatic interaction of a charge with its image: attractive for the metallic case and repulsive for insulators.

Electrodes in classical molecular simulations

To model these systems at the atomic scale, one should in principle perform quantum calculations taking into account the electronic density on both the electrode material and the electrolyte. However, this becomes computationally prohibitive for the simulation of large systems over long time scales and alternative approaches have been developed to capture the e↵ect of a metallic electrode on an electrolyte within classical simulations. The purpose of these models is not to provide an accurate description of the metal and its properties, but rather to reproduce the appropriate boundary conditions for the electrolyte. Within such a simplified description, it becomes possible to sample the configurations explored by the electrolyte at finite temperature and to investigate the properties of electrochemical interfaces. In this context, two features are particularly important: accounting for the polarization of the metal by the electrolyte, discussed in the following section, and the possibility to accumulate a net charge on the interface, e.g. in the presence of an applied voltage between two electrodes.

A first method to describe charged electrodes is to explicitly put a constant net charge on the electrode, using a surface charge on a wall [START_REF] Torrie | Electrical double layers. II. Monte Carlo and HNC studies of image e↵ects[END_REF][START_REF] Glosli | Molecular dynamics simulation of adsorption of ions from aqueous media onto charged electrodes[END_REF][START_REF] Kiyohara | Monte Carlo simulation of electrolytes in the constant voltage ensemble[END_REF][START_REF] Kiyohara | Monte Carlo Simulation of Porous Electrodes in the Constant Voltage Ensemble[END_REF] or discrete point charges [START_REF] Van Megen | The grand canonical ensemble Monte Carlo method applied to the electrical double layer[END_REF][START_REF] Crozier | Molecular-dynamics simulations of ion size e↵ects on the fluid structure of aqueous electrolyte systems between charged model electrodes[END_REF]. The electroneutrality of the system can be balanced by excess ions in the electrolyte or by an opposite charge on a second electrode. This setup would correspond to a charged pore or an isolated (open circuit) charged capacitor. From the electrostatic point of view, two oppositely and homogeneously charged walls induce a uniform electric field between them, so that it is (at least in principle) equivalent to directly apply an external electric field on the liquid confined in the electrochemical cell [START_REF] Lee | Molecular dynamics study of a dipolar fluid between charged plates[END_REF][START_REF] Hautman | Molecular dynamics simulation of water beween two ideal classical metal walls[END_REF][START_REF] Rose | Adsorption of Na + and Cl at the charged water-platinum interface[END_REF][START_REF] Smith | Simulation study of the ferrous ferric electron transfer at a metal-aqueous electrolyte interface[END_REF][START_REF] Daub | Electrowetting at the Nanoscale[END_REF]. In a real system, however, the system is rather connected to a voltage generator, which maintains a constant electric potential di↵erence between the electrodes and allows the exchange of charge between them [START_REF] Kiyohara | Monte Carlo simulation of electrolytes in the constant voltage ensemble[END_REF][START_REF] Kiyohara | Monte Carlo Simulation of Porous Electrodes in the Constant Voltage Ensemble[END_REF][START_REF] Siepmann | Influence of Surface-Topology and Electrostatic Potential on Water Electrode Systems[END_REF][START_REF] Reed | Electrochemical Interface Between an Ionic Liquid and a Model Metallic Electrode[END_REF][START_REF] Pounds | Ion adsorption at a metallic electrode: an ab initio based simulation study[END_REF][START_REF] Vatamanu | Molecular dynamics simulations of atomically flat and nanoporous electrodes with a molten salt electrolyte[END_REF][START_REF] Petersen | A computationally e cient treatment of polarizable electrochemical cells held at a constant potential[END_REF]. We will discuss the corresponding simulation setups in Section III, and now turn to the representation of polarization of the metal in classical molecular simulations.

Representing the electronic response of the metal

Beyond their net charge, a characteristic feature of metallic electrodes is their polarization by the electrolyte (this is even the main e↵ect for neutral electrodes). In the context of molecular simulations, one is primarily concerned with the e↵ect of the metal on the electrolyte and a variety of methods have been proposed. The typical slab geometry of capacitors allows a number of simplifications, including analytical expressions of the electrostatic forces acting on the electrolyte charges and of the electrostatic energy. Even though a uniform charge distribution with varying magnitude can be used to account for an applied external field or voltage, such a description lacks the lateral charge heterogeneities induced by the discrete nature of ions and molecules of the electrolyte. From the electrostatic point of view, an e cient strategy is to impose the proper electrostatic boundary conditions at the surface of the metal using the above-mentioned concept of image charges, which can be either explicit or accounted for implicitly in modified Green functions [START_REF] Torrie | Electrical double layers. II. Monte Carlo and HNC studies of image e↵ects[END_REF][START_REF] Smith | Simulation study of the ferrous ferric electron transfer at a metal-aqueous electrolyte interface[END_REF][START_REF] Petersen | A computationally e cient treatment of polarizable electrochemical cells held at a constant potential[END_REF][START_REF] Parsonage | Computer simulation of water between metal walls[END_REF][START_REF] Gardner | Water-like particles at surfaces. II. In a double layer and at a metallic surface[END_REF][START_REF] Takae | Fluctuations of local electric field and dipole moments in water between metal walls[END_REF][START_REF] Girotto | Simulations of ionic liquids confined by metal electrodes using periodic Green functions[END_REF]. The use of periodic boundary conditions (see section III) also requires special care to compute electrostatic interactions, and e cient algorithms have been developed to deal with the image charges in simulations [START_REF] Arnold | E cient algorithms for electrostatic interactions including dielectric contrasts[END_REF]. Alternatively, the electrostatic problem can also be solved numerically without resorting to image charges: The Induced Charge Computation (ICC) method treats the charge density of the solid as a dynamical variable discretized on a grid and solves the Poisson equation to obtain the induced (surface) charge, as illustrated in Figure 1b. The dielectric medium is then characterized by a space-dependent dielectric constant ✏(r) which can in principle describe arbitrarily shaped interfaces and non-homogeneous media. Based on the variational procedure of Allen et al. for solid-electrolyte interfaces [START_REF] Allen | Electrostatic potential inside ionic solutions confined by dielectrics: a variational approach[END_REF], extensions using a matrix formulation [START_REF] Boda | Computing induced charges in inhomogeneous dielectric media: Application in a Monte Carlo simulation of complex ionic systems[END_REF] or an iterative algorithm (ICC ⇤ ) [START_REF] Tyagi | An iterative, fast, linear-scaling method for computing induced charges on arbitrary dielectric boundaries[END_REF][START_REF] Breitsprecher | Electrode models for ionic liquid-based capacitors[END_REF] were proposed. Another example of this numerical approach to induced charges is the Generalized Minimal Residue (GMRES) method, which provides good performance when used in conjuction with fast Ewald solver [START_REF] Barros | E cient and accurate simulation of dynamic dielectric objects[END_REF].

N

anostructures of precious metals find many applications, for example, as catalysts, electrode materials, biomarkers, and therapeutics, including gold (Au) nanorods for photothermal cancer therapy and nanoparticles for gene delivery 1,2 . Metal nano objects have been synthesized in many shapes and sizes, however, control over nucleation, growth, and ligand interactions for nanoscale assembly remains a challenge 3-6 . Understanding the selective synthesis and structure-property relationships requires tremendous efforts by imaging, spectroscopy, and other laboratory techniques 7,8 .

Understanding and discovery can be accelerated by use of atomistic simulations up to the large nanometer scale (e.g., 100 nm) in comparison with experiment 9-15 . The Interface force field (IFF), for example, contains Lennard-Jones parameters for facecentered cubic (fcc) metals to simulate bulk solids, aqueous interfaces, and multiphase materials with polymers and biomacromolecules 9,10 . The parameters reproduce the density, surface tension, and anisotropy of surface energies of (h k l) facets, as well as the mechanical properties in excellent agreement with experiments, even better than some DFT methods 16 . Simulations using this non-polarizable model have proven helpful in understanding the adsorption mechanisms of biomolecules, as well as growth mechanisms and shape preferences of metal nanostructures using particular ligands 13,17-23 . Simulations achieved quantitative agreement with experimental observations 5,13,19,22- 25 yet mainly focused on ligands of low polarity, simple shapes, and simple surface assemblies without accounting for the effects of induced charges and external potentials. It is a shortcoming that the non-polarizable potential does not account for the contribution of the induced charges to interfacial processes during Molecular Dynamics (MD) or Monte Carlo (MC) simulations. Polarization on the metal surface was shown to affect the surface adsorption of molecules and could only be "added" a posteriori 26 . Such a posteriori calculations are significantly less accurate and impractical as they are uncoupled from the dynamics, require time-consuming post-processing of simulation outputs (coordinates and energies), and cannot be applied to corrugated metal surfaces or under external electric fields. The effects of induced charges in the metal are known to be substantial in vacuum 27-29 , under external potentials in electrodes 30,31 , in the presence of ionic liquids 21 , and at high ionic strength in solution 26 , although expected to be weaker in dilute aqueous solution 26 .

Attempts to include polarization have also been made in alternative models for metallic nanostructures and electrodes. The GolP force field adds permanent dipoles to every atom to account for effects of induced charges 11,12,32,33 . The dipoles are implemented as fixed rods and shift the image plane for positively charged vs. negatively charged species on the metal surface. Another limitation is that surface energies and mechanical properties of the metal have not been reproduced and the compatibility with biomolecular force fields requires many adjustable parameters 12 . Several further models have also been developed to describe metallic electrodes at constant potential 34-39 . Siepmann and Sprik 34 pioneered models under a constant applied potential in which variable charges are added to the electrode and their magnitude is adjusted on-the-fly according to a variational procedure. The model matches image potentials and accounts for the polarization of the electrode by the electrolyte. However, parameters for the charge distribution are then necessary and molecular dynamics simulations require on-the-fly adjustments of the charge distribution. Moreover, surface energies, as well as interfacial energies of the metals have not been validated relative to experimental data and the original energy expressions are difficult to use due to their complexity (see Supplementary Note 1 for details).

Here we introduce a simple polarizable Lennard-Jones model for metallic gold. It adds the correct amount of attractive polarization to the neat Lennard-Jones potential and retains the mobility of all atoms, as well as the other advantageous aspects of the nonpolarizable model (Fig. 1a) 10 . The model reproduces the classical image potential, lattice parameters, surface energy, and hydration energy with water in excellent agreement with experiments, and a good correlation with results from density functional theory is demonstrated. The model is compatible with

q = 0 q i q j r ij r ij r ij q = ± ±1.
E i ͑1/2͒ = E i Ј Ϯ 2ini/⑀1,

͑4͒

it is then possible to obtain a relation between i and the field E i Ј generated by all real and induced charges but the ith one. This relation, obtained by combining Eqs. ͑2͒-͑4͒, reads

i = ⑀1 2 ͩ ⑀1 -⑀2 ⑀1 + ⑀2 ͪEi Ј • ni,

͑5͒

and represents in fact the discretized form of an integral equation for i, since E i Ј can be expressed as linear combi- nation of all ͑real and induced͒ charges, but the induced one located at xi. For convenience we will denote the permittivity-dependent factor in Eq. ͑5͒ as

f = ⑀1 2 ͩ ⑀1 -⑀2 ⑀1 + ⑀2 ͪ.
͑6͒ Equation ͑5͒ is the common starting point for induced charge-based methods. 19,25,39 The strategy which is chosen to solve this equation strongly influences not only the efficiency of the methods, but also their scope, in particular, regarding the type of periodicity which can be treated.

III. THE ALGORITHM

Before introducing the ICC Ã algorithm, we will briefly review the basics of a recent similar approach, namely, ICC, 25,26 which we will use as a reference to compare the characteristics and the performances of ICC Ã . The strategy adopted in the ICC algorithm is to split the electric field E i Ј into a contribution from the real charges and one from the

induced charges E i Ј= E i real + E i ind .
Due to the linearity of Maxwell's equations, the electric field generated by the induced charges can be written as E i ind = ͚jKijj, so that Eq. ͑5͒ can always be put in the form of a system of linear equations

͚ j ͑␦ij -fKij͒j = fE i real • ni, ͑7͒
where the kernel Kij encodes the geometrical information needed to compute the contribution to the field at position xi from the induced charge located at xj. The polarization charges in ICC are then determined by solving the linear system

Au = b, ͑8͒
where Aij= ͑␦ij-fKij͒, the components of vector u are the values of the induced charges i, and the components of b are bi = fE i real • ni. This method is appealing because the matrix inversion has to be performed only at the beginning of the simulation, but it is lacks, at least in the original description of the algorithm, 25,26 treatment of long-range contributions. ͑We will turn later on the performances of ICC Ã with respect to ICC.͒

The ICC Ã algorithm basically exploits a different strategy to solve Eq. ͑5͒. No matter how the electric field is calculated ͑for example, using nearest neighbor interaction or taking into account every contribution form periodic im-ages͒, it can always be expressed as a linear function of the charges in the system, and Eq. ͑5͒ is precisely a linear system of equations that is already written in the form which is typical to iterative methods. Our aim is to employ an iterative scheme in the context of periodic systems. Among the many different algorithms 40 that can be used to implement an iterative solution for Eq. ͑5͒, the successive over relaxations ͑SORs͒ approach is a very versatile and efficient one. In a SOR scheme, the general equation

u = Bu + c, ͑9͒
is approximated with desired precision by the sequence of estimates

u ͑n+1͒ = ͑Bu ͑n͒ + c͒ + ͑1 -͒u ͑n͒ . ͑10͒
Here the superscript n denotes the iteration number and the relaxation parameter, which usually can be in the range ͑0,2͒ in order to guarantee convergence. 40 The iterative solution for the induced charges reads, explicitly,

i ͑n+1͒ = ͑fE i Ј • ni͒ + ͑1 -͒ i ͑n͒ . ͑11͒
The initial values of the induced charges can be extracted form a random distribution. At every iteration the value of the electric field is computed by means of the electrostatics algorithm of choice, and it is used to perform the next SOR approximation of the charge distribution, according to Eq. ͑11͒. One usually keeps iterating until a desired accuracy has been achieved. This can be tracked by computing

⌬ = max i=1,M ͯ i n+1 -i n i n ͯ.

͑12͒

When ⌬ is smaller than some preset tolerance value, the iterative process is stopped. This scheme is highly useful in molecular dynamics ͑MD͒ simulations, as it will converge quite fast after the first integration step. Suppose, for example, that for a given configuration of the real charges, the induced ones have been determined up to the desired accuracy. After one MD integration step the induced charge distribution has to be recalculated. However, because the mobile charges normally have moved only slightly from their previous positions, the old induced charge values already serve as good approximation. Therefore, when using these values as a starting guess for the iterative procedure at the next MD step one will obtain a sufficiently well converged solution in, at most, a handful of iterations. In practical ap- Tyagi JCP 2010 Figure 1 X X X X

N

anostructures of precious metals find many applications, for example, as catalysts, electrode materials, biomarkers, and therapeutics, including gold (Au) nanorods for photothermal cancer therapy and nanoparticles for gene delivery 1,2 . Metal nano objects have been synthesized in many shapes and sizes, however, control over nucleation, growth, and ligand interactions for nanoscale assembly remains a challenge 3-6 . Understanding the selective synthesis and structure-property relationships requires tremendous efforts by imaging, spectroscopy, and other laboratory techniques 7,8 .

Understanding and discovery can be accelerated by use of atomistic simulations up to the large nanometer scale (e.g., 100 nm) in comparison with experiment 9-15 . The Interface force field (IFF), for example, contains Lennard-Jones parameters for facecentered cubic (fcc) metals to simulate bulk solids, aqueous interfaces, and multiphase materials with polymers and biomacromolecules 9,10 . The parameters reproduce the density, surface tension, and anisotropy of surface energies of (h k l) facets, as well as the mechanical properties in excellent agreement with experiments, even better than some DFT methods 16 . Simulations using this non-polarizable model have proven helpful in understanding the adsorption mechanisms of biomolecules, as well as growth mechanisms and shape preferences of metal nanostructures using particular ligands 13,17-23 . Simulations achieved quantitative agreement with experimental observations 5,13,19,22- 25 yet mainly focused on ligands of low polarity, simple shapes, and simple surface assemblies without accounting for the effects of induced charges and external potentials. It is a shortcoming that the non-polarizable potential does not account for the contribution of the induced charges to interfacial processes during Molecular Dynamics (MD) or Monte Carlo (MC) simulations. Polarization on the metal surface was shown to affect the surface adsorption of molecules and could only be "added" a posteriori 26 . Such a posteriori calculations are significantly less accurate and impractical as they are uncoupled from the dynamics, require time-consuming post-processing of simulation outputs (coordinates and energies), and cannot be applied to corrugated metal surfaces or under external electric fields. The effects of induced charges in the metal are known to be substantial in vacuum 27-29 , under external potentials in electrodes 30,31 , in the presence of ionic liquids 21 , and at high ionic strength in solution 26 , although expected to be weaker in dilute aqueous solution 26 .

Attempts to include polarization have also been made in alternative models for metallic nanostructures and electrodes. The GolP force field adds permanent dipoles to every atom to account for effects of induced charges 11,12,32,33 . The dipoles are implemented as fixed rods and shift the image plane for positively charged vs. negatively charged species on the metal surface. Another limitation is that surface energies and mechanical properties of the metal have not been reproduced and the compatibility with biomolecular force fields requires many adjustable parameters 12 . Several further models have also been developed to describe metallic electrodes at constant potential 34-39 . Siepmann and Sprik 34 pioneered models under a constant applied potential in which variable charges are added to the electrode and their magnitude is adjusted on-the-fly according to a variational procedure. The model matches image potentials and accounts for the polarization of the electrode by the electrolyte. However, parameters for the charge distribution are then necessary and molecular dynamics simulations require on-the-fly adjustments of the charge distribution. Moreover, surface energies, as well as interfacial energies of the metals have not been validated relative to experimental data and the original energy expressions are difficult to use due to their complexity (see Supplementary Note 1 for details).

Here we introduce a simple polarizable Lennard-Jones model for metallic gold. It adds the correct amount of attractive polarization to the neat Lennard-Jones potential and retains the mobility of all atoms, as well as the other advantageous aspects of the nonpolarizable model (Fig. 1a) 10 . The model reproduces the classical image potential, lattice parameters, surface energy, and hydration energy with water in excellent agreement with experiments, and a good correlation with results from density functional theory is demonstrated. The model is compatible with In order to include the molecular details of the interface (both in the shape and its atomistic nature), one can turn to descriptions based on the electrode atoms, i.e. treating the metal at the same level of description as the electrolyte, albeit with dedicated force fields. A first class of such descriptions consists in allowing mobile charges to rearrange in response to the configuration of the electrolyte: This includes core-shell models such as the Drude oscillator [START_REF] Geada | Insight into induced charges at metal surfaces and biointerfaces using a polarizable Lennard-Jones potential[END_REF], with a charge tethered to the electrode atom via a spring (see Figure 1c), or the rod model [START_REF] Iori | Including image charge e↵ects in the molecular dynamics simulations of molecules on metal surfaces[END_REF][START_REF] Iori | GolP: An atomistic force-field to describe the interaction of proteins with Au(111) surfaces in water[END_REF], in which the auxiliary charge can rotate at a fixed distance around the atom. Such models have the advantage of being easily implemented in standard molecular simulation codes. To ensure a correct adiabatic separation of the charge dynamics and avoid instabilities and/or energy transfers, the mass of the auxiliary charge and the spring strength or rod length should be carefully chosen. A second class of models, described in the following section, considers instead fixed but fluctuating charges, illustrated in Figure 1d. One advantage of all these atomistic descriptions compared to those based on image charges is that they are not restricted to slab geometries and can be used to deal with disordered porous electrodes. They capture the polarization e↵ects by rearranging the charge distribution of the metal due to the electric potential created by the electrolyte at each step. While this is not completely in the scope of the present review, we finally briefly mention attempts at including simplified quantum mechanical treatments of the interface in classical simulations, e.g. based on the Jellium model, to represent the spilling of the electronic charge out of the electrode [START_REF] Schmickler | The interphase between jellium and a hard sphere electrolyte. A model for the electric double layer[END_REF][START_REF] Shelley | Modeling and structure of mercury-water interfaces[END_REF], or the "direct dynamics" [START_REF] Price | Molecular dynamics, density functional theory of the metal-electrolyte interface[END_REF][START_REF] Walbran | Model for the electrostatic response of the copper -water interface[END_REF]. Recent developments of mixed quantum/classical (QM/MM) simulations or tight binding approaches might also provide interesting alternatives in the near future.

q = 0 qiqj rij rij rij q = ± ±1
+ + = kr,ij (rij -r0,ij)2 Model parameters ij bonded ∑ ij nonbonded (1,2 excl) ∑ ij nonbonded (1, 2

Focus on fluctuating charge methods

We now present in some more detail the family of models describing the polarization by assigning charge distributions (typically point charges or Gaussians) to the electrode atoms, with magnitude q = {q 1 , . . . , q M } treated as additional degrees of freedom which fluctuate in response to the dynamics of ions and molecules of the electrolyte. Early models to account for the polarization of the metal by external charges in fact included both fluctuating charges and induced dipoles at the atomic sites as additional degrees of freedom [START_REF] Finnis | The interaction of a point charge with an aluminium (111) surface[END_REF][START_REF] Finnis | The interaction of a point charge with a metal surface: theory and calculations for (111), (100) and (110) aluminium surfaces[END_REF], but this idea doesn't seem to have been explored much further. The energy of the system is then expressed as a function of the positions r N = {r 1 , . . . , r N } and momenta of all electrolyte atoms (see also section II B), and the charges q = {q 1 , . . . , q M } are determined at each time step either by following an equation of motion, or by responding instantaneously to the motion of the electrolyte so as to impose the electric or electrochemical potential of the electrode atoms (see Section III).

This general idea has resulted in a variety of models, which to some extent share the same quadratic form of the energy as a function of the electrode charges q, even though the physical meaning of the parameters may di↵er, and of numerical algorithms to determine the charges at each time step. While fluctuating charge models had already been used for molecules, with the charge equilibration (QE) method or the electronegativity equalization method (EEM) [START_REF] Nalewajski | Electrostatic e↵ects in interactions between hard (soft) acids and bases[END_REF][START_REF] Mortier | Electronegativity-equalization method for the calculation of atomic charges in molecules[END_REF][START_REF] Rappe | Charge equilibration for molecular dynamics simulations[END_REF][START_REF] York | A chemical potential equalization method for molecular simulations[END_REF], their use to represent a metallic surface was first described by Siepmann and Sprik [START_REF] Siepmann | Influence of Surface-Topology and Electrostatic Potential on Water Electrode Systems[END_REF], using Gaussian charge distributions on the electrode atoms. This model allowed them to simulate a water film near the tip of a scanning tunneling microscope and was then adapted within a Born-Oppenheimer framework to simulate an electrochemical cell by Reed et al. [START_REF] Reed | Electrochemical Interface Between an Ionic Liquid and a Model Metallic Electrode[END_REF]. The electrostatic energy of the system reads

U el (r N , q) = q T Aq 2 q T B(r N ) + C(r N ) (2)
where the symmetric M ⇥ M matrix A depends on the positions of the electrode atoms and the parameters describing the charge distribution on each atom, while the components of the vector B are the electrostatic potentials due to the electrolyte on each electrode atom (see Refs [START_REF] Reed | Electrochemical Interface Between an Ionic Liquid and a Model Metallic Electrode[END_REF][START_REF] Gingrich | On the ewald summation of gaussian charges for the simulation of metallic surfaces[END_REF] for explicit expressions of A and B in the particular case of Gaussian charge distributions with 2D Ewald summation), and the scalar C corresponds to electrostatic interactions within the electrolyte. The set of electrostatic potentials on each electrode atom is given by the gradient of U el with respect to q T , @U el (r N , q) @q T = Aq B ,

and depends on the positions of the electrode atoms and on the electrolyte configuration. As a result, the set of charges satisfying the constraint of fixed electrostatic potentials = { 1 , . . . , M } for each atom, typically the same value for all atoms belonging to a given electrode and a di↵erence between the values for both electrodes, is given by q = A 1 (B + ). In practice, if the matrix A can be inverted numerically (this has to be done only once if the electrode atoms do not move), the charges are computed at each step by a simple matrix-vector multiplication. Other methods are possible, such as finding the charges by minimizing U el q T numerically, e.g. with the conjugate gradient method, or treating the fixed potential as a holonomic constraint [START_REF] Coretti | Mass-zero constrained molecular dynamics for electrode charges in simulations of electrochemical systems[END_REF].

The expression of the energy is similar in the case of the QE/EEM method, which was used to investigate electrochemical interfaces [START_REF] Streitz | Electrostatic potentials for metal-oxide surfaces and interfaces[END_REF][START_REF] Onofrio | Voltage equilibration for reactive atomistic simulations of electrochemical processes[END_REF][START_REF] Onofrio | Atomic origin of ultrafast resistance switching in nanoscale electrometallization cells[END_REF][START_REF] Liang | Applied Potentials in Variable-Charge Reactive Force Fields for Electrochemical Systems[END_REF][START_REF] Nakano | A chemical potential equalization approach to constant potential polarizable electrodes for electrochemical-cell simulations[END_REF][START_REF] Buraschi | Polarization e↵ects in dynamic interfaces of platinum electrodes and ionic liquid phases: A molecular dynamics study[END_REF], but the expressions of the matrix A and vector B involve terms such as electronegativities and chemical hardnesses, instead of the purely electrostatic picture described above. It is also related to the split charge equilibration approach, which includes bond-specific terms in the energy [START_REF] Nistor | Dielectric properties of solids in the regular and split-charge equilibration formalisms[END_REF]. Finally, Pastewka et al. pushed the concept of using fluctuating charges as a proxy for a quantum description even further, by also parametrizing the band-structure energy, which appears in the tight-binding approximation with self-consistent charges to describe the terms beyond the electrostatic interactions between them, as a function of q [START_REF] Pastewka | Charge-transfer model for carbonaceous electrodes in polar environments[END_REF]. This allowed a purely classical description of various carbon electrodes, taking into account to some extent their di↵erent band structures. The potential of such a promising strategy doesn't seem to have been much exploited so far.

B. It takes two to tango

Description of the electrolyte

On the other side of the interface, the electrolyte can be represented using various levels of sophistication. In molecular dynamics simulations, two main families of models are generally employed: all-atom models, with intramolecular interactions or rigid bonds, or coarse-grained models, o↵ering a compromise between computational cost (due to the reduced number of interacting sites and the use of larger timesteps) and accuracy [START_REF] Salanne | Simulations of Room Temperature Ionic Liquids: from Polarizable to Coarse-Grained Force Fields[END_REF]. In the case of electrochemical interfaces, coarse-grained models generally include between 4 and 10 atoms per grains [START_REF] Merlet | Influence of solvation on the structural and capacitive properties of electrical double layer capacitors[END_REF], resulting in a decrease of the number of interaction sites by one order of magnitude compared to all-atom models (note that an alternative consists in using "united-atom" models in which only the hydrogen atoms are merged with heteroatoms to which they are bonded).

In both families of models, the intermolecular interactions will then determine most of the liquid properties. Many analytical expressions are used in order to account for the short-range repulsion between the electronic clouds, the dispersion e↵ects and the electrostatic interactions. Concerning electrostatics, all the models employ partial charges which are distributed among the interaction sites of the molecules, which interact between them and with the electrode charges. In a few cases, polarization e↵ects are included by adding either charge-on-a-spring (Drude oscillator model) or induced dipoles on the atomic sites [START_REF] Tazi | Potential-induced ordering transition of the adsorbed layer at the ionic liquid / electrified metal interface[END_REF][START_REF] Bedrov | Molecular dynamics simulations of ionic liquids and electrolytes using polarizable force fields[END_REF][START_REF] Park | Interference of electrical double layers: Confinement e↵ects on structure, dynamics, and screening of ionic liquids[END_REF]. These simulations are more costly from the computational point of view since they involve additional degrees of freedom that have to be either propagated or determined by solving a set of self-consistent equations (they have a many-body character), but they should in principle be more accurate since the electrostatics interactions are of primary importance at the interface. An alternative is to use rescaled charges in order to mimic these polarization e↵ects [START_REF] Breton | Molecular modeling of aqueous electrolytes at interfaces: E↵ects of long-range dispersion forces and of ionic charge rescaling[END_REF].

Electrode-electrolyte: non-electrostatic interactions

The electrostatic interactions between the electrode and the electrolyte ions and molecules follow from their respective descriptions as charge and dipole distributions, discussed above. Another crucial aspect is to account for the repulsive and dispersive intermolecular interactions between the electrolyte atoms and the electrode surface. In the earliest studies, the surface was materialized by a structure-less wall, that took the form of a hard wall, a purely repulsive potential or di↵erent flavors of one-dimensional Lennard-Jones potentials such as the Steele potential [START_REF] Kiyohara | Monte Carlo simulation of electrolytes in the constant voltage ensemble[END_REF][START_REF] Lee | Molecular dynamics study of a dipolar fluid between charged plates[END_REF][START_REF] Hautman | Molecular dynamics simulation of water beween two ideal classical metal walls[END_REF][START_REF] Parsonage | Computer simulation of water between metal walls[END_REF]. These descriptions are computationally less demanding but lack important molecular features close to the interface, that give rise e.g. to specific adsorption sites or templating e↵ects. Atomic descriptions of the electrodes are therefore employed, using a finite number of discrete atoms, in order to account for di↵erent geometries and crystal structures [START_REF] Daub | Electrowetting at the Nanoscale[END_REF][START_REF] Siepmann | Influence of Surface-Topology and Electrostatic Potential on Water Electrode Systems[END_REF][START_REF] Geada | Insight into induced charges at metal surfaces and biointerfaces using a polarizable Lennard-Jones potential[END_REF][START_REF] Iori | GolP: An atomistic force-field to describe the interaction of proteins with Au(111) surfaces in water[END_REF]. Usual intermolecular potentials are employed such as Lennard-Jones and Born-Mayer potentials. In most of these studies, the non-electrostatic electrode-electrode interactions are neglected, because the individual electrode atoms are immobile, with a few exceptions where the electrode structure is rigid but can translate as to mimic a piston at constant pressure, as done in Ref. [START_REF] Coretti | Mass-zero constrained molecular dynamics for electrode charges in simulations of electrochemical systems[END_REF]; in only a few cases, the electrodes atoms are free to move and interact with a harmonic potential [START_REF] Spohr | Molecular dynamics simulation of a water/metal interface[END_REF]. It should be highlighted that the inclusion of explicit atoms for the intermolecular part of the interactions does not necessarily imply an atomistic treatment of the electrostatic part and vice versa. Mixing atomistic and structure-less descriptions however poses the problem of where to locate the interface (in a continuum picture) with respect to the atomic positions. Intermediate approaches have also been used, which formulate the interaction potential as a corrugated potential, avoiding the calculation of pair terms but reproducing the local roughness of the substrate [START_REF] Spohr | Computer simulation of the water/platinum interface[END_REF].

As usual, the parametrization of these non-electrostatic interactions, which may lead or not to good predictions depending on how the electrostatic part is described, is crucial. While some force fields were fitted on ab-initio calculations or to reproduce experimental data, in the absence of reliable data the choices are often based on the availability of parameters from other studies. Accurate quantum Monte-Carlo reference data were reported for interfaces between water and graphene or carbon nanotubes [START_REF] Alfè | How strongly do hydrogen and water molecules stick to carbon nanomaterials?[END_REF], that could be used to develop new classical force fields for carbon electrodes, for which no specific optimization has been proposed to date. Comparatively, more e↵orts were put in the representation of metal surfaces. For example, Heinz et al. have proposed a systematic parameterization of Lennard-Jones potentials for several face-centered cubic metals, that were shown to reproduce a few experimental data for interfaces with water, such as surface tensions [START_REF] Heinz | Accurate simulation of surfaces and interfaces of face-centered cubic metals using 12-6 and 9-6 lennard-jones potentials[END_REF]. When used in combination with a core-shell model to represent the polarization of the metal, it was necessary to reparametrize the Lennard-Jones potentials as well [START_REF] Geada | Insight into induced charges at metal surfaces and biointerfaces using a polarizable Lennard-Jones potential[END_REF]. In the case of platinum, Siepmann and Sprik showed that using a three-body function was necessary to push water molecules on top of metal atoms to represent chemisorption e↵ects [START_REF] Siepmann | Influence of Surface-Topology and Electrostatic Potential on Water Electrode Systems[END_REF]. An alternative approach was recently proposed to include these e↵ects through an attractive two-body Gaussian potential [START_REF] Clabaut | Ten facets, one force field: The GAL19 force field for water-noble metal interfaces[END_REF]. We will also come back to reactive force fields to capture the breaking and formation of bonds in Section IV D.

III. HOW TO SIMULATE POLARIZED ELECTRODE-SOLUTION INTERFACES

Beyond the choice of the force field to describe the metallic electrode, the electrolyte and the interactions between them, a molecular simulation requires the definition of the simulated system and an algorithm to sample its configurations in order to compute physical properties. In this section, we explore more specifically the available options to simulate an electrochemical cell with an applied potential di↵erence between two electrodes, as well as issues related to the sampling of the corresponding microscopic configurations. Even though the discussion below applies more generally, we illustrate the various points with a description of the electrode based on fluctuating charges. In addition, we will not discuss in detail the possibilities o↵ered by existing open source softwares for the classical molecular simulation of electrochemical interfaces with the methods described in the present review, but refer the reader to the corresponding descriptions of e.g. Metalwalls [START_REF] Marin-Laflèche | MetalWalls: A Classical Molecular Dynamics Software Dedicated to the Simulation of Electrochemical Systems[END_REF], which is dedicated to the simulation of such interfaces, or more generic simulation packages allowing such simulations (possibly with an open source modification not provided in the standard distribution) such as LAMMPS, OpenMM or ESPRESSO (see [START_REF]Some open source classical molecular dynamics simulation codes allowing the simulation of electrodes[END_REF] for links).

A. Increase tension (between electrodes) or handle (the electrolyte) with care

Simulation setup

In a typical electrochemical system, two electrodes are separated by a slab of liquid electrolyte with all dimensions along the surfaces and between them much larger than molecular ones. In order to ensure that the two interfaces do not interfere (even though the total charge accumulated on both electrodes are of course correlated, see also Section III B 1), the distance between the electrodes must be large compared to the screening length. In experiments, one also introduces a separator to prevent contact (hence short-circuit) between the electrodes. The typical dimensions of systems that can be simulated with classical molecular simulations is in the range of a few to a few tens of nanometers. This requires a simplified description of the real device and limits the range of physical systems that can be simulated in a meaningful way. In particular, one always neglects the presence of the separator, which is a safe assumption provided that the distance between the electrodes remains large compared to the screening length. This in turn sets a lower bound on the ionic concentrations that can be considered -typically ⇠0.1 mol L 1 , a limit which also emerges from the constraint of having enough ions in the simulation box to ensure a good sampling of the phase space. . (b) PBC are applied in all directions, but the simulation box also includes vacuum (here the length of the box in this direction is three times larger than the electrochemical cell) and a "slab" correction is applied when computing electrostatic interactions; as in the previous case, the two electrodes are maintained at di↵erent potentials. (c) PBC are applied in all directions, without vacuum; all electrode atoms are maintained at the same potential and the capacitor is charged by applying an electric field to the electrolyte such that = ELz, with Lz the box dimension in the direction perpendicular to the electrode.

Last but not least, the small number of atoms in the system (10 3 -10 5 ) compared to a real device requires, as always in molecular simulation of condensed matter, the use of periodic boundary conditions (PBC). Based on the above discussion, it is natural to consider PBC in the two directions along the surfaces, but not in the direction perpendicular to the electrodes, as illustrated in Figure 2a. This requires in particular computing long-range electrostatic interactions with dedicated approaches, such as 2D Ewald summation, also taking into account the description of electrode atoms by Gaussian charge distributions instead of point charges if needed [START_REF] Reed | Electrochemical Interface Between an Ionic Liquid and a Model Metallic Electrode[END_REF][START_REF] Gingrich | On the ewald summation of gaussian charges for the simulation of metallic surfaces[END_REF]. An alternative approach, more easily implemented in standard packages, is to mimick 2D PBC by considering PBC in all directions, but adding vacuum to the simulation box, as illustrated in Figure 2b. In that case, it is necessary to apply a "slab" correction when computing electrostatic interactions [START_REF] Yeh | Ewald summation for systems with slab geometry[END_REF].

Constant-voltage or finite field

The simulation of an electrochemical cell also requires a way to impose voltage between the two electrodes. The most straightforward way to achieve this, with descriptions allowing to set the potential of electrode atoms such as methods based on fluctuating charges, is to impose di↵erent values of the potential for atoms belonging to the two electrodes, i.e. one value for each electrode with a di↵erence between the two electrodes, as illustrated in panels 2a and 2b for the 2D and 3D-slab PBC cases, respectively. Recently, an alternative approach has been proposed to allow the use of 3D PBC without the need to introduce vacuum, hence decrease the computational cost, and facilitate the implementation in standard simulation packages. The method introduced by Dufils et al. [START_REF] Dufils | Simulating Electrochemical Systems by Combining the Finite Field Method with a Constant Potential Electrode[END_REF], illustrated in Figure 2c, consists in applying a finite electric field to the electrolyte via an extended Hamiltonian. Such finite field methods, developed in the framework of the modern theory of polarization [START_REF] Stengel | Electric displacement as the fundamental variable in electronic-structure calculations[END_REF], had first been adapted to investigate EDLs near charged or polar insulator-electrolyte interfaces [START_REF] Zhang | Finite field methods for the supercell modeling of charged insulator/electrolyte interfaces[END_REF][START_REF] Sayer | Charge compensation at the interface between the polar NaCl(111) surface and a NaCl aqueous solution[END_REF]. By imposing a constant potential on all the electrode atoms (there is a single electrode in this setup) using the fluctuating charge method and an electric field on the electrolyte such that = EL z , with L z the box dimension in the direction perpendicular to the electrode, one can charge the capacitor by creating two EDLs, one on each side of the electrode. If the width of the electrode is su cient, one recovers the same results (in particular for the induced charges on the electrode or the structure of the interfacial electrolyte) as in the 2D periodic case of Figure 2a, with a reduced computational cost. A further benefit of this finite field approach is that it could easily be used with other methods to describe the metal (it had already been used in 3D PBC with the ICC method [START_REF] Arnold | E cient algorithms for electrostatic interactions including dielectric contrasts[END_REF]) with explicit electrode atoms, e.g. with the core-shell model of Ref. [START_REF] Geada | Insight into induced charges at metal surfaces and biointerfaces using a polarizable Lennard-Jones potential[END_REF], which capture the polarization of the metal but does not provide a handle to control the potential of the electrodes (hence the voltage between them), or even in ab initio simulations [START_REF] Bonnet | First-Principles Molecular Dynamics at a Constant Electrode Potential[END_REF].

B. Enough is enough? Sampling configurations

While the previous section presented the di culties related to the description of the system, we now turn to the theoretical and numerical challenges pertaining to the sampling of microscopic configurations of electrode-electrolyte interfaces. Molecular simulations are inherently related to Statistical Mechanics and the fundamental problem of sampling configurations from a thermodynamic ensemble, from which physical properties can be computed as averages. Even though most molecular simulations with constant-potential electrodes are performed in the canonical ensemble (constant number of electrolyte ions and molecules, N , volume V and temperature T ), one can in principle also consider a constant normal pressure on the electrodes (at least with rigid electrodes acting as pistons), or even grandcanonical (fixed chemical potential) and Gibbs ensemble (exchange between two systems) simulations. In the latter cases, which involve the exchange of particles, it might be necessary to resort to Monte Carlo simulations, as done e.g. in Refs. [START_REF] Kiyohara | Monte Carlo simulation of electrolytes in the constant voltage ensemble[END_REF][START_REF] Kiyohara | Monte Carlo Simulation of Porous Electrodes in the Constant Voltage Ensemble[END_REF][START_REF] Punnathanam | A Gibbs-ensemble based technique for Monte Carlo simulation of electric double layer capacitors (EDLC) at constant voltage[END_REF][START_REF] Stenberg | Grand canonical simulations of ions between charged conducting surfaces using exact 3D Ewald summations[END_REF] instead of molecular dynamics more commonly employed in this field. We focus here on three separate though related aspects, which apply generally to all the above situations, but limiting ourselves to the more common canonical ensemble: the issue of global electroneutrality, the statistical mechanics of the constant-potential ensemble, and importance sampling.

Electroneutrality

When studying an electrochemical system with electrodes held at a constant potential di↵erence by a generator at a classical molecular dynamics level, the main focus and interest is rather on the electrolyte properties. Importantly, the details of the device that applies a potential di↵erence are not taken into account. Instead the models consider an open system, very much like in a grand-canonical ensemble, exchanging charge (instead of particles) with a reservoir. In order to properly separate the electrochemical cell from the electric potential generator, the interaction between them should be negligible. This requires both subsystems to be electroneutral so as to cancel the leading term in the long-range electrostatic interaction. Although in the real device there are fluctuations and the electroneutrality of the cell is only verified on average, it should be enforced at a microscopic level in molecular simulations not explicitly including the charge reservoir. This can be achieved easily in descriptions with fluctuating charges by adding a Lagrange multiplier to enforce this constraint when computing the charges of the electrode atom. This multiplier corresponds to a potential shift applied to both electrodes simultaneously, while keeping the potential di↵erence unchanged [START_REF] Coretti | Mass-zero constrained molecular dynamics for electrode charges in simulations of electrochemical systems[END_REF][START_REF] Haskins | Evaluation of molecular dynamics simulation methods for ionic liquid electric double layers[END_REF][START_REF] Scalfi | Charge fluctuations from molecular simulations in the constant-potential ensemble[END_REF].

Statistical mechanics of the constant-potential ensemble

The natural ensemble associated with the above-mentioned simulations is the canonical constant-potential ensemble, characterized by the total number of atoms N , the volume V , the temperature T and the applied potential (more precisely, the set of potentials imposed no each electrode atom, usually one value per electrode). The charge distribution on the electrode atoms is then the conjugate variable of the applied potential and its fluctuations contain useful information. The aim of molecular simulations is to sample configurations from this ensemble in order to compute meaningful statistical averages. Along with approximations in the choice of force fields and the level of description of the quantum degrees of freedom, there are approximations related to how the chosen algorithm samples this thermodynamic ensemble. Using constant charges on electrode atoms ignores any charge fluctuation around the mean. One can use external inputs or post-processing to link the total electrode charge to an equivalent applied potential, but the configurations are not sampled according to their weight in the constant-potential ensemble.

Even though in principle one could also use Monte Carlo algorithms to sample the electrolyte configurations and electrode charges, to date most simulations using fluctuating charge models used molecular dynamics. In addition, except in earlier studies using Car-Parrinello (CP) dynamics with a fictitious mass for the additional degrees of freedom [START_REF] Siepmann | Influence of Surface-Topology and Electrostatic Potential on Water Electrode Systems[END_REF], the vast majority of simulations employing fluctuating charge models rely instead on Born-Oppenheimer (BO) dynamics, by assuming a separation of time scales between the electronic and nuclear degrees of freedom [START_REF] Reed | Electrochemical Interface Between an Ionic Liquid and a Model Metallic Electrode[END_REF]. This allows using larger time steps, but reduces the full phase space to that of electrolyte configurations, since to any such configuration corresponds a single set of electrode atom charges q ⇤ . A thorough discussion of the statistical mechanics of the constant-potential ensemble can be found in Ref. [START_REF] Scalfi | Charge fluctuations from molecular simulations in the constant-potential ensemble[END_REF], where we clarified in particular which observables can be computed exactly using BO dynamics on the charges. This is e.g. the case of the average total charge of an electrode hQ tot i = hQ ⇤ tot i. In contrast, the variance of the total charge,

⌦ Q 2 tot ↵
, which is related to the di↵erential capacitance by [START_REF] Kiyohara | Monte Carlo simulation of electrolytes in the constant voltage ensemble[END_REF][START_REF] Limmer | Charge fluctuations in nanoscale capacitors[END_REF]:

C di↵ = @ hQ tot i @ = ⌦ Q 2 tot ↵ (4)
with the voltage between the two electrodes and 1/ = k B T the thermal energy, includes a contribution from the suppressed thermal charge fluctuations. The di↵erential capacitance is then given by

C di↵ = C electrolyte di↵ + C empty di↵ = ⌦ Q ⇤2 tot ↵ + C empty di↵ (5)
where

⌦ Q ⇤2 tot ↵
is the variance sampled within BO dynamics and C empty di↵ , for which an explicit expression can be found in Ref. [START_REF] Scalfi | Charge fluctuations from molecular simulations in the constant-potential ensemble[END_REF], corresponds to the empty capacitor (i.e. in the absence of thermal fluctuations of the electrolyte).

Importance sampling

As usual in molecular simulations, one may also face sampling issues due to physical processes that result in long time scales not accessible by straightforward simulations. This includes e.g. slow transport within the electrolyte, or adsorption/desorption at the interface, or even phase transitions. In that case, one can resort to dedicated approaches, such as umbrella sampling to compute free energies, or transition path sampling to compute rates and analyze mechanisms. Examples of studies in the context of electrolyte interfaces, illustrated in Section IV C, involve collective variables such as the distance of an ion to the surface or between two ions to investigate their adsorption/desorption and pair formation/dissociation, or the number of water molecules in a probe volume, to investigate the hydrophilic/phobic behavior of the interface.

An importance sampling approach more specific to electrode-electrolyte interfaces was proposed in Ref. [START_REF] Limmer | Charge fluctuations in nanoscale capacitors[END_REF]. By adapting the standard histogram reweighting approach to the case of constant-potential simulations, it is possible to combine simulations performed at di↵erent potentials in order to optimally compute the properties of the system as a function of voltage, including for voltages for which no simulations are performed. The only requirement is an overlap of the distributions of the total charge P (Q tot ), or the joint distribution P (Q tot , p), with p a property of interest such as density profiles or orientational distributions, between the simulations at various voltages. This approach has been used for example in the study illustrated in Section IV B. Even though to the best of our knowledge it has never been done to date, one could use the same idea to derive a parallel tempering algorithm, exchanging replicas between simulations at di↵erent potentials (instead of temperature in the original method) to enhance the sampling of phase space. Histogram reweighting has also been used recently in Ref. [START_REF] Takahashi | A polarizable molecular dynamics method for electrode-electrolyte interfacial electron transfer under the constant chemical-potential-di↵erence condition on the electrode electrons[END_REF] to sample the so-called vertical energy gap related to electron transfer reactions (see Section IV D). Other strategies not requiring an explicit biais on the collective variable of interest, but taking advantage of statistical tools, include indirect umbrella sampling to reweight configurations using a bias on an auxiliary collective variable, which was applied in Ref. [START_REF] Kattirtzi | Microscopic dynamics of charge separation at the aqueous electrochemical interface[END_REF] to the Madelung potential experienced by ions at an electrochemical interface (see Section IV C).

IV. WHAT CAN WE LEARN FROM MOLECULAR SIMULATIONS OF ELECTRODE-SOLUTION INTERFACES?

Once a description of the system and a sampling strategy are chosen, one can perform classical molecular simulations to compute observable properties. In this section, we provide a selection (hence by no means exhaustive list) of illustrations of such properties typical of electrode-electrolyte interfaces: capacitance, interfacial structure and dynamics and an opening towards electrochemistry. We also present in Supplemental Material an application to electrowetting. All examples have in common the use of a fluctuating charge model to describe the electrode.

A. Capacitance

We begin by considering one of the most relevant property for practical applications, namely the capacitance. More precisely, one generally applies a voltage between two metallic electrodes separated by an electrolyte, and the charge ±Q tot on both electrodes, which fluctuates in response to the thermal motion of the liquid (see Section III B 2), allows the definition of an integral capacitance C int = hQ tot i / and a di↵erential capacitance C di↵ = @ hQ tot i /@ . Both quantities coincide only when the di↵erential capacitance does not depend on voltage, i.e. when the response to applied voltage is linear. Such a linear behavior is expected from continuum electrostatics for polar liquids behaving as pure dielectric media and for dilute electrolyte solutions. However, the organization of the molecular solvent and ions at the interface do not exactly follow the assumptions leading to simple analytical expressions, and molecular simulations allow to evaluate C int and C di↵ , as well as to correlate their evolution as a function voltage with structural changes in the interfacial fluid.

In practice, one could compute the average charge hQ tot i as a function of voltage and perform a numerical derivative to obtain the di↵erential capacitance. A much more e cient and accurate approach is to use the fluctuation-dissipation relation Eq. 5, which provides C di↵ from the variance of the charge distribution in a single simulation. We show in the Supplemental Material an illustration of this method with results from Scalfi et al. for two capacitors consisting of graphite electrodes separated by pure water or a 1M aqueous NaCl solution [START_REF] Scalfi | Charge fluctuations from molecular simulations in the constant-potential ensemble[END_REF]. For pure water, the abovementioned continuum prediction for the capacitance, taking the permittivity of the SPC/E water model used in this work, results in a value 3-4 times larger than the one obtained from the simulations, because it fails to account for the e cient screening of the field by the first layers of interfacial molecules [START_REF] Willard | Water at an electrochemical interface -a simulation study[END_REF][START_REF] Bonthuis | Profile of the Static Permittivity Tensor of Water at Interfaces: Consequences for Capacitance, Hydration Interaction and Ion Adsorption[END_REF][START_REF] Jeanmairet | Study of a water-graphene capacitor with molecular density functional theory[END_REF]. Adding salt increases the capacitance compared to pure water, but this increase is moderate for the considered concentration and inter-electrode distance. As for the pure water case, the above-mentioned continuum prediction is not accurate, as expected for Debye-Hückel theory at such a high concentration. The shortcomings of such continuous descriptions to estimate the capacitance (both for pure solvent and electrolyte solutions) reflect their limitations to account for the detailed organization of the interfacial fluid, which can also be investigated by molecular simulation, as discussed in the next section.

B. Interfacial structure

One of the most straightforward yet valuable information provided by molecular simulation is the microscopic structure of the system. In the case of planar electrode-electrolyte interfaces, density profiles as a function of the position z in the direction normal to the interface can be computed from the trajectories as an ensemble average as

⇢ ↵ (z) = ⌦P i2↵ (z i z)
↵ , with the Dirac delta function and where the sum runs over atoms i of type ↵. Similar definitions can be used to compute (1D, 2D or 3D) number, charge or polarization density profiles, distributions of molecular orientation, radial distribution functions, coordination numbers, etc. From the charge density profiles, one can also obtain electrostatic potential profiles by integrating the 1D Poisson equation, to make the link with the interfacial capacitance (in the presence of ions, using the potential drop across each interface) or the dielectric permittivity of the liquid (in the absence of ions, using the slope of the potential in the bulk region). The insights from molecular simulations are particularly valuable, since they capture explicitly the e↵ects of the discrete nature of ions and solvent molecules, which lead in particular to their layering near the planar electrode, and of their complex interactions (steric, electrostatic, dispersion, ...). This allows a detailed investigation of the composition and organization of the interfacial liquids and their evolution as a function of applied voltage. The results can also serve as reference data for simpler theories, such as the ones mentioned in the introduction. As an illustration, we consider here the complex case of a water-in-salt electrolyte (WiSE), proposed recently as promising for high-voltage batteries and supercapacitors [START_REF] Suo | A New Class of Solvent-in-Salt Electrolyte for High-Energy Rechargeable Metallic Lithium Batteries[END_REF][START_REF] Wang | Superconcentrated Electrolytes for a High-voltage Lithium-ion Battery[END_REF], by reporting some results of Li et al. [START_REF] Li | Capacitive Performance of Water-in-Salt Electrolytes in Supercapacitors: A Simulation Study[END_REF]. Figure 3a shows a system consisting of a 20 mol/kg aqueous LiTFSI solution between two graphite electrodes maintained at a constant potential using the method described in section II A 4. In such an electrolyte, the number of water molecules per ion pair is lower than 3 (compared to ⇡50 for a typical 1 mol/L concentration), so that standard EDL theories are not expected to hold. The density profiles in Figure 3b clearly illustrate the layering of the interfacial fluid even in the absence of applied voltage and a bulk region far from the surfaces, as well as the change in composition and structure of both interfaces under voltage. One can note in particular the increase in water concentration near the negative electrode (left) accompanying that of Li + , which approach the surface with their limited hydration shell. The structure on the positive electrode is further illustrated in panel 3c, which shows the density profiles for oxygen and fluorine atoms from the TFSI anion for two voltages. The most dramatic structural change occurring at this surface between = 1.5 and 2 V is the reorientation of the anions, leading to a closer approach of O atoms (carrying a larger partial charge than F) to the surface, also visible on the typical snapshots of panels 3d and 3e. Even though we do not discuss this further here, it was also found in this work, using importance sampling techniques, that changes in the interfacial composition and structure as a function of voltages could be linked to peaks in the di↵erential capacitance [START_REF] Li | Capacitive Performance of Water-in-Salt Electrolytes in Supercapacitors: A Simulation Study[END_REF], as previously reported for room temperature ionic liquids [START_REF] Merlet | The electric double layer has a life of its own[END_REF]100] and concentrated electrolytes [101]. One can finally note that, even in such complex cases, molecular simulation results can be used as reference data to validate and parameterize simpler theories, as was done e.g. by McEldrew et al. in a contemporary study of a similar WiSE between uniformly charged walls [START_REF] Mceldrew | Theory of the Double Layer in Water-in-Salt Electrolytes[END_REF].

C. Interfacial dynamics

Applications of electrode-electrolyte interfaces to energy storage (batteries, supercapacitors) not only aim at improving the amount of energy that can be stored, by increasing the capacitance of the interface or the voltage, but also at reducing the charge/discharge time to increase the power that can be delivered. In the standard RC circuit picture, the charging time depends on the resistance of the bulk electrolyte and the interfacial capacitance, but not on the interfacial dynamics. The situation in realistic porous electrode materials is more complex and molecular simulation has contributed to a better understanding of the charging dynamics in these systems, also making the link with more involved equivalent circuit models used in electrochemistry [102-105]. However, even in the simpler case of a planar electrode and a solvent-based electrolyte, which is the focus of the present work, the dynamics of the interfacial fluid is modified by the presence of the solid. This in turn may modify not only the charging dynamics but also the kinetics of electron transfer reactions (which will be discussed in more detail in section IV D). We illustrate here some important findings on the dynamics at a water-platinum interface uncovered using molecular simulation.

timescales larger than nanoseconds. The characteristic time for this surface relaxation, τ s ( 1 ns, is much larger than that for typical equilibrium density fluctuations in the bulk liquid, τ b ∼ 5 ps. Therefore, although the presence of the surface introduces a static inhomogeneity, the water bound to this surface introduces a dynamic inhomogeneity. The resultant separation of timescales between bulk and surface reorganization is illustrated in Fig. 1. Fig. 1, Upper shows snapshots during slow reorganization of the surface water dipoles, whereas Fig. 1, Lower illustrates faster interfacial fluctuations.

These two features of the metal interface, the static heterogeneity of the extended interface and the slow dynamics of water at its surface, cause a decoupling of ensemble and dynamic averaging on timescales t < τ s . The decoupling implies that for a given configuration of the adlayer, liquid water swiftly equilibrates, and for t τ b the properties of the subsequent hydration layers are in dynamic equilibrium. Over intermediate timescales, however, temporal heterogeneity of the hydrogen bond network couples to the dynamically heterogeneous properties of the interface.

The next section illustrates the dependence of the adlayer structure with the exposed metal surface geometry and shows how the passivation of the hydrogen bond network on the surface creates a liquid-vapor-like interface that attracts hydrophobic particles. The subsequent section shows how frustration of the water acteristic snapshots of the adlayer obtained from our simulations, as we wetting (Fig. 2 B andD). (A liquid metal surfaces and adlayers in Fig. 2 ecules are not rendered.) For the locally fourfold coordinated and a projection of local hydrogen-bondi structure of water on the surface dipoles oriented parallel to the surf sites are occupied. At any particular exist on the surface, separating plan by 90°turns in their orientations. Fo are locally sixfold coordinated and a spacings that are commensurate wit coordination frustrates preferred bo surface has regions of local hexago rounding a vacancy, that are seen i water absorbed on the 111 surface (FCC) metals [START_REF] Kiyohara | Monte Carlo Simulation of Porous Electrodes in the Constant Voltage Ensemble[END_REF]. Because such a h tile space, this surface also has a interstitials that occupy the empty to point away from the surface on aver average coverage of about 85% of a the lattices are entirely regular, and the hydrogen-bonding network is dy order within the adlayer dictates th timescales. Similar hydrogen-bondin experimentally under ultrahigh vacu atures on Pd(111) [START_REF] Van Megen | The grand canonical ensemble Monte Carlo method applied to the electrical double layer[END_REF] and in water-h

The presence of extended interf solvated metal surface, is expected subsequent solvent layers over dis bulk correlation length. For a liqu vapor, such as water at ambien homogeneities can give rise to a de interfaces subsequently have larger graphs the mean density of water distance away from the metal surfa the adlayer depends intimately on rounding water is fairly insensitive t find for both surface geometries tha away from the interface exhibits a sh indicative of the adlayer, followed by ∼3 Å thick. Density oscillations de surface. The asymmetry between acceptors at the interface results however, the effective interaction wi as to allow the formation of capillar density oscillations seen away from

Water-Metal Interface Is Hydroph

The unbalanced attraction immedia enough to make solvation of ideal Fig. 1. Illustration of the separation of timescales between reorganizing surface configurations, which occurs on average every τs, and reorganizing the bulk density, which occurs on average every τb. Small tick marks are separated by 20 ps, which is on the order of although larger than timescales for typical density fluctuations. Large tick marks are separated by 100 ps, which is on the order of the typical relaxation times for relevant interfacial fluctuations.

2 of 6 | www.pnas.org/cgi/doi/10.1073/pnas.1301596110 CHEMISTRY formed into two separate and oppositely polarized ionic solvation shells. This solvent reorganization has been identified as the rate-limiting step for aqueous ion dissociation, and numerous efforts have been aimed at quantifying water's role in this process (1)(2)(3)(4)[START_REF] Merlet | On the Molecular Origin of Supercapacitance in Nanoporous Carbon Electrodes[END_REF]. These efforts have revealed that, although the overall process of aqueous ion dissociation is usually thermodynamically favorable, solvent reorganization leads to the emergence of a free-energy barrier that is on the order of typical thermal energies, kBT , where kB is Boltzmann's constant and T is the temperature.

Water ions, specifically hydroxide and hydronium, are the products of proton transfer to and from an individual water molecule. These ions can easily integrate into the aqueous hydrogenbonding network and leverage the Grotthuss-like shuttling of protons for delocalized and rapid transport [START_REF] Bazant | Di↵use-charge dynamics in electrochemical systems[END_REF]. This feature leads the solvation properties of water ions to differ from those of similarly sized monovalent ions [START_REF] Kiyohara | Monte Carlo simulation of electrolytes in the constant voltage ensemble[END_REF]. The solvent's role in mediating the separation of water ions is thus different from that of simple monovalent salts. It has been shown with ab initio simulation that the dissociation of bound water ions, a process known as autoionization, requires the well-timed coordination of solventinduced electric field fluctuations and the making/breaking of hydrogen bonds. These correlated fluctuations result in a nearspontaneous relocation of a proton from a neutral water molecule to a newly formed hydronium (H3O + ) [START_REF] Bazant | Di↵use-charge dynamics in electrochemical systems[END_REF][START_REF] Taylor | First principles reaction modeling of the electrochemical interface: Consideration and calculation of a tunable surface potential from atomic and electronic structure[END_REF], leaving behind a negatively charged hydroxide ion (OH ).

Because the water-mediated separation of water and classical ions occurs via different mechanisms, they can be affected differently by changes in the aqueous environment. In this work, we consider specifically the environmental changes that are associated with an electrochemical interface of an extended platinum and B. This slowly evolving water monolayer affects the structure and dynamics of the adjacent bulk liquid [START_REF] Hautman | Molecular dynamics simulation of water beween two ideal classical metal walls[END_REF]. For pure water systems, this effect is subtle in comparison with the dramatic slowdown of the monolayer itself; however, ions at this liquid interface can incorporate part of the monolayer into their solvation shell and thus couple directly to the slow dynamics of adsorbed water.

Simulating Rare Events in Heterogeneous Environments

To study ion dynamics at the aqueous electrode interface, we performed atomistic simulations of liquid water in contact with the [START_REF] Zhang | Water dynamics at electrified graphene interfaces: a jump model perspective[END_REF] surface of an extended platinum electrode. We used two different model systems: The first was designed to study the dissociation of a classical ion pair, Na + and I , and the second was designed to study the recombination of water ions, H3O + and OH . In the first model system, we describe the aqueous solution using classical force fields. Specifically, we describe water using the SPC/E model ( 27), and we describe the ions similarly, as spherically symmetric point charge particles. This efficient combination has been demonstrated to accurately reproduce the molecular structure and dynamics of liquid water as well as experimental measures of ion hydration and mobility. Although we studied this particular ion pair, previous studies on the adsorption free energies and mobilities of other alkyl halides have shown consistent qualitative changes between behaviors at and away from the electrode [START_REF] Hautman | Molecular dynamics simulation of water beween two ideal classical metal walls[END_REF]. Unfortunately, this classical nondissociative model is inadequate to describe the dynamics of water ions, whose transport is facilitated by the making and breaking of covalent OH bonds. Thus, in the second model system, we used an ab initio model of water based on density functional theory (6). Kattirtzi et al. investigated the dissociation of ion pairs near a Pt electrode using classical (for Na + -I ) and ab initio (for water ions H 3 O + -HO ) molecular dynamics [START_REF] Kattirtzi | Microscopic dynamics of charge separation at the aqueous electrochemical interface[END_REF]. Figure 4a shows configurations for the Na + -I case near a (111) Pt electrode surface, modeled using the method described in section II A 4 and the force field of Siepmann and Sprik [START_REF] Siepmann | Influence of Surface-Topology and Electrostatic Potential on Water Electrode Systems[END_REF] for the water-Pt interaction, which includes two and three body terms. The three configurations are typical for a contact ion pair, an intermediate state and a dissociated state, corresponding to the regions visible in the free energy profiles as a function of the interionic distance R, obtained by umbrella sampling (see section III B 3) and shown in panel 4b. This panel further shows that the free energy profile along this distance R is in fact only slightly a↵ected by the vicinity of the surface, with no change in the relative free energies of the states (suggesting limited changes in their solvation structure with respect to the bulk) and a small increase in the activation barrier for the dissociation/recombination, which would suggest only a slight decrease in the corresponding rates. However, as in the bulk [107], the interionic distance is not su cient to examine the mechanism. It also involves collective water fluctuations, which were described by Kattirtzi et al. using a simple collective variable, namely the Madelung potential on the ions. From the change in free energy barrier including this variable (with an increase of 1.5 k B T with respect to the bulk), one would predict a decrease in the dissociation rate near the electrode k elec ⇡ 0.2 k bulk . However the dissociation rate measured in the simulations is in fact k elec ⇡ 0.02 k bulk , i.e. the decrease is ten times larger. This is because the free energy is not su cient to predict the rate, as the dynamics of water fluctuations also plays a role [START_REF] Kattirtzi | Microscopic dynamics of charge separation at the aqueous electrochemical interface[END_REF].

The dynamics of water near a Pt surface was in fact considered in earlier work by some of the same authors [106, 108-110]. It was found that these dynamics can be highly collective, heterogeneous and slow. Figure 4c illustrates the slow time scales associated with structural rearrangements of water near a (100) Pt surface, compared to that for the relaxation of bulk density fluctuations. On this particular surface, water molecules in the first adsorbed layer form an ordered network of strong H-bonds between them. The lack of interactions with molecules beyond this layer (see the bottom parts of this panel) in fact results in an unexpected hydrophobic behavior, with e.g. a large contact angle of a water droplet on this adlayer. From the dynamical point of view, switching between distinct but equally probable H-bond patterns (see top parts of Figure 4c) requires collective rearrangements, associated with very long time scales and displaying dynamical heterogeneity [106]. While some results observed on this particular surface should apply generally to other Pt surfaces and other metals, the details of the dynamics crucially depend on the interplay between the geometry of the solid lattice and the interactions with the fluid, so that molecular simulation is ideally suited for such investigations. It was for example used recently to analyze the reorientation dynamics of water molecules at an electrified graphene interface [START_REF] Zhang | Water dynamics at electrified graphene interfaces: a jump model perspective[END_REF].

D. Towards electrochemistry

In general, simulating reactive systems is challenging in molecular simulation: Reactive force fields need to be employed, which are di cult to parameterize. Nevertheless, many electrochemical reactions do not involve bond breaking/formation, but only a change of the redox state of the species. This observation is central in Marcus theory, which was developed in order to calculate reaction free energies for electron transfer [112]. In short, this classical transition state theory introduces a key quantity, the solvent reorganization energy, associated with the solvent rearrangement when the redox state of the solute changes. In classical molecular simulation, Warshel introduced the vertical energy gap, that is the potential energy di↵erence between the oxidized and the reduced species at fixed solvent coordinates, as the key reaction coordinate for simulating electrochemical reactions [113].

Although this development allowed straightforward studies of electron transfer reactions in bulk solutions, in an electrochemical cell setup an additional complexity arises because the electrons are transferred to the electrodes instead of vacuum. Madden and co-workers proposed an approach in which the constant applied potential condition is enforced before and after the electron transfer reactions, leading to di↵erent charge organizations inside the electrodes [114]. Using this approach, they showed how the solvent reorganization energy, and hence the free energy curves, depend on the applied potential and on the distance of the redox species from the electrode in a molten salt. Surprisingly, the variations they observed where much smoother than expected from the pronounced oscillations in the mean electrical potential inside the double-layer region, a result which was confirmed in further studies [115].

The next challenge for classical simulations relies on the study of active interface. First steps in this direction were made using the Empirical Valence Bond method to describe the transfer of proton in water on silver or platinum electrodes [116,[START_REF] Wiebe | Water Structure and Mechanisms of Proton Discharge on Platinum Electrodes: Empirical Valence Bond Molecular Dynamics Trajectory Studies[END_REF]. Switching in electrometallization cells has also been studied using reactive force fields (REAXFF) together with the charge equilibration method for the electrostatic part [START_REF] Onofrio | Voltage equilibration for reactive atomistic simulations of electrochemical processes[END_REF][START_REF] Onofrio | Atomic origin of ultrafast resistance switching in nanoscale electrometallization cells[END_REF]. The electrodes, made of copper, were allowed to dissolve inside an amorphous silica and redeposit upon application/cancellation of a voltage, eventually leading to the formation of a conductive contact between the two electrodes. The extension of this approach to more conventional solvents is yet hindered by the di culty of developing accurate reactive force fields for complex systems. Another strategy was recently proposed, in which the charge transfer is modeled as a stochastic process, using methods borrowed from grand canonical Monte Carlo simulations [START_REF] Dwelle | Constant potential, electrochemically active boundary conditions for electrochemical simulation[END_REF]. We can therefore expect that more simulations of reactive electrochemical interfaces will be proposed in future years.

V. SUMMARY POINTS 1. Over the last 30 years, classical molecular simulations have emerged as an essential tool to investigate the properties of electrode-electrolyte interfaces: By providing a compromise between an atomic description and a computational cost allowing a su cient sampling of relevant electrolyte configurations, they o↵er the possibility to face the challenge posed by such interfaces bridging electrons in a solid and ions in a solvent, i.e. Quantum Chemistry and Statistical Physics.

2. We have provided an overview of the models and methods to describe the metallic character of the electrode and its interactions with the electrolyte solution, and we discussed simulation setups, algorithms and statistical tools to sample configurations from the statistical ensemble corresponding to constant-potential electrodes.

3. We have illustrated a selection of properties which can be investigated with classical molecular simulations, with examples on the capacitance, the interfacial structure and dynamics, electrowetting, as well as steps towards electrochemistry.

4. The results obtained at this classical level can be compared directly to experiments probing the interface on the molecular scale, but also serve as reference data for simpler theories of electric double layers, generally used to interpret macroscopic electrochemical experiments.

5. Even though we restricted ourselves to the simpler yet practically relevant and physically rich case of planar interfaces and solvent-based liquids, from pure solvent to water-in-salt-electrolytes, some of the methods discussed here also apply to more complex (e.g. disordered nanoporous) electrodes and to solvent-free electrolytes (room temperature ionic liquids).

VI. FUTURE ISSUES

1. The classical description of a metallic electrode remains a challenging issue. For example, most electrode materials do not behave as perfect metals. Electrostatic screening inside the solid can be captured at the continuum level using Thomas-Fermi theory, an approach which can also be exploited for classical molecular simulation. Earlier attempts included this e↵ect of screening on the electrolyte via an external potential [START_REF] Rose | Solvation of Na + and Cl at the water-platinum (100) interface[END_REF]120], but new strategies have been proposed very recently, based on fluctuating charges [START_REF] Scalfi | A semiclassical Thomas-Fermi model to tune the metallicity of electrodes in molecular simulations[END_REF] or on mobile charges [START_REF] Schlaich | Wetting transition of ionic liquids at metal surfaces: A computational approach to electronic screening using a virtual Thomas-Fermi fluid[END_REF]. This will be particularly useful to address the e↵ect of the metallic character of the electrode on the properties of the interfacial fluids, such as the nanoscale capillary freezing of ionic liquids [START_REF] Comtet | Nanoscale capillary freezing of ionic liquids confined between metallic interfaces and the role of electronic screening[END_REF].

2. The treatment of non-electrostatic interactions also o↵ers opportunities for improvement. The flexibility of the electrode has almost never been addressed, even though it might play a role in particular in porous electrodes. Due to the importance of electrochemical reactions, it is likely that reactive force fields will also gain importance, if their parametrization can be extended to more complex systems. The availability of accurate reference data from ab initio calculations, in particular quantum Monte Carlo, will allow a better parametrization of force fields based on existing models. One can also anticipate that the rapid development of Machine-Learning based approaches (e.g. force fields based on neural networks) in molecular simulation will also reach the community working on electrode-electrolyte interfaces.

3. Another promising strategy is o↵ered by hybrid approaches, coupling several levels of descriptions. As an example, a mesoscopic description of the solvent, based on classical density functional theory, was recently coupled to a fluctuating charge model of a graphene electrode [START_REF] Jeanmairet | Study of a water-graphene capacitor with molecular density functional theory[END_REF] and used to investigate electron transfer reactions in a classical solute [START_REF] Jeanmairet | A molecular density functional theory approach to electron transfer reactions[END_REF]. In the other direction, the development of classical models of the interface might also help improving the coupling between the electrolyte and the electrode in QM/MM simulations. Supplementary Figure 1 illustrates the discussion of Section 4.1 of the main text with results from Scalfi et al. for two capacitors consisting of graphite electrodes separated by pure water or a 1M aqueous NaCl solution, under a voltage = 0 V [1]. The variance of the distribution P (Q ⇤ tot ) of the total charge computed in Born-Oppenheimer dynamics simulation in the constant-potential ensemble (see Section 3.2.2), provides the contribution of the electrolyte,

C electrolyte di↵ = ⌦ Q ⇤2 tot ↵
, to the di↵erential capacitance. This contribution is large compared to that of the empty capacitor, C empty di↵ , consistently with the observation of Haskins and Lawson for ionic liquids [2] -and with the fact that it had escaped the attention of the community until recently.

II. ELECTROWETTING

Electrochemistry also goes beyond simple setups with two electrodes facing each other and an electrolyte in between. For example, electrowetting-on-dielectric setups were proposed for switchable optical devices. An example of design is made of two concentric electrodes (with di↵erent radii) separated from the liquid, which is made of a water nanodrop, by a layer of dielectric [4]. The use of constant potential simulations allowed an accurate representation of such systems at the molecular scale [3], as shown on Figure 2a. Applying a potential between the two electrodes results in a noticeable spreading of the drop, an e↵ect which is reversible. Figure 2b shows the cross-section profiles at various applied potentials, from which contact angles can be extracted. In this example the angle varies from 111 at null voltage down to 84 at 4 V. Changing the radii of the two electrodes leads to di↵erent values but the spreading e↵ect persists. At the macroscale, the variation of the contact angle with the applied potential is well predicted by the Young-Lippmann equation, but simulations showed that it does not hold anymore at the nanometric scales.

V ΔΨ Static Response. By attracting dipolar water molecules, the local electric field generally enhances materials' wetting propensity. The field effect on wetting can be characterized through nanodrop contact angle calculations. 17 Our results for a set of interelectrode potential differences in the interval 0-4 V show consistently enhanced spreading under the applied field, in analogy to experimental 34 and continuum-simulation 35,36 results. Figure 2 presents the fitted droplet contours at varied voltages for both System A and B. In Table I, we present simulated contact angles as functions of the applied voltage for both systems. We observe no hysteresis in contact angles provided at least ≃150 ps equilibration is used after each voltage change.

For constant areal capacitance, C/A sl , the Young-Lippman eq (eq 1) predicts cos θ c to rise in proportion to the voltage reminiscent of the saturation behavior obser experiments. In the absence of charge carri saturation cannot be attributed to charge r discussed in ref 25. To explore possible rel the system's capacitance C, we determin function of wetted area, A sl (U), in both sys q is the absolute value of the charge on the to maintain preset voltage U at the spec droplet spreading. Area A sl equals πR b 2 , wher R s is the radius of curvature of the drop wi V ΔΨ through nanodrop contact angle calculations. Our results for a set of interelectrode potential differences in the interval 0-4 V show consistently enhanced spreading under the applied field, in analogy to experimental 34 and continuum-simulation 35,36 results. Figure 2 presents the fitted droplet contours at varied voltages for both System A and B. In Table I, we present simulated contact angles as functions of the applied voltage for both systems. We observe no hysteresis in contact angles provided at least ≃150 ps equilibration is used after each voltage change.

= + ⎛ ⎝ ⎜ ⎞ ⎠ ⎟ R V ( ) 3 (1 cos ) (2 cos
For constant areal capacitance, C/A sl , the Young-Lippman eq (eq 1) predicts cos θ c to rise in proportion to the voltage remi expe satur discu the func q is to m drop R s is By performing additional non-equilibrium simulations where the potential was suddenly applied or released, it was observed that the retraction at zero voltage is much faster than the corresponding spreading, an e↵ect which could be explained by di↵erent liquid/solid friction properties.

Fig. 1 8 2

 18 Fig.1Polarizable Lennard-Jones model for gold. a Extension of the simple Lennard-Jones model with dummy electrons to add features of the free electron gas. The virtual electrons rest at the atom core and carry a mass of 1 au. b Visualization of the dummy electrons on the Au[START_REF] Zhang | Water dynamics at electrified graphene interfaces: a jump model perspective[END_REF] surface in the presence of an adsorbed sodium ion in vacuum. The induced charges spread across several atomic layers laterally and beneath the top atomic layer. c The energy expression contains terms for harmonic bond stretching, Coulomb energy, and van-der-Waals energy (Lennard-Jones potential). d The model uses five independent parameters (highlighted in bold) including the mass of the dummy electron me, a combination of the charge q and the bond stretching constant kr, whereby a certain ratio α = q 2 /(2kr) determines the magnitude of the image potential, as well as the Lennard-Jones parameters σ, εcore, and εe. The total mass of the gold atom (mcore + me) and the rest position of the dummy electron r0 = 0 Å are constants

FIG. 1 .

 1 FIG.1. Sketch of a section of two domains R1 and R2, characterized in the outer and inner regions by the dielectric permittivity ⑀1 and ⑀2, respectively. A discretized surface element is highlighted, and the normal vector n and the electric fields E1 and E2 at the dividing surface are shown.

Fig. 1 8 2FIG. 1 :

 181 Fig.1Polarizable Lennard-Jones model for gold. a Extension of the simple Lennard-Jones model with dummy electrons to add features of the free electron gas. The virtual electrons rest at the atom core and carry a mass of 1 au. b Visualization of the dummy electrons on the Au[START_REF] Zhang | Water dynamics at electrified graphene interfaces: a jump model perspective[END_REF] surface in the presence of an adsorbed sodium ion in vacuum. The induced charges spread across several atomic layers laterally and beneath the top atomic layer. c The energy expression contains terms for harmonic bond stretching, Coulomb energy, and van-der-Waals energy (Lennard-Jones potential). d The model uses five independent parameters (highlighted in bold) including the mass of the dummy electron me, a combination of the charge q and the bond stretching constant kr, whereby a certain ratio α = q2/(2kr) determines the magnitude of the image potential, as well as the Lennard-Jones parameters σ, εcore, and εe. The total mass of the gold atom (mcore + me) and the rest position of the dummy electron r0 = 0 Å are constants

FIG. 2 :

 2 FIG. 2: Simulation setups to simulate a capacitor consisting of graphite electrodes separated by a 1 mol L 1 aqueous NaCl solution under a voltage = 2 V; the color of the electrode atoms reflects their instantaneous charge. (a) Periodic boundary conditions (PBC) are applied only in the two directions along the surfaces, but not in the direction perpendicular to the electrodes, and the two electrodes are maintained at di↵erent potentials, corresponding to a voltage. (b) PBC are applied in all directions, but the simulation box also includes vacuum (here the length of the box in this direction is three times larger than the electrochemical cell) and a "slab" correction is applied when computing electrostatic interactions; as in the previous case, the two electrodes are maintained at di↵erent potentials. (c) PBC are applied in all directions, without vacuum; all electrode atoms are maintained at the same potential and the capacitor is charged by applying an electric field to the electrolyte such that = ELz, with Lz the box dimension in the direction perpendicular to the electrode.

FIG. 3 :

 3 FIG. 3: (a) Capacitor consisting of graphite electrodes separated by a water-in-salt-electrolyte (20 mol/kg aqueous LiTFSI solution). An all-atom model is used, with H (white), O (red), Li (green), F (orange), C (cyan), S (yellow), and N (blue) atoms, but for clarity the TFSI anions beyond interfaces are represented by isosurfaces. (b) Density profiles for Li (blue), water (orange) and TFSI (green) at = 0 V (top) and 3 V (bottom). (c) Density profiles for O and F from TFSI anions on the positive electrode, for = 1.5 and 2 V; the shaded areas indicate the first and second adsorbed layers; the key indicates the potential of the positive electrode, assumed to be + = /2). (d,e) Snapshots illustrating the structure close to the positive electrode from simulations at = 1.5 and 2 V. Adapted with permission from Ref. [97], J. Phys. Chem. C 2018, 122, 23917, Copyright (2018) American Chemical Society.

Fig. 1 .

 1 Fig.1. A and B contain typical snapshots going from a recombined state (black-bordered panel) to a dissociated state (red-bordered panel) for classical ions and water ions, respectively. The positive ion is highlighted in green, and the negatively charged ion is highlighted in yellow. C and D contain plots of the free-energy profile as a function of R, the interionic separation for classical ions and water ions, respectively. The data plotted in red correspond to ions in the bulk, and the data plotted in blue correspond to ions at the electrode interface.

FIG. 4 :

 4 FIG. 4: (a) Typical configurations for the dissociation of an ion pair (Na + in green, I in yellow) in water near a (111) Pt electrode surface, from a contact ion pair (black-bordered panel) to a dissociated state (red-bordered panel). (b) Free energy profile as a function (in units of the thermal energy kBT = 1/ ) of the distance R between the ions in the bulk (red) and near the surface (blue). (c) Illustration of the slow time scale ⌧s associated with structural rearrangements of water near a (100) Pt surface, compared to that ⌧ b for the relaxation of bulk density fluctuations (small ticks are separated by 20 ps, larger than typical values for ⌧ b ). The blue surface in the bottom panels indicates the instantaneous interface between the first adsorbed layer and the rest of the fluid, which fluctuates due to local changes in the hydrogen bond pattern. Collective fluctuations eventually leading to global rearrangements such as between the top panels occur on much longer time scales. Panels (a) and (b) adapted with permission from Ref. [93], PNAS, 2017, 114, 13374; panel (c) adapted with permission from Ref. [106], with permission of PNAS, 2013, 110, 4200.

4 .FIG. 1 :

 41 FIG. 1: Capacitors consisting of graphite electrodes separated by (a) pure water and (b) a 1M aqueous NaCl solution under a voltage = 0 V; the color bar indicates the instantaneous charge of the electrode atoms. (c) Distribution of the total charge of the electrode Q ⇤ tot for both systems, computed within Born-Oppenheimer dynamics simulations; solid lines are normalized Gaussian distributions with the corresponding standard deviation. (d) Contributions to the di↵erential capacitance per unit area C di↵ /A at 0 V: from the empty capacitor, C empty di↵ , and from the charge fluctuations induced by the electrolyte, C electrolyte di↵

  )

  Drop dimensions for different voltages ar 1 of the SI.Charge q and capacitance C depend on drop, which requires a finite time to respon of the voltage. The charge equilibration is i 4. The solid symbols in Figure5repres capacitances for the four values of the dr corresponding to U = 1, 2, 3, or 4 V. correspond to capacitances at zero voltag extrapolation to respective A sl . The capacit (hypothetical states with A sl = 0, θ c = π) diff

Figure 2 .

 2 Figure 2. Vertical cross section of the water droplet on the insulator surface, showing drop's contours (solid red line) fitted through simulation data points for different electric fields. The solid line in magenta represents the position of the surface reference plane. Left, System A; right, System B.

Figure 3 .

 3 Figure 3. Variation of the cosine of contact ang applied potential difference (U) between the e A (black) and System B (blue). Solid circles a with the dotted lines showing a likely interpol Error bars correspond to the estimated ±1°un Supporting Information). Long-dashed lines sh eq 1.

Figure 2 .

 2 Figure 2. Vertical cross section of the water droplet on the insulator surface, showing drop's contours (solid red line) fitted through simulation data points for different electric fields. The solid line in magenta represents the position of the surface reference plane. Left, System A; right, System B.

FIG. 2 :

 2 Figu appli A (b with Erro Supp eq 1

ACKNOWLEDGMENTS

The authors are grateful to many colleagues with whom they have worked or exchanged on this topic over the years, in particular Paul Madden, Michiel Sprik, David Limmer, Céline Merlet. This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement No. 771294 and grant agreement No. 863473). This work was supported by the French National Research Agency (Labex STORE-EX, Grant ANR-10-LABX-0076, and project NEPTUNE, Grant ANR-17-CE09-0046-02).

DISCLOSURE STATEMENT

The authors are not aware of any a liations, memberships, funding, or financial holdings that might be perceived as a↵ecting the objectivity of this review.