Inferring properties of macroscopic solutions from molecular simulations is complicated by the limited size of systems that can be feasibly examined with a computer. When long-ranged electrostatic interactions are involved, the resulting finite size effects can be substantial and may attenuate very slowly with increasing system size, as shown by previous work on dilute ions in bulk aqueous solution. Here we examine corrections for such effects, with an emphasis on solvation near interfaces. Our central assumption follows the perspective of Hünenberger and McCammon [J. Chem. Phys. **110**, 1856 (1999)]: Long-wavelength solvent response underlying finite size effects should be well described by reduced models like dielectric continuum theory, whose size dependence can be calculated straightforwardly. Applied to an ion in a periodic slab of liquid coexisting with vapor, this approach yields a finite size correction for solvation free energies that differs in important ways from results previously derived for bulk solution. For a model polar solvent, we show that this new correction quantitatively accounts for the variation of solvation free energy with volume and aspect ratio of the simulation cell. Correcting periodic slab results for an aqueous system requires an additional accounting for the solvent’s intrinsic charge asymmetry, which shifts electric potentials in a size-dependent manner. The accuracy of these finite size corrections establishes a simple method for *a posteriori* extrapolation to the thermodynamic limit and also underscores the realism of dielectric continuum theory down to the nanometer scale.

## I. INTRODUCTION

Molecular simulations of ions near the boundaries of liquid solutions^{1–31} have challenged and reshaped traditional microscopic perspectives on solvation thermodynamics.^{32} In the case of aqueous interfaces, they have inspired a new generation of experiments^{10,15,28,29,33–43} and theories^{4,8–10,19,24,30,44} to identify basic driving forces and their implications across a broad range of fields such as electrochemistry,^{45} aerosols,^{46,47} and protein surfaces.^{48,49} The great strength of such simulations—resolving atomistically detailed structure on ultrafast time scales—is counterbalanced, however, by the necessary use of imperfect potential energy models and by the necessary restriction to small system sizes. This paper addresses the latter problem, aiming to extract from microscopically finite simulations the behavior of macroscopically large solutions and their interfaces.

Efforts to correct for finite size effects in molecular simulation are nearly as old as the simulation methods themselves. Most notably, the use of periodic boundary conditions (PBC) mitigates the unrealistically high surface-to-volume ratios, and degrees of interfacial curvature of small systems. Even for systems that interact with relatively short-ranged interactions, however, careful treatment of interactions between distant particles is often required to obtain well-converged thermodynamic properties.^{50–64} In the case of ion solvation, the very slow decay of Coulomb interactions with distance demands an even more careful treatment. Indeed, a simulation cell bearing a net charge cannot be periodically replicated without incurring a super-extensive energetic cost. Whether implicitly or explicitly, neutralizing countercharge must be introduced to ensure a well-defined thermodynamic limit when an ion is subject to periodic boundary conditions. Here we will exclusively consider systems that are periodic in three dimensions, whose unit cells contain a single charged solute and a spatially uniform neutralizing background charge density, in addition to a collection of electroneutral solvent molecules. Due to the long range of electrostatic forces, periodic replicas of these charged constituents can generate substantial fields and potentials, which in turn induce solvent response that decays only gradually with increasing system size. These are the contributions we wish to assess and correct.

While far from macroscopic, computationally accessible systems are large enough to capture much of the complicated molecular-scale response that is missing from the classic linear response theory for solvent polarization, i.e., dielectric continuum theory (DCT). Our approach asserts that DCT accurately describes solvent response at all larger scales. As such, finite size effects in DCT should serve as useful estimates for those of more complicated molecular models. This perspective was employed by Hünenberger and McCammon^{65} to investigate the nature and magnitude of periodicity-induced artifacts for bulk solutions, building on previous studies by Hummer *et al.*^{66,67} The development described in Sec. III casts the resulting finite size corrections in a different form that facilitates analysis of interfacial systems. In Sec. III B, we confirm that our approach is equivalent to that of Ref. 65 for a charged solute in bulk solution. We also verify Hummer’s demonstration that simulations, performed using methods described in Sec. II, closely follow the predicted approach to the macroscopic limit.

Extended interfaces between coexisting phases can also be simulated using PBC. The widely used “slab geometry” includes regions of each phase (e.g., liquid and vapor) that span the periodic simulation cell in two Cartesian directions (*x* and *y*), separated by a planar interface perpendicular to the third direction (*z*). We show in Sec. II that solvation free energies computed from simulations with this geometry exhibit a strong size dependence that differs significantly from the bulk case. Results are presented both for the SPC/E model of water and also for a simple polar solvent comprising diatomic molecules with a pair of opposite charges. In Sec. III C, we derive a thermodynamic finite size correction for ions solvated by an ideal dielectric solvent in the slab geometry. This correction quantitatively reconciles the size dependence observed for the simple polar solvent.

For the case of liquid water coexisting with vapor in the slab geometry, the finite size correction obtained from DCT accounts only partly for the size dependence observed in simulations. We attribute the remaining discrepancy to the intrinsic asymmetry of water’s molecular charge distribution. In the SPC/E model, positive and negative charges are not equivalently distributed within each water molecule, in contrast to the diatomic molecules in our simple polar solvent. This asymmetry creates thermodynamic biases in the solvation of cations versus anions and also underlies the net orientation of molecules at the air-water interface. In Sec. III D, we show that charge asymmetry in the periodic slab geometry generates a size-dependent electric potential that can be easily extrapolated to the macroscopic limit. The corresponding finite size correction, added to the dielectric response correction, yields ion solvation free energies from aqueous slab simulations that are very nearly independent of unit cell size and aspect ratio. A similar charge asymmetry correction is applied to results of bulk periodic simulations in Sec. III E, yielding good agreement with the limiting behavior of slab simulations.

The size dependence of ion solvation free energies near simulated air-water interfaces thus closely follows predictions of an idealized continuum model, once effects of charge asymmetry under PBC have been reckoned. This success of DCT extends to simulation unit cells with dimensions as small as 1 nm, indicating a surprising robustness of linear response theory at nearly molecular length scales. In Sec. IV, we discuss implications of this result for microscopic mechanisms and theories of ion solvation near interfaces. In Sec. V, we conclude.

## II. SIMULATION METHODS

The general system of interest, shown schematically in Fig. 1, comprises a set of periodically replicated liquid slabs with width $w$ and dielectric constant *ϵ*, separated by regions of vacuum in the direction normal to the interface, which we take to be along the *z*-direction. The length of the primary unit cell in the *z*-direction is *L*_{z}. The liquid/vapor interface lies in the *xy*-plane and has a surface area *A* = *L*_{x}*L*_{y}, where *L*_{x} and *L*_{y} are the lengths of the simulation cell in *x* and *y*, respectively. The total volume of the simulation cell is $v$ = *L*_{x}*L*_{y}*L*_{z}. A periodically replicated solute is situated in the liquid slab at a distance *d* from one of the liquid/vapor interfaces. A bulk solution lacking liquid/vapor interfaces is obtained by setting $w$ = *L*_{z}.

We have numerically examined finite size effects on ion solvation through molecular simulations of two model solvents. The simpler of these solvents, diatomic molecules comprising oppositely charged particles separated by a short distance, closely resembles the Stockmayer fluid. The other, a simple point charge model of water, adds the complication of asymmetry under charge inversion, whose implications for finite size effects are analyzed in Sec. III D.

Our model solute is the same in both sets of simulations, namely, a Lennard-Jones particle with a point charge at its center. Parameters for the Lennard-Jones interaction, *u*(*r*) = 4*ε*[(*σ*/*r*)^{12} − (*σ*/*r*)^{6}] as a function of distance *r*, are identical to those of the SPC/E model for water, i.e., *ε* = 0.1553 kcal/mol and *σ* = 3.166 Å. This choice of *σ* corresponds roughly to the size of a Ca^{2+} ion.^{66} (The leading order correction considered in this work is, however, independent of solute size.) For consistency with standard notation, we use the symbol *σ* here to denote the Lennard-Jones diameter. The same symbol will later denote a different quantity, namely, the surface charge that results from dielectric response. The distinction should be clear from context.

In water studies, we employed the SPC/E model for interactions involving solvent molecules.^{68} This choice, made for simplicity, affects the exact values of computed solvation free energies, but our conclusions about finite size effects should apply equally well to more elaborate molecular models. Bulk liquid water simulations included a single solute and a number *N* of water molecules ranging from 64 to 756 such that the total number density was 0.033 33 Å^{−3}. Simulations of air-water coexistence included one solute and *N* = 128 water molecules in a simulation cell with lateral dimensions *L*_{x} = *L*_{y} = 12.429 Å and varying separation between periodic liquid slabs. (As the sole exception, in Fig. 6, we show results for *L*_{x} = *L*_{y} = 31.1 Å as well.) The resulting aspect ratios *L*_{z}/*L*_{x} ranged from 4 to 30. To suppress interfacial instabilities at very high aspect ratios, we introduced an external confining potential, specifically walls that effectively constrain molecules to a region slightly larger than the neat liquid slab occupies. These walls interact with the oxygen atom of each water molecule through a potential *u*_{solv-wall}(*z*), which inside the walls (|*z*| < |*z*_{wall}|) has a Weeks-Chandler-Andersen form:^{69} $usolv-wall(z)=4\epsilon solv-wall[(\delta z/\sigma solv-wall)\u221212\u2212(\delta z/\sigma solv-wall)\u22126+1/4]$ for |*δz*| < 2^{1/6}*σ*_{solv-wall}, where *δz* = *z* − *z*_{wall}, and *z*_{wall} = ±15 Å. Outside the walls (|*z*| > |*z*_{wall}|), *u*_{solv-wall}(*z*) is infinite. Otherwise (|*δz*| > 2^{1/6}*σ*_{solv-wall}), *u*_{solv-wall}(*z*) vanishes. The coordinate −*L*_{z}/2 < *z* < *L*_{z}/2 is the molecule’s position within the primary simulation cell. Away from the interfaces, the density of the liquid varied approximately between 0.032 00 Å^{−3} and 0.033 33 Å^{−3}, depending on the aspect ratio and the charge of the solute.

The simpler polar fluid we consider consists of “dumbbell” molecules, each a pair of Lennard-Jones spheres (*ε* = 0.1553 kcal/mol and *σ* = 3.166 Å) separated by a rigid bond of length *l*_{D} = 0.25 Å. Opposite charges of magnitude $$*q*_{D}$$ = 0.9*e* reside at the centers of these spheres. We consider a liquid state density of 0.025 Å^{−3}. Bulk simulations included a number *N* of dumbbells ranging from 48 to 384. Slab simulations included *N* = 96 molecules. As in aqueous simulations, lateral dimensions were fixed at *L*_{x} = *L*_{y} = 12.429 Å, while the aspect ratio *L*_{z}/*L*_{x} varied from 4 to 30. Walls were again used to confine solvent molecules to the range −15.0 Å < *z* < 15.0 Å. However, in this case, the confining potential interacts with both Lennard-Jones spheres of each dumbbell molecule.

Some of the finite size corrections detailed below involve the static dielectric constant *ϵ*. In the case of SPC/E water, the value of *ϵ* has been established by previous work and is sufficiently large that required factors of (*ϵ* − 1)/*ϵ* are nearly unity. The dumbbell fluid is considerably less polarizable so that a more precise estimate of *ϵ* is needed. From a 5 ns bulk simulation of 512 dumbbell molecules (and no solute), we estimated the dielectric constant to be *ϵ* ≈ 7.1, obtained from the fluctuation-dissipation theorem,^{70}

where **M** is the total dipole moment of the simulation cell and $\beta =(kBT)\u22121$ (*k*_{B} is Boltzmann’s constant). We append the subscript “**E** = **0**” to emphasize that this formula holds with the use of tin foil boundary conditions (see, e.g., Refs. 70–72). Similar dumbbell models have been examined in previous studies,^{73} which reported similar values of *ϵ* for comparable molecular dipole moments.

All simulations were performed with the LAMMPS simulation package.^{74} Dynamics were propagated at constant volume *V* and temperature *T* = 298 K using Langevin dynamics as implemented in LAMMPS,^{75,76} with a time step of 1 fs. The SHAKE algorithm was used to constrain all bonds and angles.^{77,78} Long range electrostatic interactions were evaluated using Ewald summation with tin foil boundary conditions. (The solvent’s net dipole would experience different forces under a different choice of boundary conditions.^{70,79–84} The solvent response of interest, however, does not directly involve this dipole, so we expect the choice of boundary conditions to have minor significance.) Ewald sums were calculated using the particle-particle particle-mesh solver,^{85} 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 Å.^{86}

The solvation free energies $Fchg(L)(q)$ we consider are precisely defined by

where *ϕ*_{solv} is the electric potential generated by the solvent, evaluated at the center of a solute; *P*^{(L)}(*ϕ*_{solv}; *q*) is its probability distribution in the presence of a solute charge *q*; and $\u27e8\u22c5\u27e9q(L)$ denotes an average over *P*^{(L)}(*ϕ*_{solv}; *q*). The superscript (**L**) specifies the dimensions **L** = (*L*_{x}, *L*_{y}, *L*_{z}) of the simulated unit cell. In the limit **L** → ∞ of large simulation cell size, $Fchg(L)(q)$ approaches the reversible work of charging a solute at infinite dilution. Evaluating the integral in Eq. (2), which is a statement of the potential distribution theorem,^{87} requires knowledge of *P*^{(L)}(*ϕ*_{solv}; 0) in its extreme wings, which we obtained from umbrella sampling. Specifically, we introduced a series of solute charges *q*/*e* = −1.0, −0.9, …, 0, …, +0.9, +1.0, which bias typical solvent potential fluctuations to a corresponding series of ranges centered on $\varphi \xafsolv(q)$, where

Samples from the probability distributions $P(L)(\varphi solv;q)=P(L)(\varphi solv;0)e\u2212\beta q\varphi solv/\u27e8e\u2212\beta q\varphi solv\u27e90$ were then reweighted and combined according to the MBAR algorithm.^{88} The resulting statistics of extreme solvent fluctuations required to compute $Fchg(L)(q)$ could as well be obtained from alternative methods such as thermodynamic integration, reversible work calculation, or irreversible transformations via Jarzynski’s identity.^{66,89,90} The biasing method we used is especially convenient for molecular dynamics simulation, as it can be achieved simply by modulating the charge of a single particle.

The value of *ϕ*_{solv}(**R**^{N}) for a given configuration **R**^{N} of solvent molecules was evaluated by inserting a solute test charge, and then subtracting contributions due to the test charge itself. In detail, we compute the total potential energy *u*(**R**^{N}; *q*_{test}) of a system including solvent, a periodic collection of test charges *q*_{test}, and a charge neutralizing background. From this, we first subtract the energy *u*(**R**^{N}; 0) of interactions among solvent molecules and of Lennard-Jones interactions with the solute. We then subtract the energy of interactions among periodic test charges and the neutralizing background, *q*_{test}*ϕ*_{wig}/2, where *ϕ*_{wig} is the canonical Wigner potential discussed in Sec. III B and Appendix B. Dividing the remainder by *q*_{test} yields the solvent potential, $\varphi solv(RN)=qtest\u22121[u(RN;qtest)\u2212u(RN;0)\u2212qtest\varphi wig/2]$.

Charging free energies $Fchg(L)(q)$ computed with these methods are presented in Figs. 2 and 3, for bulk and interfacial systems, respectively. A strong dependence on periodicity **L** is evident for both model solvents considered. For cubic simulation cells, bulk free energies vary by 10s of *k*_{B}*T* as cell dimensions grow from roughly 1 to 2 nm. Bulk results for anisotropic cells differ even more dramatically, by 100s of *k*_{B}*T* as the aspect ratio *L*_{z}/*L*_{x} grows from 2 to 12, and appear not to converge as *L*_{z} → ∞ with fixed *L*_{x}. Similarly large finite size effects are exhibited by interfacial simulations. In this case, however, charging free energies do appear to converge in the limit of infinite aspect ratio, pointing to qualitative differences in the system size dependence of periodic bulk and slab geometries. Such nontrivial convergence behavior also warns that extrapolating results naively from ever-larger systems may not achieve physically meaningfully asymptotic properties.

## III. CORRECTING FOR FINITE SIZE EFFECTS ON IONIC SOLVATION FREE ENERGIES

This work is by no means the first to address the role of finite size effects in ionic solvation free energies obtained from molecular simulation,^{65–67,91,92} and indeed, the central assumption—that DCT adequately describes long-wavelength dielectric response—follows the perspective of Ref. 65. The majority of previous work that addresses finite size effects has, however, focused on a scenario that lacks macroscopic dielectric boundaries, such as the liquid/vapor interface. Here we present a general formalism that is also applicable in such cases.

The quantity that we seek to calculate is the free energy to introduce a charge *q* to the center of a neutral cavity embedded in a solvent with dielectric constant *ϵ*, at infinite dilution,

*F*^{(∞)}(*q*) is the free energy of a macroscopic system in the presence of a single solute with charge *q*. Molecular simulations cannot directly access $Fchg(\u221e)(q)$, but can instead provide corresponding free energies under PBC,

As in Sec. II, the superscript (**L**) indicates spatial periodicity **L** = (*L*_{x}, *L*_{y}, *L*_{z}) in three Cartesian directions. Our goal is to connect the charging free energies for periodic and macroscopic systems, i.e., to calculate a correction

that can be applied to simulation results as an extrapolation to infinite dilution.

We consider two distinct contributions to Δ*F*(**L**; *q*). The first accounts for size dependence of a solvent’s structural response to solute charging, which we estimate using dielectric continuum theory. The second considers biases on ion solvation that exist even in the unperturbed solvent, specifically electric potentials experienced by a neutral solute, which are also size-dependent.

### A. Size-dependent dielectric response

DCT is defined by a simple linear response relationship between the electric field at a point **r** and the average induced polarization density **m**(**r**) at that location. Together with Poisson’s equation, this relation connects spatial variations in **m**(**r**) to the density *ρ*(**r**) of free charge,

(Throughout this formulation, we work in a unit system in which 4*πϵ*_{0} = 1, where *ϵ*_{0} is the permittivity of free space.) In our case, the free charge represents periodically arranged solutes and their compensating uniform background of charge density $\rho \xaf=q/(LxLyLz)$,

Here, *ρ*_{0}(**r**) is the charge distribution of a single solute, localized at position **r**_{0} inside the dielectric and with net charge *q* = ∫d**r***ρ*_{0}(**r**). The solute’s periodic replicas are separated by lattice vectors **b** = (*n*_{x}*L*_{x}, *n*_{y}*L*_{y}, *n*_{z}*L*_{z}), with *n*_{x}, *n*_{y}, and *n*_{z} taking on all integer values. If dielectric boundaries are present, we take them to have the same spatial periodicity as the solutes, as well as a clear dependence on **L**. For the periodic slab geometry, these boundaries are a pair of infinite parallel planes with fixed separation *w*, repeated in the *z* direction with periodicity *L*_{z} (see Fig. 1).

To calculate the charging free energy $Fchg(L)(q)$ for this theory, we exploit a generic property of linear response, namely, that the change in free energy due to an external field is half the average interaction energy between the system and field. For a dielectric material interacting with periodic solutes and their neutralizing background,

where *ϕ*_{solv}(**r**) is the electric potential generated by the polarization field **m**(**r**). Because the solvent is electrically neutral, its interaction with the uniform background charge does not contribute to Eq. (9). Determining a finite size correction Δ*F*_{DCT}(**L**; *q*) from dielectric continuum theory thus amounts to calculating the average solvent potential at the center of a solute.

We obtain the solvent potential by integrating over the volume Ω occupied by the solvent,

This volume is periodic by construction but it need not be connected. If Ω comprises disconnected regions, then a periodic series of surfaces *∂*Ω bounds these regions. The unit normal vector on these boundaries we denote by $n^$. Integrating by parts and using the divergence theorem gives

where **R** is a point on the set of boundaries *∂*Ω. The normal component of solvent polarization at a boundary, $\sigma (R;L)\u2261n^\u22c5m(R)$, serves as an effective surface charge density due to polarization. Using the basic linear response relation in Eq. (7), we finally obtain

The correction Δ*F*_{DCT}(**L**; *q*) compares response at finite **L** with the macroscopic result. In the limit **L** → ∞, all solute replicas described by *ρ*(**r**) are irrelevant except for the central charge at **r**_{0}, and the background charge density vanishes. As a result,

where Δ*σ*(**R**) = *σ*(**R**; **L**) − *σ*(**R**; ∞) is a difference between periodic and macroscopic systems. Motivated by the slab geometry of interest, we have assumed that there is a clear correspondence between boundaries for finite and infinite **L**. The **L** → ∞ slab possesses just two planar surfaces (neglecting the solute’s excluded volume^{93}), bounding the solvent from above and below; these two surfaces are also present at finite **L**. On all other boundaries included in *∂*Ω, there is no macroscopic contribution, i.e., *σ*(**R**; ∞) = 0. A more complicated set of periodic boundaries might not permit this simplification, requiring instead an explicit subtraction of surface integrals evaluated on different sets of boundaries.

For any choice of boundaries, evaluating Eq. (13) requires knowledge of the function *σ*(**R**; **L**). We have not solved a general dielectric boundary value problem, but have instead rephrased it in a convenient way. For the slab geometry, *σ*(**R**; **L**) can be worked out exactly in terms of a series of image charges, as detailed in Appendix A. More importantly, Δ*σ*(**R**) is amenable to a greatly simplifying approximation, which we will show to be very accurate in the case of periodic slabs. For a bulk system lacking boundaries, Eq. (13) does constitute a full solution. In Sec. III B, we show that this result is consistent with previous work on ion solvation in bulk solvents.

### B. Application to ion solvation in the absence of interfaces

Understanding the behavior of Δ*F*(**L**; *q*) for “bulk” periodic systems—those lacking any macroscopic dielectric boundary—has been the focus of much research. In an early study on the subject, Hummer *et al.*^{66} argued that estimates of charging free energy should include interactions of the ion with its periodic replicas and background charge. For simulations of liquid water using Ewald summation and a cubic cell of side length *L*, they showed that incorporating these “self-interactions” can yield free energies that are essentially independent of system size. The finite size correction is thus well approximated in this case by

where *ϕ*_{wig}/*q* ≈ −2.837 297/*L* (see, e.g., Ref. 94, and note that atomic units are implied in the numerical value). The Wigner potential *ϕ*_{wig}/*q* is defined as the electric potential at the site of a unit point charge due to all of its periodic replicas and a homogeneous background charge that acts to neutralize the primitive cell. (*ϕ*_{wig}/*q* is commonly referred to as *ξ*_{EW} in the literature.) Based on a dimensional analysis, Figueirido *et al.*^{91} proposed the following correction for finite dielectric constant *ϵ*:

which has the correct limiting behavior for *ϵ* → 1. Later studies extended these corrections to also account for the finite size of the ion.^{65,67} We do not consider such higher order corrections in this article.

In the absence of interfaces, our DCT expression in Eq. (13) simplifies to

For a cubic primitive cell (*L*_{x} = *L*_{y} = *L*_{z} = *L*), this result is equivalent to Figueirido’s, as expected from the DCT analysis of Ref. 65. Figure 2 exemplifies the remarkable effectiveness of this correction, echoing the conclusions of Hummer *et al.* Here we have added Δ*F*_{DCT}(**L**; *q*) to the values of $Fchg(L)(q)$ obtained from bulk molecular simulations of water (d) and the simple polar fluid (b). The DCT correction reconciles solvation free energies from simulations with significantly different periodicity, reducing discrepancies of ∼25*k*_{B}*T* down to ∼1*k*_{B}*T*. Residual finite size effects beyond Δ*F*_{DCT}(**L**; *q*) are noticeable only for the smallest simulations of the simple polar liquid. We have confirmed numerically that these differences can be largely removed with higher-order DCT corrections that account for the solute’s excluded volume.^{65,67}

For an anisotropic cell, Eq. (16) is a very straightforward generalization of the standard correction. To evaluate it numerically, we have used Ewald summation as outlined in Appendix B. In Fig. 2, we include simulation results for cuboidal simulation cells with *L*_{z}/*L*_{x} = 2 and 12, demonstrating that the DCT correction remains accurate even as the cell’s aspect ratio becomes large. The interesting aspect of this comparison is that Δ*F*_{DCT}(**L**; *q*) diverges as *L*_{z}/*L*_{x} → ∞ with *L*_{x}*L*_{y} fixed. The accuracy of our correction then suggests that $Fchg(L)(q)$ diverges in the opposite sense with growing aspect ratio. Indeed, simulation results with 756 molecules and *L*_{z}/*L*_{x} = 12 extend far beyond the plotted energy scale [$\beta Fchg(L)(\u2212e)\u2248\u2212408$], hinting at such a divergence.

### C. Application to periodic slabs

Applying the general result in Eq. (13) to systems with boundaries is considerably more difficult, as it requires calculation of the polarization surface charge *σ*(**R**; **L**) [as well as integration of *ρ*(**r**)/|**r** − **r**_{0}| over disconnected volumes]. The spatial variation of *σ*(**R**; **L**) encodes important driving forces of solvation. A point ion near a dielectric boundary experiences an electric field that grows sharply in magnitude as it approaches the surface. For a semi-infinite dielectric with a planar boundary, this field is equivalent in DCT to the force generated by an image charge opposite to the boundary, accompanied by strong spatial variation in the polarization surface charge. A poor description of this position dependence would yield a poor approximation to $\u27e8\varphi solv(r0)\u27e9q(L)$. This feature of *σ*(**R**; **L**), however, is predominantly local in space and therefore only weakly sensitive to PBC. Moreover, such local aspects of solvation are likely not well described by DCT, which lacks microscopic complexities that can dominate on small scales. The *difference* quantity Δ*σ*(**R**; **L**), we argue, should have much weaker spatial dependence and should be well described by DCT.

For the specific case of periodic slabs, a solution for *σ*(**R**; **L**) can be obtained by summing contributions from an infinite series of effective image charges, as described in Appendix A. This solution exhibits the rapid spatial variation described above, particularly as the solute ion nears the dielectric interface. But as anticipated, the image charges that contribute most strongly to this position dependence are identical in value and placement when **L** → ∞. The finite size effects of interest to our work originate instead in the long-wavelength features of polarization surface charge. We thus propose a greatly simplifying approximation, namely, that $\Delta \sigma (R)\u2248\sigma \xaf$ is constant along the planar boundaries of each dielectric slab. (This argument relies on the ion residing within the solvent’s liquid phase. As discussed in Appendix A, the approximation is inappropriate when **r**_{0} lies in the vapor phase instead.)

An appropriate value for $\sigma \xaf$ can be determined by integrating over a single periodic replica *∂ω* of the slab’s boundaries,^{95}

where *ω* denotes a single periodic replica of the slab (see Fig. 1). In using the divergence theorem, we have exploited the fact that symmetry and continuity of **m**(**r**) within the slab requires that $n^\u22c5m(R)=0$ at the lateral interfaces between periodic replicas. The use of Eqs. (7) and (8) then yields

where $w$ is the slab’s thickness and *A* = *L*_{x}*L*_{y} is the area of a single periodic replica of the upper (or lower) dielectric boundary. This value of $\sigma \xaf$ is precisely what is required for Δ*F*_{DCT}(**L**; *q*) to be finite.

The finite size correction $\Delta FDCT(unif)(L;q)$ that results from this uniform surface charge approximation, together with Eq. (18),

has a simple physical interpretation. The quantity in square brackets in Eq. (19) corresponds to an electric potential $\varphi wig*$ generated by (i) periodically arranged unit point charges (excepting the central replica); (ii) a partially neutralizing charge density that is uniform within the dielectric slabs and vanishes outside; and (iii) uniformly charged plates at the dielectric boundaries, which achieve overall neutralization. It can be regarded as a generalization of the Wigner potential, in which the uniform compensating charge outside the dielectric has been moved to the dielectric boundaries. Interestingly, in the limit *L*_{z} → ∞ with *A* and $w$ fixed, all of the compensating charge moves to these boundaries. This scenario—a single dielectric slab, with ions periodically replicated in the *x* and *y* directions, bounded above and below by uniformly charged neutralizing plates—is markedly similar to the neutralization scheme proposed by Ballenegger *et al.*^{96} for molecular dynamics simulations of charged slab systems that are aperiodic in the *z* direction.

The modification of the Wigner potential described above can be calculated exactly in closed form, as shown in Appendix B. The DCT finite size correction can then be written simply as

In the limit that *L*_{x}, *L*_{y}, and *L*_{z} all become infinite, with $w$ still fixed, $\Delta FDCT(unif)(L;q)\u21920$, an important property of the full solution for Δ*F*_{DCT} that is preserved by the uniform surface charge approximation. Note that $\Delta FDCT(unif)$ does not depend upon the solute’s position within the slab.

Applying this correction to simulation results for periodic slabs of the simple polar fluid, we can account almost completely for the substantial variation of charging free energy $Fchg(L)$ with system size. Figure 3(b) shows results for *L*_{x} = *L*_{y} ≈ 12.4 Å, slab width $w$ ≈ 27.5 Å, and various aspect ratios *L*_{z}/*L*_{x}, with the solute located equidistant from the liquid slab’s two planar boundaries. Corrected free energies for these systems vary by less than 1.05 *k*_{B}*T*.

Aqueous simulations exhibit a more complicated size dependence, which is only partially captured by $\Delta FDCT(unif)(L;q)$. Charging free energies corrected according to Eq. (20) are shown in Fig. 3(d) for *L*_{x} = *L*_{y} ≈ 12.4 Å, slab width $w$ ≈ 29.2 Å, and various aspect ratios *L*_{z}/*L*_{x}, with the solute placed equidistant from the two interfaces. The variation of $Fchg(L)(q)$ with *L*_{z} is greatly reduced by $\Delta FDCT(unif)(L;q)$, but a systematic size dependence remains, leaving differences as large as 13.4 *k*_{B}*T* unexplained. A breakdown of the uniform surface charge approximation cannot alone explain this failure since the remnant size dependence is different for cations and anions—a charge asymmetry that cannot be described by DCT. The source of this charge-asymmetric finite size effect is considered in Sec. III D.

For the sake of simplicity, Fig. 3 presents results only for the case *d* = $w$/2, i.e., with the ion situated midway between the liquid slab’s boundaries. In Sec. III D, we show that the corrections derived above are equally accurate for other values of *d*, provided the solute remains in the liquid slab.

### D. Finite size effects due to charge asymmetry

The simple polar fluid we have considered is charge symmetric, in the sense that the positive and negative charge is equivalently distributed within each dumbbell molecule. By contrast, the SPC/E model of water is charge asymmetric, with positive charges situated farther from the center of volume exclusion. Such asymmetry can lead to a physically meaningful discrimination between cations and anions. In simulations with PBC, it can also generate unphysical differences. For example, a solvent comprising volume-excluding molecules whose internal charge distribution is spherically symmetric cannot respond to the charging of a solute, yet a bulk simulation of such particles under PBC can yield a nonzero charging free energy that is different for cations and anions.^{97,98} This artifact reflects a sensitivity of *ϕ*_{solv}(**r**_{0}) to the existence of interfaces even when they are arbitrarily far away.

Bulk simulations with PBC lack interfaces regardless of periodicity so that artifacts in $Fchg(L)$ persist even in the limit **L** → ∞. The slab geometry is an intermediate case, with explicit interfaces but also artifacts and finite size effects of charge asymmetry that vanish in the **L** → ∞ limit. To demonstrate this fact, we examine the average electric potential *ϕ*_{neut} experienced by a *neutral* solute in a simulation of periodic slabs. The dependence of this “neutral cavity potential” on system size provides a way to extrapolate computed charging free energies to the macroscopic limit.

The neutral cavity potential can be written exactly as

where $v$ denotes a single unit cell of the periodic slab simulation. $\u27e8\rho solv(r)\u27e90(L)$ is the solvent’s average charge density at position **r**, with a neutral solute at position **r**_{0}. Due to volume exclusion, ⟨*ρ*_{solv}(**r**)⟩ vanishes inside the solute (and its periodic replicas). If the solute resides within the slab, this constraint requires that ⟨*ρ*_{solv}(**r**)⟩ depends on *x* and *y* as well as *z*, a significant feature for contributions to *ϕ*_{neut} from the primary simulation cell (**b** = **0**). Our interest lies, however, in contributions from more distant unit cells, which dictate the system size dependence of *ϕ*_{neut}. These distant contributions are much less sensitive to the solvent’s electrostatic inhomogeneity in *x* and *y*. We thus expect that

is a good approximation, where *const* is independent of **L**, and the quantity $\u27e8\rho solv(z)\u27e90$ has been spatially averaged over *x* and *y*. Figure 4(b) shows this laterally averaged charge density, obtained from simulations of SPC/E water with a neutral solute. Results for two different solute locations (one at the middle of the slab and the other at the liquid/vapor interface) are nearly indistinguishable. Provided the solute’s cross-sectional area *πσ*^{2}/4 is not a substantial fraction of *L*_{x}*L*_{y}, we expect the laterally averaged charge density to be insensitive to the solute’s presence at all.

Removal of the dependence on the lateral coordinates allows for great simplification. As shown in Appendix B, exploiting the symmetries of $\u27e8\rho solv(z)\u27e90$ yields

The *difference* quantity

can now be written compactly as

where we have eliminated terms that vanish due either to electroneutrality of the solvent or to the zero net dipole moment of a pure liquid slab. Numerical results for Δ*ϕ*_{neut}(**L**) are shown in Fig. 4(a), obtained directly from aqueous simulations and also from Eq. (25). They indicate that neglecting the lateral dependence of $\u27e8\rho solv(r)\u27e90(L)$ is indeed a very good approximation.

Equation (25) describes an intrinsic bias on ion solvation in a system of periodic slabs, one that vanishes in the limit *L*_{z} → ∞. We take this bias, in the absence of any solute charge, as a baseline for dielectric response, and therefore propose a full finite size correction for the slab geometry,

Equations (20) and (25) provide accurate, inexpensive, and easily implemented approximations for Δ*F*_{DCT}(**L**; *q*) and Δ*ϕ*_{neut}(**L**), respectively. Charging free energies corrected according to Eq. (26) are shown in Fig. 5(a) for an ion at the center of a periodic slab of liquid water. The remnant size dependence of $Fchg(L)+\Delta FDCT$ evident in Fig. 3(d) is captured almost exactly by the neutral cavity potential. Equation (26) accounts equally well for finite size effects observed for an ion placed at the air/water interface, as shown in Fig. 5(b).

Several features of the neutral cavity potential are noteworthy. First, within the approximation leading to Eq. (25), Δ*ϕ*_{neut}(**L**) is independent of the solute’s location, even for points well outside the slab. This invariance suggests that the bias we are correcting is dominated by very distant contributions, in effect a consequence of omitting each slab’s boundary at *x*, *y* = ∞ under PBC. Indeed, Eq. (25) resembles the conventional surface potential $\varphi surf=4\pi \u222bzvapzliq\u2009dz\u2009\u27e8\rho solv(z)\u27e9z$, where *z*_{liq} and *z*_{vap} denote locations on either side of a neat liquid/vapor interface (*z*_{vap} < *z*_{liq}). These potentials are simply related within a molecular multipole expansion up to quadrupole order, assuming that the solvent’s average polarization field is nonzero only at the interface and that its quadrupole field is uniform within the liquid phase. The neutral cavity potential then evaluates simply to Δ*ϕ*_{neut}(**L**) ≈ *ϕ*_{surf}$w$/*L*_{z}. From this relationship, our numerical results for Δ*ϕ*_{neut} suggest a surface potential of roughly −0.54 V for the SPC/E model, comparable to published values.^{12,98,100–102} Like the conventional surface potential, Δ*ϕ*_{neut}(**L**) may therefore have little to do with the molecular physics of solvation, as underscored by the possibility that Δ*ϕ*_{neut}(**L**) ≠ 0 even for spherically symmetric solvent molecules. In computing the solvent potential from a lattice sum method, we effectively adopt the so-called “P-summation scheme” (see, e.g., Refs. 73, 89, and 103) such that the surface potential contains unphysical contributions from the internal charge distribution of the solvent molecules.

For bulk liquid simulations, such artifacts of PBC are not removed by extrapolation to the macroscopic limit. For the slab geometry, Eq. (22) indicates that these artifacts are completely removed as *L*_{z} → ∞. Specifically, if the neutral solute resides in the vapor phase,^{104} **r**_{0} = **r**_{vap}, then Eq. (22) is exact, with a vanishing constant term. Integrating Eq. (23) then gives $\u27e8\varphi solv(rvap)\u27e90(L)=0$. In other words, solvent potentials are properly referenced to the vapor phase when the spacing between slabs is infinite. Based on this fact, one convenient way to evaluate Δ*ϕ*_{neut}(**L**) in slab simulations is simply to compute the average electric potential at a point within the vapor phase (and far from the interface). Also shown in Fig. 5(a) are the results from a slab simulation employing the popular Yeh-Berkowitz method^{99} for simulating slab systems, in conjunction with the neutralization scheme outlined in Ref. 96. (In this case Δ*ϕ*_{neut} = 0, and the generalized Wigner potential has been evaluated numerically in the limit *L*_{z} → ∞.) Good agreement is observed between these results and those obtained with the *a posteriori* correction derived in this work.

### E. Reconciling bulk and slab results

For the periodic slab geometry, with ions at the slabs’ centers, the limit **L** → ∞ we have considered corresponds to an individual ion in the middle of an individual slab with infinite lateral extent. This system of course still retains a finite dimension, namely, the slab’s thickness $w$. Only as $w$ grows should we expect $Fchg(\u221e)$ to approach the charging free energy for a macroscopic solution at infinite dilution. For the simple polar fluid, this remaining finite size effect is evident in Fig. 3(b), where deviations are visible between $Fchg(L)+\Delta FDCT(L)$ for periodic slab and bulk systems. A correction for finite $w$ can be estimated from DCT, as described in Appendix A, yielding an estimate of the macroscopic charging free energy away from any interface

where $Fchg(slab)$ is the result of periodic slab simulations with the corresponding finite periodicity correction Δ*F*_{DCT} + *q*Δ*ϕ*_{neut} applied. As shown in Fig. 6(a), this adjustment brings slab and bulk simulation results into good agreement for the simple polar solvent.

The same finite thickness correction should apply to the aqueous slab simulations as well. It is not sufficient, however, to reconcile differences between slab and bulk results in this case. Indeed, it is clear from Fig. 5(a) that these deviations are charge asymmetric, and thus cannot be accounted by DCT. We argue that the remaining discrepancy is due not to finite slab thickness but instead to artifacts of bulk periodic simulations. As discussed in Sec. III D, such bulk periodic systems lack interfaces even in the limit **L** → ∞, which can generate artifacts correctible only through interfacial considerations. (See Ref. 105 for a recent overview of this issue and of strategies to correct for it.) Specifically, a neutral solute at the center of a macroscopic spherical liquid droplet experiences an average electric potential,

where *r*_{liq} denotes a point in the liquid far from both the solute and the liquid/vapor interface. The first term in Eq. (29) should be well described by a bulk periodic simulation as **L** → ∞, while the latter contribution is absent. Following the reasoning of Sec. III D, we account for the interfaces missing in a bulk periodic simulation by adding *qϕ*_{surf} to the computed charging free energy. Our estimate of $Fchg(macro)$ from bulk periodic simulations is thus

where $Fchg(bulk)$ is the result of periodic bulk simulations with the corresponding finite periodicity correction Δ*F*_{DCT} applied. Numerically integrating $\varphi surf=4\pi \u222bzvapzliq\u2009dz\u2009\u27e8\rho solv(z)\u27e9z$ for all of our neutral cavity simulations (using the “P-summation” convention as described above), we obtain values of *ϕ*_{surf} ranging from −0.53 V to −0.58 V, comparable to previously reported surface potential values^{12,98,100–102} for point charge models, as well as our estimate based on Δ*ϕ*_{neut}. Figure 6(b) shows the macroscopic estimates for SPC/E water obtained from Eq. (30) with *ϕ*_{surf} = −0.58 V, a value that yields particularly close alignment of bulk and slab results. This agreement suggests that the charging free energies in Fig. 6(b) are robust assessments of ion solvation in the physically realistic scenario of a macroscopic solution with boundaries. The full macroscopic solvation free energy would require adding the thermodynamic cost to introduce a neutral solute, which should not face significant finite size effects.

## IV. DISCUSSION

Many simulation studies have examined small ions near aqueous interfaces using a slab geometry. For system sizes typical of these studies, we have demonstrated finite size effects that can shift solvation free energies by 10s of *k*_{B}*T*. While substantial in this sense, these effects are unlikely to change many qualitative conclusions from previous work. In particular, the finite size corrections of Eqs. (20) and (25) do not depend on the ion’s location within the liquid. They are therefore insensitive to a solute’s approach to the interface. Previously computed density and free energy profiles should therefore be unaffected for solute positions on the liquid side of the interface; in this case, our corrections simply amount to resetting the zero of energy. Our finite size corrections do not apply well when an ion resides outside the liquid, both because the uniform surface charge approximation of Eq. (19) fails and because non-dielectric response such as interface deformation is pronounced. Empirically, we find in this case that aligning simulation results from different system sizes requires a correction that *does* depend on the ion’s position relative to the interface. Capturing that dependence is a challenge that likely requires an approach more nuanced than DCT.

Our results for two model molecular liquids demonstrate that DCT can predict the system size dependence of interfacial ion solvation thermodynamics with striking accuracy. This success is remarkable in light of the length and energy scales involved and the complexities of interfacial response. The smallest systems we have investigated are periodic on a scale of just a few molecular diameters, well below the scale on which the liquid appears continuous. The scale of charging free energies in Fig. 6 underscores that the fluctuations governing solvation of a fully charged ion are extraordinarily rare (with relative probability <*e*^{−100}) in its absence. In the case of water, the overall response to charging a solute is markedly nonlinear. Both in bulk water and in aqueous slabs, *F*_{chg}(*q*) deviates strongly from the parabolic form characteristic of linear response. [Most notably, the curvature *d*^{2}*F*_{chg}/*dq*^{2} differs significantly at positive and negative *q* so that anions are far more favorably solvated than cations.^{66,106–109} A contribution linear in *q*, which favors cation solvation, is much weaker for appreciably charged solutes (|*q*| ≳ 0.5).^{73}] The softness of liquid/vapor interfaces and the altered hydrogen bond network structure at the liquid’s boundary contribute additional nonlinearities that shape ions’ association with the surface.^{10,18,31,110,111} And yet a continuum linear response theory describes nm-scale thermodynamic finite size effects to an accuracy of ≲*k*_{B}*T*. From the non-parabolic shape of *F*_{chg}(*q*), we know that this success of DCT cannot hold on all length scales. Our results indicate that charge-asymmetric nonlinear response in such water models is confined to the very near field of one or two solvation shells.

Our formulation of finite periodicity corrections for inhomogeneous geometries is motivated by a specific interest in how liquid water accommodates solutes near interfaces, and it builds directly on previous work from the molecular physics community. But similar finite size effects have been investigated in other areas, e.g., condensed matter physics (see, e.g., Refs. 112–116). The results and methods described here may thus be useful in a range of contexts involving polarization response, such as charged defect formation in semiconductors, and obtaining reliable reference potentials in computational electrochemistry.^{117,118}

## V. CONCLUSIONS

In this article, we have advanced a physical perspective and a mathematical framework for assessing the impact of spatial periodicity on ion solvation free energies computed from molecular simulation. Our general result Eq. (13) is a finite size correction from dielectric continuum theory, suitable for application to systems with boundaries where the local polarizability changes sharply. It requires as input long-wavelength features of induced surface polarization charge. For the case of a charged solute in a dielectric slab under periodic boundary conditions, taking this boundary charge to be uniform is an excellent approximation, which yields Eq. (20) as an easily computable correction.

For solvent molecules that are charge-asymmetric (such as water), an additional consideration outside the scope of DCT is needed to describe finite size effects in the slab geometry. Related to the surface potential, the resulting correction in Eq. (25) removes artifacts associated with the absence of distant interfaces under PBC. As *L*_{z} → ∞, this absence becomes inconsequential. By contrast, for periodic bulk simulations, the absence of interfaces is problematic even in the limit **L** → ∞. Using the surface potential to effectively restore those interfaces, we find that fully corrected bulk and slab results approach a consistent macroscopic limit.

We regard these successes as evidence that the electrostatic response of aqueous environments is well described by DCT down to length scales little larger than a molecular diameter. Theoretical work by others on the association of ions with the air/water interface has asserted that the realism of DCT extends even further in that context,^{6,8,9,19,27,119} in essence to arbitrarily small scales. The methods we have described facilitate a careful scrutiny of this hypothesis, which we pursue in ongoing work.

## ACKNOWLEDGMENTS

We would like to thank Dayton G. Thorpe and Layne B. Frechette for many fruitful discussions. The work was 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 No. DE-AC02-05CH11231.

### APPENDIX A: METHOD OF IMAGES APPLIED TO SLAB SYSTEMS

Calculating the dielectric response to a point charge near a planar boundary between two semi-infinite polarizable media is a textbook exercise.^{120} Focusing on solvation near liquid-vapor interfaces, we take one of the regions to be vacuum (*ϵ* = 1) and the other region Ω to have a large dielectric constant *ϵ*. For a point charge *q* that is located in the high-dielectric medium at a point **r**_{0}, lying a distance *d* from the interface, Poisson’s equation can be easily solved using the *method of images*. The average electric field within Ω is constructed with an effective external point charge, whose field is divergence-free inside Ω and enforces dielectric boundary conditions at the interface *∂*Ω. The electric potential *ϕ*(**r**) is then specified by the requirement that *ϕ* vanishes at infinite distance from the real charge. For a point **r** within Ω,

where

is the potential generated by an image charge *aq* at a point $r0+2dz^$ that mirrors the location of the real charge, and *a* = (*ϵ* − 1)/(*ϵ* + 1). The unit vector $z^$ is perpendicular to the interface, pointing outwards from the dielectric.

Solvation of a point charge in a dielectric slab, bounded by two planar interfaces with vapor, can be similarly solved. In this case, satisfying boundary conditions on both interfaces requires an infinite series of image charges, as depicted in Fig. 7. Here the real charge lies a distance *d*_{1} from its closest interface (the top interface in Fig. 7), and a distance *d*_{2} = $w$ − *d*_{1} from the other interface (the lower interface in Fig. 7), where $w$ is the slab thickness. Above the slab reside image charges: (i) *aq*, *a*^{3}*q*, *a*^{5}*q*, … at positions $r0+2d1z^,r0+2(w+d1)z^,r0+2(2w+d1)z^,\u2026$, and (ii) *a*^{2}*q*, *a*^{4}*q*, *a*^{6}*q*, … at positions $r0+2wz^,r0+4wz^,r0+6wz^,\u2026$, where $z^$ is an outward normal for the *upper* interface. Below the slab reside image charges: (i) *aq*, *a*^{3}*q*, *a*^{5}*q*, … at positions $r0\u22122d2z^,r0\u22122(w+d2)z^,r0\u22122(2w+d2)z^,\u2026$, and (ii) *a*^{2}*q*, *a*^{4}*q*, *a*^{6}*q*, … at positions $r0\u22122wz^,r0\u22124wz^,r0\u22126wz^,\u2026$. Since the electric potential must again vanish as **r** → ∞, *ϕ*(**r**) is still given by Eq. (A1), now with

Finite size corrections in Sec. III C involve the potential *ϕ*(**r**) evaluated at the solute’s position. For this case of a single ion in a single slab, *ϕ*(**r**_{0}) can be written exactly in terms of the Hurwitz-Lerch transcendent $\Phi nHL(z,b)$,

The potential *ϕ*_{self} due to the point ion itself is singular, but it is independent of the boundaries’ geometry and therefore cancels in all quantities of interest in the main text. In the special case that the ion is equidistant from the interfaces (*d*_{1} = *d*_{2} = $w$/2),

In Fig. 8, we plot the non-singular part of this potential, *ϕ*_{image}(**r**_{0}), for various slab widths.

The case of central interest in this paper—the periodic slab simulation geometry depicted in Fig. 1—is more complicated in two respects. First, image charges that determine the electric field $E(r)=\u2212\u2207\varphi (r)$ do not follow a simple pattern. Second, *ϕ*(**r**) need not vanish as **r** → ∞ so that additional information is needed to fully specify a boundary value problem for the electric potential. As a result, the collection of image charges used to solve $\u2207\u22c5E=4\pi \rho (r)/\u03f5$ determine *ϕ*(**r**) only up to an additive constant. Our approach circumvents this issue [see Eq. (12)] by calculating the potential directly from the polarization surface charge density *σ*(**R**), which is in turn completely specified by the electric field,

A full solution in this scheme would first determine image charges that enforce boundary conditions for $E(r)$, calculate the resulting surface charge as a function of position **R** on each boundary, and then integrate as per Eq. (12) to obtain the electric potential. This procedure would be tedious but feasible. For the physical systems we are considering, which feature an ion inside a high-dielectric solvent that coexists with vapor, an accurate approximation can be computed much more easily, as explained below.

For a periodic collection of point charges in dielectric slabs, the method of images gives a series of effective charges exemplified by Fig. 9. For a given real charge, one or more of the nearest images are identical to the non-periodic case of Fig. 7 (assuming *L*_{z} ≥ 3$w$/2). The charges and positions of more distant images are complicated by the system’s periodic structure, although simple recursion relations can be derived when *L*_{z} is an integer multiple of $w$. The infinite collection of real charges also necessitates introducing a neutralizing background charge, which generates a series of slabs of image charge density (not shown). The background charge and its images contribute importantly to $E(r)$ and thus also to *σ*(**R**). Those background contributions, however, are spatially uniform in *x* and *y*, i.e., they are independent of **R** on each boundary. We exploit this fact by splitting $\sigma (R)=\sigma \xaf+\delta \sigma (R)$ into a zero-wavevector component $\sigma \xaf$ and a component *δσ*(**R**) that varies with **R** and integrates to zero. The average surface charge density $\sigma \xaf$ can be calculated simply and without reference to image charges, as described in Sec. III C. The remainder *δσ*(**R**) receives no contribution from the background charge or its images.

Results presented in Sec. III C were obtained neglecting the lateral variation of polarization surface charge, in effect setting *δσ* to zero. The justifications for this uniform surface charge approximation are several-fold. First, the most significant contributions to *F*_{chg} through *δσ* come from the nearest image charges, which do not depend on *L*_{z} and thus cancel in the calculation of Δ*F*_{DCT}. Contributions to *σ*(**R**) from more distant images are smaller in magnitude and, due to the greater distance, significantly dominated by their impact on the average charge density, for which $\sigma \xaf$ already accounts. Finally, contributions to the surface polarization charge from nearby pairs of real and/or image charges that are equidistant from a boundary are typically of order 1/*ϵ* for large *ϵ*. For example, the physical charge *q* is largely offset by an image charge *aq* opposite to the boundary, with the residual (1 − *a*)*q* = 2*q*/(*ϵ* + 1) ∼ 1/*ϵ*. By contrast, the average charge density $\sigma \xaf\u223c1\u22121/\u03f5$ is of order unity.

The uniform surface charge approximation is less defensible when the solute ion resides instead in the low-dieletric vapor phase. Such a solute can experience significant solvation forces from multiple periodic replicas of the liquid slab; contributions from only one of these replicas will cancel in the calculation of Δ*F*_{DCT}. Contributions from nearby charges also do not offset to order 1/*ϵ* in this case. For example, the image charge nearest the interface contributes to *δσ* with magnitude (*ϵ* − 1)/(*ϵ* + 1) ∼ 1 and does not have an offsetting counterpart opposite the boundary. Empirically, we find that finite size scaling predicted from this approximation is indeed not closely followed by simulations. The discrepancy likely has multiple sources. In addition to $\Delta FDCT(unif)$ poorly approximating Δ*F*_{DCT}, we expect that linear dielectric response is a much cruder model in this case. As small ions leave the liquid phase, they can deform the interface substantially,^{110,111,121} a response that lies distinctly outside the scope of DCT.

### APPENDIX B: GENERALIZING THE WIGNER POTENTIAL

Consider a periodic charge distribution *c*(**r**) that is (i) electroneutral, $\u222bv$*d***r***c*(**r**) = 0, (ii) non-dipolar, $\u222bv$*d***r***c*(**r**) **r** = 0, and (iii) spatially uniform in two directions *x* and *y*, *c*(**r**) = *c*(*z*). The net potential *ϕ*(**r**_{0}) generated by this distribution at observation point **r**_{0} is

Exploiting the lateral uniformity of *c*(**r**), this potential may be written as a Fourier series in just one dimension,

where $\u0109(kz)=\u222b\u2212L/2L/2dz\u2009c(z)eikzz$, *k*_{z} = 2*πn*/*L*_{z}, and the sum runs over all nonzero integers *n*. This sum can be rewritten as

in terms of a function *J*(*z*) that can be expressed in closed form,

The second equality in Eq. (B4) can be shown by twice differentiating *J*(*z*), noting that $Lz\u22121\u2211kzeikzz=\delta (z)$, and solving the ordinary differential equation that results. Doing so is aided by recognizing that *d*^{2}|*z*|/*dz*^{2} = 2*δ*(*z*), that *J*(0) = *πL*_{z}/3, and that *J*(*z*) is an even function of *z*. The assumed symmetries of *c*(*z*) allow further simplification,

Because the neutralizing charge distributions encountered in Sec. III C are uniform in *x* and *y*, their contributions to Δ*F*_{DCT} can be simplified using Eq. (B5). In particular, the generalized Wigner potential $\varphi wig*$ that arises for periodic dielectric slabs differs from the conventional Wigner potential *ϕ*_{wig} only in the form of the neutralizing background. Defining

we can regard Δ*ϕ** as the potential due to periodic, uniformly charged plates with charge density $\sigma \xaf$, together with a compensating background that is uniform *outside* the dielectric slabs and zero inside. This charge distribution has precisely the form assumed for *c*(**r**) above, with

where *h*_{out}(*z*) = 1 outside the dielectric slabs (i.e., for −*L*_{z}/2 < *z* < − $w$/2 and for $w$/2 < *z* < *L*_{z}/2) and *h*_{out}(*z*) = 0 inside (i.e., for −$w$/2 < *z* < $w$/2). For an observation point inside one of the dielectric slabs (−$w$/2 < *z*_{0} < *z*/2), substituting Eq. (B7) into Eq. (B5) yields, after simple integration,

Note that this result is independent of the observation point *z*_{0}, provided it lies within the dielectric. Note also that Δ*ϕ** vanishes when *L*_{z} = $w$ so that the standard Wigner potential is recovered for a bulk periodic system.

The conventional Wigner potential, arising from periodic point charges and a homogeneous background charge $\rho \xaf=q/v$,

is most conveniently evaluated using Ewald summation,

For cubic simulation cells, the value of *ϕ*_{wig}/*q* ≈ −2.837 297/*L* is known to high precision.^{94} The finite size corrections described in this paper generally require *ϕ*_{wig} for anisotropic cells, which we computed from Eq. (B10) (with a large value of *κL*_{x} that justifies neglecting the first sum entirely). For highly anisotropic cells, the magnitude of *ϕ*_{wig} becomes very large (see Fig. 10), appearing to diverge linearly with growing aspect ratio, $e\varphi wig/q=Lx\u22121(1.0472Lz/Lx\u22123.8845)$. This divergence precisely matches that of Δ*ϕ** in Eq. (B8) so that the generalized Wigner potential $\varphi wig*$ remains finite even for very large aspect ratios.

## REFERENCES

Acknowledging that the dielectric does not penetrate the solute is essential for local polarization response, which would otherwise be singular. This local response, however, should be consistent across different system sizes; a comparison among them eliminates the singularity and justifies treating the solute as a point charge within the dielectric domain Ω.

We take the net surface charge on the slab’s two boundaries to be equal. For a single aperiodic slab, this symmetry holds regardless of the ion’s position in the slab, as can be shown by summing image charges above and below each interface. For the periodic case, we assume that symmetrically placed replicas do not upset this balance (Figs. 7 and 9).

We have in mind that the solute is sufficiently far from the interface that capillary fluctuations are not disrupted.