Dielectric response of thin water films: a thermodynamic perspective

The surface of a polar liquid presents a special environment for the solvation and organization of charged solutes, which differ from bulk behaviors in important ways. These differences have motivated many attempts to understand electrostatic response at aqueous interfaces in terms of a spatially varying dielectric permittivity, typically concluding that the dielectric constant of interfacial water is significantly lower than in the bulk liquid. Such analyses, however, are complicated by the potentially nonlocal nature of dielectric response over the short length scales of interfacial heterogeneity. Here we circumvent this problem for thin water films by adopting a thermodynamic approach. Using molecular simulations, we calculate the solvent's contribution to the reversible work of charging a parallel plate capacitor. We find good agreement with a simple dielectric continuum model that assumes bulk dielectric permittivity all the way up to the liquid's boundary, even for very thin (∼1 nm) films. This comparison requires careful attention to the placement of dielectric boundaries between liquid and vapor, which also resolves apparent discrepancies with dielectric imaging experiments.


Introduction
Interest in conned water has exploded over the last decade or so, owing principally to advances in the fabrication of devices at the nanoscale, 1-3 the potential implications for 'blue energy' and desalination, 4 and as means to understand fundamental properties of water 5,6 and its solutions. 7,8 An obvious consequence of the decreasing length scales associated with connement is an increase in the surface-to-volume ratio of liquid water, which typically amplies surface-specic effects relative to large sample geometries. The notion of nanoconned liquid water thus having properties that are inherently different to its bulk counterpart has inspired many attempts to reformulate intensive material parameters typically used to describe the bulk uid. In particular, many years of investigation along these lines [9][10][11][12][13][14][15][16][17] has focused on the static dielectric constant e liq , whose role in mediating electrostatic interactions impacts upon, e.g., solvation, capacitance and electrokinetics. Further motivation for such theoretical studies comes from recent dielectric imaging experiments 6 of water conned between two atomically at walls separated by distances as small as 0.8 nm. These imaging results were inferred to report an interfacial dielectric constant e int ¼ 2.1 (relevant to an interfacial region of thickness l int z 7.5 AE 1.5Å) that dominates the capacitance of a thin water lm. This value, typical of a bulk nonpolar liquid, signies a dramatic departure from the polarizability of bulk water, for which e liq z 80.
At the microscopic level, it is well recognized that water's interfaces exhibit local average properties that differ from the bulk liquid, varying continuously with depth within a molecular length scale l int of the surface. 18 Accordingly, many studies have aimed to rationalize conned water's electrostatic response in terms of a local dielectric constant e(z) that varies with position z along the surface normal. [9][10][11][12][13][19][20][21] Molecular simulations have estimated e(z) either from polarization uctuations, or from response to external electric elds; in either case this approach relies upon interpreting features that have been resolved at a ne scale within a theoretical framework appropriate for macroscopic dielectric materials. In this study, we pursue a different approach. Specically, we assess the ability of a simple dielectric continuum theory (DCT)-whose dielectric permittivity does not vary with depth z-to predict free energy differences when water lms are subjected to external elds. An advantage of this approach is that it is rooted in thermodynamics, which obviates the need to resolve uctuations/ response at the microscopic level. We will show that simple DCT with e(z) ¼ e bulk ¼ const. not only gives a good description of water's dielectric response under connement, but it also outperforms models that suppose a lower dielectric constant at the interface. Moreover, we also nd that for lms comprising just one or two layers of water molecules this simple DCT remains a remarkably reasonable approximation. We show that our analyses are broadly in line with the experimental observations reported in ref. 6.
The rest of the article is arranged as follows. In Sec. 2 we briey review linear response theory for dielectric uids, calling into question the notion of a permittivity that varies with position over microscopic scales. In Sec. 3 we analyze the polarization of a conned dielectric continuum under periodic boundary conditions, and derive a nite size correction for the thermodynamics of charging up a parallel plate capacitor. In Sec. 4 we use molecular simulations of simple point charge models to assess the accuracy of this correction, and compare extrapolated results with DCT predictions. We subsequently assess the performance of more complicated models in Sec. 5. In Sec. 6 we investigate the length scales at which DCT begins to fail. The sensitivity of the effective dielectric constant to the denition of lm thickness is discussed in Sec. 7. We summarize our ndings in Sec. 8.

Brief overview of dielectrics
In macroscopic DCT, the polarization P in a medium is related to the total electric eld E by the constitutive relation [22][23][24][25] 4pPðrÞ ¼ where e is the dielectric tensor, 1(r, r 0 ) ¼ Ud(r, r 0 ) with U the unit tensor, d(r, r 0 ) is Dirac's delta function, and the domain of integration is the volume occupied by the medium. Eqn (1) is a nonlocal relationship between P and E. There are two routes to arrive at the more familiar local relationship for a homogeneous, isotropic dielectric The rst is to simply assert locality i.e., e(r, r 0 ) ¼ e1(r, r 0 ). The second, more formal approach acknowledges the underlying molecular granularity, and supposes that e(r, r 0 ) is a shortranged function such that, e(r, r 0 ) x 0 for jr À r 0 j > l e . The characteristic length l e is determined by molecular correlations, and previous simulation studies suggest l e z 6Å. 26,27 If E varies slowly over distances comparable to l e , then the nonlocal relation (eqn (1)) reduces to the local one (eqn (2)), with Interfaces between different media are treated as innitely sharp boundaries within DCT. Any polarization in a medium then results in an induced surface charge s ind (x) ¼ P(x)$n(x) that occupies a region of innitesimal thickness, wheren(x) is the local surface normal. Such a scenario is, of course, an idealization of physical reality; 28 as discussed above, liquid interfaces have a nite length scale l int , which is on the order of 1 nm-or a few molecular diameters-for liquid water close to its triple point. The induced surface charge density is then understood to result from a physical charge distribution which is localized in the interfacial region, but with a thickness comparable to l int . Molecular simulations suggest that such interfacial charge distributions may vary rapidly along z. [9][10][11][12][13] Local dielectric constants obtained from simulation exhibit similar structure.
While it is reasonable to suppose that the properties of a material may differ in regions close to the interface compared to those in bulk, the notion of a local dielectric constant with variations on the molecular scale is unsettling in a couple of respects. First, in going from the nonlocal constitutive relation specied by eqn (1) to the local relation specied by eqn (2), we assumed that elds vary slowly over length scales comparable to l e z l int , so one might therefore question the appropriateness of a local dielectric constant. Second, even if one is content with the locality of e(z), DCT is a macroscopic theory, and the constitutive relations eqn (1) and (2) concern the macroscopic elds E and P. Obtaining these elds from the underlying microscopic degrees of freedom thus requires a coarse graining procedure, and it is reasonable to suppose that l e sets the minimum length scale over which any such coarse graining should be performed. Local molecular response functions that vary rapidly in space are likely important for the solvation and spatial distribution of ions, as well as electrokinetic phenomena; 10,11 it nonetheless remains challenging to reconcile variations of e(z) on the molecular scale with this viewpoint of relating coarse grained macroscopic elds (eqn (1) and (2)). By pursuing a thermodynamic perspective in this paper, which directly compares predictions of simple DCT to free energies obtained from molecular simulations, we avoid needing to compute the macroscopic elds E and P from microscopic degrees of freedom.

Using simple DCT as a finite size correction
The extent of the physical systems we have in mind are microscopic in one direction (perpendicular to the interface) but otherwise macroscopic. To represent them in computer simulations, we take the standard approach of imposing periodic boundary conditions in all three Cartesian directions. Our simulated system is thus an innite stack of thin water slabs, separated by substantial but still microscopic layers of vacuum, with an articial periodicity. Because electrostatic interactions are long in range, we anticipate nonnegligible quantitative consequences of this periodicity, particularly when its repeat length is not signicantly larger than the slab width.
To correct for such nite size effects, we adopt a strategy previously used to assess system size-dependence for ion solvation in similarly periodic slabs. Specically, we extend our work in ref. 29 to develop a nite size correction for the solvent contribution to the reversible work required to charge up a parallel plate capacitor under periodic boundary conditions (which we refer to as the 'solvation' free energy, f (L) solv ), based on the assumption that long-wavelength solvent response underlying nite size effects is well-described by DCT. 30 These predictions of DCT for charging parallel plates that bound thin water slabs serve simultaneously as a means to extrapolate computed free energies to the thermodynamic limit, and also as a test of the assumptions underlying DCT.
A representative snapshot of the system under consideration is shown in Fig. 1a. The parallel plate capacitor is approximated by two planes of N site point charges arranged on a square lattice, located at z ¼ AEw/2. The total charge of the plane at z ¼ w/2 is Q ¼ N site q site , which is equal-and-opposite to the plane at z ¼ Àw/2. The solvent water molecules are conned between these two charged planes by tightly packed volume-excluding Weeks-Chandler-Anderson 31 (WCA) particles (see Sec. 9). In most of what follows, the WCA centers and the point charges coincide, though we will also consider more general cases like those depicted in Fig. 1a. We now make two continuum approximations. First, water is treated as a dielectric slab with dielectric constant e, spanning z ¼ Àl w /2 to z ¼ +l w /2, as indicated in Fig. 1b. A value of l w appropriate to our molecular system is not a priori obvious: the WCA particles enforce very low density of oxygen atoms outside a region Àw/2 < z < w/2; given that water molecules are not point particles, however, the most realistic continuum description could involve an offset d between w and l w , i.e., l w ¼ w À d. Considerations for choosing d will be discussed later.
The two charged planes at z ¼ AEw/2 are treated in our continuum calculation as uniformly charged sheets with surface charge density q^Q/A, where A is the cross-sectional area of the simulation cell orthogonal to z. Within DCT, these charged planes enter the continuum model explicitly by introducing a discontinuity of magnitude 4pjqj in the total electric eld along z (as the planes are surrounded on either side by vacuum), irrespective of whether they are coincident with the WCA particles. In contrast, the WCA centers only enter DCT implicitly by conning the water molecules such that the thickness of the dielectric slab is l w^w À d. The continuum representation of the system is summarized in Fig. 1b. The simulation cell is periodically replicated in all three dimensions, and the periodic length along the z-direction is L.
In the ESI, † we solve the periodic continuum problem shown schematically in Fig. 1b, obtaining a total electrostatic potential in the region Àl w /2 # z # l w /2 where P is the uniform polarization of the dielectric, and we have assumed that an Ewald-style approach has been used to treat electrostatic interactions. The rst term in eqn (3) arises from the charged planes, which we denote f q . The second term arises from the polarized dielectric, and we denote this f solv . The total electric eld inside the dielectric follows directly from eqn (3): We now combine eqn (4) with the local constitutive relation (eqn (2)) to obtain an expression for P: We also show in the ESI † that the electrostatic potential at the charged plane at z ¼ Àw/2 is Similarly, for the charged plate at z ¼ +w/2 we have The solvation free energy f (L) solv ¼ q(f solv,hi À f solv,lo )/2 is the difference in reversible work (per unit area) to introduce the surface charge density q to the charged planes with and without the solvent present. Combining eqn (5)-(7) gives f ðLÞ solv ðqÞ ¼ À2pq 2 ðw À dÞ in the limit L / N we recover the expected result f ðNÞ solv ðqÞ ¼ À2pq 2 e À 1 e ðw À dÞ: the correction Df DCT (L) ¼ f (N) solv À f (L) solv we should apply for nite L is thus Eqn (10) provides a simple correction term that can be added to f (L) solv obtained from molecular simulations. The extent to which Df DCT (L) achieves consistent estimates of f (N) solv from simulations with different L is then one indicator of how well simple DCT describes the dielectric properties of water lms.

Assessing the accuracy of DCT with molecular simulations
To assess our continuum prediction of the nite size correction Df DCT (L) given by eqn (10), we will assume that e retains its bulk liquid value (e liq z 71 for SPC/E) over the entire domain Àl w /2 < z < l w /2. The only undetermined parameter in eqn (10) is then the length scale d, which determines the location of the dielectric boundaries of the solvent relative to the charged planes at z ¼ AEw/2. To determine an appropriate value of d, we note that DCT predicts an electric eld due to the solvent in the region w/2 # z < L/2 As shown in Fig. 2a, E solv can be easily obtained from simulation. (Note that, owing to the charge asymmetric distribution of individual water molecules, f solv (z) for liquid water varies across the interface even with q ¼ 0 eÅ À2 . As we are concerned with the response of the dielectric slab, in Fig. 2a we For the time being, we consider cases where the charged planes and WCA centers coincide. For a given w, we then measure E solv with q z 3 Â 10 À3 eÅ À2 for each value of L investigated, and by substituting P given by eqn (5) into eqn (11), we determine d. Results obtained with different w (see Table 1) are shown in Fig. 2b. Despite some noise, d appears to plateau as w increases; averaging results for w $ 20Å, we nd d ¼ 2.09 AE 0.17Å. This procedure is similar in spirit and effect to that of ref. 19, which locates a dielectric dividing surface based on the average potential drop across a polarized water slab. Our approach does not assign special signicance to the potential at the conning walls. More signicantly, we nd that d decreases for sub-nanometer lms, in contrast to the increase reported by ref. 19 for water between graphene sheets.
Having determined d, we are now in a position to test the appropriateness of the nite size correction given by eqn (10). To this end, in Fig. 3a and b, we show f (L) solv (q) for w ¼ 40Å and w ¼ 20Å, respectively. We focus on these values of w as they correspond to the extremal values investigated that lie in the plateau region in Fig. 2b; results for w ¼ 30Å and w ¼ 25Å are included in the ESI. † As expected, f (L) solv (q) exhibits a dependence on system size. Adding Df DCT (L) removes this dependence almost entirely, as seen in Fig. 3c and d. Also shown are results obtained by imposing vanishing electric displacement eld D ¼ 0 VÅ À1 along z, which is formally equivalent to the commonly used Yeh-Berkowitz approach for approximating 2D Ewald summation. 32,33 As results obtained with D ¼ 0 VÅ À1 should approximate L / N, they do not require a nite size correction. Importantly, excellent agreement with f (L) solv + Df DCT (L) is observed, giving us condence that eqn (10) provides a meaningful nite size correction.  The fact that the simple DCT model outlined in Sec. 3 describes the nite size behavior of f (L) solv so well suggests it is reasonable to think of thin water lms as having a uniform dielectric constant equal to that of bulk in the region they occupy. Even more tellingly, the extrapolation f (L) solv + Df DCT (L) from simulation agrees well with the continuum prediction f (N) solv in eqn (9).
To provide a physical interpretation for the length scale d, Fig. 4 shows number density proles r(z) for water's oxygen and hydrogen atoms from simulations with q ¼ 0 eÅ À2 . On these plots, we have also marked the boundary predicted by l w /2 ¼ (w À d)/2, which corresponds closely to the vanishing of average hydrogen density. Because the hydrogen atoms protrude further toward the vapor phase than the oxygen atoms, this boundary marks the outermost limit of microscopic sources of polarization uctuations. The water lm thickness we have inferred is thus the largest that could be reasonably justied based on the statistics of molecular congurations.

Assessing the validity of other models
We have shown that f (L) solv + Df DCT (L) obtained from simulation agrees well with the predictions of a simple DCT in which the dielectric constant of thin lms is identical to that of the bulk liquid (eqn (9)). If we were to decrease e, agreement with simulation data would require assigning l w a larger value than we have inferred, i.e., a value that would be difficult to justify from microscopic structure. This observation advocates against the notion that the overall dielectric permittivity of the thin lm is lower than in the homogeneous uid. By itself, however, it does not rule out a model in which the interfacial regions have a permittivity e int that is distinct from the bulk region they sandwich. For such a model, the free energy reads where l int is the width of each interfacial region. Following the dielectric imaging experiments of Fumagalli et al., 6 we take l int ¼ 7.5 AE 1.5Å and e int ¼ 2.1, and require the total width l w ¼ l bulk + 2l int to have the same value as in the uniform dielectric model: as discussed above, it is unreasonable to allow l w to increase from that value. Decreasing l w , on the other hand, offers less exibility to a model that introduces regions of low dielectric constant at the expense of those with high dielectric constant. The resulting predictions of f (N) solv,int are shown in Fig. 3c and d (labeled "bulk + interface"), where poor agreement with the simulation data is observed. Quantitatively solv,int from a dielectric continuum model, in which an interfacial layer of width l int ¼ 7.5Å is assigned a permittivity e int ¼ 2.1 much lower than in bulk liquid, computed from (eqn (12)). The shaded regions bound predictions with 6Å # l int # 9Å. Insets: snapshots from corresponding molecular dynamics simulations. different (but not signicantly improved) predictions would be obtained with different choices of l int and e int . We nd generally that l int ¼ 0 (or equivalently, e int ¼ e liq ) yields the best agreement with simulation. Evidence for this conclusion is provided in ESI. † The width l int of a notional interfacial layer differs fundamentally from the length scale d in our simple uniform DCT. They can nonetheless easily be confused. In our case z ¼ AE(w À d)/2 marks the location of a sharp interface between vapor and bulk liquid. This interface does not coincide with the location z ¼ AEw/2 of the conning charged walls because their constituent WCA particles exclude volume. d thus characterizes a region that is inaccessible to water molecules and should not be associated with the liquid. To emphasize this point, we modify the w ¼ 20Å system (Fig. 3b and d) by displacing the charged planes 5Å into vacuum, with the WCA particles xed at their original positions (i.e., the general case considered in Fig. 1a). d increases by 10Å as a result, while l w ¼ w À d is unchanged, i.e., w / 30Å and d / 12.09Å, while l w ¼ 17.91Å just as before. Changing d in this fashion clearly has nothing to do with water's interfacial dielectric properties. Fig. 5 presents results for f (L) solv and f (L) solv + Df DCT (L) for the displaced-charge system, which are virtually indistinguishable from their undisplaced counterparts in Fig. 3.
By contrast, a layer of width l int in "bulk + interface" models is clearly associated with the liquid. It is imagined to comprise water molecules whose orientational uctuations are distinct from those in bulk liquid due to the phase boundary. Multiple studies based on such models have concluded that the interfacial layer has a greatly reduced polarizability, amounting to a "dead layer" with e int z 1. 6,14,16,17 Dielectric properties of this notional dead layer may be nearly indistinguishable from vacuum, but the layer plainly belongs to the dense liquid phase within a "bulk + interface" picture.

Ultra thin films of water
We have established so far that lms of water with l w T 18Å behave quantitatively like simple dielectric continua with regard to their response to a uniform electric eld. We now investigate the behavior of 'ultra thin' lms conned between charged plates with w # 10Å. In Fig. 6a-c, we show f (N) solv obtained with simulation for w ¼ 5Å, 7.5Å, and 10Å, respectively, using d ¼ 2.09 AE 0.17Å to correct for nite size effects (eqn (10)). As the thickness of the water slab is reduced to length scales comparable to l e z 6Å, treating the water molecules as a dielectric continuum is certainly questionable. Discrepancies between f (N) solv (eqn (9)) and f (L) solv + Df DCT (L) indeed become apparent as w is decreased below 1 nm, but the relative  Fig. 2). In both cases, the dielectric boundary aligns closely with the vanishing of hydrogen atom density. error of continuum predictions is surprisingly modest. Even when w is only large enough to accommodate a single molecular layer (Fig. 6a), the continuum prediction in eqn (9) provides a reasonable ballpark estimate of the solvation free energy. For two to three molecular layers (Fig. 6b and c), quantitative agreement between simple DCT and the simulation data is recovered almost entirely.

Reconciling our results with dielectric imaging experiments
The conclusion we have drawn from computer simulationsthat the dielectric response of nanoscopically thin water lms can be anticipated from bulk properties alone-is squarely at odds with the conclusion drawn by Fumagalli et al. 6 based on dielectric imaging measurements of conned water. In this section we attempt to reconcile our results with those measurements. Assuming that our simple uniform continuum model is correct, we show how uncertainty in the thickness of a water lm can cause the apparent dielectric constant e app to depend sensitively on lm thickness. More specically, we assess the consequences of assigning a width l w + d to a lm whose actual thickness is l w , i.e., failing to account for the volume excluded by a conning substrate. Based on this assignment and the development in Sec. 3, we would expect a solvation free energy f ðNÞ solv;app ðqÞ ¼ À2pq 2 e app À 1 e app w: Equating this with the "true" free energy in eqn (9), we obtain which depends explicitly on the lm's thickness. While similar functional forms for e app (w) have been reported previously, 14,16 the physical interpretation here is different: As discussed in Sec. 5, d is not to be associated with the properties of interfacial water. Fig. 7 plots the apparent permittivity e app in eqn (14) as a function of w. Here we have set e ¼ e liq ¼ 80 and estimated d z 3.5Å for the graphene-water interface (based on density proles obtained from ab initio molecular dynamics simulations 34 ). Diminished values of e app at small w could easily, but mistakenly, be taken to signify a strong suppression of polarization uctuations and response in nanoscale water lms.
The inference of suppressed interfacial permittivity from experiments may suffer from the same issues that cause e app to depend strongly on lm thickness, much as suggested by ref. 19. To emphasize this point, in Fig. 7 we include dielectric imaging data from ref. 6, which exhibit a very similar dependence on w. As an important caveat, the samples studied by solv (q) predicted by DCT (eqn (9)), and the shaded region encompasses predictions with d ¼ 2.09 AE 0.17Å. Insets: snapshots from molecular dynamics simulations. Fig. 7 The apparent dielectric constant e app predicted by simple DCT (eqn (14)) is broadly consistent with experiment. The blue solid line is obtained with d ¼ 3.5Å, which is estimated from ab initio molecular dynamics simulations of water on graphene. The blue shaded region indicates the range of e app obtained with d ¼ 3.5 AE 2.0Å, demonstrating the sensitivity of e app to uncertainty in the film thickness. The light and dark gray dotted lines indicate e app ¼ 2.1 and e app ¼ 80, respectively.
Fumagalli et al. have a more complicated geometry than the simple "slit-pore" scenario we have considered, involving an AFM tip, multiple water channels, a graphite substrate, and hexagonal boron nitride walls. Since geometry of the dielectric boundary is precisely the issue under scrutiny here, the comparison between theory and experiment suggested by Fig. 7 should be made cautiously and only qualitatively. In our view it nonetheless suggests that the correct interpretation of measurements in ref. 6 may not in fact require invoking an interfacial dead layer.

Summary
In this article, we have probed the dielectric response of thin water lms using molecular simulations with nite size corrections from DCT. We specically calculated the solvent contribution to the reversible work of introducing charge to parallel plates that conne the lm. Our results demonstrate that response to such slowly varying external elds can be accurately captured with a DCT model whose permittivity is simply equal to that of the homogeneous liquid, in the entire region occupied by the liquid. Our analysis reveals that appropriate dielectric boundaries for these lms extend to the point where the hydrogen number density approximately vanishes, and thus incorporate all microscopic sources of polarization uctuations. This observation is consistent with our recent study, where we found that the dielectric boundary between water and spherical solutes was reasonably described by the rst peak in the solute-hydrogen radial distribution function. 35 Within this simple DCT approach, which achieves quantitative agreement with simulation for lms T1 nm, water's interfacial regions do not enter as separate domains. Loche et al. 19 have similarly concluded that water lms a few nanometers in thickness are well characterized by bulk dielectric parameters, but they report substantial deviations at smaller scales.
We also demonstrated rough consistency with experiments 6 that had previously been interpreted to imply a dielectrically dead layer of water at the liquid's boundary. This agreement is achieved by asserting that dielectric boundaries had previously not been placed appropriately. For the simple point charge model used in this study, the point where the hydrogen number density approximately vanishes coincides with the point where microscopic charge density vanishes. For polarizable models or ab initio treatments of water, it is possible that the distribution of electron density beyond the hydrogen atoms also plays a role. 36 Further investigation of this point is le for future work.
Our conclusion is further supported by results for subnanometer water lms that comprise one to three molecular layers. In these cases, one cannot sensibly discuss a bulk region, yet the simple DCT model still performs remarkably well. If anything, the apparent dielectric constant would need to increase to improve agreement with the free energy data. This result contrasts with conclusions of ref. 12 and 19, which suggest greatly suppressed polarizability at comparable scales of connement.
To be certain, applying simple DCT to these subnanometer lms cannot be carefully justied (see Sec. 2). Nonetheless, this ad hoc application of simple DCT to small length scales strongly argues against the notion of an interfacial region with low dielectric constant. Our conclusion is also in line with previous studies that have found corrections similar to Df DCT (L) for the solvation of small spherical ions in water work remarkably well down to the nanometer scale, for both bulk 30,37,38 and interfacial systems. 29 Similarly, the simple DCT model described in this article has also been found to accurately capture mean eld-like corrections for thin lms of water where electrostatic interactions are treated in a short-ranged fashion. 27 The approach we have taken is not well suited to address a free liquid-vapor interface, whose substantial topographical uctuations greatly complicate formulating and solving an appropriate dielectric boundary value problem. Previous work on the free interface has emphasized that in representative congurations the liquid phase terminates sharply at any given lateral position. 39,40 Since our results stress the importance of precisely establishing the liquid's microscopic boundary, we expect that dielectric models based on a smoothed average interface are a poor caricature of this scenario. Instead, a faithful assessment of permittivity at the free liquid-vapor interface will require attention to its undulating instantaneous structure, suggesting that a more spatially localized analysis is likely more feasible.
It would be incorrect, however, to conclude from our results that simple DCT gives a full account of polarization response at the liquid's surface. Indeed, even in bulk liquid water, the charging of small spherical solutes is characterized by shortwavelength solvent response that is not well-described by simple continuum approaches. 37 In such cases, and in contrast to the thin lms considered in this work, f (N) solv predicted by DCT is a poor estimate of that obtained from simulations. 29,35 The impact of such profound perturbations are even more pronounced for solutes near so interfaces like that between water and its vapor, where distortions of the interface result in nonlinearities beyond the scope of current theoretical treatments. 41 But for perturbations that vary slowly in space, like the uniform elds considered here, the results of this study add to a growing body of work that supports a surprisingly simple view of water's surface (and, by extension, thin lms) as a dielectric medium: its local permittivity is equivalent to the bulk dielectric constant, all the way down to nanometer length scales.

Methods
All simulations followed the basic setup shown in Fig. 1a. Two planes of N site ¼ 100 point charges were placed on a square lattice at z ¼ AEw/2. Water molecules, modeled with the SPC/E potential, 42 were conned to the region Àw/2 # z # w/2 by volume-excluding WCA particles, 31 whose centers, for the most part, coincided with the point charges at z ¼ AEw/2. The interaction between an individual WCA particle and a water molecule is dened by where s ¼ 2.5Å, 3 ¼ 0.1 kcal mol À1 , and r is the distance between the WCA particle and the oxygen atom of the water molecule. As described in Sec. 5, for w ¼ 20Å, we also performed simulations where the planes of point charges were displaced 5Å into vacuum, but leaving the rest of the system unchanged. For each value of w, simulations of 5 ns (following at least 100 ps of equilibration) were performed with q site /e ¼ 0, 1 Â 10 À3 , ., 5 Â 10 À3 . The total volume of the simulation cells was 12.75Å Â 12.75Å Â L, where L takes values as indicated throughout the article. All simulations were performed with the LAMMPS simulations package. 43 Lennard-Jones interactions between water molecules were truncated and shied at 10Å. Long range electrostatic interactions were evaluated using particle-particle particle-mesh Ewald summation, 44 with parameters chosen such that the RMS error in the forces were a factor 10 5 smaller than the force between two unit charges separated by a distance of 1.0Å. 45 Where indicated in the text, the electric displacement eld along z was set to zero, using the implementation given in ref. 46. The geometry of the water molecules was constrained using the RATTLE algorithm. 47 Temperature was maintained at 298 K with a Nosé-Hoover chain thermostat 48,49 with a damping constant of 100 fs. A time step of 2 fs was used throughout. The number of water molecules used in the simulations is given in Table 1.
The free energy of charging up parallel plate capacitors was computed by averaging electric potentials appropriately. Let 4 (i) solv,hi (R N ) and 4 (j) solv,lo (R N ) denote the instantaneous electric potentials due to the solvent with conguration R N at site i of one of the point charges in the plane at z ¼ w/2, and site j in the plane at z ¼ Àw/2, respectively. The total solvation free energy F (L) solv (Q) is then dened by P (L) (D4 solv ; Q) is the probability distribution of D4 solv in the presence of two charged planes with total charges AE Q, in a simulation box of length L; and h$i (L) 0 denotes an average over P (L) (D4 solv ; 0). Similar to our previous studies, 29,35 F solv (Q) was computed using the MBAR algorithm. 50 The solvation free energies per unit area, f (L) solv that we consider are then obtained by dividing F solv (Q) by the cross-sectional area of the simulation cell.

Data availability
The data that supports the ndings of this study and input les for the simulations are openly available at the University of Cambridge Data Repository, https://doi.org/10.17863/ CAM.81959.

Author contributions
S. J. C. and P. L. G. conceived and conducted research. Both authors wrote the paper, at all stages.

Conflicts of interest
There are no conicts to declare.