Effects of surface rigidity and metallicity on dielectric properties and ion interactions at aqueous hydrophobic interfaces

Electrostatic forces play a major role in the interactions between matter on the nano-scale and consequently determine the macroscopic behavior of many materials. Electrostatic interactions between charged objects such as lipid membranes, proteins, molecules, and ions are profoundly influenced by the surrounding water.¹ On the linear response level, the time-averaged effect of water on electrostatic interactions is quantified through the static dielectric constant ɛ, which is spatially constant in homogeneous bulk systems and reduces the force between two charges by a factor of 80 compared to vacuum. Close to an interface, however, the effect of water on electrostatic interactions is more intricate. The water density near an interface deviates from its bulk value and the proximity of a surface restricts the molecular conformations. Therefore, it is not surprising that interfacial water has a different structure and a different dielectric response compared to bulk, as was shown in molecular simulations as well as in experiments.^2–11 The change of the dielectric properties has a pronounced influence on the behavior of ions at interfaces¹² and, for example, influences their distribution at solid^13,14 as well as fluid (i.e., liquid–liquid or liquid–vapor) interfaces,^15,16 but also fundamentally influences the transport of ions across interfaces, in particular between an aqueous and an apolar phase.¹⁷ Studying the distribution and transport of ions and charged molecules is of fundamental importance in various fields of physics, chemistry, and biology, from heterogeneous catalysis to ion extraction and from electrochemical applications like energy storage to drug delivery.^18–22 For understanding the influence of interfacial water effects on charged objects, the determination of their solvation free energy profile across interfaces is of great importance, because the free energy profile represents the thermodynamic driving force for transport. The free energy itself is influenced by non-electrostatic as well as electrostatic contributions, where the latter are fundamentally determined by the dielectric water properties at the interface.

To shed light on the water properties at hydrophobic planar interfaces, we here use molecular dynamics simulations and investigate the interfacial dielectric behavior as well as the free energy profile of chloride ions at different aqueous interfaces. In order to disentangle the effects of rigidity and metallicity of the surface, we investigate metallic as well as nonmetallic versions of graphite, which is a very rigid surface, and treat vapor–water (short for vapor–liquid water) and hexane–water interfaces, which are two examples of soft fluid interfaces, as well as a (nonmetallic) graphene–water interface, which again is rigid. Using previously developed methods for extracting the spatially resolved dielectric tensor,^6,9,23 we find pronounced differences in the dielectric water response at the different surfaces, which reflect varying dielectric interfacial excesses that can be explained by shifted positions of the dielectric dividing surface (DDS). By splitting the ionic free energy profiles into electrostatic and non-electrostatic Lennard-Jones (LJ) contributions, we find that the repulsion of ions from an interface is dominated by the electrostatic contribution for fluid interfaces and by the non-electrostatic contribution at solid interfaces, in agreement with previous findings for a thiocyanate ion.²⁴ Not surprisingly, a sharp dielectric-interface model for the ionic free energy profile fails for the fluid interfaces, which are characterized by broad and smeared-out dielectric response profiles for both perpendicular and parallel directions. Furthermore, by expanding the Coulomb part of the ionic free energy profile in powers of the ionic charge, we demonstrate that for all interfaces, nonlinear dielectric effects, as quantified by higher-order expansion coefficients, are modified significantly as the ion approaches the surface. However, the magnitudes of the different nonlinear components are different for each interface type, which can be understood by symmetry arguments: The cubic contribution is similar for all surface types and accounts for the nonlinear coupling between ionic hydration and water orientation at the planar interface, the quartic contribution differs between different surface types and accounts for interface deformation, in particular at the soft fluid–fluid interfaces. The sharp-interface dielectric model for the quadratic part fails for vapor and hexane because the interfaces are smeared out (not because they are flexible); interface flexibility shows up mostly in the quartic contribution. Surprisingly, the metallic nature of graphite only plays a marginal role in determining the interfacial water structure and the interfacial free energy profile of single ions. Deep in the vapor or the hexane phase, we find that several solvating water molecules stay bound to a chloride and sodium ion in equilibrium, with the average number of hydration waters being different for chloride and sodium and for vapor and hexane.

We simulate interfaces between liquid water and water vapor, hexane, graphene, and metallic and nonmetallic graphite. For water, the nonpolarizable SPC/E model is used.²⁵ Production runs are executed in the NVT ensemble at a temperature of T = 300 K. The first set of simulations is without ions. In a second set, we add a chloride or a sodium ion with a charge q = ∓e at variable positions as well as a sodium or a fluoride counterion, respectively, fixed in the center of the water slab, which avoids artifacts in charged inhomogeneous systems,²⁶ the interaction of the ion with its counterion and all periodic images is subtracted using analytical expressions. Ion LJ parameters are taken from Ref. 27.

The simulation box of the vapor–water system has an extension of 5.1 × 5.4 × 70 nm³ and contains 9024 water molecules that spontaneously form a liquid-water slab with a thickness of roughly 10 nm and a vapor phase with a thickness of roughly 60 nm. To avoid motion of the water slab toward the ion, we apply a harmonic force on the center of mass of the system in the z direction with a force constant of 1000 kJ mol⁻¹ nm⁻². It was shown previously²⁸ that such a force has no detectable effect on the structure and the dynamics of the water slab.

The hexane–water system has an extension of 4 × 4 × 12.5 nm³ and contains 2000 water molecules and 630 hexane molecules. The hexane is modeled with the all-atom OPLS (Optimized Potentials for Liquid Simulations) force field²⁹ and, thus, carries constant partial charges. Before production runs, the phase-separated system is equilibrated for 1 ns in the NPT ensemble at 1 bar pressure. The water slab thickness is roughly 5 nm. The nonaqueous phases of the vapor–water and the hexane–water systems are saturated with water molecules, as discussed below.

The graphene–water system has an extension of 5.1 × 5.4 × 70 nm³, where the water is confined between two graphene sheets with a separation of 9.9 nm. The carbon atoms are frozen in space and nonpolarizable and do not carry any charge, the graphene sheets are, therefore, nonmetallic. We adjust the water number between the graphene slabs to a value of 8701 molecules, for which the water chemical potential equals the bulk value of (−11.706 ± 0.003) k_BT.³⁰ 

All simulations above are performed using the GROMACS molecular dynamics simulation package,³¹ and the systems are periodically replicated in all directions. In the case of the vapor–water and the graphene–water interfaces, a slab-correction is used.³² 

We also simulate metallic and nonmetallic versions of graphite using the MetalWalls molecular dynamics simulation package,³³ employing periodic boundary conditions in lateral directions only. For the metallic system, a constant electric-potential boundary condition is used on one of the two graphite slabs, while the second graphite slab is treated as nonmetallic. In these simulations, the charge on the carbon atoms is allowed to fluctuate with the constraints of a constant electric potential at the atomic positions and global electroneutrality.³⁴ The aqueous phase is treated with constant point charges and creates a fluctuating electric field, which, in turn, determines the polarization charges on the metallic surface atoms. The metallic graphite–water system has an extension of 3.41 × 3.69 nm², with a separation between the graphite slabs of 5.51 nm. Each graphite slab is composed of three graphene sheets with a separation of 0.3354 nm consisting of 480 carbon atoms per sheet and 2071 water molecules in between. Graphite is used instead of graphene to ensure that the induced charges in the metallic slab decay to zero on the third sheet. As for the graphene systems, all carbon atoms are frozen in space and the water number is chosen so that the water density in bulk equals that for the graphene system. All force field parameters are identical to the graphene systems. For direct comparison, the same system is also run treating both graphite slabs as nonmetallic but keeping all other parameters the same.

Further information about the methods and simulation details are given in Sec. S1 in the supplementary material.

Figure 1 shows profiles of the total mass density ρ_m(z), the parallel dielectric response ɛ_‖(z), and the inverse perpendicular dielectric response $ε⊥−1(z)$ as well as the laterally averaged electrostatic potential φ(z) for the water–vapor, water–hexane, and water–graphene systems. All profiles are shifted according to the Gibbs dividing surface (GDS) that is based on the water mass density profiles and shown by the vertical black dotted line in Fig. 1(d), that means, z = 0 denotes the GDS position for all systems. The water phase is located in the right half-space (z > 0) whereas the nonaqueous phase is located in the left half-space (z < 0). Figures 1(a)–1(c) show simulation snapshots of the three systems and illustrate that the interfacial water layer is highly corrugated at the vapor and hexane interfaces, while it is perfectly flat at the rigid graphene sheet. ρ_m(z) shown in Fig. 1(d) exhibits the expected sigmoidal shape for the liquid vapor–water and hexane–water systems with a broadening that originates mostly from the intrinsic interfacial water roughness and to a lesser extent from capillary waves.³⁵ In contrast, ρ_m(z) of the graphene system shows the typical oscillations due to water layering at rigid flat interfaces and a peak at z = −0.19 nm due to the graphene sheet. In Sec. S3 of the supplementary material, we show that the water–vapor interfacial profiles are reproduced rather well by convoluting the graphene profiles with an interface position distribution. This suggests that local details at the water–vapor interface are to some degree smeared out by interface corrugation and capillary fluctuation effects.³⁵ 

The parallel ɛ_‖(z) and the inverse perpendicular $ε⊥−1(z)$ dielectric profiles for the fluid vapor–water and hexane–water interfaces show sigmoidal shapes that are qualitatively similar to the shapes of the corresponding mass density profiles. The parallel profile ɛ_‖(z) for graphene exhibits pronounced maxima and is roughly proportional to the corresponding mass density profile, but the perpendicular profile $ε⊥−1(z)$ shows multiple sign changes and is not related to the water mass density profile near the graphene wall in any simple manner.^6,9 Note that the sign changes of $ε⊥−1(z)$ reflect divergencies of ɛ_⊥(z), a clear indication of overscreening at rigid surfaces, which is only observed for the perpendicular component.^6,9 The vertical dashed lines in Figs. 1(e) and 1(f) depict the DDS positions, which are constructed analogously to the Gibbs dividing surface [see Eq. (S4) in the supplementary material] and quantify the dielectric excess of each interface. We see that the DDS for all systems and for both parallel and perpendicular directions is located to the left of the GDS, meaning that the dielectric interfacial excess is generally positive; in other words, water at hydrophobic surfaces has a higher dielectric response per water molecule at the interface than in bulk. The location of the DDS is different for each system and the ordering of the three systems does not agree between the parallel and perpendicular directions. Whereas for the parallel component, the ordering from left to right is hexane, graphene, vapor, meaning that hexane has the highest and vapor the smallest surface excess, the ordering for the perpendicular component is hexane, vapor, graphene, meaning that graphene has the lowest excess. A high dielectric excess means that the dielectric screening of a charge close to the interface is strong, which, in turn, leads to reduced electrostatic interactions. It transpires that interfacial electrostatics is expected to be weakest for hexane for a fixed location relative to the GDS, but due to the different ordering of the parallel and perpendicular components of graphene and vapor, it is not clear whether the graphene–water or the vapor–water system shows the strongest electrostatic interactions at the interface (a point we will come back to in further detail below). In summary, we find pronounced anisotropy of the dielectric interfacial behavior: Not only are the parallel and perpendicular dielectric profiles very different, but also the parallel and perpendicular DDS positions are different and there does not seem to be a simple correlation of the DDS positions with the rigidity of the interface. The only trend we can discern is that for soft interfaces, the parallel DDS is closer to the GDS than the perpendicular DDS, whereas the opposite behavior is found for the rigid graphene surface.

Figure 1(g) shows the laterally averaged electrostatic potential $φ(z)=−∫−∞zdz′∫−∞z′ρ(z′′)dz′′/ε0$⁠, obtained by integrating twice over the charge density profile ρ(z) due to the water (and hexane where applicable) partial charges. We see that the potential for the graphene system shows pronounced oscillations, which seem correlated to the density and the parallel dielectric profiles but not to the perpendicular dielectric profile. In contrast, the potential in the fluid systems changes monotonically. The potential consists of a dipolar contribution due to interfacial water ordering and a quadrupolar contribution that is only weakly dependent on the interfacial water ordering.⁶ We see that the potential in the water bulk is negative for all systems but differs considerably between the different systems; the vapor–water system reaches a potential of φ(z → ∞) = −0.6 V far from the interface, while for the graphene and the hexane systems, we find reduced values of φ(z → ∞) = −0.38 V and φ(z → ∞) = −0.37 V, respectively. In fact, the negative sign is produced by the dominant negative quadrupolar water contribution, which is an intrinsic water property and depends sensitively on the water force field but is independent of the interfacial structure; in ab initio simulations, the quadrupolar water contribution has a positive sign, as has been amply discussed in the literature.^36,37 The dipolar contribution due to water interfacial orientation is less dependent on the water force field; therefore, the difference in the total potential seen in Fig. 1(g) is expected to be a robust feature of these systems, indicating that the orientational interfacial water ordering is more pronounced for vapor compared to graphene and hexane. Note that φ(z) is the electrostatic potential experienced by an infinitesimally small test charge that laterally averages over the entire water phase (including the interior of the water molecules) and is, therefore, unrelated to the electrochemical potential experienced by an ion, as will be discussed further below.

Obviously, the interfacial dielectric properties influence the interfacial ion solvation and the ionic interaction with an interface, but since an ion significantly perturbs the interfacial water structure and possibly induces nonlinear dielectric effects, it is not clear to what extent dielectric profiles are sufficient to predict the ion free energy close to an interface. We, thus, place a chloride ion at variable positions relative to the interface and extract its free energy profile using thermodynamic integration (TI), as done previously in Ref. 38. Figure 2 shows the free energy profiles F(z) for the vapor, hexane, and graphene systems, where the reference point F = 0 is located in the center of the water slab. Figure 2(c) shows the total free energy profile F(z) = F_LJ(z) + F_Coul(z), while in Figs. 2(d) and 2(e), the separate Lennard-Jones and Coulomb contributions are shown, which are subsequently obtained from thermodynamic integration (see the supplementary material, Sec. S1 E). In other words, F_LJ(z) is the free energy profile of an uncharged LJ sphere, while F_Coul(z) is the free energy needed for charging the sphere; thus, both F_LJ(z) and F_Coul(z) contain contributions from LJ as well as Coulomb interactions. We find that F(z) is repulsive for all three systems. The decomposition demonstrates that the interplay between the LJ and Coulomb contributions is dramatically different for the systems. While the total free energy is dominated by repulsive Coulomb interactions for the vapor–water and hexane–water systems, it is dominated by repulsive LJ interactions for the graphene–water system.

As the ion approaches the vapor–water and hexane–water interfaces, the Coulomb contribution in Fig. 2(e) increases monotonically, which is due to interfacial polarization effects and corresponds to image-charge repulsion in a simple continuum dielectric model.³⁹ In contrast, the Lennard-Jones contribution in Fig. 2(d) becomes attractive if the ion approaches the vapor or hexane phase, which is caused by the positive (i.e., unfavorable) solvation free energy of a Lennard-Jones sphere in water.^40–42 Note that for negative z values, we observe the formation of a water finger around the ion in the vapor and hexane phases, as shown in Figs. 2(a) and 2(b) and reported previously for vapor^28,43–46 and organic solvents.^43,44 The hexane phase produces a lower total free energy for z < 0, because of a less repulsive Coulomb contribution. This can only partly be rationalized by the slightly higher bulk dielectric constant of our nonpolarizable hexane force-field model (ɛ_hexane = 1.01 from Fig. 1) compared to water vapor. We also find that F_LJ(z) in the hexane phase is slightly more repulsive than in vapor, which one could be tempted to associate with the different LJ interaction parameters between chloride and hexane, σ_Cl,Hex = 0.4 nm and ϵ_Cl,Hex = 0.34 kJ mol⁻¹, and between hexane and hexane, σ_Hex,Hex = 0.35 nm and ϵ_Hex,Hex = 0.28 kJ mol⁻¹. However, based on the larger chloride–hexane LJ interaction parameter ϵ_Cl,Hex compared with the hexane–hexane parameter ϵ_Hex,Hex, one would expect the net LJ solvation free energy to be more favorable in hexane compared to vapor. We conclude that the different total solvation free energies of a chloride ion in hexane and vapor must be related to details of the microscopic solvation structure in hexane and vapor, as will be discussed further below. In contrast to the hexane and vapor systems, and as previously shown in Ref. 38, the total free energy for graphene is dominated by the Lennard-Jones repulsion, as seen in Fig. 2(d).

Finally, comparing the Coulomb contributions for each surface type in Fig. 2(e), we find that at the graphene interface the Coulomb contribution increases most steeply, which is related to the different water mass density profiles at the different surfaces, being very smooth for the fluid interfaces and rather sharp for the graphene layer, as shown in Fig. 1(d). A second difference between the systems is that the ion loses solvating water molecules as it approaches the graphene layer, while for the vapor and hexane systems, the formation of a water finger maintains a complete hydration shell around the ion even when the ion crosses the GDS, which results in a reduced Coulomb free energy repulsion in the nonaqueous phase.

If the ion moves further away from the interface into the vapor or the hexane phase, the water finger breaks but the ion stays hydrated by a water solvation shell that consists of a few water molecules.^28,45,46 The solvating water molecules are in equilibrium with the water-saturated vapor or hexane phase and with the water bulk phase. From simulations without an ion, we obtain for the vapor phase a water number density of c_vap = 7 · 10⁻⁴ nm⁻³, which corresponds to an ideal gas pressure of P_vap = 30 mbar, slightly lower than the experimental water-vapor pressure of P_vap = 42 mbar at 300 K.⁴⁷ For the hexane–water system we find a slightly lower but comparable water concentration of c_hex = 5.9 · 10⁻⁴ nm⁻³ in the hexane phase, in good agreement with the experimental value c_hex = 6.3 · 10⁻⁴ nm⁻³⁴⁸ and leading to a corresponding partial water pressure in the hexane phase of P_hex = 24 mbar based on the ideal gas equation. See the supplementary material, Fig. S2, for a time series of the number of water molecules in the vapor phase.

The equilibrium properties of the water solvation shell around an ion in the nonaqueous phase is investigated in separate simulations by fixing an ion sufficiently far from the interface at a distance of z = −7.0 nm in the vapor or hexane phase. The simulation snapshots in Figs. 3(a)–3(d) demonstrate the accumulation of hydration water around an ion over time, starting initially with an unhydrated ion configuration. Figures 3(e) and 3(f) show the number of water molecules $NH2O$ in a spherical shell of radius 1 nm around a chloride and a sodium ion as a function of time, again starting with an initially unhydrated ion. The data demonstrate that chloride and sodium attract roughly the same number of water molecules in their hydration shell, which suggests a similar hydration free energy of the two ions. In Sec. S1 C of the supplementary material, we demonstrate that all hydrating water molecules bind within a single hydration shell around each ion. Additionally, we show in Sec. S6 of the supplementary material that the ion hydration structure at the water–vapor interface using the TIP4P/ɛ water force field⁴⁹ is comparable to the structure obtained for the SPC/E force field, which is used for all our results.

We rationalize these results by an analysis of the binding free energy of a single water molecule onto an ion.²⁸ The hydrating water molecule is described by the partial oxygen charge, q_O, located at a distance B_O,i from the ion i, and the partial hydrogen charges q_H = −q_O/2 at a distance of B_H,i from the ion center; see insets in Fig. 4 for a schematic description of the hydration model. The hydration free energy is the sum of the Coulomb free energy and the translational and rotational entropy losses upon water binding to an ion,

where v_naq,j is the molecular water volume either in fully saturated vapor, v_{naq, vapor} = 1/c_vap = 1415 nm⁻³, or in hexane, v_{naq, hexane} = 1/c_hex = 1698 nm⁻³ (see the supplementary material, Sec. S1 B), and $vhyd, i=4πBO,i3/3$ is the approximate hydration volume, i.e., the volume available for a hydrating water molecule. We estimate the chloride–oxygen distance as B_{O, Cl} = 0.32 nm and the sodium–oxygen distance as B_{O, Na} = 0.24 nm from the radial distribution functions g(r) in liquid water shown in Fig. 4. Additionally, B_{H, Cl} = 0.27 nm and B_{H, Na} = 0.40 nm are calculated using the SPC/E geometric parameters²⁵ (see Sec. S4 for details on the calculations). Using the SPC/E partial charges, we find U_Coul/k_BT ≈ −24 for chloride and $≈−44$ for sodium. The (rescaled) translational entropy is S_trans ≈ 11 regardless of whether the ion is chloride or sodium and the surrounding medium is vapor or hexane, see the supplementary material, Table S2, for the exact values. The (rescaled) rotational entropy loss upon water binding has been estimated as S_rot ≈ 2 for a similar system.⁵⁰ In conclusion, the Coulomb free energy in Eq. (1) outweighs the entropy loss and we estimate a favorable hydration free energy of about U_bind ≈ −11 k_BT for the first water molecule that adsorbs onto a chloride ion. For the sodium ion, we find a much larger binding energy of U_bind ≈ −31 k_BT. The negative energy corroborates our simulation result that chloride and sodium are hydrated in saturated water vapor and hexane. Note that the binding free energy of subsequently adsorbing water molecules will be reduced due to Coulomb repulsion between water molecules. That the number of adsorbing water molecules is roughly the same for chloride and sodium ions indicates that the first hydration shell around both ions has roughly the same size and that a second hydration is free-energetically unstable.

We find from our simulations in Figs. 3(e) and 3(f) that the number of hydration waters in hexane is slightly larger compared to vapor. This cannot be explained within our simple binding model for a single water molecule and is most likely due to the different structure of the hydration shell in vapor and hexane and related to the reduced interfacial tension between water and hexane compared to vapor. Additionally, we note that the equilibration time for the hydration shell in Figs. 3(e) and 3(f) is much longer in hexane compared to vapor, which is due to the low diffusivity of individual water molecules in hexane compared to vapor.

From the distribution of $NH2O$ shown in Figs. 3(g) and 3(h) and Gaussian fits according to

we extrapolate P(⁠$NH2O$⁠) to $NH2O$ = 0 and thereby extract the hydration free energy of the ions according to F₀ = −k_BT ln P(0) in the hexane and in the vapor phase. The results are given in Table I.

The fit results are subject to large errors due to statistical noise in the simulated distributions and also because of possible non-Gaussian contributions to the distribution for $NH2O$ → 0. Nevertheless, we find as a trend that (i) the solvation free energy F₀ in the hexane phase is substantially larger (i.e., more favorable) compared to the vapor phase and (ii) F₀ is slightly larger for sodium compared to chloride. The second trend is related to the smaller size of the sodium ion, which leads to a stronger hydration of sodium compared to chloride, as also follows from our simple ion–water binding model above. The first trend is in line with the larger number of hydration waters in hexane compared to vapor and, as mentioned above, most likely caused by a reduced interfacial tension between the hydration shell and the surrounding hexane molecules compared to vapor.

Many surfaces are metallic, meaning that charges in the environment induce surface charges such that the potential is constant inside the metal. This also holds for some carbon-based materials such as graphene and graphite. One would expect that surface metallicity has a strong influence on the interfacial water structure and the interactions with ions, since the water partial charges and ionic charges will interact with polarization charges on the metal. In order to check for such effects, we perform constant potential simulations for metallic graphite, where the charges on the metallic atoms are allowed to fluctuate in response to the fluctuating electric field created by the aqueous phase (for more details, see the Methods section and Sec. S1 in the supplementary material). Note that the water and ions have constant charges and, thus, no charge transfer is allowed between the electrolyte and the metallic surface. For comparison, we also perform simulations for nonmetallic graphite.

In Fig. 5, we compare our results for metallic and nonmetallic graphite, where we also add results for nonmetallic graphene for comparison that we discussed before. We note that the LJ force-field parameters and the lateral carbon lattice structure of graphene and both graphite versions are exactly the same. In Figs. 5(a) and 5(b), we show simulation snapshots for the nonmetallic and metallic graphite slab versions, respectively, with the charges on the carbon atoms indicated by color. It is clearly seen that for the metallic graphite, the negative chloride ion induces a positive local surface charge, while a water hydrogen atom close to the surface induces a negative local surface charge. Panels (c)–(f) show structural properties at the interface: the water mass density ρ_m(z), the number density ρ_N(z) of the oxygen and hydrogen atoms (where the hydrogen density is divided by two for comparison with the oxygen density), the charge density ρ(z) as well as the mean cosine of the angle between the water dipole moment and the surface normal, cos(θ). The profiles for the three different surfaces look very similar and are characterized by a strongly enhanced water density and pronounced orientational ordering in the first hydration layer. Water molecules right at the surface point their hydrogen atoms toward the surface, cos(θ) < 0, followed by water molecules within the density maximum that point their oxygens toward the surface, cos(θ) > 0. In particular, we observe only marginal differences in the profiles between metallic and nonmetallic graphite, in agreement with previous reports:^51–54 A slightly stronger orientation of the first water layer on metallic graphite can be deduced from the larger hydrogen number density for z < 0.05 nm in Fig. 5(d) and from the more negative average cos(θ) close to the GDS for z < 0.15 nm in Fig. 5(f).

In Fig. 5(g), we compare the free energy profile for a chloride ion that approaches metallic graphite (red solid line) and nonmetallic graphite (cyan solid line), obtained from Umbrella Sampling (see the supplementary material, Sec. S1 E), note that the red solid line is almost completely hidden behind the cyan solid line. The free energy profile for the nonmetallic graphene obtained from TI is also shown (data points). A very large free-energetic repulsion is found close to the surface. In line with our results for the interfacial water properties in Figs. 5(c)–5(f), no difference can be discerned between the free energy profiles in Fig. 5(g), which is somewhat surprising since one would, based on continuum dielectric theory, expect image-charge repulsion from the nonmetallic graphite and image-charge attraction to the metallic graphite surface.

To bring out small deviations between the chloride free energy at metallic and nonmetallic graphite, we show in Fig. 5(h) the free energy difference given by

As expected, this difference is positive, meaning that the ion is repelled more from the nonmetallic surface than from the metallic surface, but the free energy difference is rather small compared to the total free energy and only amounts to roughly 5%. At the short distances where ΔF(z) is sizable, for z < 0.1 nm, the water is depleted from the surface and the water density in Fig. 5(c) is practically zero. However, this does not mean that the dielectric effect of the interfacial water at graphite is absent, since dielectric effects are rather long-ranged and the DDS positions of both parallel and perpendicular dielectric profiles at graphene in Fig. 1 are negative and, thus, the interfacial water exhibits a positive dielectric excess.

Since the LJ parameters of the metallic and nonmetallic graphite surfaces are the same, ΔF(z) must be caused by electrostatic polarization effects. In order to understand the results in Fig. 5(h) in more depth, we compare with a dielectric continuum model for the interaction of a charged sphere with a sharp dielectric interface.³⁹ The radius of the sphere a = 0.254 nm used in the model reflects the effective linear dielectric radius of a chloride ion and is defined by the linear Born model.³⁸ We neglect the tensorial character of the interfacial water dielectric response here and use an isotropic model with two dielectric constants: ɛ₁ for the surface and ɛ₂ for the water. We set ɛ₁ = ∞ for metallic graphite and ɛ₁ = 1 for nonmetallic graphite and use ɛ₂ = 70, which corresponds to the dielectric constant of the SPC/E water model.⁵⁵ This simple isotropic dielectric model predicts the free energy difference as

where U_iso(λ, z) is the free energy profile of a charged spherical shell at a sharp dielectric interface, given in Eq. (S5) in the supplementary material and first derived in Ref. 39, and λ = (ϵ₂ − ϵ₁)/(ϵ₁ + ϵ₂) is the dielectric contrast. $zDDSnon$ and $zDDSmet$ are the positions of the DDS for the nonmetallic and metallic graphite–water interfaces. Our neglect of the interfacial tensorial dielectric properties is justified, since it was shown for nonmetallic surfaces that the polarization free energy for a tensorial dielectric three-region model, with an anisotropic dielectric layer between z_DDS,⊥ and z_DDS,‖, agrees well with a much simpler isotropic dielectric two-region model with the DDS position given by z_DDS = (z_DDS,⊥ + z_DDS,‖)/2 (see Fig. S5 in the supplementary material).³⁸ We, thus, take for the DDS at the nonmetallic graphite the value for nonmetallic graphene from Figs. 1(e) and 1(f), approximately given by $zDDSnon=−0.11nm$ (taken as the mean of the parallel and perpendicular DDS positions). For the metallic graphite, we cannot extract the dielectric interfacial properties from our simulations, because of the limited simulation length of the constant potential simulations caused by their high additional computational cost. Therefore, we treat $zDDSmet$ as a fitting parameter. The best nonlinear least-square fit of Eq. (4) over the whole z range yields $zDDSmet=−0.1nm$ and is shown in Fig. 5(h) as a red solid line. The pronounced deviations from the simulation data are presumably due to nonlinear dielectric effects, which are neglected in the model but are present in the simulations, as we show in the next section. Most strikingly, we find that $zDDSmet≈zDDSnon$⁠, meaning that the dielectric properties of the interfacial water are suggested to be the same on the metallic and nonmetallic graphite surfaces. This is surprising at first sight but totally in line with the fact that all water structural properties in Figs. 5(c)–5(f) are almost indistinguishable on the metallic and nonmetallic graphite surfaces. As a matter of fact, the precise position of the DDS $zDDSmet$ has only a minor influence on the polarization free energy difference ΔF(z), as shown by the additional curves in Fig. 5(h) where we choose $zDDSmet=±0.2$ nm.

We note that, different from the behavior we here observe for a chloride ion on graphite, sodium ions in a 1 molar sodium chloride solution have been shown to adsorb quite strongly onto metallic gold surfaces.^54,56 This different behavior presumably is linked to the sodium ion being smaller than chloride, but it is most likely also to the different parametrization of the metallic surface in those studies: In Ref. 54, a Drude-oscillator model was used, while in Ref. 56, sodium adsorption was only observed for a large Gaussian width used to characterize the spatial polarization charge distribution of the metallic surface atoms, corresponding to the rather small hardness of gold atoms. The hardness of carbon atoms is higher than that of gold atoms, and the parametrization of the present work was previously shown to reproduce experimental data.⁵⁷ In line with our results, Ref. 56 also finds only marginal differences in the water structure on metallic and nonmetallic surfaces, which supports our finding that metallicity is only relevant for ions very close to the surface.⁵⁸ 

It is well known that ion hydration is asymmetric and nonlinear in the ion charge.^38,59–63 To resolve this asymmetry and nonlinearity as an ion approaches an interface, we determine the linear and nonlinear parts of the electrostatic free energy by expanding F_Coul(z) in powers of the ion charge q,

where we extract the z-dependent coefficients by polynomial fits to our simulation data.³⁸ Details on the method and exemplary fits are shown in Sec. S1 E and Fig. S4 in the supplementary material. In Fig. 6, we show the resulting coefficients for the graphene, vapor, and hexane systems as a function of the distance z of the ion from the GDS. Obviously, the coefficients for the different surfaces decay to the same values at a distance less than a nanometer from the interface and agree with results from independent bulk simulations, indicated by red dashed lines. Figure 6(d) shows the linear coefficient, which corresponds to the electrochemical potential Φ(z) and is defined with respect to a reference position in vacuum, which exhibits only mild deviations between the different interfaces if the ion is in the water phase, i.e., for z > 0. Note that Φ(z) in Fig. 6(d) differs substantially from φ(z) in Fig. 1, since φ(z) corresponds to the laterally averaged electrostatic potential in the water phase while Φ(z) is the potential inside a LJ sphere, which receives substantial contributions from the perturbed water structure around it. A more detailed comparison between φ(z) and Φ(z) is presented in Fig. S7 in the supplementary material.

Figure 6(e) shows the prefactor of the electrostatic part that is proportional to the squared ion charge q², A(z). Here, we observe rather pronounced deviations between the different surfaces even within the water phase. For z > 0, A(z) is substantially higher for the vapor–liquid interface than for hexane and graphene, pointing to more pronounced linear dielectric repulsion from the vapor phase compared to hexane or graphene. Since A(z) reflects the linear dielectric contribution to the polarization free energy, it is this part that can be compared with predictions from continuum linear dielectric theory. We, therefore, add results from the analytical solution of Poisson’s equation using a three-region anisotropic dielectric model that accounts for different positions of the perpendicular and parallel DDS,³⁸ which are taken from Fig. 1. As demonstrated before,³⁸ the three-region anisotropic dielectric model predictions agree nicely with the simulations for the graphene system, but the model is not able to quantitatively reproduce the simulation results for the hexane–water system and even less so for the vapor–water system. The disagreement presumably originates from interfacial broadening of the dielectric profiles for the fluid–fluid interfaces, shown in Fig. 1, which is not covered in the analytical model that assumes sharp dielectric interfaces.

Figure 6(f) shows the prefactor B(z) of the cubic dielectric contribution that is proportional to q³, which is remarkably similar for all three interfaces [similar to the potential Φ(z) in Fig. 6(d)]. By symmetry, this term reflects dielectric nonlinear effects that are odd with respect to the sign of the ion charge. We, thus, conclude that the position-dependence of the B(z) profile results from the nonlinear coupling of the local water orientation around the ion with the water ordering at the planar interface, which seems to be rather similar in magnitude for the three different interfacial systems and does not seem to be affected much by the different broadening of the interfacial dielectric profiles shown in Fig. 1 (which should produce an effect that is even in the ion charge). In contrast, the prefactor C(z) of the quartic contribution that is proportional to q⁴ shows pronounced deviations between the different interface types, similar to the quadratic term A(z). In fact, because C(z) is very different for graphene on the one hand and the two soft interfaces on the other, we conclude that C(z) reflects the capability of the water interface to deform in response to the presence of an ion. We note that our fits become unstable for ion positions far in the nonaqueous phase for z < −1.0 nm, presumably due to the formation of a highly volatile and fluctuating water finger around the ion²⁸ (see Fig. S4 in the supplementary material for details). Thus, when a long water finger is present, a fourth-order polynomial is not sufficient to accurately describe nonlinear interfacial dielectric effects.

We investigate the dielectric properties of interfacial water and the free energy profile of a chloride ion at different hydrophobic surfaces, ranging from soft, such as the vapor–liquid water interface and the hexane–water interface, to rigid, like graphene and graphite. We also consider two different versions of graphite: a metallic and a nonmetallic one. The consideration of such a wide range of systems allows us to separately study the effects of surface rigidity and metallicity on the dielectric and electrostatic properties at hydrophobic–water interfaces.

We find that at soft interfaces, the interfacial water roughness and (to a lesser extent) capillary waves broaden the interfacial mass, charge density, dielectric, and electrostatic potential profiles. As a consequence, effective continuum models that employ a sharp dielectric boundary do not work well for soft interfaces, even when exclusively trying to model the linear dielectric response. In contrast, at the rigid graphene surface, all these profiles exhibit pronounced layering effects and the dielectric profiles are strongly anisotropic. The dielectric excess, quantified by the shift between the dielectric dividing surface position relative to the Gibbs dividing surface, is positive for all considered surfaces but differs significantly between the different surfaces. This means that water at hydrophobic surfaces has a higher dielectric response than in bulk, but this excess is highly surface-type specific. This dielectric surface specificity can be accounted for in continuum dielectric models by proper positioning of the dielectric dividing surface position.

The free energy profile of a chloride ion is characterized by strong repulsion from the low-dielectric surfaces, which is dominated by electrostatics for the soft surface but by Lennard-Jones repulsion for the rigid surfaces. This has far-reaching consequences for the correct modeling of ion–surface interactions in coarse-grained models. Chloride and sodium ions in vapor and hexane are hydrated and surrounded by roughly 5–10 water molecules in equilibrium, which drastically lowers their free energy in the nonaqueous phase.

Comparing nonmetallic and metallic graphite, we find only minute changes in the interfacial water density and orientation profiles. Also, there are only insignificant differences of the chloride free energy profiles between nonmetallic and metallic graphite, which is rationalized by the fact that the repulsion is dominated by non-electrostatic interactions at a rigid surface; the small free energy difference can be modeled by a dielectric model that correctly accounts for the dielectric-dividing surface positions at the two different surfaces. Nonlinear dielectric effects are pronounced at all surfaces and describe nonlinear hydration effects as well as nonlinear interface-deformation effects at the soft interfaces. Linear tensorial dielectric theory with sharp dielectric interfaces works well for the graphene–water interface but not for the diffuse vapor–water and hexane–water interfaces.

The supplementary material contains a detailed summary of the simulation methods, including a discussion of the graphene–water interaction force field used in the simulations, the expression for the electrostatic energy of a spherical charged shell at a dielectric interface, a discussion of the effect of capillary waves and roughness on interfacial profiles, details on the model for the binding of a single water molecule onto an ion, a comparison of the laterally averaged interfacial electrostatic potential and the electrostatic potential acting on an ion, and results for the TIP4P/ɛ water force field.

We acknowledge support provided by Deutsche Forschungsgemeinschaft Grant No. IRTG-2662 Project No. 434130070 “Charging into the future,” the MaxWater initiative from the Max-Planck Society, and by the European Research Council under the European Union’s Horizon 2020 research and innovation program (Grant Agreement Nos. 835117 and 863473). We gratefully acknowledge computing time on the HPC clusters at the physics department and ZEDAT, FU Berlin.

The authors have no conflicts to disclose.

Philip Loche: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Laura Scalfi: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Supervision (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Mustakim Ali Amu: Data curation (supporting); Investigation (supporting); Methodology (supporting). Otto Schullian: Validation (equal); Writing – review & editing (supporting). Douwe J. Bonthuis: Conceptualization (equal); Investigation (equal); Methodology (equal); Supervision (equal); Writing – original draft (supporting); Writing – review & editing (supporting). Benjamin Rotenberg: Conceptualization (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Supervision (equal); Writing – original draft (supporting); Writing – review & editing (supporting). Roland R. Netz: Conceptualization (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Project administration (equal); Supervision (equal); Writing – original draft (equal); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request. The code used for the constant potential simulations is openly available in the repository https://gitlab.com/ampere2/metalwalls.