Improving Continuum Models to Deﬁne Practical Limits for Molecular Models of Electriﬁed Interfaces

We develop a continuum theory of electrolyte solutions in contact with a metal electrode, based on a generalized free energy functional, and use it to explore the structure of the electric double layer at different electrochemical conditions. The model captures the effects of speciﬁc adsorption of ions and solvent polarization, and can be applied on the same footing to cases with non-zero faradaic current, beyond the classical double-layer regime. These advances permit the prediction of peculiar character in the ion proﬁles at electrode potentials near the redox level, exploration of the electrochemical stability of the interface, and differentiation between the mechanisms of electron and ion transport and associated time scales. The developed methodology enables us to self-consistently determine the fundamental limits for a microscopic description of biased interfaces, in terms characteristic sizes and time scales of relevant processes, within atomistic and ab initio molecular dynamics simulations. Mathematical of Electrochemical Systems at Multiple Scales

As the efficiency of electrochemical conversion is dictated by the energetics of elementary charge transfer reactions occurring at the electrode/electrolyte interface, the formulation of methods of control over these reactions has become the pivotal goal of fundamental material science for advanced energy storage systems. The ultimate goal of theoretical electrochemistry then is to explore the energetics of the relevant processes from an atomistic perspective by means of ab initio molecular dynamics (AIMD) where ionic and electronic components of the system are described explicitly and are constrained by various realistic conditions, such as the external bias, the temperature and the electrolyte concentration. However, despite more than a decade of efforts this goal is still far from being fulfilled.
There are several interconnected fundamental problems to resolve before an atomistic approach could be routinely used to describe electrochemistry, if only at the level of a half-cell (i.e., a single electrodeelectrolyte interface). First, there are at least three distinct regimes of performance of the electrochemical interface and each of these is characterized by different sets of the time and space scales: a) the classical electric double layer (EDL) regime when no charge transfer is assumed, b) the regime of an electrode at equilibrium with the electrolyte (open circuit conditions, (OCP)) when only an exchange current in a form of mutually compensating fluxes of charged particles is possible, and c) the "true" electrochemical regime with a non-zero net faradaic current. As long as the goal of atomistic simulations is to reach a state when one would observe the actual electrochemical reactions driven by the bias, the computational electrode potential scanning protocol should carry the system through a set of (quasi-)equilibrium points at which the atomistic structure is consistent with the targeted external conditions. The key requirement for such simulations is thus thermodynamic stability of the interface that encompasses all degrees of freedom. Currently, there is no such uniform atomistic methodology for the electrode potentials in a range of ±1 V around the OCP since it would have to cover the processes whose characteristic times range from attoseconds (electron transfer (ET) time) to hundreds of picoseconds (solvent dielectric relaxation) at length scales ranging from angstroms to tens of nanometers. Instead, we currently have methodologies that either elaborate the structure of the classical EDL [1][2][3] or study the charge transfer per se assuming consistency between the atomistic structure of the EDL and the external conditions. 4 The crossover between the two approaches is not well defined, such * Electrochemical Society Member. z E-mail: abaskin@lbl.gov that, for an unknown system, the simulated electrochemical activity may be an artifact of the methodology, which relies on an initial guess for the structure of the EDL that might prove inconsistent with the aforementioned external conditions or some basic assumptions. The second problem is the experimental relevance of the computational electrode potential that is a long-standing generic issue for theoretical electrochemistry intimately related to the boundary conditions of the half-cell setup. In conventional electrochemistry, 5,6 the absolute single electrode potential, E, is defined by two contributions -the Galvani potential drop, ϕ, across the interface, which reflects the response of the electrolyte, and the chemical potential of electrons in the electrode μ e . The relation between these two is not trivial and strongly depends on the regime at which the electrode operates. Therefore, the computational electrode potential that is comparable to experimental values is not known a priori, since at each potential the interface is characterized by a (unique) thermodynamic (steady) state which is unknown a priori and which must instead be determined a posteriori. This, in turn, implies that in contrast to the currently adapted potentiostat schemes 7-10 the electrode potential cannot be a free external parameter and should rather be the output of a selfconsistent procedure.
The necessity to account for the response of the electrolyte in evaluating the electrode potential itself has a profound effect on the atomistic modeling of biased interfaces and beyond. A distinct spatial inhomogeneity of the EDL 11 at different regimes challenges the very idea of defining the redox potentials through overall estimates of the free energy for bulk regions of the cell. Likewise, one needs extra measures of caution when computationally assessing the electrochemical stability of an electrolyte. Currently, the AIMD approach suffers from unrealistic configurations of the interfaces due to the limited number of solvent molecules explicitly included in calculations which leads to ambiguity in the definition of the electrode potential, 12 insufficient accounting for the dielectric response of the liquid environment, [13][14][15][16] erroneous assessments of the potential of zero charge (PZC), 15,17,[18][19][20][21][22] the spatial profile of the average electrostatic potential in the EDL and its capacitance. Moreover, in practically realizable unit cells for AIMD even an elementary fluctuation of the ion distribution or polar solvent molecules near the interface, caused by either thermal motion or by the loss of inventory accompanying reduction, causes a drastic change, on the order of ±1.0 V, in the calculated electrode potential. 23 Such large fluctuations of the potential might then in turn trigger an artificial spontaneous cascade of ion reductions that might lead to a run-away regime. To circumvent this problem, one needs to limit stud-Journal of The Electrochemical Society, 164 (11) E3438-E3447 (2017) E3439 ies to those with a large portion of concentrated electrolyte 16 where both fluctuations in the charge profile of solvent and ions are negligible with respect to the thermodynamic average. However, calculations at this scale may not be currently achievable due to computational overheads.
Lastly, the lack of a computationally tractable unified description of the continuous electronic and discrete solute/solvent degrees of freedom constitutes another fundamental problem for the atomistic description of biased interfaces. Customarily, the typical condensedphase AIMD setup, which ideally should be open at either end away from the interface, for electrons at the electrode and for dissolved ions in the electrolyte, is constrained by periodic boundary conditions (PBC) in conjunction with the requirement of global charge neutrality needed to define the potential drop across the interface. A specified concentration is realizable within the atomistic picture only as integer ratios of solute species to solvent molecules. However, in the conventional picture of continuous densities of electronic states in metals, the variation of chemical potential of electrons is realizable via essentially non-integer electron population. Obviously, such a treatment is not possible for discrete ions within a finite simulation. In order to represent realistic situations, the boundary conditions (applied bias) and the molecular content of the AIMD setup must be mutually compatible in terms of system size as well as sampling time. For a given upper bound on system size, only a discrete and limited set of boundary conditions is then possible and vice versa. Therefore, in practice, the AIMD setup for a biased interface must be subjected to additional serious limitations suiting the goals of simulations. For instance, the steady state regime at (or beyond) the redox potential may become inaccessible to AIMD as it would require simulations at constant chemical potential of ions 24 to replenish the lost ones, while the actual goal usually is to study the nascent "just forming" states resulting from irreversible ET.
As a conceptual alternative to the AIMD schemes, continuum approaches represent the electrolyte as a structureless medium and its response to the external polarization reflects a thermodynamically averaged estimate of the ion distribution with respect to distance from the interface that is consistent with externally imposed boundary conditions (e.g. the electrolyte concentration and the potential difference) at the expense of the loss of information about local molecular-scale spatial and dynamic variations. A class of such hybrid continuum models that combine the continuum description of an electrolyte and the atomistic DFT-based treatment of the electrode is based on the generalized Poisson equation (GPE) 25 with a smooth spatially dependent dielectric function ε( r) 23,26-28 and various polarizable continuum models. [29][30][31] These mean-field models are generalized in the framework of classical DFT [32][33][34][35] in terms of free energy functionals [ϕ; ρ + , ρ − ] whose variation with different boundary conditions for the electrostatic potential provides the thermodynamically optimal spatial characteristics of the interface. The accounting for steric effects, 36-41 short-range coulomb correlations, and many-body interactions introduced by the self-energy corrections to the mean-field approximation 42-45 enables going beyond the standard Poisson-Boltzmann approximation (PBA) 46 and explains overscreening and crowding phenomena as well as the "camel-bell" transition in differential capacitance for ionic liquids. 47,48 Depending on a particular implementation, the corresponding ad hoc modifications are introduced either directly into an action-like free energy functional resulting in the Helmholtz-like form of GPE 49,50 or lead to specific changes in the boundary conditions for ϕ( r ). 8,51,52 Despite these advances, these continuum models do not include the effects of specific adsorption, thus predicting unrealistic electrostatic potential profiles known to be dominated by specifically adsorbed ions. More importantly, by virtue of the general theorems of electrostatics, the boundary conditions of the GPE for the electrostatic potential at the electrode, ϕ 0 , and the electrode surface charge density, σ, are uniquely connected, which makes it impossible to account for the effects of ET on the structure of the EDL. The same holds true for any aspect of dynamics at the interface. This limits the scope of continuum models to the classical EDL regime. However, the theory does not provide estimates of its limits in terms of the applied external bias, which may lead to an erroneous picture if the applied potential for unknown system extends beyond the safe potential (or surface charge) range of the EDL, when only reorganization of the electrolyte is assumed. Another problem of continuum models is the incomplete connection to the equilibrium thermodynamics of the metal-liquid interface, whose description is usually done in terms of the unique (for given conditions) redox potential as a starting point. In contrast, for continuum models, the starting point is the state of zero charge of the surface (zero potential with respect to its value in the bulk) characterized by flat ion density profiles which, from a thermodynamic stand point, is a non-equilibrium state defined by the so-called potential of zero charge (PZC). The proper evaluation of the offset between these starting points for the electrode potential requires the introduction of the equilibration mechanism, whether it is due to dipolar screening that induces the accumulation of extra surface charge in the pure solvent or the actual equilibration of the metal electrode with the electrolyte containing its own ions. In the framework of the hybrid continuum models that treat metal electrodes exposed to dielectric media quantum mechanically 47 this problem is solved only heuristically and the ambiguity of mapping the electrode potentials of continuum models onto the thermodynamic scale remains.
To go beyond the classical EDL regime, a combination of ET Tafel-like kinetic equations 9,53 with ion transport kinetics of continuous media 54 is used. Despite some progress in thermodynamically consistent continuum transport theory, 55,56 this approach lacks selfconsistency and does not resolve the effect of ET on the structure of the double layer due to the application of a linearized version of Ohm's law to the steady-state regime for the homogeneous system. At this point, the input from AIMD that provides immediate access to configurational and even chemical/electronic degrees of freedom would be very valuable, but within current limits on length and time scales, such atomistic scale ET events may cause changes which would push the system far from thermodynamic equilibrium (e.g., significantly altering concentration) thus devaluating further estimates of its characteristics from AIMD.
Due to the obvious complementary character of the AIMD approach and continuum models, a self-consistent combination of these two would permit a more realistic description of the electrodeelectrolyte interface. In particular, instead of relying on the spontaneous equilibration of the slow molecular subsystem, from an atomistic computational stand point it would be much more efficient to have it pre-equilibrated by means of a continuum model for each value of the external electrode potential. Likewise, by extending the applicability of a continuum model to the regime when the electrode potential is close to (or beyond) the redox potential would allow the formulation of a robust computational protocol for AIMD simulations that would avoid the pitfalls of having the molecular system far from (quasi-)equilibrium. At the same time, using continuum models with correct boundary conditions for ϕ( r), one can evaluate useful criteria relevant to AIMD, such as the minimum required size of the system and the realistic range of the electrode potential for a given electrolyte concentration.
In pursuit of developing such a self-consistent methodology for modeling biased interfaces, in this work we propose a continuum phenomenological model to describe the effects of ET on the structure of the EDL as well as the effects of specific adsorption of ions. With properly tuned input parameters, this advanced model should permit both design and validation of thermodynamically consistent AIMD simulations and enable exploration of the limits of a microscopic description of biased interfaces.

Theory
We start with the PBA of the free energy functional Journal of The Electrochemical Society, 164 (11) E3438-E3447 (2017) between monovalent ions is incorporated in the energy term and core-core repulsions are neglected. With a self-consistent mean-field electrostatic potential ϕ( r ), the entropy term associated with solvent molecules, 39,57 and a mass conservation term defined by the chemical potentials of positive and negative ions μ ± , 36,37 the functional [ϕ( r ); ρ + , ρ − ], that from now on we will call the Gouy-Chapman with the exclusion volume (GCEV) model, gives rise to a modified Poisson-Boltzmann equation (MPBE) with Neumann or Dirichlet boundary conditions: ∂ϕ Here we propose to generalize this GCEV model by adding terms describing non-electrostatic specific interactions 1 ( r ) and 2 ( r ) of monovalent (z = 1) ions of effective size a with the electrode surface: By definition, 1 ( r ) and 2 ( r ) adsorption potentials, as a part of the free energy functional, correspond to the interactions of ions with the neutral electrode in the infinitely diluted limit. Therefore, these potentials do not include any dependencies on surface coverage or bias in contrast to the adsorption energies used in the phenomenological electrochemistry (e.g. Frumkin isotherm 58,59 ). Indirectly, 1 ( r ) and 2 ( r ) should depend on the bias as it tunes the overlap of the ion frontier orbitals with the electrode surface states. However, since we disregard the electronic structure of the electrode, these quantum effects will be ignored here. However, the proposed model keeps the dependencies of the free energy of adsorption of ions on the bias and surface coverage since it operates with the self-consistent electrostatic potential that takes into account finite size effects and all interactions between ions except non-Coulombic ones. We also ignore repulsion between solvent molecules and thus ignore the layering in the double layer. 60 Here, in contrast to other studies where the dielectric function is explicitly dependent on either the electric field ε(E) 61 or the electronic density of solutes 26 or metal electrodes ε(ρ) 49,50 causing dielectric saturation in regions of high electric field or electronic density, we use a spatially dependent dielectric function ε( r). As has been demonstrated, 62 a smooth "step-like" function ε(z) which decays from its bulk value to the high frequency limit in the vicinity of the electrode (at sufficiently high surface electron density σ) can accurately reproduce a response consistent with Booth's ansatz for ε(E). In this way, we ignore the dependence of ε(z) on the surface charge density, which will, however, be reflected in an underestimation of the differential capacitance in the small regime of electrode potentials around PZC.
The variation δ δϕ = 0 and δ δρ± = 0 leads to the modified Poisson-Boltzmann equation (MPBE) that reads as follows: where n b is the bulk concentration of the ions and a is the characteristic size of the solvated ion. Its solution defines the extreme saddle point of the functional and spatial profiles ρ ± for counterions that become dependent on each other. As both ρ + and ρ − are usually complex functions of coordinates parametrized by the electrode potential, it is more convenient to explore the properties of the model by analyzing the "phase relationship" f (ρ + , ρ − ) = 0. For the Gouy-Chapman model (as a → 0), the phase relation is well known ρ + ρ − = n 2 b . 46,47 As we found, in the GCEV model 36 ρ + and ρ − are dependent through the relation ρ+ρ− Due to the fundamental symmetry of the coulomb interactions, the common feature of this class of models is that f (ρ + , ρ − ) = 0 is independent on the boundary conditions, as it is exemplified in Fig. 1a. The inversion symmetry ρ + ↔ ρ − is guaranteed by the equal size of the ions. The introduction of non-electrostatic interactions 1 ( r ) and 2 ( r ) results in the concomitant modifications of f (ρ + , ρ − ) = 0 which now reflects the interplay between adsorption potentials, electrostatic forces and the effects of volume exclusion. As a consequence, the phase relation becomes dependent on the electrode potential ϕ 0 (see Fig. 1b).
Eq. 2 describes the double layer operating as an ideal nonlinear capacitor. As the electrochemical potential of electrons , at a certain overpotential one would expect non-zero faradaic current and significant changes in the structure of the double layer associated with the "leakage" of the capacitor. However, with a single self-consistent MPBE it is impossible to capture these effects as explained above as a general property of equations of electrostatics.
To resolve this issue, we introduce a simplified kinetic equation which describes the time evolution of the concentration of the reducing species, ρ + , that is characterized by two different time scales. The kinetics of the forward reduction process ρ + → ρ 0 is characterized by τ → whereas the kinetics of the backward oxidation processes ρ 0 → ρ + is associated with τ ← . Here, ρ 0 is the time and spatially dependent concentration of the reduced species. Finally, T D introduces the kinetics related to the migration/diffusion of ions in the double layer. For simplicity's sake, here we consider only the effect of cation reduction. The equation reads as follows: wherec(z) describes the portion of ρ + (z, t) or ρ 0 (z, t) that gets reduced/oxidized at coordinate z. Here we assume a tunneling regime whenc(z) = e − z a 0 , and a 0 specifies the "depth" of the layer where ρ + (z) effectively reduces/oxidizes. The functionc(z) reflects the mechanism of the electrochemical process. ρ + (z, 0) is the initial cation concentration profile corresponding to the ideal capacitor regime (as defined by the solution of the MPBE above). Thus, the first two terms in Eq. 4 describe the competition between the spatially dependent ET and the thermodynamic force bringing new cations from the bulk electrolyte to restore the initial profile that corresponds to the extremum for the free energy functional of Eq. 2. The last term in Eq. 4 describes the oxidation rate of newly reduced ρ 0 (z, t). Here we also neglect the effect of the accumulated ρ + (z, t) and do not consider any mechanism of diffusion from the surface. Note, that the PBA does not take into account any aspects of dynamics. The combination of Eqs. 3 and 4 is not self-consistent, yet provides an approximate way to include the effects of ET on the structure of the EDL.
At steady state, Finally, the ET times, τ → and τ ← depend on the difference between the reduction level eφ red and the electrochemical potentials of electron. We assume that In order to find the self-consistent solution for the electrostatic potential as well as for the corresponding density profiles for anions, we rewrite the free energy functional * [ϕ( r); ρ − (r )] with the fixed cation spatial profile ρ * + (z) whose variation gives a new MPBE: with the same boundary conditions as for Eq. 3. Eqs. 3, 5, 6 and 8 constitute the formalism which describes the effects of ET on the structure of the double layer.

Results and Discussion
Effects of specific adsorption and electron transfer.-Eq. 3 captures several important features. The double-layer capacitance is dominated by the effects of adsorption of ions. 63 In Fig. 2 we show the progression of differential capacitance, C di f = dσ dϕ , for two different adsorption energies of anions with E ads = 3k B T (blue solid line) and E ads = 6k B T (red solid curves) for concentrations n b = 0.1, 0.2, 0.4 and 0.8 M. Here we assume that the electrode is "the sticky surface" for anions with the adsorption energy profile defined by 1 = E ads (1 − e −k(z−z 0 ) ) 2 − E ads where k = 2Å −1 and z 0 = 0, and no specific adsorption for cations 2 = 0. The effective size of ions is set as a = 3Å far enough from the saturation limit ∼10Å at n b = 0.8 M which should guarantee the minimum of C dif , at the PZC in the absence of ion adsorption, dielectric constant ε = 80 and temperature T = 300K . For sufficiently negative electrode potentials ϕ < −0.15 eV, the effects of specific adsorption vanish, and both pictures give the same results. However, at higher adsorption energy of anions, one sees not only the shift of the minima of C dif , but also a peculiar "camel-bell" transition. 37,38,40 Similar trends were experimentally observed for aqueous solutions of simple anions for polycrystalline Ag-electrodes. 63,64 The values of C dif are overestimated. However, Eq. 3 with advanced adsorption and steric repulsion potentials 65,66 may easily produce, without introducing the ad hoc term ε r 0 , the effect of the compact Helmholtz layer with a linear profile of the electrostatic potential resulting in more realistic values of C dif .
Another useful feature of Eq. 3 is that it allows us, in principle, to estimate the Esin-Markov coefficient and the electrosorption valency l B = −z( ∂σ M ∂σ i ) E , where σ M is the charge density on the metal and σ i is the charge density connected with the specifically adsorbed species, 64,67 and to compare that with experimental values of l B . That opens the way to theoretically evaluate the charge transfer coefficient, λ, and refine the model based on experimental results. 68,69 With both the effects of lattice saturation and specific adsorption accounted for, the concentration profiles ρ − (z) and ρ + (z) can be viewed as a generalization of a Frumkin isotherm. The model presented above also naturally includes the interactions with the image charge induced on the electrode surface but does not account for the contribution from the "quantum capacitance". That affects the estimates of the PZC accordingly. Lastly, a smoothed step-like dielectric function ε(z) = 1 + ε−1 1+ex p(−α(z−z 0 )) 26,49 similar to that found in approximate solutions of the Poisson equation with the Booth's-like ε(E), 61,62 where α = 5 and z 0 = 0.5Å in Eq. 3 renormalizes the overall capacitance making it more realistic as compared to that from the GCEV model, as shown in Fig. 2 by green and orange curves corresponding to the anion adsorption energies E ads = 3k B T and 6k B T , respectively.
The analysis of the model that combines the simple ET kinetic Eq. 4 and the self-consistent Eq. 3 allows us to draw the following observations. At μ e φ red (ideal capacitor regime), τ → τ ← but at this electrode potential not enough ρ 0 is created to drive the current in the anodic direction. The second term in Eq. 5 induces asymmetry between ρ + (z) and ρ 0 (z) that leads to non-zero current at μ e (ϕ) = φ red , where τ → = τ ← = τ 0 . Further decrease of the electrode potential promotes a cathodic current by decreasing τ → but at the same  Black dashed curves -differential capacitance of the GCEV model. 36 time it depletes the cation concentration next to the electrode, slowing down reduction. On the other hand, the reduced ρ 0 will be less and less likely oxidized due to the increase of τ ← . This competition of rival factors results in a steady state current which depends on the initial mutual alignment of the reduction level φ red and initial E f . However, at the finite scanning rate of voltammetry, the current may appreciably deviate from the steady magnitude and should exhibit hysteresis.
In Fig. 3 we show the effects of ET on the ion concentration profiles and the electrostatic potential as compared to the regular EDL predicted by the GCEV model. No effects of specific adsorption are included. Different electrochemical conditions are defined by parameter 0 (see the inset of Fig. 3) which specifies the offset between E f of the neutral electrode and the reduction potential φ red and thus its electrochemical activity. As τ 0 and T D are set as 1 and 100, respectively, τ → = τ ← = T D at the potential corresponding to 0 = 0.115eV. At the potentials φ red − μ e ≤ 0.115 eV, τ → ≤ T D and the EDL behavior starts deviating from that of the ideal capacitor. A peculiar non-monotonic concentration profile of anions (Fig. 2c) is a result of the competition between attraction to cations, repulsion from the negatively charged electrode and steric factors which preclude high local concentration of both types of ions at the same distance. If the backward oxidation of ρ 0 is suppressed (τ ← τ → at μ e = φ red ), the model predicts the nearly full depletion of ρ + at z = 0 which at small 0 may lead to the depolarization of the electrode surface (σ > 0 at ϕ < 0). The model also predicts the regimes at which the electrode surface is covered by anions even at ϕ < 0.
The effects of ET impact the differential capacitance of the double layer. Though the differential capacitance is strictly speaking well defined only in the range of potentials within the classical double layer, here we show the "renormalization" of surface charge density in the presence of non-zero current. In Fig. 4 we show the corresponding modifications of C dif as the electrode potential approaches φ red . Focusing first on the effect of varying solution concentration (blue curves) at constant offset 0 = 0.22 and at φ red − μ e ≈ k B T (below −0.15 V ), we notice C dif drops dramatically (compare with Fig. 1 or the dashed curves). Next, as the electrochemical activity, 0 , of the electrode increases progressively (red curves), we see changes in C dif in the same voltage range. The inclusion of the effects of specific adsorption and solvent polarization mixed with ET will have a profound effect on the EDL structure affecting the PZC. In principle, the parameters for the presented model can be extracted from EIS experiments mapping the characteristics of ET resistance and ion diffusion with corresponding characteristic time scales T D , τ → and τ ← .

AIMD simulations of the biased interface: fundamental limits
and implications from continuum models.-Thermodynamic consistency with interfacial electrochemical boundary conditions requires that AIMD simulations meet several interconnected requirements. Let us make some common assumptions on typical AIMD calculations relevant to condensed phase interfacial studies of half-cells: (1) We employ periodic boundary conditions (PBCs), where the semi-infinite solid electrode surface is modeled as a sufficiently thick slab of atoms extended periodically in the plane of the surface -typically less than 10 atomic layers may suffice for most simple metals. These PBCs define two relevant characteristics of the associated orthorhombic supercell -its overall length in the direction of the surface normal, l, and its area in the plane of the surface, S; (2) The liquid/electrolyte region is defined by another "slab" sufficiently thick to preserve a particular volume with bulk behavior (material density, molecular diffusivity, average concentration and spatial correlation of dissolved species, etc.) at some distance from the electrode surface; (3) overall the length of the AIMD simulation cell is sufficient to also include interfacial phenomena different from the bulk behavior in both the electrode and electrolyte, for example, to encompass the electric double layer. If we assume the electrode to be metallic, then these interfacial differences will be limited (excluding the outermost layer of metal atoms) to the electrolyte -this would not be so for any insulating or semiconducting layers formed on the electrode surface, perhaps as models of a solid electrolyte interphase. So, the volume of the AIMD supercell can be obtained as the product S (l solid +l interface +l bulk ). Within this volume we In all that follows, we will refer to one specific interface as the active electrochemical interface. Clearly, a slab of finite thickness has two sides that meet the electrolyte under PBCs, but we assume that at finite bias, inversion symmetry is broken and only one surface becomes electrified. This can be achieved by the initial arrangement of the ion distribution, which should mimic the profile determined from the continuum model. Specifically, within the "bulk" electrolyte volume, S l bulk , there are an equal number of cations and anion, while the mismatch defining the number of surface electrons, n e , is limited to the interfacial region, Sl interface . We assume that thermal fluctuations are insufficient to drive the system out of this metastable state toward the lower energy state with the excess charge split evenly on either side of the slab. By virtue of the necessarily discrete nature of the atomistic description of the ions, the charge neutrality requirement limits n e to integer values and sets some minimum limits on the surface change density at the electrode (described below). In the absence of an explicit reference electrode, an internal reference for the half-cell AIMD setup is provided by the calculated, self-consistent, thermodynamically-averaged, electrostatic potential in the bulk region of the explicitly included electrolyte. Vital characteristics of the supercell are: its ability to screen electric fields, such as that generated by the electronic surface charge density, within a certain length l scr and a particular time scale. At thermal equilibrium, material specifics define these length and time scales as the Debye length, δ D , and the slowest relaxation time τ x 70 (e.g. the inverse rotational diffusion coefficient), respectively. Under an external bias, the Gouy-Chapman length l GC 71 dictates the variation of the distributions of ions and defines the number of ionic species within a layer next to the electrode as compared to its value in the bulk. Thus l GC via the link with the surface charge density limits the minimum electrolyte concentration and the electrode surface area, S, connecting the efficiency of the time averaging over the simulation length and the "static" (or relaxed) picture given by continuum models. Though the applicability of these lengths and time scales is, strictly speaking, limited by the idealized diluted electrolytes (ionic strength below 0.2 M, electrode potentials below 80 mV 72 ), they can still be used to estimate the critical parameters of AIMD setups.
The combination of the aforementioned conditions leads to the following limitations for the AIMD simulations of biased interfaces (here we consider the charging of the electrode with electrons).
1. Due to the discrete nature of the ionic charge imbalance, both the surface charge density, σ, and the electrode potential, ϕ, are quantized values with the minimum "quanta" σ = ± 1 / S and ϕ = L( σ) defined by the capacitance of the EDL. For the AIMD setups with a fixed number of ions, the maximum imbalance n e ≤ N T − − 1, where N T − is the total number of anions defined by the bulk electrolyte concentration, n b , and the size of the interfacial region. 2. The surface area has to be within S min ≤ S ≤ S max range. S max is limited by the computational feasibility of the simulation (maximum number of solvent molecules). For aqueous solution, S max [Å 2 ] ≤ N water

L[Å]
, where L is the length of the solution part of the cell. The lowest limit for S min is dictated by the minimum electrode potential ϕ min (in electrochemical scale) that has to be in the safe EDL region at least several k B T away from reduction potential φ red for the minimum ionic imbalance (n e = 1). However, for an AIMD practitioner S min may rather be limited by other factors such as incompatibility of the bulk electrolyte concentration with the target bias (see below). 3. In either conditions (constant surface charge density 73 or constant chemical potential of electrons 74 ), the size of the cell must be such that thermal fluctuations of polar solvent molecules and ionic species do not exceed a certain limit. As has been shown, 73 for a pure polar solvent, the probability distribution of the electrode potential derived from Marcus theory for linear solvent polarization has a Gaussian shape with a standard deviation independent of the actual electrode potential. In the realistic AIMD setup with water (L = 10Å and S = 100Å 2 , 32 water molecules), the standard deviation of the electrode potential ϕ 2 = 0.51 V , which is much larger than the thermal energy. Clearly, for such a small AIMD setup, correspondingly large fluctuations might easily induce electron transfer events across the solid-liquid interface with little correspondence to reality. Explicit inclusion of ionic species 8,9 may introduce even more serious problems, since the interfacial fluctuations of discrete ion distributions may induce even larger fluctuations of the electrode potential or may shift the energy of a LUMO orbital for a molecular species in the electrolyte into resonance with the electrode Fermi level. The rigorous estimations of the effects of fluctuations is a non-trivial problem that requires one to use not only the fluctuation-dissipation theorem 70 which relates the differential capacitance and the variance of the total charge q 2 = C di f k B T or the electrode potential ϕ 2 = k B T C di f but some other kinetic characteristics of the solution, e.g. diffusion constants, viscosity, etc. To the best of our knowledge the general expression for the spatially dependent standard deviation of ion densities ρ 2 ± (z; σ) is unknown for the advanced continuum models and will be a subject of future studies. Note that the necessary constants in the ET model described above may be determined from consistent AIMD simulations, especially the position in energy of available electronic levels in the electrolyte and estimates of the time constants for charge transfer based on constrained electronic structure calculations. 75 4. To properly describe the interface, the cell has to comprise the interfacial region characterized by l scr (screening length) and the bulk region whose dimensions allow for the proper description of ion-ion correlations at finite concentration ( √ S ≥ 2l c , l bulk ≥ 2l c , where l c is the ion-ion correlation length and l bulk is the length of the bulk electrolyte region, defined already). With a fixed number of ions in the simulation cell, the bulk region will effectively work as a buffer to damp oscillation in the ion profile at the interface which is in fact what may define its dimension. To withstand a certain potential drop across the interface, there will be corresponding ion density profiles with a spatial ionic imbalance. For a monovalent binary aqueous electrolyte (ε = 80, T = 300 K , n b = 1M, z = 1, δ D ≈ 3Å) contained in the simulation cell with a size L = 10Å and S = 150Å 2 , the minimum surface charge density σ min = 10.7 μC/ cm 2 . However, if the concentration is n b = 1M, this unit cell contains only one ion pair which along with the estimates of l GC = 3.4Å at σ min = 10.7 μC/ cm 2 suggests that this AIMD setup is a very rough representation of the physical system where the screening of the surface charge may be obtained only at the expense of very long simulation time. Furthermore, the strong implicit correlation between charged species in neighboring periodic replicas of the supercell will generate unphysical long-ranged electrostatic errors that may also prevent fair comparisons with measurements. This example also shows that for a certain concentration the surface area S may not be too small as it would require an excessive length of the simulation cell L to accommodate the appropriate number of ions to support the interfacial ionic imbalance which is again tantamount to a very long simulation time.
Here, we propose to use our continuum model to assess the critical limits for the AIMD setup to meet the requirements of thermodynamic stability from the standpoint of the bias. Specifically, to determine practicable AIMD simulations, we focus on the size of the AIMD setup in terms of number of atoms, which, for condensed phase systems, is directly linked to spatial dimensions via material density. For the moment, we omit discussion of the necessary time scales for converged thermodynamic averaging, with the expectation that this is merely an overall multiplicative factor in determining computational cost and tractability.
For comparison, we use both the GCEV model and our model that includes the effects of the specific adsorption of anions and spatial dependence of the dielectric function. Other effects, such as ET can be included in a similar way. The input parameters (bulk concentration, dielectric function and surface charge density) are used to obtain the information about the self-consistent electrode potential and the ion density distributions that are further used as an input for the AIMD simulations. In this way, we can try to reproduce the preequilibrated state compatible with the boundary conditions and skip the slow equilibration of the molecular subsystem. The AIMD is the then used to obtain anew the ion profiles as an updated input for the continuum model. In Fig. 5 we show the conceptual scheme of input/output parameters.
To assess the screening length, l scr , we require that at each induced surface charge density, σ, the deviation of the average electrostatic potential at that distance from the overall potential drop meets the condition | ϕ(lscr ) ϕ | ≤ 0.02 which guarantees the approximate neutrality of the interfacial region. To define the length of the bulk region, l b , we require it to be long enough and contain enough ions to support the ionic imbalance and preserve approximately the bulk solution concentration in the overall cell: N T where N T ± -total number of ions of each sign in the cell L = l scr + l b , N c ± = S ∫ lscr 0 ρ c ± (z)dz -is the number of ions in the interfacial region predicted by the continuum model. For a given surface area S we evaluate the length of the cell L as a function of the concentration for each value of the charge imbalance n e and thus evaluate the minimum size of the cell for each value of the electrode potential. For the model that includes the effect of specific adsorption of anions and spatial variance of dielectric function we use the parameters described below.
1+ex p(−α(z−z 0 )) , where α = 5, and z 0 = 0.5Å. The effective size of ions is set as a = 4Å, dielectric constant ε = 80 and temperature T = 300K . Empirically we found that l scr ∼ 1 √ n b and has a very weak dependence on the boundary conditions similarly to what one would expect from the Debye screening.
In Fig. 6a we show the results obtained with the GCEV model for the surface charge density 8 μC cm 2 ≤ σ ≤ 40 μC cm 2 keeping the overall neutrality of the cell. The estimates show that for concentrations n b ≤ 0.5 M AIMD simulations would require excessively long cells that for the higher charge imbalance (e.g. n e = 5) can easily exceed computational feasibility (150Å → 1000 water molecules, black dashed line). Alternatively, these limitations can be seen in terms of maximum applied potential (see Fig. 5b) that cell can support. For concentration 0.5 M, 1 M and 2 M, ϕ max is 0.12 V, 0.2 V and 0.38 V , respectively. As long as the GCEV continuum model captures well the electrolyte response, these estimates seriously challenge the attempts to simulate atomistically thermodynamically stable interfaces. To meet the stability requirements, the short AIMD cells (≤ 20Å) can be used to simulate interfaces only at the potentials smaller or   comparable to kT = 0.025 V which may not be relevant for any practical applications.
The consideration of more realistic continuum models helps to alleviate these severe size limitations. In Figs. 6c, 6d we show the estimates of the cell sizes when specific adsorption of ions and spatial variance of dielectric function are accounted for. The strong specific adsorption (E ads = 0.5 eV ) of anions reduces the variation of cell lengths within the same range of the surface charge density (see Fig.  5c) since at these conditions the capacity is dominated by the adsorption. Therefore, the biggest effect of it is expected to be well seen at lower potentials as it is shown in Fig. 5d. Indeed, for concentration 0.5 M, 1 M and 2 M, ϕ max is 0.21 V, 0.42 V and 0.65 V , respectively. This result indicates that AIMD simulations of stable interfaces at ϕ ∼ 0.5 V can be practicable at realistic concentrations as they naturally include such effects.
The attempts to evaluate the size limits for interfaces at higher biases encounter serious problems inherent to continuum models per se. In the absence of ion-ion correlations, this simplified description predicts the over-packing of cations and the full depletion of anions in the interfacial region. The problems become especially evident at lower concentrations at which N c − < 1 in the large region (l scr ∼ 15 − 20Å). At this point, AIMD outputs should serve to refine the continuum model in a self-consistent fashion to construct eventually more realistic interfaces.

Conclusions
In the presented work we develop the model that combines the classical continuum theory of the electric double layer with kinetic equations to incorporate the effects of specific adsorption and electron transfer within the formalism based on the free energy functional. Specifically, we achieved the following results.
1) Introduction of non-electrostatic (e.g., specific adsorption) interactions 1 and 2 of ions with the electrode permits treating the system at the level of the differential modified Poisson-Boltzmann equation and leads to significant modifications of the "phase relationship" f(ρ + , ρ − ) = 0 and the differential capacitance of the double layer as compared to the GCEV model. 2) Accounting for the non-linear solvent polarization effect results in the appreciable renormalization of differential capacitance.
These advances permit direct comparison of our model to the experimental measurements of the electrode surface coverage by specifically adsorbed ions (e.g. Hurwitz-Parson method). This in turn enables an access to the details of ion adsorption potential.
3) We show that the description of the electric double layer at electrode potentials close to the redox potentials can be efficiently performed with a formalism based on a pair of coupled modified Poisson-Boltzmann equations. This scheme predicts the peculiar structure of the double layer in the presence of non-zero faradaic current and significant changes in the differential capacitance at potentials near the redox levels. Our model is equally suitable for either const − ϕ or const − σ boundary conditions.
When the formalism is used with the parameters obtained from ab initio calculations and/or relevant experiments, e.g. EIS, it allows us to merge the continuum models self-consistently with AIMD. Using continuum models, we perform an analysis of the thermodynamic consistency and the limiting characteristics of an AIMD setup in terms of the cell sizes and the applied surface electron density/electrode potential that leads us to the following results. 4) We provide estimates of AIMD cell sizes as functions of applied potentials that meet the requirements of thermodynamic stability of biased interfaces and compatibility of the content of the electrolyte (ion concentration, number of solvent molecules) with the target conditions. 5) These estimates indicate that for realistic electrolyte concentrations and currently used AIMD cell sizes, the validity of ab initio approaches is limited to very small electrode potentials (ϕ max < 0.05 V). Simulations at the higher biases would imply strongly non-equilibrium configurations of the solution and incompatibility with the targeted external conditions. 6) Estimates of the critical parameters for AIMD simulations based on continuum models that include the effects of specific adsorption and solvent polarization show that simulations of practical voltages (ϕ ∼ 0.1 V) are within reach of atomistic simulations but require special computational protocols to avoid time-consuming pre-equilibration stages.
We believe that our findings pave the way toward a self-consistent methodology for modelling biased interfaces with the continuum level description as an indispensable part.