Interfacial water asymmetry at ideal electrochemical interfaces

Chemistry at the electrochemical interface underpins a wide range of energy conversion,¹ energy storage,^2,3 and chemical synthesis technologies.⁴ An important feature of electrochemical processes is their sensitivity to the electrode potential, which provides an additional mode of control not available for reactions in the liquid phase and at solid–gas interfaces. Designing electrodes and electrolytes to fully exploit this, and target specific chemical reactions requires a comprehensive understanding of the thermodynamics and atomic-scale charge response of the electrochemical interface.

Measurements of electrochemical capacitance as a function of the electrode potential provide a sensitive experimental probe of the overall charge response of the electrochemical double layer. However, the capacitance depends on the solvent dielectric response in the inner layer, ionic response in the diffuse layer, and specific adsorption of ions. It is not straightforward to disentangle these effects, especially near the potential of zero charge (PZC). The capacitance near the PZC is dominated by the low capacitance of the diffuse layer at the low ionic concentrations typically used in experiments. The capacitance at potentials far from the PZC is determined by a combination of dielectric saturation in the inner layer and ion packing effects. It remains unclear from experiments whether the inner layer capacitance (after removing the diffuse layer capacitance dip) exhibits a potential of maximum capacitance (PMC) coincident with the PZC or away from the PZC, indicating an asymmetric charge response of the inner layer. This inner layer response is particularly important for the energetics of chemical reactions at the surface, and it is, therefore, critical to disentangle it from the diffuse layer and ion effects.

A complementary probe of electrochemical interfaces that is more sensitive to the inner layer response is the potential of maximum entropy (PME),⁵ measured from the temperature dependence of the electrode–electrolyte potential difference, V, at a fixed charge density.⁶ Specifically, the electrode charge density σ at which the potential does not change with temperature corresponds to the charge of maximum entropy (CME) because the thermodynamic relation ∂V/∂T|_σ = −∂S/∂σ|_T implies that ∂S/∂σ = 0 at this point. (The corresponding potential is the PME.) This is most cleanly realized by measuring potential transients following laser-induced heating of the interface as this minimizes other effects of a temperature change.^5,7–10 The PME is typically close to and slightly below the PZC for metal electrodes in aqueous electrolytes,^6,7,11,12 e.g., 0.1 V below PZC for Ir(111),¹³ allowing it to be used as an approximate measure of the PZC.¹⁴ The corresponding charge (CME) is ∼−5 μC/cm² for Au(111)⁷ and between −4 and −6 μC/cm² for mercury.⁶ The apparently universal negative CME in metal–water interfaces is attributed to the oxygen end of water facing the electrode at the PZC, requiring a negative electrode charge to counter this preferential orientation and increase the entropy.⁷ This asymmetry of interfacial water also agrees with recent spectroscopic evidence¹⁵ and molecular dynamics (MD) studies on interfacial water dynamics,¹⁶ prompting the question whether the inner layer capacitance is similarly asymmetric.

In particular, do the charges of maximum entropy (CME) and capacitance (CMC) coincide? These are thermodynamically distinct quantities with CME corresponding to ∂V/∂T|_σ = 0 and the CMC to ∂²V/∂σ²|_T = 0 and could be different, in general. In the simplest formulations of an asymmetric continuum dielectric, these quantities may be expected to coincide. For example, if the nonlinear dielectric response $ϵ(E)$ is assumed to peak at a non-zero electric field $E$ at the interface, instead of at a zero field like the bulk response, then, the dipole-orientation entropy will also peak at the same field. This is because the nonlinearity of the dielectric response stems primarily from the competition between the potential energy of the molecular dipoles in the electric field and the entropy of dipole reorientation.¹⁷ In this case, the capacitance and entropy of the solvent layer at the interface will peak at the same interfacial field $E$ and, hence, at the same surface charge density σ. However, this simplified picture does not account for capacitance and entropy contributions from adsorption or electron transfer effects or from beyond the solvent layer (e.g., ions), which could lead to differences between the CMC and CME. Consequently, evaluating the relationship between CMC and CME will be invaluable in developing simplified models of the electrochemical interface, such as continuum solvation models for first-principles electrochemistry,¹⁸ which correctly account for the asymmetry in the charge response of the interface.

Computational prediction of capacitance and entropy of electrochemical interfaces under identical conditions would provide great insight into the relation between the charge response and thermodynamics of the interface. However, this has been challenging due to limitations of MD simulations that can address both these properties. Ab initio MD (AIMD) simulations can capture all relevant physical effects by density-functional theory (DFT) treatment of the electrons, in principle, but are limited in time and length scales required to accurately model the capacitance of the interface. Classical MD simulations can achieve the required scales to model electrochemical interfaces^19–24 but require special care for capacitance predictions.^25,26 In particular, such simulations must account for the electronic response of the electrode not being localized to the surface atoms,²⁷ which can be done by shifting the “effective electron-response plane” based on separate DFT calculations or by incorporating simplified electron-response models such as Thomas–Fermi screening.^28,29 Such techniques have been applied to electrochemical capacitance with ionic liquid electrolytes^28,30–32 but infrequently for metallic electrodes with aqueous electrolyte.³³ Predicting the entropy of the electrochemical interface to evaluate PME or CME has remained even more challenging. Pioneering attempts based on analyzing fluctuations of the work function in ab initio MD simulations³⁴ have been limited in their quantitative comparison to experimental PME due to the computational cost and the difficulty in referencing the electrochemical potential.³⁵ 

In this article, we combine classical MD simulations of ideal aqueous metal electrode interfaces with both the effective electron-response plane approach,^28,30 and the electronic response from DFT calculations. Prediction of absolute capacitance for a specific aqueous metal electrode interface still remains a challenge, so we instead predict a family of potential-dependent capacitance curves corresponding to interfaces with different peak values of capacitance. We find that the potential of maximum capacitance (PMC) relative to PZC is always negative and depends on the peak capacitance value but that the charge of maximum capacitance (CMC) is constant across the family and depends only slightly on the technique used for incorporating the electronic response. We then compare the CMC to the charge of maximum entropy (CME), predicted by directly simulating the heating of the same electrochemical interfaces in MD. We show that the CMC and CME are both negative, indicating asymmetric response of water with the same sign for charge response and thermodynamics but that their magnitudes are distinct at a level well above the accuracy of our predictions.

The capacitance and entropy of real electrochemical interfaces can be strongly affected by the interactions between specific electrodes and electrolytes. In particular, adsorption at the interface can occur to a varying degree for electrolyte ions and water molecules. Experimentally, it can be hard to eliminate ion adsorption pseudocapacitance contributions to the measured capacitance in situations of high surface coverage of adsorbates.³⁶ Additionally, strongly hydrophilic surfaces may adsorb water and alter the surface dipole.^35,37 Here, we target an ideal electrochemical interface without any such effects that add complexity to the interpretation of the CMC and CME.

Figures 1(a) and 1(b) show a typical snapshot of our overall classical MD simulation cell (details in Sec. II B in the following) and a close-up of the interfacial region. With charges on the surface metal atoms and the closest H atoms in water separated by d > 1 Å, these configurations contain a “vacuum” capacitance contribution in series of ϵ₀/d ≈ 8.8 μF/cm², where ϵ₀ is the permittivity of vacuum. The direct capacitance prediction from classical MD must therefore be smaller than this value, as indeed shown in the lowest curve in Fig. 1(c). This is much smaller than the typical ∼50 μF/cm² peak double-layer capacitance of aqueous metal electrodes.^36,37

In a real electrochemical interface, the electronic response of both the metal and water extends significantly past the planes of the corresponding atoms, substantially narrowing the effective vacuum gap d. This effect can be captured, in principle, by ab initio MD simulations with full quantum-mechanical treatment of all electrons, but the nanosecond time scale of ion equilibration is computationally prohibitive, especially for unit cells with thousands of atoms that include a statistically significant number of ions. DFT calculations of metal–water interfaces with frozen water geometries are feasible to calculate electron spill-over effects but miss the dielectric response from solvent-dipole reorientation. Such calculations are therefore difficult to combine with other models to predict the overall interfacial capacitance.

Consequently, we adopt the effective electron-response plane approach³⁰ to narrow the vacuum gap d by virtually moving the electrode charge location toward the electrolyte. As long as the electrode and electrolyte charge densities do not overlap, this does not change the electric field due to the electrode in the electrolyte region by Gauss’s law. Therefore, this change does not alter the MD trajectories and amounts to measuring the electrostatic potential from the MD simulation at a modified location in the final analysis. In principle, DFT can predict the effective charge density response location, and this works reasonably for graphene–water interfaces.²⁸ However, for metal–water interfaces with significantly higher capacitance, the effective gap is nearly zero, leading to extreme sensitivity of the predicted capacitance to the charge location.

To circumvent this issue, we focus not on the absolute capacitance of a specific metal–water interface but the behavior of the potential-dependent capacitance for a family of ideal metal–water interfaces with different charge response locations. Further, the location of the metal atoms in the MD simulation is no longer relevant in the analysis of the capacitance: they serve primarily to set up an interfacial potential and electric field for the electrolyte. Similarly, there is no particularly meaningful or stable spatial location in the liquid profile to reference the plane location. Therefore, we reference the electron response plane location based on the predicted capacitance curves, allowing us to more conveniently compare the trends between different methods in the following. Specifically, we pick the reference Δz = 0 for the effective electron response plane to be the location that yields a capacitance C_PZC = 40 μF/cm² at the PZC. We pick this value because it is within the range of typical metal–water interface capacitance and does not lead to appreciable overlap with the electrolyte charge density. With Δz = 0 as the maximum capacitance curve, we calculate the capacitance for several values of Δz > 0 moving the electrode charge away from the electrolyte. Note that at just Δz = 0.3 Å, the peak capacitance already drops to 20 μF/cm² [Fig. 1(c)], smaller than that for most metal–water interfaces, while the unmodified MD charge density corresponds to a Δz ≈ 2.4 Å.

In summary, we predict a family of capacitance curves for ideal metal–water interfaces by varying the location of the effective electrode charge response location indexed by Δz based on a specific value of capacitance at PZC. Note that the water molecules are free to move closer to the electrode with increasing electric field magnitude (in both directions away from PZC), therefore accounting for any electrostriction effects that are particularly important for highly compressible fluids such as ionic liquids³⁸ and that have also been observed in confined aqueous electrolytes.³⁹ We also investigate the effect of the charge-dependent metal electronic response by directly combining electron density profiles from DFT with the MD charge density (Sec. II C).

We use the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS⁴⁰) to perform classical MD simulations of aqueous 1 mol/l NaF electrolyte between Ag(100) electrodes in a 45 × 45 × 60 ų unit cell with periodic boundary conditions. The chosen ionic concentration is deliberately high compared to typical electrochemical experiments for two reasons: (1) the diffuse layer capacitance is less significant at this concentration, allowing us to focus on inner layer properties, and (2) lower concentrations will require larger unit cells and present greater statistical sampling challenges. We pick NaF as the simplest of nominally non-adsorbing electrolytes, although F⁻ may exhibit greater specific adsorption than more complex compound ions such as Cl$O4−$ or $KPF6−$⁠. The space between the electrodes is 45 × 45 × 46 ų of which Fig. 1(a) shows a section with ∼10 Å depth visible into the plane of the page.

Each simulation includes two back-to-back half cells with more than enough distance to separate the two diffuse layers, given that the screening length is $∼3$ Å for this electrolyte at 298 K (using the Debye screening length as a rough estimate). We set up the same charge in both half cells to create nominally inversion-symmetric unit cells, avoiding issues with long-range dipole interactions and ion equilibration between the two cells. The electrodes are treated with a single layer of charged atoms with effective potentials capturing the interaction of the electrolyte with an Ag(100) slab as discussed in the following and are separated by 14 Å vacuum, found to be sufficient to suppress the interaction between periodic images of the electrolyte. For each randomly initialized configuration discussed in the following, we perform energy minimization followed by NVT simulations at 298 K with a 2 fs time step, discard the first 1 ns for equilibration, and capture statistics over 9 ns.

We use the rigid extended simple point charge (SPC/E) water model⁴¹ with molecule geometry constrained by the SHAKE algorithm^42,43 and with Lorentz–Berthelot mixing of Lennard-Jones parameters (arithmetic for σ and geometric for ϵ) except between the Na⁺ and F⁻ ions for accurate treatment of the ion-pair interactions.⁴⁴ The SPC/E water model underestimates the dielectric constant of water and the air/water surface tension by 10%–20%¹⁷ but captures their trends with temperature and should suffice for our exploration of asymmetric charge response.⁴⁵ We parameterize a Morse potential for the short-ranged interaction between the electrode atoms and water/ions from electronic DFT calculations of a single molecule/atom next to a neutral Ag(100) surface, as detailed in the supplementary material. This effective interaction parameterized to DFT includes the image-charge attraction between water molecules and the wall. We therefore do not explicitly require a fixed-potential treatment of the metal electrode to capture this effect⁴⁶ but neglect dependence of the electrode–electrolyte interaction on the potential by using fixed metal atom charges. Note that while we require specific metal atom parameters for the MD simulation, we analyze the results more generally for the properties of ideal metal–water interfaces as discussed previously.

We initialize the simulation cell with a random distribution of water molecules created using a simple Monte Carlo insertion method that enforces a minimum O–O separation of 2 Å. We then randomly replace a selected number of water molecules with Na⁺ or F⁻ ions, one each from equally sized bins in the z direction; this ensures a uniform initial spatial distribution of ions to mitigate ion equilibration times. Finally, we picked the number of water molecules and ions iteratively based on trial MD simulations to ensure a bulk density of 1 g/cm³ for water and 1 mol/l for NaF far in the center of the simulation cell. For the neutral electrode, we end up with 2809 water molecules and 46 Na⁺ and F⁻ ions each.

We charge the electrodes in steps of one additional ion of the same type (either Na⁺ or F⁻) in each half cell, amounting to a step of $1e−/(45Å)2≈0.79μ$C/cm² with a compensating charge distributed equally among all the surface metal atoms. We extend these simulations up to 20 extra ions of each type, thereby spanning electrode charges from σ ≈ (−16 to +16) μC/cm². While the difference in ion numbers is fixed by charge neutrality of the simulation cell, the total number of ions may vary with the electrode charge. Grand canonical simulations could ensure that the density of ions approaches the target bulk value (1 mol/l each) in the center of the simulation cell but are challenging to perform at the scale required here. Instead, we evaluate two simplified schemes of setting the ion numbers. The first scheme fixes the sum of cation and anion numbers (at 2 × 46 = 92 in our case). The second scheme fixes the number of minority ions with the same charge as the electrode (at 46 in our case) and increases the number of majority ions of opposite charge as the electrode. Table S3 in the supplementary material lists the ion numbers for all electrode charges in both schemes, and we contrast their predictions in the following.

For each of the 41 electrode charge values (including neutral) in both ion-number schemes, we perform five independent MD simulations starting from different random configurations, yielding ten half cells for each charge point. We then compute the planarly averaged charge density profile (as a function of z) with the electrode charge density offset by different amounts as discussed in Sec. II A. For each offset, we solve a 1D Poisson equation to get the electrostatic potential profile for each electrode charge and measure the potential V between the electrode and the bulk region of the electrolyte. Note that this is exactly equivalent to planarly averaging the 3D Poisson equation because the Coulomb kernel and planar average are both diagonal operators in reciprocal space and, hence, commute with each other. For capacitance calculations, we only need differences in the potential far from the interface and therefore do not need to worry about short-ranged contributions to the local potential seen by an ion near the interface (Madelung potential).⁴⁷ Additionally, the interface potential difference, V, that we calculate from a point charge model will differ by an overall constant compared to more realistic charge distributions of the water molecules and ions,⁴⁸ but this does not impact the differences in the potential used in V—PZC and the differential capacitance evaluation. Finally, we compute C = ∂σ/∂V from a cubic spline fit of the (σ, V) data to obtain a differential capacitance curve. (See Fig. S2 in the supplementary material for details.) Note that we perform each charge simulation from completely independent random configurations, and therefore, the smoothness of the obtained charging curves confirms adequate equilibration of our simulations.

Figure 2 contrasts the two schemes for setting ion numbers for a given electrode charge. Figure 2(a) shows that at the highest electrode charge of +16 μC/cm², the bulk density of both ions is depleted to ∼60% of the initial 1 mol/l value when the sum of ion numbers is fixed. In contrast, fixing the minority ion number keeps the bulk ionic concentration much closer to the targeted value. Figure 2(b) shows the family of capacitance curves obtained from simulations using both schemes. Interestingly, the capacitance is mostly unchanged between these schemes, except for a slightly larger capacitance in the fixed-minority scheme for the most strongly charged electrodes. We can understand this by estimating the minimum diffuse layer capacitance, C_d = ϵ/λ_D, where ϵ is the bulk permittivity and λ_D is the Debye screening length of the electrolyte (ignoring activity factors of ions for this order of magnitude estimation). C_d changes from $≈230μ$F/cm² at 1 mol/l concentration to $≈180μ$F/cm² at 0.6 mol/l. In series, this can introduce a change of at most 3 μF/cm² when the total capacitance is around 50 μF/cm². Notice that the impact of ionic concentration changes is negligible here because of the high 1 mol/l concentration that leads to a high diffuse capacitance and would be much more significant for simulations with lower ionic concentrations.

Finally, to evaluate the charge and potential of maximum entropy, we repeat the entire set of classical MD simulations mentioned above at 318 K and compute dV/dT for each electrode charge density σ from finite-difference derivatives between 298 and 318 K. Further details on the analysis of these results are discussed below in Sec. III C.

So far, we discussed capacitance prediction from the classical MD simulations by offsetting the electrode charge location to account for electronic response. We also present results in the following that include the electronic response of the electrode from DFT. We perform electronic DFT calculations of a seven-layer Ag(100) surface in the JDFTx code⁴⁹ with the Perdew–Burke–Ernzerhof (PBE) exchange-correlation functional,⁵⁰ ultrasoft pseudopotentials⁵¹ at kinetic energy cutoffs of 20 hartrees for the wavefunction and 100 hartrees for the charge density, a 12 × 12 × 1 k-point mesh, and Fermi smearing of 0.01 hartrees. The slabs are separated by 16 Å vacuum, and truncated coulomb potentials are used to remove periodic interaction between the slabs.⁵² We apply electric fields perpendicular to the surface from 0 to 17.8 V/nm in steps of 0.89 V/nm (21 calculations) to obtain equal and opposite surface charge densities in steps of 0.79 μC/cm², matching the MD simulations mentioned above.

We extract the difference in the electron density from the zero-field simulation as the DFT electrode charge density profile for each of the 41 surface charge densities simulated in MD (zero and 20 charge magnitudes of each sign). The above-mentioned DFT simulations implicitly include the nonlinear response of the electrode to surface charge densities. To test the effect of the nonlinearity, we also calculate the linear-response change in the electron density to an infinitesimal electric field using density-functional perturbation theory (DFPT) with all other parameters identical to the DFT calculations mentioned above. We then replace the planarly averaged charge density contribution from the classical MD electrode with these electronic charge density profiles (both from DFT and DFPT) before solving the Poisson equation for the potential to analyze the impact of electronic charge response on the capacitance.

We begin with an analysis of the variation of electrolyte charge distributions with electrode charge as predicted from ensembles of long-time classical MD simulations. Figure 3 compares ion density and total charge density profiles averaged over ten half cells each for various electrode charges with negative electrode charges on the left and positive electrode charges on the right. Note that the results on the left and right are from separate simulations: each MD simulation contains electrodes of the same sign to avoid overall unit cell dipoles, as discussed in Sec. II B.

Each ion density in Fig. 3(a) increases substantially from the bulk values in the vicinity of the oppositely charged electrode—Na⁺ near the negative electrodes on the left and F⁻ near the positive electrode on the right, as expected. Correspondingly, the ion densities are suppressed near the like-charged electrode. The peak ion densities are ∼4 Å away from the electrode surfaces (at ±23 Å in z). In another 2–3 Å further, the ion profiles transition to a rapid decay toward the bulk density, as expected.

Importantly, the ion profiles are not equal for the neutral electrode: there is a small excess of F⁻ closer to the electrode with a small peak of Na⁺ further out. We would intuitively expect such profiles for a slightly positively charged electrode but instead find it for a neutral electrode. Consequently, the ion profiles we expect for a neutral electrode would instead appear for slightly negatively charged electrodes. This gives a first indication that properties we expect to be symmetric about zero charge from classical continuum models, such as dielectric response and capacitance, might be centered at a negative charge instead.

The ion densities discussed above are further from the electrode than the first layer of water, and the net charge density profiles of the electrolyte shown in Fig. 3(b) are dominated by water for the first 3 Å from the electrode. In particular, the charge density adjacent to a neutral electrode starts with a nearer positive H peak, followed by a negative O peak further away, in agreement with AIMD predictions for metallic surfaces.⁵³ When the electrode is charged positively, the nearer H peak is suppressed in magnitude and the further-away O peak is enhanced. In contrast, for negative electrode charges, the strength of both the H and O peaks is enhanced, leading to a larger charge density response than the positive case. Most importantly, the charge response is closer to the electrode on the negative side: this should lead to a smaller potential difference for the same charge magnitude and, hence, a larger capacitance on the negative side.

As discussed in Sec. III A, we expect the capacitance of the interface to peak at negative electrode charges, consistent with the predicted capacitance curves shown in Figs. 1(c) and 2(b). Next, we turn to a quantitative analysis of the capacitance asymmetry and the location of its maximum. Figure 4 shows the family of capacitance curves for ideal aqueous electrochemical interfaces with different peak capacitances obtained with different models of the electrode charge density and different schemes of positioning the electrode charge density relative to the electrolyte charge density from MD with the corresponding charge densities at the reference position, Δz, shown in Fig. 5.

First, Fig. 4(a) shows the direct prediction from MD after offsetting the electrode charge density location by different Δz, as detailed in Sec. II A. The potential of maximum capacitance (PMC) is always negative, but its magnitude is inversely proportional to the peak capacitance value (see Fig. S3 in the supplementary material). Instead, if we look at the electrode charge density at the PMC, it is constant across the family, resulting in a charge of maximum capacitance (CMC) of −2.9 μC/cm² (with an uncertainty $∼0.1μ$C/cm² as estimated in Fig. S2 in the supplementary material).

We can understand this behavior by noting that the asymmetry in the response of the water is a built-in polarization for the neutral electrode. For a specific value of the interfacial electric field E, this built-in polarization is neutralized, and we can expect this point to be the location of maximum capacitance. This interfacial electric field is E = σ/ϵ₀ by Gauss’s law, directly determined by Gauss’s law, while the electrode potential depends on the overall electrostatic potential profile and the electrode charge response location (Δz in particular). Hence, the CMC should be invariant across the family of capacitance curves, while the PMC depends sensitively on the overall capacitance. Finally, note that the capacitance increases slightly for large potentials: this is due to an increased density of the electrolyte in response to high electric fields at the interface (electrostriction).¹⁷ The magnitude of this effect will be sensitive to the electrode–electrolyte interaction potential and may be obscured by ion adsorption at potentials far from the PZC in the experiment anyway. Consequently, we focus on the asymmetric behavior of the capacitance close to the PZC in the following.

We next consider the impact of deviations from the simple point charge model of the electrode charge density (sheet charge in the planar average) considered so far. First, we use the linear-response charge density profile from DFPT (see Sec. II C) with a constant shape scaled to each electrode charge density value [Fig. 5(b)]. Figure 4(b) shows the result of using this charge density profile instead of the MD charge but with the same placement scheme as mentioned above: at specific locations in z (with various Δz offsets). We find absolutely no difference in the capacitance curves and CMC, which can be explained by Gauss’s law as long as the electrode and electrolyte charge densities do not overlap. The exponential tail of the electronic charge response of the electrode does, in fact, overlap partially with the electrolyte charge density for the highest-capacitance Δz = 0 case [Fig. 5(b)], but this impacts the capacitance negligibly. The overall potential difference for a given electrode charge density shape would match that for the sheet charge case when the sheet charge is placed at the center-of-charge of the charge distribution. This center-of-charge location is absorbed into our definition of Δz based on the C_PZC = 40 μC/cm² criterion (Sec. II A), and so, the family of capacitance curves remains unchanged.

Next, let us account for the actual variation of the electron density of the electrode. First, consider the effect of the electron density on the short-ranged potential of the liquid: an increased electron density would lead to higher repulsion that pushes non-bonded liquid atoms away. Figure 4(c) includes this effect by placing the electrode and electrolyte charges based on the overlap $n̄(r⃗)=∫dr⃗′n(r⃗′)nwater(r⃗−r⃗′)$ of a water molecule’s electron density $nwater(r⃗)$ and the electrode electron density $n(r⃗)$⁠, as parameterized in the SaLSA solvation model.^18,54 Specifically, the separation at which $n̄$ crosses $nc=1.42×10−3a0−3$ correlates with the non-bonded distance of nearest approach.⁵⁴ We place the DFPT electrode charge density relative to the MD electrolyte density profiles based on this condition for each electrode charge and the offset by various Δz to obtain the family of capacitance curves. Interestingly, we find that the charge density profiles in Fig. 5(c) are unchanged from the previous case. Correspondingly, the capacitance curves and CMC still do not change [Fig. 4(c)], indicating a negligible effect of the change in short-ranged repulsion with electrode charge density.

Finally, we use the full nonlinear variation of electrode charge density profile from DFT in Fig. 4(d), which leads to a small but noticeable difference in the capacitance curves and CMC. In particular, the magnitude of the CMC and the overall asymmetry are reduced compared to all previous cases. Essentially, in the nonlinear response of the DFT, electron repulsion makes it is easier to positively charge the electrode by removing electrons than to negatively charge it by adding electrons. This effect favors higher response on the positively charged side [Fig. 5(d)] and therefore reduces the overall asymmetry toward negative charges due to the water at the interface. The net result is still a negative CMC but reduced slightly in magnitude to −2.2 μC/cm².

We showed that the capacitance peak occurs for negatively charged electrodes with a CMC of −2.9 μC/cm² based on the MD charge densities, which reduces in magnitude to −2.2 μC/cm², accounting for nonlinearities in the electronic response of the electrode from DFT. We found that the charge is a better measure of the maximum point compared to the potential because the charge directly determines the interfacial electric field seen by the water surface, while the potential depends more globally on the overall electrostatic potential. Similarly, we expect the potential of maximum entropy (PME) relative to PZC to be inversely proportional to the peak capacitance for a family of ideal electrochemical interfaces with different capacitances, while the charge of maximum entropy (CME) will be constant. Consequently, we will focus on comparing the CME to the CMC predicted above.

We predict CME by directly mimicking the experimental approach: heat the electrochemical interface and measure the change in the electrode potential ∂V/∂T at fixed charge. Specifically, we repeat the entire set of simulations used to generate the above-mentioned results (which were for 298 K) at 318 K and compute ∂V/∂T as a finite-difference derivative for each electrode charge density, σ (Fig. 6). We find that ∂V/∂T crosses zero with a positive slope near −5 μC/cm², and since ∂S/∂σ|_T = −∂V/∂T|_σ, this implies ∂S/∂σ will cross zero at this point with a negative slope. Therefore, σ ≈ −5 μC/cm² should be a local maximum of entropy as a function of electrode charge.

The error bars shown in Fig. 6 are calculated as the standard deviation of the ∂V/∂T calculated from each of the ten half cells simulated for each charge. Note the extremely small magnitude of the results: the potentials only change by ∼10 mV over our 20 K baseline; the CME prediction therefore requires large ensembles with sufficient statistics and a careful analysis. To precisely pin down the zero-crossing of ∂V/∂T, we fit a general Kernel ridge regression model so as to not bias the result by choosing a more restrictive function form such as a polynomial. We repeat the fit for 1000 re-samplings of the data by selecting the result from different half cells at each charge. The band shown in Fig. 6 is the 90% confidence interval obtained from this ensemble of fits. Finally, from the zero-crossing of this band, we can quantify the 90% confidence interval of the CME to be −5.1 ± 0.6 μC/cm².

Our predicted CME for ideal metal–water interfaces is in excellent agreement with experimental measurements of −5 μC/cm² for gold⁷ and −4 to −6 μC/cm² for mercury.⁶ Most interestingly, it differs from the capacitance peak location (CMC) of −2.9 to −2.2 μC/cm² from Sec. III B. In addition, note that the 20 K difference used in the calculation of ∂V/∂T does not affect the interface properties appreciably and cannot be the reason for the difference between CMC and CME: the CMC shifts by at most 0.1 μC/cm² between 298 and 318 K (Fig. S4 of the supplementary material). [Figure S5 additionally confirms this by repeating the ∂V/∂T calculation from 298 to 308 K simulations, finding a confidence interval of (−6.1 ± 1.6) μC/cm² for the CME, in agreement with the 20 K difference result mentioned above.] The CMC and CME indicate asymmetric charge and thermodynamic response in the same direction: higher for negatively charged electrodes. However, by carefully calculating both quantities from the same set of MD simulations, we can unambiguously conclude that the CMC and CME do not coincide even for ideal electrochemical interfaces. This provides a renewed incentive to experimentally measure the CMC, which as discussed in the Introduction, is challenging because the low ionic concentrations typically used to avoid ion adsorption lead to a capacitance dip that obscures the precise location of the capacitance maximum.

We have performed classical molecular dynamics simulations and evaluated the capacitance and potential of maximum entropy for aqueous, charged metallic interfaces. We find distinct, non-coincident values for the CMC and CME. For surfaces with large capacitance, the potentials of minimum entropy and maximum capacitance will be very similar. Our findings of the asymmetric response of interfacial water open new questions about the electrochemical interface and the response properties themselves.

Future work is necessary to understand why the CME and CMC do not coincide even for ideal interfaces. This could stem from the different spatial regions that contribute to each effect. The entropy is sensitive to the entire polarized region of each half cell in the interface, while the interfacial (series) capacitance is most sensitive to the lowest-capacitance region in space: closest to the metal. These spatial dependencies could be explored further by varying the ionic concentration, and the generality of CME–CMC differences can be tested using MD simulations of other asymmetric solvents such as acetonitrile.¹⁸ 

Going beyond ideal interfaces, extensions of this approach can systematically quantify the impact of specific electrode and electrolyte properties on the charge response and thermodynamics of the double layer. In particular, hydrophobicity of the electrode and ion sizes are known to impact the structure of interfacial water.^55,56 Finally, the detailed capacitance and entropy predictions of charged interfaces presented here will facilitate future development of more accurate solvation models for electrochemistry.

See supplementary material for details on force field para-meters, extraction of capacitance, and capacitance and entropy calculations at additional temperatures.

A.S. and R.S. acknowledge support by the U.S. Department of Energy, Office of Science, Basic Energy Sciences, under Award No. DE-SC0022247. All calculations were carried out at the Center for Computational Innovations at Rensselaer Polytechnic Institute. We thank Dr. Thomas P. Moffat for his useful suggestions.

The authors have no conflicts of interest to disclose.

The data that support the findings of this study are available from the corresponding author upon reasonable request.