# Dielectric response of thin water films: A thermodynamic perspective

Stephen J. Cox  Yusuf Hamied Department of Chemistry, University of Cambridge,
Lensfield Road, Cambridge CB2 1EW, United Kingdom<EMAIL_ADDRESS>Phillip
L. Geissler  Chemical Sciences Division, Lawrence Berkeley National
Laboratory, Berkeley, CA 94720, United States.  Department of Chemistry,
University of California, Berkeley, CA 94720, United States.
<EMAIL_ADDRESS>
###### Abstract

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 ($\sim 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.

##  I Introduction

Interest in confined water has exploded over the last decade or so, owing
principally to advances in the fabrication of devices at the nanoscale Geim
and Grigorieva (2013); Celebi  _et al._ (2014); Jain  _et al._ (2015), the
potential implications for ‘blue energy’ and desalination Bocquet (2020), and
as means to understand fundamental properties of water Algara-Siller  _et al._
(2015); Fumagalli  _et al._ (2018) and its solutions Nakamuro  _et al._
(2021); Wang  _et al._ (2021). An obvious consequence of the decreasing length
scales associated with confinement is an increase in the surface-to-volume
ratio of liquid water, which typically amplifies surface-specific effects
relative to large sample geometries. The notion of nanoconfined 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 fluid. In particular, many years of
investigation along these lines Ballenegger and Hansen (2005); Bonthuis,
Gekle, and Netz (2011, 2012); Schlaich, Knapp, and Netz (2016); Loche  _et
al._ (2018); Zhang (2018); Olivieri, Hynes, and Laage (2021); Matyushov
(2021); Mondal and Bagchi (2021) has focused on the static dielectric constant
$\epsilon_{\rm 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 Fumagalli  _et al._ (2018) of water confined between two
atomically flat walls separated by distances as small as 0.8 nm. These imaging
results were inferred to report an interfacial dielectric constant
$\epsilon_{\rm int}=2.1$ (relevant to an interfacial region of thickness
$\ell_{\rm int}\approx 7.5\pm 1.5$ Å) that dominates the capacitance of a thin
water film. This value, typical of a bulk nonpolar liquid, signifies a
dramatic departure from the polarizability of bulk water, for which
$\epsilon_{\rm liq}\approx 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 $\ell_{\rm int}$ of
the surface Rowlinson and Widom (2002). Accordingly, many studies have aimed
to rationalize confined water’s electrostatic response in terms of a local
dielectric constant $\epsilon(z)$ that varies with position $z$ along the
surface normal Ballenegger and Hansen (2005); Bonthuis, Gekle, and Netz (2011,
2012); Schlaich, Knapp, and Netz (2016); Loche  _et al._ (2018, 2020); Jalali
_et al._ (2020); Hamid  _et al._ (2021). Molecular simulations have estimated
$\epsilon(z)$ either from polarization fluctuations, or from response to
external electric fields; in either case this approach relies upon
interpreting features that have been resolved at a fine scale within a
theoretical framework appropriate for macroscopic dielectric materials. In
this study, we pursue a different approach. Specifically, 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 films are subjected to external fields. An advantage of this
approach is that it is rooted in thermodynamics, which obviates the need to
resolve fluctuations/response at the microscopic level. We will show that
simple DCT with $\epsilon(z)=\epsilon_{\rm bulk}=\mathrm{const.}$ not only
gives a good description of water’s dielectric response under confinement, but
it also outperforms models that suppose a lower dielectric constant at the
interface. Moreover, we also find that for films 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. Fumagalli  _et al._ , 2018.

The rest of the article is arranged as follows. In Sec. II we briefly review
linear response theory for dielectric fluids, calling into question the notion
of a permittivity that varies with position over microscopic scales. In Sec.
III we analyze the polarization of a confined dielectric continuum under
periodic boundary conditions, and derive a finite size correction for the
thermodynamics of charging up a parallel plate capacitor. In Sec. IV 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. V. In
Sec. VI we investigate the length scales at which DCT begins to fail. The
sensitivity of the effective dielectric constant to the definition of film
thickness is discussed in Sec. VII. We summarize our findings in Sec. VIII.

##  II Brief overview of dielectrics

In macroscopic DCT, the polarization $\mathbf{P}$ in a medium is related to
the total electric field $\mathbf{E}$ by the constitutive relation Fulton
(1978a, b, 1983); Caillol (1992)

$4\pi\mathbf{P}(\mathbf{r})=\int\\!\mathrm{d}\mathbf{r}^{\prime}\,\left[{\bm{\epsilon}}(\mathbf{r},\mathbf{r}^{\prime})-\mathbf{1}(\mathbf{r},\mathbf{r}^{\prime})\right]\cdot\mathbf{E}(\mathbf{r}^{\prime}),$
(1)

where ${\bm{\epsilon}}$ is the dielectric tensor,
$\mathbf{1}(\mathbf{r},\mathbf{r}^{\prime})=\mathbf{U}\delta(\mathbf{r},\mathbf{r}^{\prime})$
with $\mathbf{U}$ the unit tensor, $\delta(\mathbf{r},\mathbf{r}^{\prime})$ is
Dirac’s delta function, and the domain of integration is the volume occupied
by the medium. Equation 1 is a nonlocal relationship between $\mathbf{P}$ and
$\mathbf{E}$. There are two routes to arrive at the more familiar local
relationship for a homogeneous, isotropic dielectric

$4\pi\mathbf{P}(\mathbf{r})=(\epsilon-1)\mathbf{E}(\mathbf{r}).$ (2)

The first is to simply assert locality i.e.,
${\bm{\epsilon}}(\mathbf{r},\mathbf{r}^{\prime})=\epsilon\mathbf{1}(\mathbf{r},\mathbf{r}^{\prime})$.
The second, more formal approach acknowledges the underlying molecular
granularity, and supposes that
${\bm{\epsilon}}(\mathbf{r},\mathbf{r}^{\prime})$ is a short-ranged function
such that, ${\bm{\epsilon}}(\mathbf{r},\mathbf{r}^{\prime})\simeq 0$ for
$|\mathbf{r}-\mathbf{r}^{\prime}|>\ell_{\epsilon}$. The characteristic length
$\ell_{\epsilon}$ is determined by molecular correlations, and previous
simulation studies suggest $\ell_{\epsilon}\approx 6$ ÅZhang, Hutter, and
Sprik (2016); Cox (2020). If $\mathbf{E}$ varies slowly over distances
comparable to $\ell_{\epsilon}$, then the nonlocal relation (Eq. 1) reduces to
the local one (Eq. 2), with

${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\epsilon\mathbf{U}=\int\\!\mathrm{d}\mathbf{r}^{\prime}\,{\bm{\epsilon}}(\mathbf{r},\mathbf{r}^{\prime}).}$

Interfaces between different media are treated as infinitely sharp boundaries
within DCT. Any polarization in a medium then results in an induced surface
charge $\sigma_{\rm
ind}(\mathbf{x})=\mathbf{P}(\mathbf{x})\cdot\hat{\mathbf{n}}(\mathbf{x})$ that
occupies a region of infinitesimal thickness, where
$\hat{\mathbf{n}}(\mathbf{x})$ is the local surface normal. Such a scenario
is, of course, an idealization of physical reality Jackson (1999); as
discussed above, liquid interfaces have a finite length scale $\ell_{\rm
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 $\ell_{\rm int}$.
Molecular simulations suggest that such interfacial charge distributions may
vary rapidly along $z$ Ballenegger and Hansen (2005); Bonthuis, Gekle, and
Netz (2011, 2012); Schlaich, Knapp, and Netz (2016); Loche  _et al._ (2018).
Local dielectric constants obtained from simulation exhibit similar structure.

Figure 1: Molecular and continuum representations of the system considered.
(a) Water molecules (oxygen atoms in blue) are confined between volume-
excluding WCA particles (light gray). Dark gray circles represent point
charges: negative on the left, positive on the right, separated by a distance
$w$. (b) In the continuum representation, these planes of point charges are
approximated as uniformly charged sheets, as indicated by the dashed dark gray
lines. The effect of the WCA particles enters implicitly by bounding the
solvent, itself represented as a continuum with dielectric constant
$\epsilon$, within a slab of thickness $\ell_{\rm w}=w-\delta$, where
$\delta/2$ indicates the distance between the solvent-vapor dielectric
boundary, and the charged planes. In both (a) and (b), the simulation cell is
periodically replicated in all three dimensions, and its length in the
direction normal to the charged planes is $L$.

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 specified by Eq. 1 to the local relation specified by Eq. 2, we
assumed that fields vary slowly over length scales comparable to
$\ell_{\epsilon}\approx\ell_{\rm int}$, so one might therefore question the
appropriateness of a local dielectric constant. Second, even if one is content
with the locality of $\epsilon(z)$, DCT is a _macroscopic_ theory, and the
constitutive relations Eqs. 1 and 2 concern the macroscopic fields
$\mathbf{E}$ and $\mathbf{P}$. Obtaining these fields from the underlying
microscopic degrees of freedom thus requires a coarse graining procedure, and
it is reasonable to suppose that $\ell_{\epsilon}$ 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 Bonthuis, Gekle, and Netz (2011, 2012); it nonetheless remains
challenging to reconcile variations of $\epsilon(z)$ on the molecular scale
with this viewpoint of relating coarse grained macroscopic fields (Eqs. 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 fields $\mathbf{E}$
and $\mathbf{P}$ from microscopic degrees of freedom.

##  III 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 infinite stack of thin water slabs, separated by
substantial but still microscopic layers of vacuum, with an artificial
periodicity. Because electrostatic interactions are long in range, we
anticipate nonnegligible quantitative consequences of this periodicity,
particularly when its repeat length is not significantly larger than the slab
width.

To correct for such finite size effects, we adopt a strategy previously used
to assess system size-dependence for ion solvation in similarly periodic
slabs. Specifically, we extend our work in Ref. Cox and Geissler, 2018 to
develop a finite 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_{\rm solv}^{(L)}$), based on the assumption that long-wavelength
solvent response underlying finite size effects is well-described by DCT
Hünenberger and McCammon (1999). 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_{\rm
site}$ point charges arranged on a square lattice, located at $z=\pm w/2$. The
total charge of the plane at $z=w/2$ is $Q=N_{\rm site}q_{\rm site}$, which is
equal-and-opposite to the plane at $z=-w/2$. The solvent water molecules are
confined between these two charged planes by tightly packed volume-excluding
Weeks-Chandler-Anderson Weeks, Chandler, and Andersen (1971) (WCA) particles
(see Sec. IX). 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 $\epsilon$, spanning
$z=-\ell_{\rm w}/2$ to $z=+\ell_{\rm w}/2$, as indicated in Fig. 1b. A value
of $\ell_{\rm 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
$\delta$ between $w$ and $\ell_{\rm w}$, i.e., $\ell_{\rm w}=w-\delta$.
Considerations for choosing $\delta$ will be discussed later.

The two charged planes at $z=\pm w/2$ are treated in our continuum calculation
as uniformly charged sheets with surface charge density $q\equiv 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 $4\pi|q|$ in the total electric field
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 confining the water
molecules such that the thickness of the dielectric slab is $\ell_{\rm
w}\equiv w-\delta$. 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 Supplementary Material (SM), we solve the periodic continuum problem
shown schematically in Fig. 1b, obtaining a total electrostatic potential in
the region $-\ell_{\rm w}/2\leq z\leq\ell_{\rm w}/2$

$\phi(z)=4\pi q\bigg{(}-\frac{zw}{L}+z\bigg{)}+4\pi
P\bigg{(}-\frac{z(w-\delta)}{L}+z\bigg{)},$ (3)

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 first term in Eq. 3 arises from the charged planes, which we
denote $\phi_{q}$. The second term arises from the polarized dielectric, and
we denote this $\phi_{\rm solv}$. The total electric field inside the
dielectric follows directly from Eq. 3:

$E=-4\pi q\bigg{(}1-\frac{w}{L}\bigg{)}-4\pi
P\bigg{(}1-\frac{w-\delta}{L}\bigg{)}.$ (4)

We now combine Eq. 4 with the local constitutive relation (Eq. 2) to obtain an
expression for $P$:

$P=-\frac{(\epsilon-1)(1-\frac{w}{L})q}{1+(\epsilon-1)(1-\frac{w-\delta}{L})}.$
(5)

We also show in the SM that the electrostatic potential at the charged plane
at $z=-w/2$ is

$\phi_{\rm solv,lo}=2\pi P\bigg{(}\frac{w(w-\delta)}{L}-(w-\delta)\bigg{)}.$
(6)

Similarly, for the charged plate at $z=+w/2$ we have

$\phi_{\rm solv,hi}=2\pi P\bigg{[}-\frac{w(w-\delta)}{L}+(w-\delta)\bigg{]}.$
(7)

The solvation free energy $f^{(L)}_{\rm solv}=q(\phi_{\rm solv,hi}-\phi_{\rm
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 Eqs. 5, 6 and 7 gives

$f^{(L)}_{\rm solv}(q)=-2\pi
q^{2}(w-\delta)\frac{(\epsilon-1)(1-\frac{w}{L})^{2}}{1+(\epsilon-1)(1-\frac{w-\delta}{L})}.$
(8)

In the limit $L\to\infty$ we recover the expected result

$f^{(\infty)}_{\rm solv}(q)=-2\pi q^{2}\frac{\epsilon-1}{\epsilon}(w-\delta).$
(9)

The correction $\Delta f_{\rm DCT}(L)=f^{(\infty)}_{\rm solv}-f^{(L)}_{\rm
solv}$ we should apply for finite $L$ is thus

$\Delta f_{\rm DCT}(L)=2\pi
q^{2}(w-\delta)(\epsilon-1)\Bigg{[}\frac{(1-\frac{w}{L})^{2}}{1+(\epsilon-1)(1-\frac{w-\delta}{L})}-\frac{1}{\epsilon}\Bigg{]}.$
(10)

Equation 10 provides a simple correction term that can be added to
$f^{(L)}_{\rm solv}$ obtained from molecular simulations. The extent to which
$\Delta f_{\rm DCT}(L)$ achieves consistent estimates of $f^{(\infty)}_{\rm
solv}$ from simulations with different $L$ is then one indicator of how well
simple DCT describes the dielectric properties of water films.

##  IV Assessing the accuracy of DCT with molecular simulations

To assess our continuum prediction of the finite size correction $\Delta
f_{\rm DCT}(L)$ given by Eq. 10, we will assume that $\epsilon$ retains its
bulk liquid value ($\epsilon_{\rm liq}\approx 71$ for SPC/E) over the entire
domain $-\ell_{\rm w}/2<z<\ell_{\rm w}/2$. The only undetermined parameter in
Eq. 10 is then the length scale $\delta$, which determines the location of the
dielectric boundaries of the solvent relative to the charged planes at $z=\pm
w/2$. To determine an appropriate value of $\delta$, we note that DCT predicts
an electric field due to the solvent in the region $w/2\leq z<L/2$

$E_{\rm solv}=4\pi P\frac{(w-\delta)}{L}.$ (11)

As shown in Fig. 2a, $E_{\rm solv}$ can be easily obtained from simulation.
(Note that, owing to the charge asymmetric distribution of individual water
molecules, $\phi_{\rm 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 have plotted $\Delta_{0}\phi_{\rm
solv}(z)\equiv\phi_{\rm solv}(z)-\phi_{\rm solv,0}(z)$, where $\phi_{\rm
solv,0}(z)$ is the average electric potential profile with $q=0$ $e$/Å2.)

For the time being, we consider cases where the charged planes and WCA centers
coincide. For a given $w$, we then measure $E_{\rm solv}$ with $q\approx
3\times 10^{-3}\,e/{\rm\AA}^{2}$ for each value of $L$ investigated, and by
substituting $P$ given by Eq. 5 into Eq. 11, we determine $\delta$. Results
obtained with different $w$ (see Table 1) are shown in Fig. 2b. Despite some
noise, $\delta$ appears to plateau as $w$ increases; averaging results for
$w\geq 20$ Å, we find $\delta=2.09\pm 0.17$ Å. This procedure is similar in
spirit and effect to that of Ref. Loche  _et al._ , 2020, which locates a
dielectric dividing surface based on the average potential drop across a
polarized water slab. Our approach does not assign special significance to the
potential at the confining walls. More significantly, we find that $\delta$
decreases for sub-nanometer films, in contrast to the increase reported by
Ref. Loche  _et al._ , 2020 for water between graphene sheets.

Figure 2: (a) Average electrostatic potential due to the solvent (solid blue)
and charged planes (dashed orange) with $q\approx 3\times 10^{-3}$ $e$/Å2,
$w=30$ Å and $L=120$ Å. The vertical dotted lines indicate the positions of
the WCA particles. For the solvent, we plot $\Delta_{0}\phi_{\rm
solv}(z)=\phi_{\rm solv}(z)-\phi_{\rm solv,0}(z)$, where $\phi_{\rm solv,0}$
is the average potential with $q=0$ $e$/Å2. The average electric field due to
the solvent is used to determine $\delta$. (b) The inferred displacement
$\delta/2$ between WCA particles and the dielectric boundary depends only
weakly on the width of the liquid slab. Each point in the plot of $\delta$ vs.
$w$ is the average of 5 simulations with different values of $L$. Averaging
results for $w\geq 20$ Å gives $\delta=2.09\pm 0.17$ Å, which is used
throughout. The shaded orange region indicates a 95% confidence interval. Note
that in both (a) and (b), results have been obtained from simulations where
the WCA particles and charged planes coincide. Figure 3: Dependence of
solvation free energy $f_{\rm solv}^{(L)}(q)$ on system size $L$, shown in (a)
and (b) for $w=40$ Å and $w=20$ Å, respectively. The values of $L$ for $w=40$
Å are are indicated in the legend of panel (a); those for the thinner liquid
slab are shown in (b). In both cases the WCA particles coincide with the
charged planes. Adding $\Delta f_{\rm DCT}(L)$ given by Eq. 10 largely removes
this sensitivity, as seen in (c) and (d) for $w=40$ Å and $w=20$ Å,
respectively. DCT predictions for $f_{\rm solv}^{(\infty)}(q)$ (Eq. 9) are
plotted as dashed gray lines. Black squares and gray triangles show results
obtained with $D=0$ V/Å for the smallest and largest values of $L$,
respectively. The pink dotted lines show predictions of $f^{(\infty)}_{\rm
solv,int}$ from a dielectric continuum model, in which an interfacial layer of
width $\ell_{\rm int}=7.5$ Å is assigned a permittivity $\epsilon_{\rm
int}=2.1$ much lower than in bulk liquid, computed from (Eq. 12). The shaded
regions bound predictions with $6\,{\rm\AA}\leq\ell_{\rm int}\leq
9\,{\rm\AA}$. Insets: Snapshots from corresponding molecular dynamics
simulations.

Having determined $\delta$, we are now in a position to test the
appropriateness of the finite size correction given by Eq. 10. To this end, in
Figs. 3a and 3b, we show $f_{\rm solv}^{(L)}(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 SM. As expected, $f_{\rm
solv}^{(L)}(q)$ exhibits a dependence on system size. Adding $\Delta f_{\rm
DCT}(L)$ removes this dependence almost entirely, as seen in Figs. 3c and 3d.
Also shown are results obtained by imposing vanishing electric displacement
field $D=0$ V/Å along $z$, which is formally equivalent to the commonly used
Yeh-Berkowitz approach for approximating 2D Ewald summation Yeh and Berkowitz
(1999); Zhang and Sprik (2016). As results obtained with $D=0$ V/Å should
approximate $L\to\infty$, they do not require a finite size correction.
Importantly, excellent agreement with $f^{(L)}_{\rm solv}+\Delta f_{\rm
DCT}(L)$ is observed, giving us confidence that Eq. 10 provides a meaningful
finite size correction.

The fact that the simple DCT model outlined in Sec. III describes the finite
size behavior of $f_{\rm solv}^{(L)}$ so well suggests it is reasonable to
think of thin water films as having a uniform dielectric constant equal to
that of bulk in the region they occupy. Even more tellingly, the extrapolation
$f_{\rm solv}^{(L)}+\Delta f_{\rm DCT}(L)$ from simulation agrees well with
the continuum prediction $f_{\rm solv}^{(\infty)}$ in Eq. 9.

To provide a physical interpretation for the length scale $\delta$, Fig. 4
shows number density profiles $\rho(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 $\ell_{\rm w}/2=(w-\delta)/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 fluctuations.
The water film thickness we have inferred is thus the largest that could be
reasonably justified based on the statistics of molecular configurations.

Figure 4: Number density profiles $\rho(z)$ for hydrogen (dashed blue) and
oxygen (solid blue) atoms of water, with $q=0$ $e$/Å2 for (a) $w=40$ Å and (b)
$w=20$ Å. In both cases the WCA particles coincide with the charged planes.
The vertical dashed line shows the location $z=w/2$ of WCA particles, and the
vertical dotted line indicates the dielectric boundary at $z=(w-\delta)/2$.
(The shaded region indicates the same 95 % confidence interval as in Fig. 2.)
In both cases, the dielectric boundary aligns closely with the vanishing of
hydrogen atom density.

##  V Assessing the validity of other models

We have shown that $f_{\rm solv}^{(L)}+\Delta f_{\rm DCT}(L)$ obtained from
simulation agrees well with the predictions of a simple DCT in which the
dielectric constant of thin films is identical to that of the bulk liquid (Eq.
9). If we were to decrease $\epsilon$, agreement with simulation data would
require assigning $\ell_{\rm 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 film is lower than in the homogeneous fluid. By
itself, however, it does not rule out a model in which the interfacial regions
have a permittivity $\epsilon_{\rm int}$ that is distinct from the bulk region
they sandwich. For such a model, the free energy reads

$f^{(\infty)}_{\rm solv,int}=-2\pi q^{2}\bigg{[}\ell_{\rm
bulk}\left(\frac{\epsilon-1}{\epsilon}\right)+2\ell_{\rm
int}\left(\frac{\epsilon_{\rm int}-1}{\epsilon_{\rm int}}\right)\bigg{]}.$
(12)

where $\ell_{\rm int}$ is the width of each interfacial region.

Following the dielectric imaging experiments of Fumagalli _et al._ Fumagalli
_et al._ (2018), we take $\ell_{\rm int}=7.5\pm 1.5$ Å and $\epsilon_{\rm
int}=2.1$, and require the total width $\ell_{\rm w}=\ell_{\rm
bulk}+2\ell_{\rm int}$ to have the same value as in the uniform dielectric
model: As discussed above, it is unreasonable to allow $\ell_{\rm w}$ to
increase from that value. Decreasing $\ell_{\rm w}$, on the other hand, offers
less flexibility to a model that introduces regions of low dielectric constant
at the expense of those with high dielectric constant. The resulting
predictions of $f^{(\infty)}_{\rm solv,int}$ are shown in Figs. 3c and 3d
(labeled “bulk+interface”), where poor agreement with the simulation data is
observed. Quantitatively different (but not significantly improved)
predictions would be obtained with different choices of $\ell_{\rm int}$ and
$\epsilon_{\rm int}$. We find generally that $\ell_{\rm int}=0$ (or
equivalently, $\epsilon_{\rm int}=\epsilon_{\rm liq}$) yields the best
agreement with simulation. Evidence for this conclusion is provided in SM.

Figure 5: Solvation free energies with the charged planes moved 5 Å into
vacuum. (a) $f_{\rm solv}^{(L)}(q)$ exhibits the same finite size dependence
as in Fig. 3b, implying the same value of $w-\delta$ and thus demonstrating
that the layer of width $\delta/2$ should not be associated with the liquid
phase. (b) Adding $\Delta f_{\rm DCT}(L)$ to the results in (a), with $w=30$ Å
and $\delta=12.09$ Å, essentially removes dependence on $L$ entirely. Inset:
Snapshot from a molecular dynamics simulation showing the position of the
charged planes relative to the WCA centers (see Fig. 1). The double headed
arrow indicates $w$.

The width $\ell_{\rm int}$ of a notional interfacial layer differs
fundamentally from the length scale $\delta$ in our simple uniform DCT. They
can nonetheless easily be confused. In our case $z=\pm(w-\delta)/2$ marks the
location of a sharp interface between vapor and bulk liquid. This interface
does not coincide with the location $z=\pm w/2$ of the confining charged walls
because their constituent WCA particles exclude volume. $\delta$ 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 (Figs. 3b and 3d) by displacing the charged planes 5 Å into vacuum,
with the WCA particles fixed at their original positions (i.e., the general
case considered in Fig. 1a). $\delta$ increases by 5 Å as a result, while
$\ell_{\rm w}=w-\delta$ is unchanged, i.e., $w\to 30$ Å and $\delta\to 12.09$
Å, while $\ell_{\rm w}=17.91$ Å just as before. Changing $\delta$ in this
fashion clearly has nothing to do with water’s interfacial dielectric
properties. Fig. 5 presents results for $f_{\rm solv}^{(L)}$ and $f_{\rm
solv}^{(L)}+\Delta f_{\rm DCT}(L)$ for the displaced-charge system, which are
virtually indistinguishable from their undisplaced counterparts in Fig. 3.

By contrast, a layer of width $\ell_{\rm int}$ in “bulk+interface” models is
clearly associated with the liquid. It is imagined to comprise water molecules
whose orientational fluctuations 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 $\epsilon_{\rm int}\approx 1$ Zhang (2018, 2018); Matyushov
(2021); Mondal and Bagchi (2021); Fumagalli  _et al._ (2018). 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.

##  VI Ultra thin films of water

We have established so far that films of water with $\ell_{\rm w}\gtrsim 18$ Å
behave quantitavely like simple dielectric continua with regard to their
response to a uniform electric field. We now investigate the behavior of
‘ultra thin’ films confined between charged plates with $w\leq 10$ Å. In Figs.
6a, 6b and 6c, we show $f^{(\infty)}_{\rm solv}$ obtained with simulation for
$w=5$ Å, 7.5 Å, and 10 Å, respectively, using $\delta=2.09\pm 0.17$ Å to
correct for finite size effects (Eq. 10). As the thickness of the water slab
is reduced to length scales comparable to $\ell_{\epsilon}\approx 6$ Å,
treating the water molecules as a dielectric continuum is certainly
questionable. Discrepancies between $f_{\rm solv}^{(\infty)}$ (Eq. 9) and
$f_{\rm solv}^{(L)}+\Delta f_{\rm DCT}(L)$ indeed become apparent as $w$ is
decreased below 1 nm, but the relative 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 Eq. 9 provides a
reasonable ballpark estimate of the solvation free energy. For two to three
molecular layers (Fig. 6b and 6c), quantitative agreement between simple DCT
and the simulation data is recovered almost entirely.

Figure 6: Solvation free energy $f_{\rm solv}^{(L)}(q)+\Delta f_{\rm DCT}(L)$
in ultra thin water films, for (a) $w=5$ Å, (b) $w=7.5$ Å, and (c) $w=10$ Å.
For each value of $w$, simulations with $L=2w,3w,\ldots,6w$ have been
performed. In all cases the WCA particles coincide with the charged planes.
The dashed line shows $f_{\rm solv}^{(\infty)}(q)$ predicted by DCT (Eq. 9),
and the shaded region encompasses predictions with $\delta=2.09\pm 0.17$ Å.
Insets: snapshots from molecular dynamics simulations.

##  VII Reconciling our results with dielectric imaging experiments

Figure 7: The apparent dielectric constant $\epsilon_{\rm app}$ predicted by
simple DCT (Eq. 14) is broadly consistent with experiment. The blue solid line
is obtained with $\delta=3.5$ Å, which is estimated from _ab initio_ molecular
dynamics simulations of water on graphene. The blue shaded region indicates
the range of $\epsilon_{\rm app}$ obtained with $\delta=3.5\pm 2.0$ Å,
demonstrating the sensitivity of $\epsilon_{\rm app}$ to uncertainty in the
film thickness. The light and dark gray dotted lines indicate $\epsilon_{\rm
app}=2.1$ and $\epsilon_{\rm app}=80$, respectively.

The conclusion we have drawn from computer simulations—that the dielectric
response of nanoscopically thin water films can be anticipated from bulk
properties alone—is squarely at odds with the conclusion drawn by Fumagalli
_et al._ Fumagalli  _et al._ (2018) based on dielectric imaging measurements
of confined 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 film can cause
the apparent dielectric constant $\epsilon_{\rm app}$ to depend sensitively on
film thickness. More specifically, we assess the consequences of assigning a
width $l_{\rm w}+\delta$ to a film whose actual thickness is $l_{\rm w}$,
i.e., failing to account for the volume excluded by a confining substrate.
Based on this assignment and the development in Sec. III, we would expect a
solvation free energy

$f_{\rm solv,app}^{(\infty)}(q)=-2\pi q^{2}\frac{\epsilon_{\rm
app}-1}{\epsilon_{\rm app}}w.$ (13)

Equating this with the “true” free energy in Eq. 9, we obtain

$\epsilon_{\rm
app}(w)=\left(1-\frac{\epsilon-1}{\epsilon}\frac{w-\delta}{w}\right)^{-1},$
(14)

which depends explicitly on the film’s thickness. While similar functional
forms for $\epsilon_{\rm app}(w)$ have been reported previously Zhang (2018);
Matyushov (2021), the physical interpretation here is different: As discussed
in Sec. V, $\delta$ is not to be associated with the properties of interfacial
water.

Fig. 7 plots the apparent permittivity $\epsilon_{\rm app}$ in Eq. 14 as a
function of $w$. Here we have set $\epsilon=\epsilon_{\rm liq}=80$ and
estimated $\delta\approx 3.5$ Å for the graphene-water interface (based on
density profiles obtained from _ab initio_ molecular dynamics simulations
Tocci  _et al._ (2020)). Diminished values of $\epsilon_{\rm app}$ at small
$w$ could easily, but mistakenly, be taken to signify a strong suppression of
polarization fluctuations and response in nanoscale water films.

The inference of suppressed interfacial permittivity from experiments may
suffer from the same issues that cause $\epsilon_{\rm app}$ to depend strongly
on film thickness, much as suggested by Ref. Loche  _et al._ , 2020. To
emphasize this point, in Fig. 7 we include dielectric imaging data from Ref.
Fumagalli  _et al._ , 2018, which exhibit a very similar dependence on $w$. As
an important caveat, the samples studied by 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. Fumagalli  _et al._ , 2018 may not in
fact require invoking an interfacial dead layer.

##  VIII Summary

In this article, we have probed the dielectric response of thin water films
using molecular simulations with finite size corrections from DCT. We
specifically calculated the solvent contribution to the reversible work of
introducing charge to parallel plates that confine the film. Our results
demonstrate that response to such slowly varying external fields 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 films
extend to the point where the hydrogen number density approximately vanishes,
and thus incorporate all microscopic sources of polarization fluctuations.
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 first peak in the solute-hydrogen radial distribution
function Cox, Mandadapu, and Geissler (2021). Within this simple DCT approach,
which achieves quantitative agreement with simulation for films $\gtrsim 1$
nm, water’s interfacial regions do not enter as separate domains. Loche _et
al._ Loche  _et al._ (2020) have similarly concluded that water films 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 Fumagalli  _et al._
(2018) 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 Kathmann  _et al._ (2011).
Further investigation of this point is left for future work.

Our conclusion is further supported by results for subnanometer water films
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 Refs. Schlaich, Knapp, and Netz, 2016; Loche
_et al._ , 2020, which suggest greatly suppressed polarizability at comparable
scales of confinement.

To be certain, applying simple DCT to these subnanometer films cannot be
carefully justified (see Sec. II). 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 $\Delta
f_{\rm DCT}(L)$ for the solvation of small spherical ions in water work
remarkably well down to the nanometer scale, for both bulk Hummer, Pratt, and
García (1996); Figueirido, Del Buono, and Levy (1995); Hünenberger and
McCammon (1999) and interfacial systems Cox and Geissler (2018). Similarly,
the simple DCT model described in this article has also been found to
accurately capture mean field-like corrections for thin films of water where
electrostatic interactions are treated in a short-ranged fashionCox (2020).

The approach we have taken is not well suited to address a free liquid-vapor
interface, whose substantial topographical fluctuations greatly complicate
formulating and solving an appropriate dielectric boundary value problem.
Previous work on the free interface has emphasized that in representative
configurations the liquid phase terminates sharply at any given lateral
position Willard and Chandler (2010, 2014). 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 short-wavelength solvent response that is not well-described
by simple continuum approaches Hummer, Pratt, and García (1996). In such
cases, and in contrast to the thin films considered in this work,
$f^{(\infty)}_{\rm solv}$ predicted by DCT is a poor estimate of that obtained
from simulations Cox and Geissler (2018); Cox, Mandadapu, and Geissler (2021).
The impact of such profound perturbations are even more pronounced for solutes
near soft interfaces like that between water and its vapor, where distortions
of the interface result in nonlinearities beyond the scope of current
theoretical treatments Cox  _et al._ (2020). But for perturbations that vary
slowly in space, like the uniform fields 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 films) as a dielectric medium: Its
local permittivity is equivalent to the bulk dielectric constant, all the way
down to nanometer length scales.

##  IX Methods

All simulations followed the basic setup shown in Fig. 1a. Two planes of
$N_{\rm site}=100$ point charges were placed on a square lattice at $z=\pm
w/2$. Water molecules, modeled with the SPC/E potential Berendsen, Grigera,
and Straatsma (1987), were confined to the region $-w/2\leq z\leq w/2$ by
volume-excluding WCA particles Weeks, Chandler, and Andersen (1971), whose
centers, for the most part, coincided with the point charges at $z=\pm w/2$.
The interaction between an individual WCA particle and a water molecule is
defined by

$u_{\rm
WCA}(r)=\begin{cases}4\varepsilon\left[\left(\frac{\sigma}{r}\right)^{12}-\left(\frac{\sigma}{r}\right)^{6}\right],&r\leq
2^{1/6}\sigma,\\\ 0,&r>2^{1/6}\sigma,\end{cases}$ (15)

where $\sigma=2.5$ Å, $\varepsilon=0.1$ kcal/mol, and $r$ is the distance
between the WCA particle and the oxygen atom of the water molecule. As
described in Sec. V, 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_{\rm
site}/e=0,1\times 10^{-3},\ldots,5\times 10^{-3}$. The total volume of the
simulation cells was $12.75\,{\rm\AA}\times 12.75\,{\rm\AA}\times L$, where
$L$ takes values as indicated throughout the article.

All simulations were performed with the LAMMPS simulations package Plimpton
(1995). Lennard-Jones interactions between water molecules were truncated and
shifted at 10 Å. Long range electrostatic interactions were evaluated using
particle-particle particle-mesh Ewald summation Hockney and Eastwood (1988),
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 ÅKolafa and Perram (1992). Where indicated in the text, the
electric displacement field along $z$ was set to zero, using the
implementation given in Ref. Cox and Sprik, 2019. The geometry of the water
molecules was constrained using the RATTLE algorithm Andersen (1983).
Temperature was maintained at 298 K with a Nosé-Hoover chain thermostat
Shinoda, Shiga, and Mikami (2004); Tuckerman  _et al._ (2006) 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.

Table 1: Number of water molecules $N_{\rm wat}$ for each value of $w$ investigated (WCA centers coincide with point charges). $w$/Å | 5 | 7.5 | 10 | 20 | 25 | 30 | 40
---|---|---|---|---|---|---|---
$N_{\rm wat}$ | 14 | 27 | 41 | 93 | 125 | 143 | 206

The free energy of charging up parallel plate capacitors was computed by
averaging electric potentials appropriately. Let $\varphi^{(i)}_{\rm
solv,hi}(\mathbf{R}^{N})$ and $\varphi^{(j)}_{\rm solv,lo}(\mathbf{R}^{N})$
denote the instantaneous electric potentials due to the solvent with
configuration $\mathbf{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_{\rm solv}^{(L)}(Q)$ is then defined by

$\displaystyle\exp\left[-\beta F_{\rm solv}(Q)\right]=\langle\exp\left[-\beta
Q\Delta\varphi_{\rm solv}\right]\rangle_{0}^{(L)}$
$\displaystyle=\int\\!\mathrm{d}(\Delta\varphi_{\rm
solv})\,P^{(L)}(\Delta\varphi_{\rm solv};0)\exp\left[-\beta
Q\Delta\varphi_{\rm solv}\right],$ (16)

where

$\displaystyle\Delta\varphi_{\rm solv}(\mathbf{R}^{N})$
$\displaystyle=\frac{1}{N_{\rm site}}\left(\sum_{i\in\rm hi}\varphi^{(i)}_{\rm
solv,hi}(\mathbf{R}^{N})-\sum_{j\in\rm lo}\varphi^{(j)}_{\rm
solv,lo}(\mathbf{R}^{N})\right)$ $\displaystyle=\varphi_{\rm
solv,hi}(\mathbf{R}^{N})-\varphi_{\rm solv,lo}(\mathbf{R}^{N});$ (17)

$P^{(L)}(\Delta\varphi_{\rm solv};Q)$ is the probability distribution of
$\Delta\varphi_{\rm solv}$ in the presence of two charged planes with total
charges $\pm Q$, in a simulation box of length $L$; and
$\langle\cdot\rangle_{0}^{(L)}$ denotes an average over
$P^{(L)}(\Delta\varphi_{\rm solv};0)$.

Similar to our previous studiesCox and Geissler (2018); Cox, Mandadapu, and
Geissler (2021), $F_{\rm solv}(Q)$ was computed using the MBAR algorithmShirts
and Chodera (2008). The solvation free energies per unit area, $f_{\rm
solv}^{(L)}$ that we consider are then obtained by dividing $F_{\rm solv}(Q)$
by the cross-sectional area of the simulation cell.

## Supplementary Material

Supplementary Material includes a detailed derivation of results obtained for
the periodic continuum model presented in Sec. III. Results for the
“bulk+interface” model with different parameters are also presented, along
with those obtained with $w=30$ Å and $w=25$ Å.

## 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 conflicts to declare.

###### Acknowledgements.

We are grateful to Laura Fumagalli for sharing her experimental results. We
also thank Gabriele Tocci for sharing his results from _ab initio_ molecular
dynamics. S.J.C is a Royal Society University Research Fellow (URF\R1\211144)
at the University of Cambridge. P.L.G is supported by the U.S. Department of
Energy, Office of Basic Energy Sciences, through the Chemical Sciences
Division (CSD) of Lawrence Berkeley National Laboratory (LBNL), under Contract
DE-AC02-05CH11231. For the purposes of open access, the authors have applied a
Creative Commons Attribution (CC BY) licence to any Author Accepted Manuscript
version arising.

## Data Availability Statement

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

## References

  * Geim and Grigorieva (2013) A. K. Geim and I. V. Grigorieva, Nature 499, 419 (2013).
  * Celebi  _et al._ (2014) K. Celebi, J. Buchheim, R. M. Wyss, A. Droudian, P. Gasser, I. Shorubalko, J.-I. Kye, C. Lee, and H. G. Park, Science 344, 289 (2014).
  * Jain  _et al._ (2015) T. Jain, B. C. Rasera, R. J. S. Guerrero, M. S. Boutilier, S. C. O’hern, J.-C. Idrobo, and R. Karnik, Nat, Nanotechnol. 10, 1053 (2015).
  * Bocquet (2020) L. Bocquet, Nat. Mater. 19, 254 (2020).
  * Algara-Siller  _et al._ (2015) G. Algara-Siller, O. Lehtinen, F. Wang, R. R. Nair, U. Kaiser, H. Wu, A. K. Geim, and I. V. Grigorieva, Nature 519, 443 (2015).
  * Fumagalli  _et al._ (2018) L. Fumagalli, A. Esfandiar, R. Fabregas, S. Hu, P. Ares, A. Janardanan, Q. Yang, B. Radha, T. Taniguchi, K. Watanabe, G. Gomila, K. S. Novoselov, and A. K. Geim, Science 360, 1339 (2018).
  * Nakamuro  _et al._ (2021) T. Nakamuro, M. Sakakibara, H. Nada, K. Harano, and E. Nakamura, J. Am. Chem. Soc. 143, 1763 (2021).
  * Wang  _et al._ (2021) L. Wang, J. Chen, S. J. Cox, L. Liu, G. C. Sosso, N. Li, P. Gao, A. Michaelides, E. Wang, and X. Bai, Phys. Rev. Lett. 126, 136001 (2021).
  * Ballenegger and Hansen (2005) V. Ballenegger and J.-P. Hansen, J. Chem. Phys 122, 114711 (2005).
  * Bonthuis, Gekle, and Netz (2011) D. J. Bonthuis, S. Gekle, and R. R. Netz, Phys. Rev. Lett. 107, 166102 (2011).
  * Bonthuis, Gekle, and Netz (2012) D. J. Bonthuis, S. Gekle, and R. R. Netz, Langmuir 28, 7679 (2012).
  * Schlaich, Knapp, and Netz (2016) A. Schlaich, E. W. Knapp, and R. R. Netz, Phys. Rev. Lett. 117, 048001 (2016).
  * Loche  _et al._ (2018) P. Loche, C. Ayaz, A. Schlaich, D. J. Bonthuis, and R. R. Netz, J. Phys. Chem. Lett. 9, 6463 (2018).
  * Zhang (2018) C. Zhang, J. Chem. Phys. 148, 156101 (2018).
  * Olivieri, Hynes, and Laage (2021) J.-F. Olivieri, J. T. Hynes, and D. Laage, J. Phys. Chem. Lett. 12, 4319 (2021).
  * Matyushov (2021) D. V. Matyushov, J. Phys. Chem. B 125, 8282 (2021).
  * Mondal and Bagchi (2021) S. Mondal and B. Bagchi, J. Chem. Phys. 154, 044501 (2021).
  * Rowlinson and Widom (2002) J. J. S. Rowlinson and B. Widom,  _Molecular Theory of Capillarity_ (Courier Dover, Minola, NY, USA, 2002).
  * Loche  _et al._ (2020) P. Loche, C. Ayaz, A. Wolde-Kidan, A. Schlaich, and R. R. Netz, J. Phys. Chem. B 124, 4365 (2020).
  * Jalali  _et al._ (2020) H. Jalali, H. Ghorbanfekr, I. Hamid, M. Neek-Amal, R. Rashidi, and F. Peeters, Phys. Rev. E 102, 022803 (2020).
  * Hamid  _et al._ (2021) I. Hamid, H. Jalali, F. M. Peeters, and M. Neek-Amal, J. Chem. Phys. 154, 114503 (2021).
  * Fulton (1978a) R. L. Fulton, J. Chem. Phys. 68, 3089 (1978a).
  * Fulton (1978b) R. L. Fulton, J. Chem. Phys. 68, 3095 (1978b).
  * Fulton (1983) R. L. Fulton, J. Chem. Phys. 78, 6865 (1983).
  * Caillol (1992) J. Caillol, J. Chem. Phys. 96, 7039 (1992).
  * Zhang, Hutter, and Sprik (2016) C. Zhang, J. Hutter, and M. Sprik, J. Phys. Chem. Lett. 7, 2696 (2016).
  * Cox (2020) S. J. Cox, Proc. Natl. Acad. Sci. 117, 19746 (2020).
  * Jackson (1999) J. D. Jackson,  _Classical electrodynamics_ , 3rd ed. (Wiley, New York, 1999).
  * Cox and Geissler (2018) S. J. Cox and P. L. Geissler, J. Chem. Phys. 148, 222823 (2018).
  * Hünenberger and McCammon (1999) P. H. Hünenberger and J. A. McCammon, J. Chem. Phys. 110, 1856 (1999).
  * Weeks, Chandler, and Andersen (1971) J. D. Weeks, D. Chandler, and H. C. Andersen, J. Chem. Phys. 54, 5237 (1971).
  * Yeh and Berkowitz (1999) I.-C. Yeh and M. L. Berkowitz, J. Chem. Phys. 111, 3155 (1999).
  * Zhang and Sprik (2016) C. Zhang and M. Sprik, Physical Review B 94, 245309 (2016).
  * Tocci  _et al._ (2020) G. Tocci, M. Bilichenko, L. Joly, and M. Iannuzzi, Nanoscale 12, 10994 (2020).
  * Cox, Mandadapu, and Geissler (2021) S. J. Cox, K. K. Mandadapu, and P. L. Geissler, J. Chem. Phys. 154, 244502 (2021).
  * Kathmann  _et al._ (2011) S. M. Kathmann, I.-F. W. Kuo, C. J. Mundy, and G. K. Schenter, J. Phys. Chem. B 115, 4369 (2011).
  * Hummer, Pratt, and García (1996) G. Hummer, L. R. Pratt, and A. E. García, J. Phys. Chem. 100, 1206 (1996).
  * Figueirido, Del Buono, and Levy (1995) F. Figueirido, G. S. Del Buono, and R. M. Levy, J. Chem. Phys. 103, 6133 (1995).
  * Willard and Chandler (2010) A. P. Willard and D. Chandler, J. Phys. Chem. B 114, 1954 (2010).
  * Willard and Chandler (2014) A. P. Willard and D. Chandler, J. Chem. Phys. 141, 18C519 (2014).
  * Cox  _et al._ (2020) S. J. Cox, D. G. Thorpe, P. R. Shaffer, and P. L. Geissler, Chem. Sci. 11, 11791 (2020).
  * Berendsen, Grigera, and Straatsma (1987) H. J. C. Berendsen, J. R. Grigera, and T. P. Straatsma, J. Phys. Chem. 91, 6269 (1987).
  * Plimpton (1995) S. Plimpton, J. Comput. Phys. 117, 1 (1995).
  * Hockney and Eastwood (1988) R. W. Hockney and J. W. Eastwood,  _Computer simulation using particles_ (CRC Press, 1988).
  * Kolafa and Perram (1992) J. Kolafa and J. W. Perram, Mol. Sim. 9, 351 (1992).
  * Cox and Sprik (2019) S. J. Cox and M. Sprik, J. Chem. Phys. 151, 064506 (2019).
  * Andersen (1983) H. C. Andersen, J. Comput. Phys. 52, 24 (1983).
  * Shinoda, Shiga, and Mikami (2004) W. Shinoda, M. Shiga, and M. Mikami, Phys. Rev. B 69, 134103 (2004).
  * Tuckerman  _et al._ (2006) M. E. Tuckerman, J. Alejandre, R. López-Rendón, A. L. Jochim, and G. J. Martyna, J. Phys. A 39, 5629 (2006).
  * Shirts and Chodera (2008) M. R. Shirts and J. D. Chodera, J. Chem. Phys. 129, 124105 (2008).
  * Wirnsberger  _et al._ (2016) P. Wirnsberger, D. Fijan, A. Šarić, M. Neumann, C. Dellago, and D. Frenkel, J. Chem. Phys. 144, 224102 (2016).
  * Pan, Yi, and Hu (2017) C. Pan, S. Yi, and Z. Hu, Phys. Chem. Chem. Phys. 19, 4861 (2017).

Supplementary Material

##  S1 Derivation

The system we consider is shown schematically in Fig. 1. Two planes with
charge density $\pm q$ are located at $z=\pm w/2$. We will consider a more
general case than in the main article. Here, a linear dielectric occupies the
region $-(w/2-\delta_{\rm lo})\leq z\leq(w/2-\delta_{\rm hi})$, such that
$\delta_{\rm lo}+\delta_{\rm hi}=\delta$; while the electrostatic potential is
sensitive to the values of $\delta_{\rm lo}$ and $\delta_{\rm hi}$, we will
show that $\Delta f_{\rm DCT}(L)$ only depends on their sum. The boundaries of
the dielectric are situated at $\xi_{\rm hi}=w/2-\delta_{\rm hi}$ and
$\xi_{\rm lo}=w/2-\delta_{\rm lo}$. The polarization of the medium is $P$.

### Potential due to the charged plates

The potential due to the charged plates is,

$\phi_{q}(z)=4\pi\int_{\rm
cell}\\!\mathrm{d}z^{\prime}\rho_{q}(z^{\prime})J(z-z^{\prime}),$ (S1)

with

$\rho_{q}(z)=q\big{[}\delta_{\rm D}(z-w/2)-\delta_{\rm D}(z+w/2)\big{]},$ (S2)

where $\delta_{\rm D}(x)$ is the Dirac delta-function, andWirnsberger  _et
al._ (2016); Pan, Yi, and Hu (2017); Cox and Geissler (2018)

$J(z)=\text{const.}+\frac{z^{2}}{2L}-\frac{|z|}{2}.$ (S3)

Inside the region occupied by the charged sheets, $-w/2\leq z\leq w/2$, we
have

$\phi_{q}(z)=4\pi q\bigg{(}-\frac{zw}{L}+z\bigg{)}.$ (S4)

Similarly, for $w/2<z\leq L/2$,

$\phi_{q}(z)=4\pi q\bigg{(}-\frac{zw}{L}+\frac{w}{2}\bigg{)},$ (S5)

while for $-L/2\leq z<-w/2$,

$\phi_{q}(z)=4\pi q\bigg{(}-\frac{zw}{L}-\frac{w}{2}\bigg{)}.$ (S6)

### Potential due to a uniformly polarized dielectric

A uniformly polarized dielectric generates the same electric potential as a
charge distribution comprising two uniformly charged planes,

$\rho_{\rm solv}(z)=P\big{[}\delta_{\rm D}(z-\xi_{\rm hi})-\delta_{\rm
D}(z+\xi_{\rm lo})\big{]}.$ (S7)

This leads to the following potential,

$\phi_{\rm solv}(z)=4\pi P\bigg{[}-\frac{z(\xi_{\rm hi}+\xi_{\rm
lo})}{L}+\frac{\xi^{2}_{\rm hi}-\xi^{2}_{\rm
lo}}{2L}+\frac{1}{2}\big{(}|z+\xi_{\rm lo}|-|z-\xi_{\rm hi}|\big{)}\bigg{]}$
(S8)

For the region occupied by the dielectric we have ($-\xi_{\rm lo}\leq
z\leq\xi_{\rm hi}$),

$\phi_{\rm solv}(z)=4\pi P\bigg{[}-\frac{z(\xi_{\rm hi}+\xi_{\rm
lo})}{L}+z+\frac{\xi^{2}_{\rm hi}-\xi^{2}_{\rm lo}}{2L}+\frac{\xi_{\rm
lo}-\xi_{\rm hi}}{2}\bigg{]},$ (S9)

while for $\xi_{\rm hi}<z\leq L/2$

$\phi_{\rm solv}(z)=4\pi P\bigg{[}-\frac{z(\xi_{\rm hi}+\xi_{\rm
lo})}{L}+\frac{\xi^{2}_{\rm hi}-\xi^{2}_{\rm lo}}{2L}+\frac{\xi_{\rm
hi}+\xi_{\rm lo}}{2}\bigg{]},$ (S10)

and for $-L/2\leq z\leq-\xi_{\rm lo}$ we have,

$\phi_{\rm solv}(z)=4\pi P\bigg{[}-\frac{z(\xi_{\rm hi}+\xi_{\rm
lo})}{L}+\frac{\xi^{2}_{\rm hi}-\xi^{2}_{\rm lo}}{2L}-\frac{\xi_{\rm
hi}+\xi_{\rm lo}}{2}\bigg{]}.$ (S11)

### The total potential

The total potential is simply the linear superposition of potentials due to
the charged planes and the solvent, $\phi(z)=\phi_{q}(z)+\phi_{\rm solv}(z)$.
Most important for the derivation is the region $-\xi_{\rm lo}\leq
z\leq\xi_{\rm hi}$,

$\phi(z)=4\pi q\bigg{(}-\frac{zw}{L}+z\bigg{)}+4\pi P\bigg{[}-\frac{z(\xi_{\rm
hi}+\xi_{\rm lo})}{L}+z+\frac{\xi^{2}_{\rm hi}-\xi^{2}_{\rm
lo}}{2L}+\frac{\xi_{\rm lo}-\xi_{\rm hi}}{2}\bigg{]}.$ (S12)

The potential in each of the remaining regions is listed below.

For $-L/2\leq z<-w/2$:

$\phi(z)=4\pi q\bigg{(}-\frac{zw}{L}-\frac{w}{2}\bigg{)}+4\pi
P\bigg{[}-\frac{z(\xi_{\rm hi}+\xi_{\rm lo})}{L}+\frac{\xi^{2}_{\rm
hi}-\xi^{2}_{\rm lo}}{2L}-\frac{\xi_{\rm hi}+\xi_{\rm lo}}{2}\bigg{]}.$ (S13)

For $-w/2\leq z<-\xi_{\rm lo}$:

$\phi(z)=4\pi q\bigg{(}-\frac{zw}{L}+z\bigg{)}+4\pi P\bigg{[}-\frac{z(\xi_{\rm
hi}+\xi_{\rm lo})}{L}+\frac{\xi^{2}_{\rm hi}-\xi^{2}_{\rm
lo}}{2L}-\frac{\xi_{\rm hi}+\xi_{\rm lo}}{2}\bigg{]}.$ (S14)

For $\xi_{\rm hi}<z\leq w/2$:

$\phi(z)=4\pi q\bigg{(}-\frac{zw}{L}+z\bigg{)}+4\pi P\bigg{[}-\frac{z(\xi_{\rm
hi}+\xi_{\rm lo})}{L}+\frac{\xi^{2}_{\rm hi}-\xi^{2}_{\rm
lo}}{2L}+\frac{\xi_{\rm hi}+\xi_{\rm lo}}{2}\bigg{]}.$ (S15)

For $w/2<z\leq L/2$:

$\phi(z)=4\pi q\bigg{(}-\frac{zw}{L}+\frac{w}{2}\bigg{)}+4\pi
P\bigg{[}-\frac{z(\xi_{\rm hi}+\xi_{\rm lo})}{L}+\frac{\xi^{2}_{\rm
hi}-\xi^{2}_{\rm lo}}{2L}+\frac{\xi_{\rm hi}+\xi_{\rm lo}}{2}\bigg{]}.$ (S16)

Note that $\xi_{\rm hi}+\xi_{\rm lo}=w-\delta$, where $\delta=\delta_{\rm
hi}+\delta_{\rm lo}$.

### Linear response

Equations S12–S16 provide general expressions for the total electrostatic
potential for the periodic continuum model considered in Fig. 1. As $P$
depends upon the electric field, a self-consistent solution is required. In
the case that the dielectric medium is linearly responding, however, the
solution is analytically tractable. Consider the electric field inside the
dielectric. From Eq. S12 we find for $-\xi_{\rm lo}\leq z\leq\xi_{\rm hi}$,

$E=-4\pi q\bigg{(}1-\frac{w}{L}\bigg{)}-4\pi
P\bigg{(}1-\frac{w-\delta}{L}\bigg{)}.$ (S17)

Applying the local constitutive relation, $4\pi P=(\epsilon-1)E$, we find

$P=-(\epsilon-1)\bigg{[}q\bigg{(}1-\frac{w}{L}\bigg{)}+P\bigg{(}1-\frac{w-\delta}{L}\bigg{)}\bigg{]},$
(S18)

or rearranging,

$P=-\frac{(\epsilon-1)(1-\frac{w}{L})q}{1+(\epsilon-1)(1-\frac{w-\delta}{L})}.$
(S19)

From Eq. S13, it is clear that the potential at the charged plate at $z=-w/2$,
due to the polarized dielectric is

$\phi_{\rm solv,lo}=2\pi P\bigg{[}\frac{w(w-\delta)}{L}+\frac{\xi^{2}_{\rm
hi}-\xi^{2}_{\rm lo}}{L}-(w-\delta)\bigg{]}.$ (S20)

Similarly, for the charged plate at $z=+w/2$ we have,

$\phi_{\rm solv,hi}=2\pi P\bigg{[}-\frac{w(w-\delta)}{L}+\frac{\xi^{2}_{\rm
hi}-\xi^{2}_{\rm lo}}{L}+(w-\delta)\bigg{]}.$ (S21)

The solvation free energy is $f^{(L)}_{\rm solv}=q(\phi_{\rm
solv,hi}-\phi_{\rm solv,lo})/2$. Combining Eqs. S20, S21 and S19 gives,

$f^{(L)}_{\rm solv}=-2\pi
q^{2}(w-\delta)\frac{(\epsilon-1)(1-\frac{w}{L})^{2}}{1+(\epsilon-1)(1-\frac{w-\delta}{L})}.$
(S22)

In the limit $L\to\infty$ this gives,

$f^{(\infty)}_{\rm solv}=-2\pi q^{2}\frac{(w-\delta)(\epsilon-1)}{\epsilon}.$
(S23)

The finite size correction we must apply is $\Delta f_{\rm
DCT}(L)=f^{(\infty)}_{\rm solv}-f^{(L)}_{\rm solv}$. Thus,

$\Delta f_{\rm DCT}(L)=2\pi
q^{2}(w-\delta)(\epsilon-1)\Bigg{[}\frac{(1-\frac{w}{L})^{2}}{1+(\epsilon-1)(1-\frac{w-\delta}{L})}-\frac{1}{\epsilon}\Bigg{]}.$
(S24)

##  S2 Sensitivity of $f_{\rm solv,int}^{(\infty)}$ to $\epsilon_{\rm int}$
and $\ell_{\rm int}$

In Fig. S1 we plot $f_{\rm solv}^{(L)}(q)+\Delta f_{\rm DCT}(L)$ for $w=40$ Å
and $w=20$ Å (see Fig. 3), but with $f_{\rm solv,int}^{(\infty)}$ (Eq. 12)
parameterized with $\epsilon_{\rm int}=10$ and $\ell_{\rm int}=6.0\pm 1.5$ Å.
We argue that $\ell_{\rm int}=\ell_{\epsilon}\approx 6$ Å sets a lower bound
on reasonable values of $\ell_{\rm int}$. As discussed in the main article,
increasing $\epsilon_{\rm int}$ and decreasing $\ell_{\rm int}$, while
imposing the constraint $\ell_{\rm w}=\ell_{\rm bulk}+2\ell_{\rm int}$ will
obviously reduce discrepancies between $f_{\rm solv,int}^{(\infty)}$ and
$f_{\rm solv}^{(L)}(q)+\Delta f_{\rm DCT}(L)$, as evidenced by Fig. S1.
Nonetheless, it is clear that $f^{(\infty)}_{\rm solv}$ given by Eq. 9 still
provides a superior description of the simulation data.

Figure S1: $f_{\rm solv}^{(L)}(q)+\Delta f_{\rm DCT}(L)$ with (a) $w=40$ Å and
(b) $w=20$ Å. These data are the same as shown Figs. 3c and 3d, except
$f^{(\infty)}_{\rm solv,int}$ (pink dotted lines) is plotted with
$\epsilon_{\rm int}=10$ and $\ell_{\rm int}=6.0\pm 1.5$ Å. While discrepancies
between $f^{(\infty)}_{\rm solv,int}$ and $f_{\rm solv}^{(L)}(q)+\Delta f_{\rm
DCT}(L)$ are reduced compared to Figs. 3c and 3d, $f^{(\infty)}_{\rm solv}$
given by Eq. 9 (gray dashed lines) still gives a superior description of the
simulation data.

##  S3 Results with $w=30$ Å and $w=25$ Å (WCA centers coincide with the
charged planes)

In Fig. S2 we present results for $f_{\rm solv}^{(L)}(q)$ and $f_{\rm
solv}^{(L)}(q)+\Delta f_{\rm DCT}(L)$ obtained with $w=30$ Å and $w=25$ Å,
where in both cases, the positions of the WCA particles coincide with the
charged planes. We draw the same conclusions as from Fig. 3 in the main
article.

Figure S2: Dependence of solvation free energy $f_{\rm solv}^{(L)}(q)$ on
system size $L$, shown in (a) and (b) for $w=30$ Å and $w=25$ Å, respectively.
The values of $L$ for $w=30$ Å are are indicated in the legend of panel (a);
those for the thinner liquid slab are shown in (b). In both cases the WCA
particles coincide with the charged planes. Adding $\Delta f_{\rm DCT}(L)$
given by Eq. 10 largely removes this sensitivity, as seen in (c) and (d) for
$w=30$ Å and $w=25$ Å, respectively. DCT predictions for $f_{\rm
solv}^{(\infty)}(q)$ (Eq. 9) are plotted as dashed gray lines. Black squares
and gray triangles show results obtained with $D=0$ V/Å for the smallest and
largest values of $L$, respectively. The pink dotted lines show predictions of
$f^{(\infty)}_{\rm solv,int}$ from a dielectric continuum model, in which an
interfacial layer of width $\ell_{\rm int}=7.5$ Å is assigned a permittivity
$\epsilon_{\rm int}=2.1$ much lower than in bulk liquid, computed from (Eq.
12). The shaded regions bound predictions with $6\,{\rm\AA}\leq\ell_{\rm
int}\leq 9\,{\rm\AA}$. Figure S3: Number density profiles $\rho(z)$ for
hydrogen (dashed blue) and oxygen (solid blue) atoms of water, with $q=0$
$e$/Å2 for (a) $w=30$ Å and (b) $w=25$ Å. In both cases the WCA particles
coincide with the charged planes. The vertical dashed line shows the location
$z=w/2$ of WCA particles, and the vertical dotted line indicates the
dielectric boundary at $z=(w-\delta)/2$. (The shaded region indicates the same
95 % confidence interval as in Fig. 2.) In both cases, the dielectric boundary
aligns closely with the vanishing of hydrogen atom density.