Dielectric constant of water in the interface

The Maxwell field has played a prominent role in the theories of dielectrics for two mostly disconnected reasons. First, in the case of a homogeneous field produced by a planar capacitor, one gets the direct experimental access to the Maxwell field E through the voltage on the plates V and the distance between them d: E = V/d. The response of the dielectric to the applied voltage is then conveniently characterized in terms of the susceptibility to the Maxwell field. The second reason is that the scalar electrostatic potential ϕ is the solution of the dielectric boundary value problem.¹ The field E follows directly from that solution as E = − ∇ϕ.

The first reason for the importance of E can be viewed as both an advantage and disadvantage since E itself is never accessible experimentally and only the line integral ∫E ⋅ dℓ = V, producing the voltage V, can be measured.² In the case of an inhomogeneous field there is no way to extract the field from the integral and experiment generally does not have direct access to inhomogeneous Maxwell fields. The problem was realized already at the time of birth of the electromagnetic theory. Since inhomogeneous fields cannot be accessed directly, Thompson suggested to use small cavities to measure internal fields inside dielectrics.³ This approach has in fact been realized by modern-day spectroscopy, which allows one to evaluate the local field acting on a dye molecule through the field-induced shift of its spectral line.^4,5 However, the connection between such a local field and the macroscopic Maxwell field has been elusive beyond the standard prescriptions of the dielectric theory.⁶ In addition, the ability to spatially resolve the distribution of the electric field and inhomogeneous polarization within molecular systems of nanometer scale has been limited.⁷ 

From the theoretical perspective, the Maxwell field is well defined by the Coulomb law. The starting point is the overall microscopic electric field E_m, combining the field E₀ of the external charges (vacuum field) with the electric field of all molecular bound charges distributed with the charge density ρ_b (“b” stands for the bound charge). The result is obviously

where

The Maxwell field is produced from this equation as a result of two steps: (i) statistical average 〈E〉 of the instantaneous E_m over the configurations of a statistical ensemble and (ii) coarse graining of 〈E〉 over a “physically small” volume averaging out the microscopic correlations between the molecules of the material.⁸ This volume is not precisely defined and, in fact, is never explicitly involved. The theory, as it is formulated for bulk dielectrics and interfaces, instead introduces coarse graining through constitutive relations as we discuss next.

By taking the divergence of E_m and substituting ∇ ⋅ E₀ = 4πρ₀ for the density of the external charge ρ₀, one arrives at ∇ ⋅ E_m = 4π(ρ₀ + ρ_b). Further, due to the conservation of charge, the instantaneous density of bound charge can be replaced with the divergence of the polarization vector field P_m, such as ρ_b = − ∇ ⋅ P_m.⁸ One arrives at the equation for instantaneous fields

which looks very much like the standard Maxwell equation, except that the fields in this equation refer to an arbitrary statistical configuration of the system. Of course, this equation is just a different form of the Coulomb law, which applies to microscopic dimensions and arbitrary configurations of charges. The two-step averaging and coarse graining procedure mentioned above will produce the averaged and smoothed-out fields E and P and the corresponding electric displacement vector D = E + 4πP. The Maxwell equation for this coarse grained displacement vector follows from Eq. (3) as ∇ ⋅ D = 4πρ₀.

Equation (3) and its coarse grained version still cannot be solved without applying a closure, or constitutive, relation between P_m and E_m or between P and E. The connection between microscopic fields P_m and E_m is a complex problem of statistical mechanics of liquids.⁹ It is therefore assumed that coarse graining helps in eliminating this complexity and leads to local constitutive relations between coarse grained fields

This constitutive relation thus establishes the direct proportionality between the vector fields P and E through the susceptibility χ, which is a scalar parameter for isotropic materials. Empirical evidence suggests that this approximation, when used for macroscopic uniformly polarized dielectrics, yields the bulk dielectric susceptibility χ_s, which is a material property, i.e., a parameter characterizing bulk dielectric and independent of the sample shape (the surface effects die off in the macroscopic limit). Correspondingly, the dielectric constant of bulk dielectric, ϵ_s = 1 + 4πχ_s, is a material property as well.

This result is quite non-trivial since even for coarse grained vector fields the susceptibility χ₀ to the field of external charges E₀ does not share insensitivity to the surface effects (boundary conditions). χ₀ is not a material property, and it depends on the shape of the sample through the dielectric boundary value problem. Given that the inhomogeneous Maxwell field E is not accessible experimentally, most problems of interest for applications involving inhomogeneous fields (solvation of molecules, solvent-induced shifts of spectral lines, interfacial problems, etc.) are formulated¹⁰ in terms of the response to an inhomogeneous external electric field E₀. Nevertheless, the Maxwell field has to be introduced in order to solve the problem since only this field is believed to provide local constitutive relations between P and E required to arrive at the Poisson equation. The locality of the Maxwell field for inhomogeneous external fields does not have firm experimental support and is likely an approximation. This fundamental difficulty is responsible for many complications arising in the general problem of electric polarization of interfaces.^11–14 

The problem of interfacial polarization is solved in dielectric theories by replacing the microscopic fields E_m and P_m at each point in the interface with the corresponding coarse grained fields and then applying the local constitutive relation (4) to each point of the interface. When substituted to the differential form of the Coulomb law in Eq. (3), this formulation leads to the Poisson equation for E fully specified in terms of the external charges. However, there is no factual coarse graining when this procedure is applied to microscopic problems, and it is nearly impossible even to define an algorithm of volume coarse graining when fields are changing on the scale of molecular dimensions. The Poisson equation is obtained in such cases by direct substitution E_m → E and P_m → P and the subsequent use of the constitutive relation. As already mentioned, coarse graining of microscopic fields is not achieved directly by averaging over a judiciously chosen volume, but is produced by applying a specific local form of the constitutive relation. The smooth function E, obtained from the solution of the Poisson equation, then leads to a smooth P, instead of a highly oscillatory function characteristic of interfaces.^14–16 It is the constitutive relation that replaces coarse graining over a small volume in converting the microscopic fields into macroscopic fields.

Since coarse graining is in fact not performed, one can adopt a somewhat different, and less restrictive, form of the constitutive relation involving only the statistically averaged fields in the interface

Of course, Eq. (5) is an approximation. The question we address here is how to build a consistent theory of interfacial polarization when this approximation is applied.

The advantage of Eq. (5) over Eq. (4) is that, in contrast to volume coarse graining, statistical averages are well-defined even on the microscopic length-scale and one can proceed with ensemble-based algorithms of defining susceptibilities. In other words, in contrast to smoothly varied functions P and E in Eq. (4), the corresponding fields in Eq. (5) are highly oscillatory, as usually produced by liquid-state theories and numerical simulations. We do not intend to apply Eq. (5) to the entire interface, but only to the dividing surface separating the solute from the solvent. Furthermore, this relation is applied to projections on the surface normal only and thus becomes a scalar relation, in contrast to the vector form in Eq. (4). The formalism then connects the microscopic calculations in the interface to the electrostatic boundary value problem.

If the constitutive relation is the only step separating the microscopic differential Coulomb law in Eq. (3) from the dielectric boundary value problem, one wonders if this procedure can be supplemented with susceptibilities reflecting the microscopic structure of the polarized interface, i.e., the susceptibility χ in Eq. (5). The standard Maxwell dielectric boundary value problem in fact implements one additional approximation of replacing χ in Eq. (5) with the susceptibility χ_s of bulk dielectric.¹ This approximation is not required if one recognizes the bulk and surface susceptibilities as two separate parameters. Once this separation is achieved, any scalar interface susceptibility can be used in the boundary value problem. Not surprising, the idea of an effective susceptibility or, in other words, of interface dielectric constant has been actively discussed in the literature.^{11,15,17–20} As we point out in more detail below, any such definition requires disconnecting the properties of the interface from the properties of bulk dielectric. The interface susceptibility then becomes a surface property, not directly linked to the bulk dielectric constant.

Some of the phenomenological recipes proposed to deal with microscopic interfaces, such as the popular distance-dependent dielectric constant for solvation problems,²¹ do not withstand the scrutiny of microscopic formulations.²² The problem with such formulations is that spatial correlations of interfacial polarization are typically highly oscillatory functions^14–16 and do not allow defining simple distance-dependent susceptibilities. If any meaningful microscopic susceptibility has a chance to enter the standard boundary value problem, it should be consistently derived from the microscopic Coulomb law in Eq. (3) and not introduced as an ad hoc phenomenological recipe justified by fitting to experimental data or results of numerical simulations. Providing such a consistent approach is the goal of this article. In other words, the main question addressed here is what is the effective (interface) dielectric constant, absorbing into itself the microscopic properties of the interface, that should enter the standard dielectric boundary value problem? We provide a general formulation of the problem, followed by specific calculations of the dielectric response of water interfacing a spherical solute.

When one takes the statistical average in Eq. (3), one arrives at ∇ ⋅ 〈D〉 = 0 inside the dielectric where there are no external charges. This relation translates, through Gauss’ theorem, into the condition of continuity of the projection of 〈D〉 on the unit vector $nˆ$ normal to the interface. This condition can be written as⁸ 

where σ is the surface charge density determined by the normal projections of the polarization density $Pni=nˆ⋅Pi$ (i = 1, 2) in two media in contact in the interface,

In standard dielectric theories, the surface charge density σ screens the external charge. It means that if a probe charge is placed at a large distance from the interface, the effective force between the external charge q and the probe charge is reduced by the opposite charge of the interface polarization and an effective charge q_eff, instead of q, is measured by the force. This is illustrated in Fig. 1 for the simple case of a spherical cavity of radius a with charge q placed at its center, as discussed in the numerical simulations of aqueous solutions below. The positive charge q will, in dielectric theories, create the opposite in sign surface charge density $σ=−(1−ϵs−1)(q/S)$⁠, where S = πa² is the area of the cavity. The effective charge producing the measurable force on an external probe charge, q_eff = q + σS = q/ϵ_s, is then reduced by the dielectric constant of the dielectric ϵ_s.

When the constitutive relation (5) is applied to statistically averaged fields 〈E〉 and 〈P〉 in Eq. (3), 〈E〉 = − ∇〈ϕ〉 satisfies the Laplace equation Δ〈ϕ〉 = 0 inside the dielectric where there are no external charges. The properties of the interface enter the problem through the boundary condition in Eq. (7). Therefore, the goal of reformulating the standard Maxwell boundary value problem needs to focus on introducing microscopic properties of the interface into the boundary conditions of the Laplace equation.

Equation (7) suggests that the only property of the interface one needs to supply to the solution of the Laplace equation is the surface charge density or the normal projection of the polarization density. The linear response approximation⁹ provides the desired property in terms of a non-local susceptibility function χ₀(r, r′) (generally a tensor) depending on two coordinates in the interface,^10,23

where the integral is over the entire space and the 2-rank tensor of susceptibility is

Here, β = 1/(k_BT) and δP = P − 〈P〉.

The susceptibility χ₀ in Eq. (9) is a second rank tensor defined by the corresponding Cartesian components.¹⁰ For some geometries of the interface, it is convenient to consider specific projections of χ₀. For instance, for the planar interface, one defines parallel (⁠$nˆ∥$⁠) and perpendicular (⁠$nˆ⊥$⁠) projections^15,18 as the scalar functions $χ∥=nˆ∥⋅χ0⋅nˆ∥$ and $χ⊥=nˆ⊥⋅χ0⋅nˆ⊥$⁠. Similarly, for spherical solutes which we consider below, one can define the scalar projection on the radial direction,¹⁵ $χ0rr=rˆ⋅χ0⋅rˆ$⁠. Such definitions become less useful for interfaces of arbitrary shape. Using longitudinal and transverse symmetries of the polarization field provides a more general formulation.^13,23 Our goal here does not involve calculating distance-dependent projections of the susceptibility. We focus instead on the normal projection of the polarization density field in Eq. (7), taken at the dividing surface, which can be defined for an arbitrary interface. Once defined, it completes the boundary value problem for the interface electrostatic potential 〈ϕ〉.

The two-point tensor χ₀(r, r′) depends on two positions, r and r′, separately to reflect its interface character and the involvement of three-body solute-solvent-solvent correlations. This needs to be contrasted with the non-local susceptibility of bulk dielectrics depending only on r − r′. To simplify the problem, one can assume that the length of polarization correlations in the interface is much shorter than the characteristic dimension of the interfacial region and apply the local approximation neglecting such correlations altogether,¹¹ 

This approximation obviously eliminates the integral in Eq. (8) shifting the focus to the inhomogeneous susceptibility χ₀(r). It can be obtained by integrating Eq. (9) over r′,

where M is the total dipole moment of the dielectric, δM = M − 〈M〉. Analogs of this equation for different symmetries of the interface have been proposed by Stern and Feller¹⁸ and by Ballenegger and Hansen^11,15,17 and extensively used in a number of recent simulations of interfacial polarization.^{19,20,24–26} We note that the local approximation becomes exact in the limit of a uniform external field considered by Stern and Feller.¹⁸ 

Before we proceed to the exact formula for the susceptibility tensor, not involving the local approximation of Eq. (10), it is useful to connect χ₀(r) to ϵ_s. One obtains for an isotropic dielectric

where the integration is performed over the volume Ω of the dielectric and $Tr[χ0]=∑αχ0αα$⁠. The fluctuation expression on the left-hand side of this equation enters the Kirkwood-Onsager equation for the dielectric constant²⁷ and thus provides the connection between the volume integrated susceptibility and the bulk dielectric constant ϵ_s. There is, however, no direct connection between χ₀(r) and ϵ_s. In other words, polarization fluctuations δP(r) still carry microscopic information, no matter how far from the interface. These fluctuations are coarse-grained, with the microscopic information lost, by volume integration.

Even though the local approximation in Eq. (10) provides a simple resolution of the problem, it is not required.¹² One can take into account the longitudinal character of the field of external charges E₀ and the fact that ∫P_T(r′) ⋅ E₀(r′) dr′ = 0 for the transverse projection of the polarization field P_T (Helmholtz theorem²). Therefore, only the longitudinal projection of the polarization P(r′) = P_L(r′) enters χ₀ in Eqs. (8) and (9). The longitudinal projection of the polarization is in turn connected to the electrostatic field of the bound charge as¹² 4πP_L = − E_b, where E_b is given by Eq. (2). Given that E_b = − ∇ϕ_b, one can apply the Gauss theorem to eliminate the integral in Eq. (8). The result is the exact relation for 〈P_n〉 not requiring the use of the local approximation

Here,

is the fluctuation of the Coulomb interaction energy of the dielectric with the external charges q_i; ϕ_bi is the electrostatic potential of the bound charges of the dielectric at the location of the charge q_i, δϕ_bi = ϕ_bi − 〈ϕ_bi〉.

It is important to emphasize that no specific assumptions regarding either the origin of the polarization density P_n or the electrostatic energy U^C have been introduced in deriving Eq. (13). Both parameters can be microscopic quantities sampled by numerical computer simulations. For instance, a polar liquid with molecular dipoles m_j with coordinates r_j will have the polarization density $Pn(r)=nˆ⋅∑jmjδr−rj$⁠. Correspondingly, U^C can be viewed as the energy of Coulomb interactions of all atomic partial charges of the dielectric with the external charges. This property is routinely provided by numerical simulations. Although P_n is calculated in the interface and thus can be viewed as a local property, the non-local Coulomb interactions of bound charges are incorporated into U^C. The right-hand side in Eq. (13) is fundamentally a three-body, solute-solvent-solvent correlation function, incorporating all chain diagrams of dipolar solvent-solvent interactions responsible for dielectric screening.²⁸ 

Equation (13) is the exact solution for the problem of surface charge density in the interface assuming linear response to the field of external charges. Deriving it does not require constitutive relations. If the constitutive relation (5) is adopted, one can find the statistically averaged electrostatic potential of the interface 〈ϕ〉 from the solution of the Laplace equation and, in addition, ask the question of what susceptibility, or dielectric constant, can be assigned to the interface. Such scalar interface susceptibility can be defined by the equation

This constitutive relation, even though it looks similar to Eq. (4), is in fact a weaker condition. It applies to one projection only, instead of all three Cartesian projections in Eq. (4), and only in the interface, instead of the entire dielectric in Eq. (4).

Equation (15) can be substituted back to Eqs. (6) and (7) to produce the boundary conditions for the Poisson equation

where ϕ₀ is the electrostatic potential of external charges remaining continuous at the dividing surface. This is the only place where the susceptibility of the interface enters the boundary value problem. The local constitutive relations, Eqs. (4) and (5), applied globally to the entire dielectric in dielectric theories, are replaced here with the constitutive relation in Eq. (16) applied to the dividing surface only.

The constitutive equation (15) might be a reasonable approximation for a few molecular layers in the interface, but is not expected to hold globally. Likewise, the susceptibility χ_0n, and the interface dielectric constant ϵ_int defined for spherical solutes below, are parameters characterizing the interface. We, therefore, do not expect them to approach the dielectric susceptibility or the dielectric constant of the bulk material in any specific limit. Even for a macroscopic interface, χ_0n is still an interfacial parameter (like the surface tension), which should not be expected to be simply related to the bulk susceptibility χ_s.

Solving the boundary value problem, i.e., the Poisson equation supplemented by the boundary conditions in Eqs. (6) and (7), yields the electrostatic potential in the interface 〈ϕ〉. Given the external charge q one also obtains the electrostatic energy q〈ϕ〉 and the electrostatic free energy (1/2) q〈ϕ〉 (the Born equation discussed below). However, in contrast to the standard dielectric boundary value problem,¹ this solution does not provide a direct access to the macroscopic polarization, or the macroscopic dipole, of the sample measured in the plane capacitor dielectric experiment. The reason is mentioned in the previous paragraph and is additionally illustrated in Fig. 2 drawn for the simple case of a polar liquid in contact with a uniformly charged plane producing a uniform electric field. When this simplified geometry is used in dielectric theories, the knowledge of the polarization at the dividing surface provides the knowledge of the polarization in the bulk, which is equal to the boundary polarization for the planar geometry. For non-planar geometry, the field of external charges decreases from the surface to the bulk, but the rules of calculating the polarization from the field are uniform for the entire dielectric. The result is the same as for the planar geometry: knowing how to calculate the polarization at the dividing surface implies the ability to apply the same rule to the entire dielectric sample.

This universality is not maintained in the microscopic formulation presented here. While the susceptibility at the dividing surface is the only property of the liquid needed for the boundary value problem, it is not the same as the susceptibility in the bulk. Therefore, finding the macroscopic dipole, connected to the bulk dielectric constant, requires volume integration of $nˆ⋅χ0⋅nˆ$ for the slab geometry or integration of Tr[χ₀] for the spherical geometry, as in Eq. (12). The interface susceptibility χ_0n at the dividing surface becomes a surface parameter and does not provide a direct link to the dielectric properties of the bulk.

Here we apply the arguments presented above to the problem of water polarized at the interface of a spherical solute. The normal to the interface is defined outward from the dielectric,⁸ $nˆ=−rˆ$⁠, $rˆ=r/r$ (Fig. 1). A further simplification of the geometry is achieved by locating the external charges at the center of the solute. All susceptibility tensors become scalars, with the only non-zero diagonal radial component $χ0rr=rˆ⋅χ0⋅rˆ$⁠. We will drop the indexes for brevity with the notation $χ0(r)=χ0rr(r)$⁠. For this specific type of interface, the local approximation leads to

where P_r = − P_n and M_r = − M_n denote the radial projections of the corresponding vectors (Fig. 1).

Since the electric field of the central charge q is E_0n = − q/r², one can define the linear distance-dependent susceptibility χ_0n(r) analogous to the one in Eq. (15)

Here, ϕ_b is the electrostatic potential produced by the dielectric at the center of the spherical solute where the external charge is placed. The interface susceptibility χ_0n follows from this function by adopting r = a, i.e., the radius of the spherical surface separating the solute from the surrounding dielectric. We note that our Eq. (18) is equivalent to its integral form earlier derived by Ballenegger and Hansen,¹⁵ 

where the integral over r′ produces δϕ_b in our Eq. (18).

The definition of the position of the dividing dielectric surface presents a major difficulty for all dielectric theories, and it is not going to go away in our formulation recasting the problem of a microscopic polarized interface as the dielectric boundary problem. The question we are addressing here is what is the susceptibility or the surface charge density that needs to enter the boundary conditions once such a dividing surface is defined. A question relevant to this goal is how sensitive such a definition would be to possible variations of the position of the dividing surface. One ideally wants a robust definition, little sensitive to changes in the cavity radius a.

Figure 3 shows χ₀(r) in the local approximation [Eq. (17)] and the exact χ_0n(r) [Eq. (18)] calculated from molecular dynamics (MD) simulations performed in this study. The simulations are done for TIP3P water²⁹ interfacing spherical solutes of varying diameter and interacting with the oxygen of water by the Kihara potential (a hard-sphere repulsion core with the radius R_HS combined with a surface layer of soft Lennard-Jones potential).^30,31 One has to keep in mind that correlation functions in Eqs. (17) and (18) fundamentally reflect three-particle solute-solvent-solvent correlations. Relatively long MD simulations, ∼200 ns, were therefore required to converge them for each solute studied here. More details on the simulation protocol are given in the supplementary material,³² and here we discuss the results.

It is clear from the calculations that the local [Eq. (17)] and exact [Eq. (18)] formulations for the radial interface susceptibility generally agree with each other. The exact formulation is obviously preferable since it is free of the locality assumption. Both results show an oscillatory behavior of the interface susceptibility, leading to potential uncertainties when the cavity radius is altered. Some type of averaging over the oscillations, or coarse graining, is needed to arrive at a robust definition of interface susceptibility and the corresponding dielectric constant. An approach developed previously¹² and adopted here is to represent χ_0n(r) in Eq. (18) as the derivative of the correlation function based on the integrated dipole moment M_r(r) of water within the sphere of radius r

If differential in the above equation is taken at each point, one recovers χ_0n(r), with its oscillatory behavior shown in Fig. 3. Alternatively, instead of taking the differential at each point, we determine the linear slope with respect to r to average out the oscillations of χ_0n(r) caused by molecular granularity. This linear slope then provides us with the scalar coarse grained susceptibility of the interface with oscillations averaged out. Figure 4 shows that indeed the slope can be well defined from the correlation function 〈δM_r(r) δϕ_b〉 calculated in a few molecular layers.

The susceptibility χ_0n to a radial external field E_0n can be associated with the interface dielectric constant ϵ_int according to the relation (ϵ_int − 1)/(4πϵ_int) = χ_0n,^12,15 which leads to

The values, obtained from the slopes of the radial correlation functions shown in Fig. 4, are presented in Fig. 5. We find that ϵ_int decreases slowly from ∼9 to 4 as the effective size of the solute increases from ∼5.5 Å to 18.5 Å. Overall, the value of the interface dielectric constant is much smaller than the bulk value for TIP3P water, ϵ_s ≃ 97.³³ The definition of ϵ_int by Eq. (21) is prone to numerical instabilities when 4πχ_0n becomes greater than unity due to calculation errors. ϵ_int is not required for the solution of the boundary value problem in Eq. (16) and χ_0n is sufficient. It is presented here solely because of the history of the subject casting the dielectric boundary value problem in terms of the dielectric constant instead of susceptibility.

The results shown by circles in Fig. 5 are obtained for neutral Kihara solutes. Even though Eq. (13) contains the electrostatic interaction energy with the external charge of the ion, which is proportional to the charge magnitude, charge cancels out when the surface susceptibility is defined by dividing the surface polarization by the ion field in Eq. (15). We therefore operate in the linear response domain, when one can assume that the presence of an external charge does not alter the structure of the interface used to perform the statistical averages. That this is indeed the case is demonstrated by simulating Kihara solutes with positive and negative charges placed at their centers. These results are shown by diamonds in Fig. 5 and are indistinguishable from the results obtained for neutral solutes (see the supplementary material for more details³²). Our simulations are indeed consistent with the linear response approximation.

The small value of the interface dielectric constant of water has potentially dramatic consequences for the problem of hydration. Our formalism anticipates that ϵ_int is used in the dielectric boundary value problem. Therefore, the solvation free energy of a spherical ion carrying charge q and assigned the cavity radius a is given by the standard Born equation^34,35 F = − (1/2) χ_Bq², where the Born solvation susceptibility is

Here, χ_0n(a) and ϵ_int(a) indicate that the dependence of the Born solvation susceptibility on the cavity radius can be more complex than the traditionally anticipated a⁻¹ scaling.²² 

The reasons for the relative success of the Born equation in predicting the free energy of solvation and its dramatic failure in describing entropy of solvation have long been known.^36–39 Both are related to the low sensitivity of the Born formula to the solvent properties when the bulk dielectric constant ϵ_s ≫ 1 is used instead of ϵ_int in Eq. (22). The traditional form of the Born equation significantly underestimates the entropy of solvation³⁶ since the term $ϵs−2∂ϵs/∂T$⁠, appearing in the entropy, is too small. This deficiency can be potentially remedied if, according to our calculations, ϵ_int ≪ ϵ_s. The final verdict requires knowledge of ϵ_int(T). Our estimate for TIP3P water gives $∂ϵint/∂TV≃−0.8×10−4K−1$ (supplementary material³²), which is significantly lower than $∂ϵs/∂TP=−0.36K−1$ of bulk water. Whether this low value is shared by more realistic force fields of water is not clear at the moment.

We discuss here a formalism connecting the Maxwell boundary value problem with the microscopic structure of the interface. In other words, the paper asks the question: What is the dielectric constant that should enter the boundary conditions in the Laplace equation describing a polarized dielectric interface? The problem is formulated in terms of the interface susceptibility or, alternatively, the interface dielectric constant. This property is calculated from an exact equation statistically averaging correlated fluctuations of the interface polarization density and the electrostatic energy of external charges interacting with the polarized dielectric. Evaluated by MD simulations of water interfacing spherical solutes, the interface dielectric constant is found to be significantly lower than the corresponding bulk value.

This research was supported by the National Science Foundation (Grant No. CHE-1464810) and through XSEDE (Grant No. TG-MCB080116N). The authors are grateful to Daniel Martin for his help with the setup of simulations.