Electrosorption at metal surfaces from first principles

Electrosorption of solvated species at metal electrodes is a most fundamental class of processes in interfacial electrochemistry. Here, we use its sensitive dependence on the electric double layer to assess the performance of ab initio thermodynamics approaches increasingly used for the first-principles description of electrocatalysis. We show analytically that computational hydrogen electrode calculations at zero net-charge can be understood as a first-order approximation to a fully grand canonical approach. Notably, higher-order terms in the applied potential caused by the charging of the double layer include contributions from adsorbate-induced changes in the work function and in the interfacial capacitance. These contributions are essential to yield prominent electrochemical phenomena such as non-Nernstian shifts of electrosorption peaks and non-integer electrosorption valencies. We illustrate this by calculating peak shifts for H on Pt electrodes and electrosorption valencies of halide ions on Ag electrodes, obtaining qualitative agreement with experimental data already when considering only second order terms. The results demonstrate the agreement between classical electrochemistry concepts and a first-principles fully grand canonical description of electrified interfaces and shed new light on the widespread computational hydrogen electrode approach.


INTRODUCTION
In recent years, calculations based on ab initio thermodynamics have increasingly contributed to unraveling key processes in interfacial electrochemistry; e.g., in batteries, fuel cells, and other electrocatalytic systems. In such studies, the electrochemical environment and operation conditions are suitably represented in the form of thermodynamic reservoirs [1][2][3][4][5] . Within a grandcanonical setup these reservoirs are then coupled to predictivequality (typically density-functional theory (DFT)) total energy calculations for the electrode, to deduce electrochemical stabilities and activities; see e.g., refs. 2,[6][7][8][9][10][11] . In the application to interfacial electrochemistry a key challenge to this general ab initio thermodynamics concept is the necessity to exchange electrons with a reservoir representing the electrode potential. In principle, this requires to perform DFT calculations in various charge states. This clashes with the common representation of the electrode as finite slab using periodic-boundary conditions, for which straightforward calculations can only be performed at zero total charge of the cell. For this reason, in most practical application, the early computational hydrogen electrode (CHE) approach 2,6,12 relies solely on the energetics of charge-neutral electrode calculations, in the absence of the electrochemical double layer. The dependence on the electrode potential is then included in the analysis as an a posteriori shift of the electrochemical potential of the electrons taken from the reservoir, whose number is a priori fixed according to the charge-neutrality condition.
Conceptually, this restriction can be readily lifted through the introduction of meaningful counter-charge models in the vacuum region between the periodic slabs 1,13 . This then allows for a surface polarization and charging of the electrode slab while still maintaining the overall charge neutrality of the supercell. In the recent fully grand-canonical (FGC) approaches this is realized in practice by the use of an implicit solvation model in the DFT calculations 5,9,11, , which captures the effects of both electrolyte polarization and solvent screening, and renders the calculations largely insensitive to details of the counter-charge model employed 35 at conventional thicknesses of the vacuum/ dielectric region 32 . Within this approach, variable-charge calculations can nowadays be performed at high computational efficiency in a number of major DFT codes, e.g., ENVIRON 18,36 , VASPsol 21 , FHI-aims 25 , Jaguar 29 , BigDFT 24 , GPAW 31,37 , or JDFTx 38 .
The capabilities of FGC approaches have been highlighted in several publications and include a potential-dependence of chemical reaction steps 8,[29][30][31]39,40 , potential-induced surface reconstructions or shifts of stable adsorption sites 5,11,13,41 . However, there is presently still limited understanding of the methodological differences between the still prevalent zerocharge CHE and the variable-charge FGC approaches. In this work we therefore present a consistent thermodynamic framework to compare both types of approaches and highlight the simplifications made in practical calculations. In particular, we show by analysis of a generic second-order Taylor expansion of the interface energies with respect to surface charge 7,11,20,34,39,[42][43][44][45] that the CHE approach can be understood as a first-order approximation to the FGC energetics, while to second order, double-layer charging is represented by changes in work function and interfacial capacitance.
We illustrate these differences between CHE, second-order expansion and FGC in the application to two classic probes of interfacial electrochemistry at metal electrodes, namely non-Nernstian shifts of proton electrosorption peaks at low-index Pt crystals 5,[46][47][48][49][50] and the non-integer electrosorption valency of halide ions at Ag(111) [51][52][53][54] . While neither of these observations can be reproduced at the CHE level, already the second-order expansion introduces them qualitatively, if not even in quantitative agreement with existing experimental data. We find that with current simple implicit solvation models higher-order terms captured in the FGC approach do not generally lead to further significant improvements. Apart from the intrinsic shortcomings of the continuum solvation itself 55 , we also discuss alternative sources for the remaining discrepancies with experiments. As an important corollary, the present work also demonstrates analytically that the FGC energetics does not tend toward CHE results in the limit of single adsorbates in infinitely large surface unit-cells. Instead, the discussed applications suggest that second-order calculations are a computationally efficient approximation to FGC calculations.

RESULTS
Ab initio thermodynamics and interfacial free energy The central quantity in our ab initio thermodynamics approach to surface electrochemistry is the Gibbs excess energy of interface configuration α 4,5,13 , The specific interface configuration α is hereby characterized by its surface geometry (e.g., the position of adsorbates) and chemical composition. In this work we will focus on adsorbates at metal electrodes and we will specify this composition in terms of the number of neutral substrate atoms N α s , the number of (possibly charged) adsorbate species N α a , and the number of excess electrons N abs;α e (conveniently measured as absolute surface charge, cf. below). The additional dependence on temperature T and pressure p in the grand-canonical setup is henceforth dropped for ease of notation. Appropriately normalized to the surface area A, the intensive interface excess free energy is and denotes the cost of creating the electrode surface α in the electrochemical environment at the applied potential Φ E , as it compares the Gibbs free energy of the electrode G α surf with the sum of the free energies of its constituents in their corresponding (bulk) reservoirs. These free energies are defined by the bulk chemical potential μ s of the (neutral) substrate atoms and the electrochemical potentials (highlighted with a tilde) of the (charged) adsorbate speciesμ a and of the excess electronsμ e . The latter relates to the applied electrode potential viaμ e ¼ ÀeΦ E . Without loss of generality, we choose the extensive Gibbs excess energies for the entire electrode as central base quantities in this work, and will therefore explicitly include a surface area dependence in the equations wherever required.
Note that the sign convention followed in Eq. (1) considers negative Gibbs free energies and (electro)chemical potentials, with more negative values corresponding to more stable or more dilute configurations. The electrode potential Φ E is measured according to electrochemistry conventions, with increasing values away from the zero-reference vacuum level, such that e.g., the experimental standard hydrogen electrode (SHE) lies at +4.44 V on this absolute scale 56 . In multicomponent systems different neutral and charged species may exist and appropriate sums over the indices s and a with appropriate stoichiometries need to be introduced in Eq. (1). We stress that this need extends also generally to all interfacial species in the double layer, such as water or solvated ions [57][58][59] . However, in this work we restrict ourselves to implicit solvation models in which only the specifically bound ionic adsorbates a are treated explicitly in Eq. (1).
It is important to point out that the so-defined excess energies have often no direct thermodynamic relevance, as they still depend on the specific interface configuration α. The true thermodynamic excess energy G exc corresponds to the minimum G α exc in the space of all possible interface configurations α 5 , and this then also defines the true ground-state electrode configuration in terms of N s , N a , and N abs e . We note in passing that the two mathematical steps of formulating the excess energies G α exc as in Eq. (1) and then minimizing with respect to configurations α correspond to a Legendre transform, which is why G exc approximates in fact the grand canonical free energy of the system. In a comprehensive all-explicit simulation the optimum interface configuration that defines G exc would eventually emerge directly from an appropriate, ideally exhaustive sampling of all corresponding degrees of freedom. In practice, G exc is instead often approximately determined by computing G α exc for a set of candidate configurations α with fixed N α s ; N α a ; N abs;α e , and the configuration with minimal G α exc is declared to be the most stable one [3][4][5]60,61 . It is also possible to compare candidate configurations and excess energies that already minimize G α exc in a sub-space of all interface configurations α. In the electrochemical context, this refers notably to charge-equilibrated excess energies G α exc that minimize G α exc with respect to the number of electronic charges N abs;α e for a configuration α with defined N α s and N α G α exc then defines the cost of creating an interface with composition N α s and N α a at a given applied potential Φ E , as would be obtained equivalently in a constant potential calculation 27,30,62 .
As common in ab initio thermodynamics approaches to surface systems 4 , the difference in the solid-state terms in Eq. (1), i.e., G α surf and any bulk-like reservoir μ s , is generally approximated as Here, E DFT;α surf is the total energy of candidate structure α and E DFT bulk;s is the total energy per atom of the bulk-like reservoir of neutral species s, both total energies being typically calculated with DFT. ΔF corr,α is a correction term due to the changes in the vibrational and configurational degrees of freedom of the adsorbates in candidate structure α as compared to their corresponding reservoir. This free energy term is commonly assumed to be independent of the charge state N abs;α e of the candidate structure. One should note that even though E DFT;α surf is denoted here as a total energy, it does contain free energy contributions due to solvent screening and electrolyte response when an implicit solvation model is employed in the DFT calculation 5,32 . is thus the number of electrons that need to be exchanged with the external electron reservoir. In the present application to adsorbates at metal electrodes i.e., ÀeN abs;α e is the sum of charges to compensate for the ion charges q a N α a of all adsorbates of charge q a and the net electronic surface excess charge ÀeN net;α e that is itself compensated exactly by the electrolyte counter charges in the diffuse layer. Without meaningful counter-charge model, any net charge of a metallic surface poses a significant challenge in calculations with periodic boundary conditions 63 . In its practical realization, without explicitly represented electrolyte ions (c.f. ref. 64 ) the CHE 2,6,61 approach therefore considers only charge-neutral candidate structures with N net;α e ¼ 0. In order to do so, it thus assumes that each ion of charge q a also drags q a /e electrons onto the surface upon adsorption. Any electrochemical reaction involving a proton transfer would then, for instance, necessarily become a protoncoupled electron transfer (PCET). For such charge-neutral N.G. Hörmann et al. structures, solvation effects, at least on the level of implicit solvation models, are often small. This is why corresponding CHE work is often also based on DFT calculations performed without any solvation treatment at all.
FGC simulations in implicit solvation models instead bypass the overall supercell charge neutrality restriction by balancing any N net;α e ≠ 0 through added electrolyte counter charges in the dielectric region between the periodically repeating slabs. At sufficiently high electrolyte concentrations, the details of the electrolyte model are thereby largely irrelevant for the overall energetics 5,19,32,35,39 . Practical calculations then perform the minimization of G α exc in Eq. (1) with respect to N abs;α e for every candidate structure by explicitly calculating a number of finite charge states. Within the approximation in Eq. (4) and assuming a symmetric slab setup the minimization of Eq. (3) yields the condition where the underbraced equality arises from Janak's theorem and the employed reference for the electrostatic potential (see Methods). Equation (6) thus states that the minimization is equivalent to choosing a system such that its work function eΦ α is equal to the externally applied (target) potential Φ E times the elementary charge e. Thus, the charge-equilibrated G α exc can be directly obtained by interpolating the dependence of the energy and the charge on the work function and evaluating G α exc for charge values where the work function corresponds to Φ E 5 , similar to what is done in the generalized CHE and related constantpotential approaches 27,30,62 .
Obviously, this procedure increases the computational burden, as each candidate structure has to be calculated in different charge states. On the other hand, certain electrochemical observables can be simulated that would otherwise be inaccessible within the CHE approximation. Noteworthy, these include non-linear variations of the interface energy with applied potential (cf. Lippmann equation 5,47,65-67 ) potential-induced surface reconstructions 5,13,41 that maintain the overall stoichiometry in terms of N a and N s , and pH-shifts of electrosorption peaks on the RHE scale as well as non-integer electrosorption valencies 51,52,54 , as will be explicitly illustrated here.
Grand canonical energetics within a quadratic approximation Further insight into the differences between FGC and CHE results can be gained by analyzing a second-order Taylor expansion of E DFT;α surf with respect to the net electronic surface excess charge N net;α |fflfflfflfflfflfflfflfflffl ffl{zfflfflfflfflfflfflfflfflffl ffl} where E DFT;α surf;0 is the total energy of the uncharged system, Φ α 0 the work function and c α 0 the extensive interfacial capacitance at the potential of zero charge (PZC). c α 0 relates to the more common, area-normalized capacitance C α 0 measured in μF/cm 2 via Here and in the following, a subscript 0 will refer to quantities determined at the PZC. The explicit dependence of the total and excess energies on N α s and N α a is dropped for ease of notation.
Within this second-order expansion, the minimization of the excess free energy of candidate structure α with respect to the number of electrons can be carried out analytically following Eq. (6) to obtain Inserting Eqs. (7) and (8) into Eq. (1) then yields the chargeequilibrated excess energy G α exc up to second-order as : The first bracket in Eq. (10) is identical to the excess free energy G α exc;CHE ðΦ E Þ as it would result within the CHE approach. The second term G α exc;DL corresponds to the energy cost of charging the double layer.
This derivation highlights that the CHE approximation is equivalent to a first-order approximation in Φ E of the FGC energetics. For each candidate structure α the CHE energy expression is therefore accurate around the PZC of α where the term G α exc;DL becomes small. On the other hand, for potentials much different than the PZC can become relevant. Within the CHE model, the only potentialdependence of G α exc ðΦ E Þ arises from the term N α a ðq a Φ E Àμ a Þ. This dependence is thus always linear in Φ E , at variance with the typical parabolic potential dependence of electrochemical interface energies 5,47,68 .
From the previous derivation we furthermore see that the relative stability of two electrode configurations α, β differing only in the number of adsorbates a and at applied potential Φ E is and ΔG αÀβ exc;DL ðΦ E Þ ¼ À Within the CHE approximation the relative stability of systems with the same number of adsorbates N α a ¼ N β a is thus potentialindependent, excluding the possibility to simulate potentialdependent diffusion barriers, potential-induced shifts of stable adsorption sites, or potential-induced surface reconstructions that conserve the number of adsorbates (and generally also substrate atoms) at the electrode, as e.g., observed for the quasi hexagonal reconstruction of Au(100) 5 . Intriguingly, all these shortcomings are lifted when including the next higher order term in Φ E , namely the capacitive charging contribution of the double layer G αÀβ exc;DL -an energy contribution that is widespread in the classical electrochemical literature [51][52][53][68][69][70][71][72][73][74][75][76][77][78] . We will explicitly demonstrate below how this contribution leads to a shift of the electrosorption peaks with pH, and how it relates to the so-called electrosorption valency, another classic concept in phenomenological electrochemistry [51][52][53] .
The double-layer correction ΔG αÀβ exc;DL ðΦ E Þ can be further simplified by a change of variables according to In this notation, ΔG αÀβ exc;DL becomes very suggestive. The first term corresponds to the (average) applied potential drop across the electric double layer ðΦ E À Φ 0 Þ times the charge difference Δq net ¼ c β 0 ΔΦ 0 necessary to bring both systems α and β to an equal work function. We denote it with the superscript net , as the term corresponds to the difference in net surface charges (8) of both systems. As shall be shown below it can also be directly related to the adsorbate dipole 79 and the electrosorption valency. The second term results from a change in the capacitive energy due to different interfacial capacitances Δc 0 of the two electrode configurations. Note that identical interface energy contributions were already proposed by Frumkin 68 to understand electrosorption of charge-neutral species.
pH shifts of proton electrosorption peaks at platinum Electrosorption peaks in cyclic voltammetry (CV) measurements are known to be sensitive probes of the electric double layer 5,46,47 . They therefore form ideal observables to assess critically the capabilities and quantitative performance of the CHE and the FGC approaches. Here, we concentrate in particular on proton electrosorption at two low-index single crystal surfaces of platinum, Pt(111) and Pt (100). For both systems, small but significant shifts of the electrosorption peaks with pH even on the RHE scale have been observed experimentally [48][49][50] . Such shifts go beyond what would simply be predicted by the Nernst equation and are certainly an intriguing object of study.
Within ab initio thermodynamics we approximate the electrosorption peak position 5 where the first-order CHE term is and the double-layer correction is where λ 1 ðpHÞ ¼ 4:44V À Φ 0 À kBT e ln ð10ÞpH and λ 2 ðpHÞ ¼ λ 1 ðpHÞ À ΔΦ0 2 : (18) Here, we have specified the general variables in Eqs. (11) and (14) for proton electrosorption as q a = +e andμ a ¼μ H , and have introduced the charge difference between the clean and Hcovered surface Δq net ¼ c clean Þ of the work functions of the H-covered and clean surfaces, respectively.
Within the CHE approximation (ΔG HÀclean exc;DL ¼ 0), the condition is readily resolved for U sorp;CHE RHE by exploiting the electrochemical potential of a proton, which is given by the SHE conditions as 5 where μ(H 2 (g)) is the chemical potential of hydrogen gas at normal conditions (298 K, 1 bar). This then yields where in the second step we have followed the general approximation of Eq. (4) and introduced the average H adsorption energy E H ads plus the (negligibly small) free energy correction ΔF corr,H due to the change of vibrational and configurational degrees of freedom of the adsorbed hydrogen. The CHE approximation thus essentially predicts the proton electrosorption peak on the RHE scale simply at the potential corresponding to minus the average hydrogen adsorption energy. In fact, this holds in general, i.e., the CHE approximation would always equate electrosorption and adsorption for any kind of adsorbate. This is why in the CHE literature electrosorption and adsorption are often used interchangeably.
Most importantly, Eq. (20) does not depend anymore on pH, i.e., within the CHE approximation the proton electrosorption peak does not exhibit any shifts on the RHE scale. If we were to transform U sorp;CHE RHE to an absolute scale, the peak would thus shift simply as predicted by the Nernst equation, at variance with the experimental observations (cf. further below). To this end, we acknowledge that lateral interactions among the adsorbed H atoms and configurational entropy terms might introduce a variation of the average adsorption energy and free energy correction term with surface coverage 9 . However, this will rather lead to a broadened electrosorption peak than natively introduce a pH-dependence in U sorp;CHE RHE . Similarly, different interface models in the underlying DFT calculations, e.g., the surface in vacuum, in various implicit solvation models or with differing adsorbed icelike water layers, will primarily affect only the absolute value of the calculated adsorption energy and therewith the predicted pHindependent electrosorption peak position. The only way non-Nernstian shifts can be rationalized within the CHE model is by a pH-induced change in the adsorption energy, assuming e.g., a pH dependence of the interfacial water structure, or secondary effects such as the adsorption of co-ions, as further discussed below.
In contrast, already the second-order double-layer correction ΔG HÀclean exc;DL ðU sorp RHE Þ introduces intrinsic pH-dependencies, and leads to a quadratic equation defining U sorp;2nd RHE (cf. Eqs. (16)-(18)). This quadratic equation and even more so the higher-order terms contained in the FGC approach yield therefore complex non-Nernstian shifts of electrosorption peaks on the RHE scale. As an example, we point out that already the linear term in ΔG HÀclean exc;DL ðU sorp RHE Þ (∝Δq net ) effectively renormalizes the CHE peak position, while introducing an additional, linear dependence on pH as where we have introduced the per-adsorbate quantity ΔQ net ¼ Δq net NH . ΔQ net is largely independent of the surface coverage, given that the extensive quantity Δq net is expected to scale approximately linearly with N H (/ ΔΦ α 0 ðN H Þ, see e.g., Supplementary Fig. 10). Thus, whenever ΔQ net /e is not vanishingly small, the corrections and linear pH shift of Eq. (21) become significant. N.G. Hörmann et al. Figure 1 illustrates the performance of the various approaches, by collecting the calculated proton electrosorption peak positions on Pt(111) and Pt (100). Concentrating first on the data obtained with an electrode surface model containing 1/4 ML hydrogen coverage in the energetically most stable fcc (Pt(111)) and bridge sites (Pt(100)), respectively, we see qualitatively the same trends among the theoretical approaches on both surfaces. In contrast to the constant U sorp;CHE RHE predicted by the CHE approximation, already the second-order approximation yields non-Nernstian shifts of U sorp;2nd RHE at both surfaces. We note that a significant contribution to these shifts arises from the second, quadratic term in Eq. (18), i.e., from the large change in the interfacial capacitance ΔC 0 upon proton electrosorption. In fact, for Pt(100), this contribution is large enough that it actually changes the shift direction from an upwards shift to higher U sorp RHE from the linear Δq net term only (cf. Eqs. (18) and (21) and the signs of the work function change for H on Pt in Table 1) to a small shift downwards. However, this holds, of course, only, if an implicit solvation model is used in the underlying DFT calculations as done in this work.
Using a surface slab in vacuum would yield highly underestimated differential interfacial capacitances and would correspondingly then also yield only minor shifts of the electrosorption potential even within the second-order approximation.
At both surfaces, the FGC approach yields shifts of the electrosorption potential U sorp;FGC RHE that almost coincide with the second-order approximation at low pH. Towards higher pH, the two approaches then exhibit increasingly different results, with the FGC approach leading to smaller absolute shifts. These findings are consistently obtained also for other H adsorption sites and at other surface coverages (see Supplementary Figs. 5-9) . Figure 1 illustrates this explicitly by also showing the results calculated for full ML coverage. The primary effect of changing coverage is an offset of the various U sorp RHE , as expected from the different average H adsorption energies. U sorp;2nd RHE and U sorp;FGC RHE exhibit only small changes of their slope and curvature. This is in agreement with the often rigid shift of experimental CV peaks at overall unaltered shape 5 . From the above derivations, we see that this will occur whenever the change of work function and interfacial capacitance scales approximately linearly with the surface coverage. Then the double-layer correction G HÀclean exc;DL ðU sorp RHE Þ will scale similarly with N H as the CHE term G HÀclean exc;CHE ðU sorp RHE Þ, and the pH terms in the overall condition for U sorp;2nd RHE , Eq. (16), become coverage independent. Such a behavior can be expected for most adsorbates, whenever adsorbate-induced work function changes derive from local dipoles and capacitance variations from local alterations of the capacitance.
While the comparison of the various theoretical approaches sketches a rather consistent picture for both Pt(111) and Pt(100), clear differences arise in Fig. 1 in the comparison to the existing experimental data [48][49][50] . For Pt(111), the non-Nernstian shift of the proton electrosorption peak is reproduced very nicely by the FGC approach. Even though this holds also for the absolute positions, we note that this comparison is less well-defined. First, we infer only approximately the absolute experimental positions from the better resolved half peak-currents of the CVs in refs. 48,49 . (see also Supplementary Fig. 11). Second, absolute theoretical positions are sensitively affected by the employed approximate DFT exchangecorrelation functional, through the change in the H adsorption energy. The good agreement reached with the PBE functional could thus be partially fortuitous as discussed further below. More important is the fact that the FGC approach yields non-Nernstian shifts of the correct sign and order of magnitude without having to invoke a pH-dependent change of the adsorption energy. To this end, it is important to point out that a previous theoretical CHE study by Karlberg et al. 46 found no impact on the adsorption energy upon application of an electrostatic field and inclusion of explicit interfacial water molecules. In contrast, most recent results from ab initio molecular dynamics suggest a small influence of explicit water on average adsorption energies 81 . Nevertheless, even such refined dynamical simulations on neutral slabs, i.e., remaining at the CHE-level in terms of electrode charge equilibration, could not explain the non-Nernstian pH shifts of U sorp RHE at this surface. We therefore believe that the rationalization of these shifts within the FGC approach is indeed correct for the right reason, i.e., that the appropriate account of charge equilibration in the FGC scheme is key to describe proton electrosorption here.
On Pt(100) the performance of the FGC approach as measured by the existing experimental data is not as good. The peak position is significantly misaligned and the non-Nernstian shift with pH is predicted in the wrong direction and with wrong magnitude. The FGC approach predicts an almost negligible shift, in contrast to experiments where the magnitude of the shift is of comparable size as at Pt(111), just in the opposite direction. As already pointed out, the misalignment is most likely related to the  50 for Pt (100). The linear fits are guide to the eye. Table 1. Work function (Φ 0 ) and area-normalized interfacial capacitances (C 0 ) at the PZC for the systems studied in the main text at the lowest (clean) and highest studied coverage θ. More data in the SI.

System
Adsorption site  82,83 yields, for instance, a peak position lowered by 140 meV to 0.33 V, which agrees very well with the average experimental position. At the same time though, this functional would also lower the predicted Pt(111) peak position by 130 meV and thus worsen the apparent agreement with experiment for this surface. Accepting a typical ±200 meV uncertainty of GGA-level absolute adsorption energies, the truly critical point assessing the description of the electrochemistry are therefore not the absolute positions, but the non-Nernstian shifts of the electrosorption peaks. To this end, the inability of our calculations to describe the experimental pH shifts at Pt(100) are most likely related to the absence of the explicit water structure in the present calculations. As shown convincingly in a recent study by Cheng et al. 84 , an increasing distance between metallic surface and interfacial water for higher pH values lead to a decrease in the adsorption energy and thereby to a peak shift of correct magnitude on the RHE scale. In general, any such influence of surface-specific interfacial water 58,59,[85][86][87][88][89][90][91][92] can not be captured in calculations based only on implicit solvation. This highlights that present-day FGC calculations still need to be seen as a numerically efficient approximation-that as a next step might need to be refined by appropriately extending them to mixed explicit/implicit solvation models 39,58 . We also note in passing that Pt(100) has no proper double layer region (co-adsorbed H, OH) 49 ; co-ion effects might thus also contribute to the overall pH-shift on the RHE scale 93 . As such, proton electrosorption will remain an intriguing test system for future more explicit calculations to address the complexity of the interface.
Electrosorption valencies of halide ions at Ag(111) Electrosorption valencies are another fundamental and, as we shall see, closely related concept in the context of electrosorption in interfacial electrochemistry 51,52 . The electrosorption valency l a (or formal charge number 77 ) denotes the number of electrons that flows onto the electrode upon electrosorption of one adsorbate a at applied potential Φ E . While some debate exists of how to actually measure l a or relate it to other concepts like partial charge transfer [51][52][53][54][70][71][72][73][74][75][76][77][78] , l a is a well-defined thermodynamic quantity. According to the Lippmann equation 47,65 , the electronic surface charges are given by the derivative of the interface energy with respect to the applied potential. Within the present ab initio thermodynamics framework, this definition translates to obtaining the absolute number of electrons N α;abs e ðΦ E Þ of an electrode configuration α from the derivative of the charge-equilibrated excess energy G α exc ðΦ E Þ with respect to Φ E , Without loss of generality, we will focus here again on a symmetric slab setup, where on each of the two identical surfaces one adsorbate a electrosorbs to an empty site of an electrode configuration α characterized by N a adsorbates of the same type and at the same site type. This results in a new electrode configuration with 2(N a + 1) adsorbates (and, as before, without changing the number of substrate atoms N s ). The electrosorption valency of this process is then where the factor 1/2 renormalizes the symmetric slab setup. If we insert the general second-order expansion (Eqs. (11) and (14)), the electrosorption valency can again be separated into a CHE and a DL part At the CHE level we recover the assumption dictated by chargeneutrality that each adsorbed ion of charge q a drags q a /e electrons onto the surface. The electrosorption valency is thus purely given by the formal ionic charge: it is a constant integer number, independent of the applied potential and surface structure or composition. This simplification is lifted by the twoterm DL correction. The first potential-independent term renormalizes the formal ionic charge by the charge difference Δq net ¼ c 2Na 0 ΔΦ 0 and therewith allows for non-integer electrosorption valencies. The second term introduces an explicit potential dependence. To this end, it is important to point out that in this derivation both DL correction terms formally contain extensive differential capacitances. However, the intensive electrosorption valency is still well defined, as we consider the idealized situation where only one adsorbate electrosorbs on each side of the symmetric electrode slab. This is further clarified when switching to a mean-field picture that characterizes the adsorbate layer only in terms of the surface coverage θ = N a *A site /A, where A site is the crystallographically determined surface area per adsorption site. This allows to describe the finite change in adsorbate numbers in terms of changing coverage, e.g., This expression then contains only intensive quantities: A site , the differential capacitance C θ 0 , and the derivative of the work function with respect to θ. In this formulation, the DL-correction to the electrosorption valency becomes l DL a ðΦ E ; θÞ ¼ i.e., in essence a charge per adsorbate necessary to keep the work function constant upon electrosorption. ΔQ net is identical to the term introduced before in the discussion of linear pH shifts for proton electrosorption in Eq. (21), where the derivative was approximated by the average change per adsorbate. Considering that dΦ θ 0 dθ ¼ À 1 ϵ0Asite d z can be related to the normal component of the adsorbate dipole d z 79 (with electrochemical sign convention), we recognize that the terms in Eq. (27) are closely related to the dipolar and capacitive contributions to the electrosorption valency that are recurrently discussed in classical electrochemistry works 47,51,52,94,95 . In fact, the relation between d z and l a is identically reported by Schmickler 95 .
To illustrate the effects of the DL corrections, we consider the  27) and therefore the double-layer correction to the electrosorption valencies. Figure 2 provides the corresponding valencies, computed with the CHE approximation, with the full second-order double-layer correction and only considering the second-order ΔQ net term. By construction, the CHE approximation yields the formal charge −1 for the entire halide ion series. Instead, the second-order terms introduce varying non-integer values which vary with coverage. To this end, it is important to note that the extrapolations of the corresponding lines to the limit θ → 0 do not coincide with the CHE results, supporting the fundamental differences between FGC and CHE calculations also for individual adsorbates in the limit of infinite cell size. In this respect, it is reassuring to see that already the consideration of the second-order ΔQ net term reproduces the expected trend of the computed electrosorption valencies with the Pauling electronegativity over the entire low-coverage range. This also extends to the full DL correction, but the actual values obtained depend then on the chosen applied potential (cf. Eq. (27)). Deferring a detailed comparison to experimental data to a forthcoming publication, we include in Fig. 2 measured "integral" electrosorption valencies for low coverages 54 . Also in light of the large uncertainties of different measurement methods 53,54 , the agreement obtained with the DL corrected calculations is very encouraging, and confirms the importance of FGC approaches toward a reliable first-principles description of electrified interfaces.

DISCUSSION
Ab initio thermodynamics offers a computationally efficient and ideally predictive access to questions in interfacial electrochemistry. In the CHE approximation it enjoys high popularity, in particular in applications to electrocatalysis. In this work, we used electrosorption as a fundamental elementary step in electrocatalysis to compare this approach to FGC ab initio thermodynamics. The present analysis underscores the importance of double-layer charging, missing in the zero net-charge CHE approach. In the context of electrosorption this manifests itself in the ability to describe non-Nernstian pH-shifts of electrosorption peaks and non-integer electrosorption valencies; both effects can be captured at the variable-charge FGC level. As evaluated against existing experimental data for proton electrosorption at Pt electrodes, and halide ion electrosorption at Ag electrodes, the present first-principles description yields in general results closer to experiment than pure CHE calculations and is partly already semi-quantitative. Analysis of a second-order approximation to the FGC energetics allows to rationalize the relevance of the potential of zero charge, the interfacial capacitance, and work function shifts induced by adsorbates e.g., through their dipoles. All the expressions derived here for the electrochemical stability of adsorbates and their electrosorption valencies are in agreement with classical electrochemistry works 47,[51][52][53][68][69][70][71][72][73][74][75][76][77][78]94,95 , closing the gap between phenomenological and first-principles treatments of electrified interfaces. On a more practical level, the analysis with the second-order model provides clear insight into the limitations of zero-charge calculations and guidelines for their usage: • Increased deviations of the CHE energetics from FGC results are expected whenever the electrode potentials Φ E considered are significantly different from the potential of zero charge Φ 0 , as the correction terms scale with (Φ E − Φ 0 ).
• Energetic differences between CHE and FGC calculations are proportional to the magnitude of the interfacial capacitance C 0 ; CHE energetics thus becomes problematic for high electrolyte concentrations with large C 0 values. Furthermore, a non-negligible influence of the electrolyte chemistry can be expected, due to ion-specific interfacial capacitances 96,97 . DFT slab calculations in vacuum underestimate C 0 significantly and do therefore not allow to assess the effects of an applied potential.
• Major differences between CHE and FGC calculations will occur for adsorbates with significant dipole moments, whenever the work function change per adsorbate is large. This is in line with other studies 30,40 and was noted early 7,98 .
The improved description at the FGC level comes at a somewhat increased computational cost. This does not so much refer to the actual DFT calculations themselves, thanks to efficient implementations of current implicit solvation models. Instead, it refers to the larger number of DFT calculations of the system in different charge states. Depending on the specific application this may typically mean~10 times more calculations. While this will typically not be prohibitive, we note that the number of calculations required to determine the quantities for the secondorder model is generally smaller, as it only includes parameters obtained at the PZC, namely adsorption energies, Φ α 0 and the interfacial capacitance C α 0 . Considering the success of this model in describing the electrosorption phenomena studied in this work, calculations at this level might therefore be an appealing intermediate option for systems where computational cost is truly limiting. Having said this, we emphasize that the remaining inaccuracies of FGC calculations can still be significant. Notable factors are the approximate DFT exchange-correlation functional and their errors in adsorption energies, the neglect of explicit water 58,59,86 , or missing co-ions 93 . In addition, the present analytical derivations indicate that discrepancies between allimplicit and all-explicit interfacial capacitances 39,55,92,99 and work functions 58 can have an impact as well. We expect many of these limitations to be overcome by FGC-type schemes with implicit/ explicit hybrid descriptions of interfacial water, an approach we will pursue in the future.

DFT calculations
All DFT calculations reported below are performed with the Quantum ESPRESSO package 100 (PWscf), the PBE exchange-correlation functional 101  , see text). The valencies l θ are calculated with the CHE approximation (solid black lines), the full secondorder double-layer (DL) correction (dotted blue lines) and the second-order ΔQ net term (dashed red lines) (see Eqs. (25) and (27)). Experimental "integral" electrosorption valencies for the low-coverage regime are from Foresti et al. 54 (orange stars). and pseudopotentials from the SSSP library 102 (v0.7, PBE, efficiency) with density and wave function cutoffs of 360 and 45 Ry, respectively. The two systems, H on Pt and halide ions on Ag are studied in a periodic, symmetric slab setup with substrate and adsorbate degrees of freedom described explicitly. For Pt we use (2 × 2) cells with 7 and 8 Pt layers for the (111) and (100) surfaces, respectively; for Ag(111) we use a ( ffiffiffiffiffi 12 p ffiffiffiffiffi 12 p ) supercell consisting of 6 Ag layers. Slabs are separated by ≈17 Å. Brillouin zone integrations are performed using Γ-centered Monkhorst-Pack meshes (Pt: (10 × 10 × 1), Ag: (4 × 4 × 1)) and a cold smearing 103 of 0.01 Ry (Pt) or 0.02 Ry (Ag), yielding a numerical accuracy better than 1 meV/Å 2 in the interface energies.
As implicit solvation model we use the SCCS implementation of ENVIRON 18,36,104,105 with optimized interfacial parameters (ρ min ¼ 0:0013; ρ max ¼ 0:01025, α = β = γ = 0) and a Helmholtz-layer representation of the electrolyte via planar counter charges. For this setup, good agreement is found between experimentally measured and predicted PZCs, interfacial capacitances, and electrosorption behavior for Pt electrodes, as reported in earlier work 5 . In particular, the work function-when determined as the energy difference between the Fermi level and the flat electrostatic potential in the implicit region-and the electrode potential Φ E on the absolute scale are identical, which is why all electrochemical potentials (ionic and electronic) in this work are expressed on the absolute scale.
As reference energies for the adsorbates we use the Gibbs free energies of the diatomic molecules (p = 1 bar, T = 298 K) in vacuum including the vibrational and rotational degrees of freedom. The free energy corrections ΔF corr are determined as the Helmholtz free energy contributions due to coverage-independent vibrations of the adsorbed species (calculated via finite-differences), neglecting small configurational entropy contributions. Specifically, in the case of H@Pt we thus discriminate only between hollow and bridge sites with the (111)-fcc and (100)-bridge vibrations for 0.25 monolayer (ML) coverage taken as per-adsorbate values.
For the determination of the charge-equilibrated excess energies, we optimize the geometries for the interfacial candidate structures for a set of finite surface charges N net e (≥9 values centered around 0) and evaluate their work functions Φ α and total energies E DFT;α surf , as exemplified in Fig. 3a for the case of H at fcc sites on Pt(111). Values are chosen to sample a potential range of~4-5 V around the PZC (N net e ¼ 2AC 0 ðΦ E À Φ 0 Þ with C 0 ≈30 μF/ cm 2 ) for H on Pt, and in a smaller range of~2 V for the halide systems. As described above and more thoroughly in ref. 5 , charge-equilibrated excess energies G α exc (dots in Fig. 3c) can then be obtained directly for the set of electrode potentials Φ E , which correspond to the obtained set of work functions Φ α ðN net;α e Þ. This involves evaluation of Eq. (1) with the described approximations and for the DFT energies with corresponding surface charges E DFT;α surf N net;α e ðΦ α Þ À Á (dots in Fig. 3a). The lines in Fig. 3 correspond to polynomial interpolations using the data of the studied set of charges. The close agreement between a second-order model and the FGC results for this system are supported by the reported residuals in Fig. 3b, and the accurate parabolic interpolation in Fig. 3c (error~1 meV/Å 2 ) (see also Supplementary Figs. 1-4). The interfacial capacitances follow directly from the second derivative of the interpolant (cf. Eq. (9)). Table 1 reports work functions Φ 0 and area-normalized interfacial capacitances C 0 for the systems studied in the main text, at the lowest and highest coverage. Values for other coverages and adsorption sites are reported in the Supplementary Tables 1-4.

DATA AVAILABILITY
All relevant computational results are provided in the reported tables. including a second-order polynomial fit. b Fit residues for a secondand fifth-order fit. c Charge-equilibrated excess energies G α exc ðΦ E Þ=2A for clean and H(fcc)-covered Pt(111) at applied electrode potential Φ E for pH = 0, including second-order fits.