Assessing long-range contributions to the charge asymmetry of ion adsorption at the air-water interface
Stephen J. Cox, Dayton G. Thorpe, Patrick R. Shaffer, Phillip L. Geissler
AAssessing long-range contributions to the charge asymmetry of ionadsorption at the air-water interface
Stephen J. Cox, Dayton G. Thorpe,
2, 3
Patrick R. Shaffer, and Phillip L. Geissler
2, 4 Department of Chemistry, University of Cambridge, Lensfield Road, Cambridge CB2 1EW,United Kingdom Chemical Sciences Division, Lawrence Berkeley National Laboratory, Berkeley, CA 94720,United States. Department of Physics, University of California, Berkeley, CA 94720, United States. Department of Chemistry, University of California, Berkeley, CA 94720, United States. (Dated: 2 October 2020)
Anions generally associate more favorably with the air-water interface than cations. In addition to solute sizeand polarizability, the intrinsic structure of the unperturbed interface has been discussed as an importantcontributor to this bias. Here we assess quantitatively the role that intrinsic charge asymmetry of water’ssurface plays in ion adsorption, using computer simulations to compare model solutes of various size andcharge. In doing so, we also evaluate the degree to which linear response theory for solvent polarization isa reasonable approach for comparing the thermodynamics of bulk and interfacial ion solvation. Consistentwith previous works on bulk ion solvation, we find that the average electrostatic potential at the centerof a neutral, sub-nanometer solute at the air-water interface depends sensitively on its radius, and thatthis potential changes quite nonlinearly as the solute’s charge is introduced. The nonlinear response closelyresembles that of the bulk. As a result, the net nonlinearity of ion adsorption is weaker than in bulk, butstill substantial, comparable to the apparent magnitude of macroscopically nonlocal contributions from theundisturbed interface. For the simple-point-charge model of water we study, these results argue distinctlyagainst rationalizing ion adsorption in terms of surface potentials inherent to molecular structure of theliquid’s boundary.Counter to expectations from conventional theories ofsolvation, there is a large body of both computationaland experimental evidence indicating that small ions canadsorb to the air-water interface.
Implications acrossthe biological, atmospheric and physical sciences haveinspired efforts to understand the microscopic drivingforces for ions associating with hydrophobic interfaces ingeneral.
A particular emphasis has been placed onunderstanding ion specificity, i.e., why some ions exhibitstrong interfacial affinity while others do not. Empir-ical trends indicate that ion size and polarizability areimportant factors, as could be anticipated from conven-tional theory. More surprisingly, the sign of a solute’scharge can effect a significant bias, with anions tendingto adsorb more favorably than cations.Here we examine the microscopic origin of this chargeasymmetry in interfacial ion adsorption. We specificallyassess whether the thermodynamic preference can be sim-ply and generally understood in terms of long-range bi-ases that are intrinsic to an aqueous system surroundedby vapor. By “long-range” and “nonlocal” we refer tomacroscopically large scales, i.e., collective forces that arefelt at arbitrarily long distance. Such a macroscopicallylong-range bias is expected from the air-water interfacedue to its average polarization, and by some measures thebias is quite strong. By contrast, “local” contributionscomprise the entire influence of a solute’s microscopic en-vironment, including electrostatic forces from moleculesthat are many solvation shells away – any influence thatdecays over a sub-macroscopic length scale.The importance of macroscopically nonlocal contribu- tions has been discussed extensively in the context of ionsolvation in bulk liquid water, which we review in Sec. Ias a backdrop for interfacial solvation. The notion thatsuch contributions strongly influence charge asymmetryof solvation at the air-water interface has informed the-oretical approaches and inspired criticism of widely usedforce fields for molecular simulation.
A full under-standing of their role in interfacial adsorption, however,is lacking.In the course of this study, we will also evaluate thesuitability of dielectric continuum theory (DCT) to de-scribe the adsorption process. DCT has provided anessential conceptual framework for rationalizing water’sresponse to electrostatic perturbations. But a more pre-cise understanding of its applicability is needed, particu-larly for the construction of more elaborate models (e.g.,with heterogeneous polarizability near interfaces )and for the application of DCT to evermore complex(e.g., nanoconfined ) environments.
I. CHARGE ASYMMETRY IN BULK LIQUID WATER
Our study of interfacial charge asymmetry is stronglyinformed by previous work on the solvation of ions in bulk liquid water. In this section we review importantperspectives and conclusions from that body of work, asa backdrop for new results concerning ions at the air-water interface. a r X i v : . [ c ond - m a t . s o f t ] O c t A. Distant interfaces and the neutral cavity potential
A difference in adsorption behaviors of anions andcations is foreshadowed by the fact that ion solvation inmodels of bulk liquid water is also substantially chargeasymmetric. Born’s classic model for the charging of asolute captures the basic scale of solvation free energies,as well as their rough dependence on a solute’s size. We will characterize the size of a solute by its radius R of volume exclusion, the closest distance that a watermolecule’s oxygen atom can approach without incurringa large energetic penalty. Contrary to Born’s result, com-puter simulations indicate that the sign of the charge ofsmall ions can significantly influence their charging freeenergy F chg ( q, R ) i.e., the work involved in reversibly in-troducing the solute’s charge q . This dependence ismost easily scrutinized for simple point charge (SPC)models of molecular interactions, where an ion’s chargecan be varied independently of its other properties. InSPC/E water, for instance, charging a solute roughlythe size of fluoride ( R F ≈ .
317 nm) has an asymme-try, F chg ( e, R F ) − F chg ( − e, R F ) ≈
16 kcal/mol, almost 30times larger than thermal energy k B T . Here, e is themagnitude of an electron’s charge.The ultimate origin of charge asymmetry in liquid wa-ter is of course the inequivalent distribution of positiveand negative charge in a water molecule itself. On av-erage, the spatial distribution of positive and negativecharge is uniform in the bulk liquid, but any breaking oftranslational symmetry will manifest the distinct statis-tics of their fluctuating arrangements. A neutral, solute-sized cavity in water, for example, experiences an im-mediate environment in which solvent molecules have anonvanishing and spatially varying net orientation. Theinternal charge distributions of these oriented solventmolecules generate a nonzero electric potential at thecenter of the cavity, whose sign and magnitude are notsimple to anticipate. By our characterization, this elec-trostatic bias is local in origin – the total contribution ofmolecules beyond a distance r from the cavity decays tozero as r increases.The inequivalent spatial distribution of positive andnegative charge in water can generate spatially nonlocal biases as well, effects that extend over arbitrarily largedistances. Any point in the bulk liquid is macroscopicallyremoved from the physical boundaries of the liquid phase(e.g., interfaces with a coexisting vapor phase), but thosedistant boundaries may nonetheless impact the thermo-dynamics of bulk ion solvation. This expectation stemsfrom a textbook result of electrostatics: an infinitely ex-tended (or completely enclosing) dipolar surface, withpolarization pointing along the surface normal, generatesa discontinuity in electric potential. This voltage offsetdoes not decay with distance from the interface, and thusmeets our criterion for macroscopic nonlocality. A two-dimensional manifold of polarization density is certainlya crude caricature of a liquid-vapor interface, but for apolar solvent whose orientational symmetry is broken at . . . . . R [nm] − − − − − β e φ n e u t bulkinterface FIG. 1. The average electric potential φ neut at the center ofa neutral cavity varies considerably with the cavity’s radius R . Moreover, this dependence differs for the solute at z = z liq (“bulk”) and z = z int (“interface”). The error bars indicate95% confidence intervals. its boundaries, a similarly long-range potential from theinterface is expected to bias the solvation of charged so-lutes, even macroscopically deep inside the liquid phase.The average electric potential φ neut at the center ofa neutral cavity, which we call the “neutral cavity po-tential”, sums these local and extremely nonlocal con-tributions. The former depends on the cavity’s size(or more generally on the geometry of the solute repre-sented by the cavity). The latter, interfacial contributionshould, by contrast, be insensitive to such microscopicdetails, since the distant surface is unperturbed by thesolute. The net electrostatic bias from these two sourcescan be straightforwardly calculated in computer simu-lations, not only for SPC models but also with ab initio approaches. Fig. 1 shows φ neut ( R ) for cavities inbulk liquid SPC/E water (properly referenced to vaporfollowing Ref. 37). Negative potentials of a few hundredmV, varying by nearly a factor of two as R grows from0.2 nm to 1 nm, echo results of previous studies. Distin-guishing quantitatively between local and nonlocal con-tributions to φ neut , however, is surprisingly confounding,even for the exceedingly strict definition of nonlocalityconsidered here.One strategy to remove local contributions from φ neut is to consider the limit R = 0. In this extreme case theprobe – in effect a neutral, non-volume excluding solute– does not break translational symmetry and induces nostructural response. Given the lack of local structure,the presumably nonlocal quantity φ neut (0) = φ surf is of-ten called the “surface potential”. Lacking volume exclu-sion, however, this probe explores the liquid phase uni-formly, including even the interior of solvent moleculeswhere electrostatic potentials can be very large. A dis-turbing ambiguity results: The value of φ surf can be sen-sitive to modifications of a solvent model that have noimpact on the solvation of any volume-excluding solute.Refs. 39 and 41 illustrate this issue vividly, constructing‘smeared shell’ variants of SPC models with identical sol-vation properties but very different values of φ surf . Thisvariation in surface potential corresponds to differencesin the so-called Bethe potential, which is discussed fur-ther in the Supporting Information (SI).A related, and somewhat more molecular, approachto isolating the electrostatic bias from a distant phaseboundary is to sum contributions to φ neut only frommolecules that reside in the interfacial region. For amacroscopic droplet of liquid water, one could classifyeach molecule in a given configuration as either interfa-cial or bulk based on its position relative to the interface.The restricted sum φ d = (cid:42) N (cid:88) j ∈ interface (cid:88) α q α | r jα | (cid:43) (1)could then be considered as a macroscopically long-ranged, surface-specific component of φ surf that is appro-priately insensitive to a solvent molecule’s internal struc-ture. Here r jα denotes the position of site α in molecule j , whose charge is q α , relative to the center of the droplet. φ d depends significantly, however, on the way moleculesare notionally divided between surface and bulk. This de-pendence, which has been demonstrated previously, we calculate explicitly and generally in the SI. Writtenin the form φ d = − π (cid:90) z vap z liq dz P ( z )where P ( z ) is the solvent dipole density at a dis-placement z from the interface, it reveals φ d asthe well-known “dipole component” of the surfacepotential. Here, z liq and z vap indicatepoints within the bulk liquid and bulk vapor, respectively.For SPC/E water, a surface/bulk classification in Eq. 1based on the position of a water molecule’s center ofcharge gives a value φ centerd = −
40 mV that differs froman oxygen atom-based classification, φ Od = 240 mV, evenin sign. Because water molecules are not point parti-cles, there is no unique way to define an interfacial pop-ulation, and as a result no unique value of φ d , thoughattempts have been made to define an optimal choice. And because molecules near the liquid’s boundary arenot strongly oriented on average, the range of plausiblevalues for φ d is as large as their mean.The ambiguities plaguing interpretations of φ surf and φ d are one and the same. Indeed, if we consider an in-terfacial population of charged sites rather than intactmolecules, then φ surf and φ d become equal. (When defin-ing an interface of intact molecules, φ surf and φ d differby the so-called Bethe potential, whose analogous am-biguity is described in SI.) φ neut has been characterizedas a two-interface quantity, combining thebias φ d from the distant solvent-vapor interface togetherwith the remaining “cavity” bias φ c = φ neut − φ d fromthe local solute-solvent interface. From the perspectivewe have described, these two interfaces are not truly sep-arable, even if a macroscopic amount of isotropic bulk liquid intervenes between them – they must be definedconsistently, and the manner of definition substantiallyinfluences the change in electrostatic potential at eachinterface. This is not to say that such a decompositioncannot be useful. Indeed, for computationally demand-ing ab initio approaches it can be convenient to considerlocal and nonlocal contributions to φ neut such that, ina first step, φ c can be obtained from relatively smallsimulations of the bulk under periodic boundary condi-tions. The effects of φ d can then be accounted for ina subsequent step involving simulations of the neat air-water interface. Such an approach was used to good effectin Ref. 30 to calculate the solvation free energy of LiF.Nonetheless, this still amounts to an arbitrary choice ofdividing surface, making it challenging to assign aphysical interpretation to φ d and φ c individually. Differ-ent, and equally plausible, ways of partitioning moleculescan give different impressions of the two interfaces. Onlythe sum φ neut = φ c + φ d is unambiguous.Establishing an absolute electrostatic bias on the bulkliquid environment due to a distant interface is thushighly problematic for water. A direct scrutiny of thisnonlocal contribution, based on the fundamentally am-biguous potential φ d , is untenable. Instead, we assessthe relative importance of local and nonlocal biases bycomparing the solvation properties of different ions. Lo-cal contributions can depend sensitively on features likesolute size R and charge q , while macroscopically non-local contributions cannot. Long-range influence of theinterface might therefore be clarified by dependence ofthe neutral cavity potential on R . In particular, dom-inance by the distant liquid-vapor interface would im-ply weak variation of φ neut with solute size, which in-fluences only microscopically local structure. The solutesize-dependence shown in Fig. 1 does not support sucha dominance. Growing the cavity from R = 0 .
24 nm to0 . φ neut by roughly 100 mV, followed by anincreasing trend for larger cavities. As emphasized inRefs. 30 and 41, the role of local charge asymmetry is farfrom negligible over this range of solute size.It is tempting to expect the large- R behavior of φ neut to reveal a strictly interfacial component, since localforces attenuate in magnitude when solvent moleculescannot approach the probe position closely. As oth-ers have noted, however, neutral cavities larger than R = 1 nm induce a solvent environment with the basiccharacter of the air-water interface. In the limit of large R , drying at the solute-solvent interface will generate acavity potential that cancels the oppositely oriented dis-tant interface with the vapor phase, yielding φ neut ≈ This asymptotic cancellation should begin for nanoscalecavities, though effects of local interface curvature maycause φ neut to decay slowly towards zero. Judging fromour results, there is no intermediate plateau value of φ neut that could reasonably be assigned to a single liquid-vaporinterface. B. Solvation thermodynamics and the asymmetrypotential
The difficulty of uniquely identifying a surfacedipole component of φ neut notwithstanding, the rel-evance of such neutral probe quantities for ionsolvation thermodynamics has also been thoroughlyexamined. As an essentialthermodynamic measure of solvation, we examine thefree energy change F solv ( q, R ) when a solute ion is re-moved from dilute vapor and added to the liquid phase.This change could be evaluated along any reversible paththat transfers the solute between phases, and differentpaths can highlight different aspects of solvent response.For studying charge asymmetry, a particularly appeal-ing path first creates a neutral, solute-sized cavity in theliquid, with reversible work F cav ( R ). The second step,whose free energy change F chg ( q, R ) was discussed above,introduces the solute’s charge. The charge asymmetryof interest compares solvating a cation and anion of thesame size; since F cav is insensitive to the solute’s charge,its contribution to F solv = F cav + F chg cancels in the dif-ference F solv ( q, R ) − F solv ( − q, R ) = F chg ( q, R ) − F chg ( − q, R )(2) ≡ qψ ( q, R ) (3)Eq. 3 defines an asymmetry potential ψ , an analogue of φ neut that accounts for solvent response.The connection between ψ ( q, R ) and φ neut can be madeprecise through a cumulant expansion of F chg in powersof q , F chg ( q, R ) = q (cid:104) φ solv (cid:105) − βq (cid:104) ( δφ solv ) (cid:105) + O ( q ) , (4)where (cid:104)· · · (cid:105) denotes a canonical average in the pres-ence of a neutral solute-sized cavity, φ solv is the fluc-tuating electric potential at the center of the cavity dueto the surrounding solvent (so that φ neut = (cid:104) φ solv (cid:105) ),and δφ solv = φ solv − φ neut . The O ( q ) term in Eq. 4describes linear response of the solvent potential φ solv to the solute’s charging. This response, which could becaptured by a Gaussian field theory `a la DCT, is chargesymmetric by construction. The asymmetry potential ψ ( q, R ) = φ neut ( R ) + O ( q ) is therefore equivalent to φ neut within linear response.Previous work has demonstrated that water’s re-sponse to charging sub-nanometer cavities is significantlynonlinear. In ψ ( q, R ) the breakdownof linear dielectric behavior is evidenced by deviationsaway from the limiting value ψ (0 , R ) = φ neut ( R ). Fig. 2ashows our numerical results for the asymmetry potentialas a function of q for solutes in bulk liquid SPC/E wa-ter. For large solutes ( R > ∼ . ψ ismodest as q increases from 0 to e . For smaller cavities,linear response theory fails dramatically, in that chargeasymmetry changes many-fold as the solute is charged.In the case of a fluoride-sized solute, the asymmetry at full charge ( eψ ( e, R F ) ≈ k B T ) is qualitatively differ-ent than in linear response ( eφ neut ( R F ) ≈ − k B T ). ForSPC models of bulk liquid water, the ultimate electro-static bias in solvating cations and anions of this sizeclearly cannot be attributed to the innate environmentof a neutral cavity, much less to the structure of a dis-tant interface. Ab initio molecular dynamics studies havereached a similar conclusion. SPC simulations of bulk liquid water indicate that thenonlinearity of solvent response to solute charging has astep-like character:
For one range of solute charge( q < q c ), the susceptibility d (cid:104) φ solv (cid:105) q /dq is approximatelyconstant. In the remaining range ( q ≥ q c ), d (cid:104) φ solv (cid:105) q /dq is also nearly constant, but with a different value. Piece-wise linear response (PLR) models inspired by this ob-servation give a broadly reasonable description of bulksolvation thermodynamics throughout the entire range − e < q < + e . In our discussion of ion adsorption be-low, we will assess the suitability of a PLR model forinterfacial solvation as well. II. CHARGE ASYMMETRY IN ION ADSORPTION
In bulk liquid water, an electric potential from itsbounding interfaces cannot be unambiguously identified.Even the sign of the bias generated by a liquid-vapor in-terface is unclear. Moreover, the nonlinear local responseto solute charging can exert a bias on ion solvation thatsignificantly outweighs the charge asymmetry due to dis-tant interfaces.Solvation within the interfacial environment is hardlyless complex, juxtaposing the fluctuating intermoleculararrangements of bulk water together with broken symme-try and the microscopic shape variations of a soft bound-ary. It is thus unlikely that complications described inSec. I for bulk liquid are much eased in the interfacialscenario. We should not expect, for example, that theneutral cavity potential for a solute positioned near theinterface will be dominated by a simple nonlocal contri-bution. Nor should we expect the accuracy of linear re-sponse approximations to be greatly improved, such that φ neut is predictive of charge asymmetric solvation.The adsorption of an ion to the interface, however,concerns the difference in solvation properties of bulkand interfacial environments. To the extent that non-linear response and local structuring at the interface aresimilar to those in bulk liquid, their effects may cancel,or at least significantly offset, in the thermodynamics ofadsorption. Our main results concern this possibility ofcancellation, which would justify regarding macroscopi-cally nonlocal contributions to φ neut as the basic originof charge asymmetry in ion adsorption.We begin by establishing that biases on solvation atthe interface are complicated in ways that qualitativelyresemble biases in bulk. As before, we consider soluteswith a range of sizes and charges, now positioned at theliquid’s boundary (illustrated in Fig. 3a). The free ener- (a) bulkinterface (b) FIG. 2. Ion solvation in water is both asymmetric and non-linear, as quantified by the asymmetry potential ψ ( q, R ; z ).Results are shown for solutes (a) in the bulk liquid, and(b) near the air-water interface, spanning ranges of charge0 < q ≤ e and solute size 0 . ≤ R ≤ . ψ < q , in-dicating that weakly charged cations are more favorably sol-vated than anions. For the smaller solutes, ψ increases with q , a signature of non-linear response. Anions consequentlybecome more favorably solvated at large q . For the larger so-lutes ( R = 0 .
75 nm and R = 1 . ψ on q . gies and potentials defined in Sec. I for bulk solution nowacquire dependence on the Cartesian coordinate z thatpoints perpendicular to the mean surface. SI shows thedetailed location z int we designate as adsorbed for eachion. In all cases z int lies near the Gibbs dividing surface,where the solvent density falls to half its bulk value. Thelarger solutes occupy considerable volume, so that thesolvent density profile in our finite simulation cell changesnoticeably with their height z . A precise interfacial solute location is therefore difficult to justify. When neutral andlocated near z int , however, these nanometer-size solutestend to deform the instantaneous phase boundary, just as they induce local drying in bulk solution. Thisresponse essentially fixes their location relative to the in-stantaneous interface, so that their solvation propertiesshould be fairly insensitive to the choice of z int .The neutral cavity potential for interfacial solutes isshown in Fig. 1. As was observed for the bulk liquid, φ neut is consistently negative over the range R = 0 .
24 nmto R = 1 nm but varies significantly with solute size.In this case the potential increases nearly monotoni-cally with R , though the values of φ neut (0 .
75 nm) and φ neut (1 nm) are statistically indistinguishable within oursampling. Just as for bulk liquid, we expect φ neut to van-ish in the limit R → ∞ . Here, drying at the surface ofvery large solutes effects a distortion of the liquid-vaporinterface that places the probe (located at the cavity’scenter) distinctly in the vapor phase. Judging from ourresults, the asymptotic approach to this limit is quiteslow for interfacial solutes. Nonetheless, φ neut changes bynearly 40% over the range of R considered, emphasizingthe importance of local, solute-dependent contributions.As concluded for the bulk solvent, macroscopically nonlo-cal potentials arising from orientational structure of theair-water interface do not dominate the charge asymme-try experienced by neutral solutes at z int .The response to charging a solute at the air-water in-terface is strongly nonlinear, to a degree comparable withbulk response. A similarly important role of nonlinearresponse at interfaces has been reported previously. The resulting q -dependent charge asymmetry closely re-sembles bulk behavior, as quantified by the asymmetrypotential ψ ( q, R ; z ), whose dependence on solute positionwe now make explicit. Fig. 2b shows simulation resultsfor ψ ( q, R ; z int ) for SPC/E water. On the scale that ψ changes as q increases from 0 to e , the charging responsein bulk liquid and at the interface are nearly indistin-guishable by eye. This close similarity suggests that thepredominant source of nonlinearity lies in aspects of localresponse which are not so different in the two environ-ments.Comparing ψ ( q, R ; z int ) with ψ ( q, R ; z liq ), and φ neut ( R ; z int ) with φ neut ( R ; z liq ), gives a sense forfeatures of solvation that most strongly shape ionadsorption. Similarities point to aspects of solventstructure and response which are largely unchangedwhen an ion moves to the interface. These contributionsmay be important for solvation in an absolute sense,but their cancellation indicates a weak net influence onadsorption thermodynamics.For all values of R we considered, φ neut is less negativeat z int than at z liq . In the simplest conception of the liq-uid’s boundary as a layer of nonzero dipole density, onewould expect the nonlocal component of φ neut to atten-uate steadily in magnitude as a solute moves from theliquid phase into the interfacial region, and then vanishas the solute enters vapor. Whether this rough picture isconsistent with the observed shift in φ neut depends on thesign of the nonlocal potential φ d . Unfortunately this signis uncertain, as described in Sec. I, due to the intrinsicambiguity in dividing molecules between bulk and surfaceregions. Ref. 41 calculated a positive dipole component ofthe surface potential, φ d = +260 mV. Within the simplecontinuum picture, this value suggests a downward shiftin φ neut as z increases from z liq to z int , in contrast toour simulation results. A different partitioning scheme,however, can give φ d <
0, suggesting an upward shift, aswe observe in simulation.Although the direction of change in φ neut might be an-ticipated from the sign of φ d , the magnitude of this shiftvaries considerably with solute size. For R = 0 .
24 nm, | φ neut | is reduced by about 15% when the cavity is placedat the interface. For R = 0 .
415 nm the reduction isgreater than 50%. This variation cannot arise from non-local biases, which are insensitive to the size or chargeof a solute. A distinct, macroscopically nonlocal contri-bution could manifest as a nonzero asymptotic value of∆ ads φ neut = φ neut ( R ; z int ) − φ neut ( R ; z liq ) at intermediate R ; according to our data, if such a limit exists it occursfor solutes larger than 1 nm.The similarity between the asymmetry potentials ψ ( q, R ) for solutes in the bulk and at the interface of-fers some hope that complicating factors of nonlinear re-sponse cancel out in the adsorption process. The extentof this cancellation is quantified by an adsorption asym-metry potential∆ ads ψ ( q, R ) = ψ ( q, R ; z int ) − ψ ( q, R ; z liq ) , (5)= (2 βq ) − ln (cid:20) ρ int ( − q, R ; ρ bulk ) ρ int (+ q, R ; ρ bulk ) (cid:21) , (6)where ρ int is the average number density of a solute at z = z int , given its concentration ρ bulk in bulk solution. Eq. 6highlights the direct relationship between ∆ ads ψ ( q, R )and the relative adsorption propensities of cations andanions: For dilute solutes with opposite charge, equalsize, and equal bulk concentration, exp [2 βq ∆ ads ψ ( q, R )]directly indicates the enhancement of anions over cationsat the interface, as shown in Figs. 3 and 4. From the pre-ceding discussion of the asymmetry potential itself, it isclear that ∆ ads ψ ( q → , R ) = ∆ ads φ neut . The full de-pendence of ∆ ads ψ on q thus incorporates the adsorptionbehavior of the neutral cavity potential as well as thecorresponding solvent response to charging. Our numer-ical results for ∆ ads ψ ( q, R ) are the central contributionof this paper.The adsorption asymmetry potential ∆ ads ψ ( q, R ), asdetermined from simulations of the SPC/E model, areplotted as a function of q in Fig. 4 for several values of R .For the smaller solutes, the scale on which ∆ ads ψ changesupon charging is dramatically smaller than the asymme-try potentials themselves. Nonlinear solvent response inthese cases cancels substantially in the process of adsorp-tion, but by no means completely. Despite the partialcancellation, ∆ ads ψ still varies by more than 100 mV as (a)(b) FIG. 3. The propensity for an ion to adsorb to the air-waterinterface depends strongly on the sign of its charge. (a) Snap-shot of an iodide-sized anion ( R = 0 .
415 nm) at the interface.The system comprises a free-standing slab of liquid water sur-rounded on either side by its vapor. (Only one of the twointerfaces is shown.) The z direction is indicated by the ar-row. The size of the solute is depicted schematically by thedashed circle. (b) Potential of mean force ∆ F as a functionof ion position z , for a solute charge q = +0 . e (solid orange)and q = − . e (dashed blue). The anion adsorbs much morestrongly to the interface than the cation for this solute size.The dotted green line indicates the connection between thesefree energy profiles and the adsorption asymmetry potentialin Eq. 6. q increases from 0 to e , comparable in magnitude to φ d and φ neut . For R = 0 .
24 nm and R = 0 .
317 nm, this vari-ation is sufficient to change even the sign of ∆ ads ψ , andtherefore to change the sense of charge bias: Small mono-valent cations “adsorb” more favorably to the air-waterinterface than do anions of the same size. In this sizerange, however, the adsorbed state is unstable relativeto the fully solvated ion in bulk solvent unless q is verysmall in magnitude.As was previously observed for bulk solvation, we findthat the response to charging a solute at the air-waterinterface, while nonlinear on the whole, is roughly piece-wise linear. Deviations from piecewise linearity are gen-erally stronger in the interfacial case. It is therefore lessstraightforward to parameterize an interfacial piecewiselinear response model, i.e., to identify a crossover charge q c at which the susceptibility d (cid:104) φ solv (cid:105) q /dq changes dis-continuously. The SI presents plausible choices for q c and these limiting susceptibilities for our three small-est solutes, from which adsorption asymmetry poten-tials ∆ ads ψ (PLR) can be readily computed. The resultingPLR predictions are plotted in Fig. 4b. Two basic fea-tures of our simulation results are accurately capturedby this phenomenological description. Specifically, (i)for small solute charge, ∆ ads ψ is an approximately con-stant or modestly increasing function of q , and (ii) a morestrongly decreasing trend of ∆ ads ψ follows for larger q .Nearly quantitative agreement can be obtained for aniodide-sized solute, R = 0 .
415 nm. Smaller solutes ex-hibit a more complicated charge dependence that lies be-yond a simple PLR description. We note that this test ofPLR is a demanding one, given the small scale of ∆ ads ψ relative to ψ ( q, R ; z int ) and ψ ( q, R ; z liq ) individually. Tothe extent that PLR is a successful caricature, these re-sults suggest that the adsorption charge asymmetry atfull charging ( q = e ) derives from a combination of fea-tures of solvent response, including an interface-inducedshift in the crossover charge q c at which the characterof linear response changes. The neutral cavity potential φ neut figures into this combination as well, but by nomeans does it dominate for these solute sizes.For the larger solutes we examined, the nonlinearityof solvent response to charging is not pronounced, ei-ther in bulk liquid or at the interface. The differencein nonlinearity of these environments is necessarily alsonot large, with e ∆ ads ψ changing by less than k B T overthe range q = 0 to q = e . This small variation is com-parable in scale to those of ψ ( q, R ; z int ) and ψ ( q, R ; z liq )themselves. Judged on that scale, the cancellation of non-linear response is in fact less complete for R = 0 .
75 nmthan for smaller solutes. As we have discussed, cavitieswith
R > ∼ z coincides with the Gibbs dividing surface. Whensuch a solute is endowed with sufficient charge, wettingwill occur at its surface, eventually raising the interfaceto effectively move the solute into the liquid phase. Thissolvent response, which originates in the physics of phaseseparation, is intrinsically nonlinear. For large R , a so-lute charge well in excess of e is required to fully inducethis structural change, but at the nanoscale it may man-ifest as an incipient nonlinearity for q ≈ e .In summary, the adsorption asymmetry potential∆ ads ψ depends significantly on solute size R and charge q . Neither of these sensitivies can arise from intrinsicorientational bias at the neat air-water interface. Long-range electrostatic forces from oriented molecules at theliquid’s boundary, which contribute importantly to sur-face potentials like φ surf and φ d , are inherently unaffectedby the presence, size, or charge of a sufficiently distant solute. These results highlight the importance of localsolvent structure and response for charge asymmetry ininterfacial ion adsorption, and they highlight the dan-ger of inferring solvation thermodynamics from ion-freequantities such as φ surf and φ d . . n m . n m . n m (a)(b) FIG. 4. Linear response theory cannot faithfully describe thedifferences between adsorption profiles of sub-nanometer an-ions and cations, as demonstrated in (a) by variations in ad-sorption asymmetry potential ∆ ads ψ with both R and q . Forthe smallest solutes ( R < ∼ . ads ψ even changes signas q increases. In this size range, fully charged cations aremore abundant at the interface than anions (with the samebulk concentration). At larger R , solutes with q = − e adsorbmore strongly than those with q = + e . As the solute diameterapproaches R = 1 nm, nonlinear response during the chargingprocess becomes much less pronounced. Values of R are indi-cated in the legend. (b) A PLR model (heavy lines) predicts∆ ads ψ is initially flat, followed by a steady decrease as q in-creases. This qualitatively captures the simulation data (lightlines), although it fails to capture the leveling off at large q seen for R = 0 .
240 nm and 0 .
317 nm.
III. DISCUSSION AND CONCLUSIONS
In this study, we set out to understand whether ornot charge asymmetry in interfacial ion adsorption couldbe understood in terms of macroscopically long-ranged,collective forces intrinsic to water. For ion solvationin bulk, difficulties in unambiguously determining suchlong-ranged contributions were already apparent fromprevious results. Building on that work, our results showthat for SPC models of water such a simple mechanisticpicture is inadequate for interfacial solvation as well. Inaddition to the difficulties in partitioning molecules be-tween ‘near’ and ‘distant’ interfaces, complex nonlinearresponse also underlies substantial shortcomings of tryingto rationalize ion adsorption from surface potentials thatcharacterize biases of the undisturbed air-water interface.The nonlinearities in F chg ( q, R ) for bulk and interfacialenvironments, while similar, are sufficiently different thatthe process of adsorption is also substantially nonlinear.A compelling inference of adsorption tendencies from in-trinsic properties of the undisturbed liquid and its inter-face with vapor requires information that is more subtlethan an average electric potential and macroscopic di-electric susceptibility. As highlighted by the potentialdistribution theorem, this information can in principlebe gleaned from equilibrium statistics of the undisturbedsolvent. But in terms of fluctuations in electric potential,it involves high-order correlations whose physical mean-ing is not transparent.In previous work we developed and tested finite sizecorrections for computer simulations of interfacial ionsolvation. Based on DCT, these corrections proved tobe quite accurate even for simulation unit cells withnanometer dimensions. Our conclusion that DCT is afaithful representation of aqueous polarization responsedown to nanometer length scales is reinforced by the re-sults of this paper. In particular, when charging a so-lute of diameter R = 1 nm, solvent response on an abso-lute scale is linear to a very good approximation, bothin bulk liquid and at the interface. The results of Figs. 2and 4, however, also indicate that 1 nm marks the valid-ity limit of linear response. When charging a cavity with R = 0 .
75 nm, nonlinear contributions to charge asym-metry are quantitatively important; for smaller solutessuch nonlinear contributions become not just importantbut instead dominant. In passing, we note that even forthe larger solutes, a significant charge asymmetry per-sists, both for bulk solvation and for adsorption to theinterface. This persistent bias weighs against the basis ofthe tetra-phenyl arsonium/tetra-phenyl borate (‘TATB’)extrathermodynamic assumption, an issue that has alsobeen raised by others.
The highly simplified description of molecular interac-tions in SPC models is certainly a crude approximationto real microscopic forces. But the specific ion effectsit exhibits cannot be ascribed simply to an errant sur-face potential. Indeed, discrepancies between models inpotentials such as φ d (whose definition requires an arbi- trary convention), φ surf (which pertains to a solute thatdoes not exclude volume), or even φ neut (which for sub-nanometer solutes does not account for the strong asym-metry of solvent response) are not greatly alarming. φ surf and φ d can vary significantly among different models, butthey do not weigh on ion solvation thermodynamics in adirect way, either in bulk liquid or at the air-water inter-face. (This does not contradict their use for computing F chg once a choice for partitioning molecules between theinterface and bulk has been made.) By contrast, trendsin F solv and ∆ ads ψ at full charging reflect on essential mi-croscopic mechanisms that underlie specific ion adsorp-tion. SPC models may be best viewed as caricatures ofa disordered tetrahedral network, with intrinsic chargeasymmetry due to the distinct geometric requirementsof donating and accepting hydrogen bonds. These es-sential features of liquid water are often associated withnonlinear response in solvation. By implicating non-linearities of precisely this kind as sources of ion-specificadsorption properties, our results support the use of SPCmodels as a physically motivated test bed for exploringthe microscopic basis of surprising trends in interfacialsolvation. Conversely, our results underscore the limita-tions of DCT and notions of long-ranged contributionsfrom unperturbed interfaces, which do not describe es-sential local aspects of the chemical physics underlyingion adsorption and its charge asymmetry. The conse-quences of this shortcoming are likely to be exacerbatedin confined geometries. Work to move beyond standardDCT approaches is an active area of research (e.g. Refs.24–26,69) and it is hoped that the results presented inthis study will help to guide future theoretical develop-ments.
IV. METHODS
All simulations used the SPC/E water model; so-lutes were represented as Lennard-Jones (LJ) sphereswith a central charge q . The SHAKE algorithm was usedto maintain a rigid water geometry. Periodic bound-ary conditions were imposed in all three Cartesian di-rections, with the liquid phase spanning two directionsin a slab geometry. Long-range Coulomb interactionswere summed using the particle-particle particle-meshEwald method.
A spatially homogeneous backgroundcharge was included to maintain electroneutrality andthus guarantee finite electrostatic energy. For solute sizes
R < .
75 nm, the system comprised 266 water moleculeswith simulation cell dimensions 2 × × . . For R ≥ .
75 nm the simulation cell size was 3 . × . × . and we used 1429 water molecules. Solvent density pro-files that indicate the interfacial location z int for eachsolute are given in the SI. A time step of 1 fs was usedfor all simulations. A temperature of 298 K was main-tained using Langevin dynamics, as implementedin the LAMMPS simulation package, which was usedthroughout.Due to the long range of Coulomb interactions, ionsolvation in polar solvents has important contributionseven from distant solvent molecules. Thermodynamicestimates from molecular simulations are thus subject tosubstantial finite size effects, which have been the focusof many studies. In Ref. 37 we showed for liquidwater in a periodic slab geometry that values of φ neut depend on simulation box size in a slowly decaying butpredictable way. The limit of infinitely separated peri-odic images can thus be obtained with a simple finitesize correction, which amounts to referencing electric po-tential values to the vapor phase. We have applied thiscorrection to all potentials reported in this paper. Thepotential of mean force ∆ F ( q, R ; z ) for ions in periodicliquid slabs are, by contrast, nearly independent of sim-ulation cell size for z ≤ z int . To compute ∆ F ( q, R ; z ), we followed the same proce-dure as outlined in Ref. 18, namely umbrella samplingand histogram reweighting with MBAR. To calculate ψ ( q, R ; z ) for a given choice of R and z , simulations wereperformed with q/e = − . , − . , . . . , +0 . , and + 1 . q values allows for ample overlap amongprobability distributions P q (cid:0) φ solv (cid:1) of the electrostaticpotential at the center of the solute (Fig. S10). For R ≤ .
415 nm statistics were obtained from trajectories5 ns in duration. For R = 0 .
75 nm and 1 . . . P ( φ solv ) over a correspondingly broad range of φ solv . F chg was computed by averaging exp( − βqφ solv )according to the distribution P ( φ solv ), as prescribed byWidom’s potential distribution theorem, e − βF chg = (cid:90) d φ solv P (cid:0) φ solv (cid:1) e − βqφ solv (7)The integral in Eq. 7 was performed numerically. ACKNOWLEDGMENTS
S.J.C (02/15 to 09/17), D.G.T, P.R.S and P.L.G weresupported by the U.S. Department of Energy, Office ofBasic Energy Sciences, through the Chemical SciencesDivision (CSD) of Lawrence Berkeley National Labo-ratory (LBNL), under Contract DE-AC02-05CH11231.Since 10/17, S.J.C. has been supported by a Royal Com-mission for the Exhibition of 1851 Research Fellowship. P. Jungwirth and D. J. Tobias, Chem. Rev. , 1259 (2006). R. R. Netz and D. Horinek, Annu. Rev. Phys. Chem. , 401(2012). D. E. Otten, P. R. Shaffer, P. L. Geissler, and R. J. Saykally,Proc. Natl. Acad. Sci. USA , 701 (2012). M. Mucha, T. Frigato, L. M. Levering, H. C. Allen, D. J. Tobias,L. X. Dang, and P. Jungwirth, J. Phys. Chem. B , 7617(2005). P. B. Petersen, R. J. Saykally, M. Mucha, and P. Jungwirth, J.Phys. Chem. B , 10915 (2005). L. Piatkowski, Z. Zhang, E. H. Backus, H. J. Bakker, andM. Bonn, Nature Commun. , 4083 (2014). D. Verreault, W. Hua, and H. C. Allen, J. Phys. Chem. Lett. ,3012 (2012). D. Liu, G. Ma, L. M. Levering, and H. C. Allen, J. Phys. Chem.B , 2252 (2004). M. D. Baer and C. J. Mundy, J. Phys. Chem. Lett. , 1088 (2011). D. Ben-Amotz, J. Phys.: Condens. Matter , 414013 (2016). J. Noah-Vanhoucke and P. L. Geissler, Proc. Natl. Acad. Sci.USA , 15125 (2009). A. Arslanargin and T. L. Beck, J. Chem. Phys. , 104503(2012). M. D. Baer, I.-F. W. Kuo, D. J. Tobias, and C. J. Mundy, J.Phys. Chem. B , 8364 (2014). T. L. Beck, Chem. Phys. Lett. , 1 (2013). A. P. dos Santos and Y. Levin, Faraday Discuss. , 75 (2013). Y. Levin, Phys. Rev. Lett. , 147803 (2009). Y. Levin, A. P. dos Santos, and A. Diehl, Phys. Rev. Lett. ,257802 (2009). D. L. McCaffrey, S. C. Nguyen, S. J. Cox, H. Weller, A. P.Alivisatos, P. L. Geissler, and R. J. Saykally, Proc. Natl. Acad.Sci. USA , 13369 (2017). S. Ou, Y. Hu, S. Patel, and H. Wan, J. Phys. Chem. B ,11732 (2013). S. Ou and S. Patel, J. Phys. Chem. B , 6512 (2013). C. Caleman, J. S. Hub, P. J. van Maaren, and D. van der Spoel,Proc. Natl. Acad. Sci. USA , 6838 (2011). M. D. Baer, A. C. Stern, Y. Levin, D. J. Tobias, and C. J.Mundy, J. Phys. Chem. Lett. , 1565 (2012). Y. Levin and A. P. dos Santos, J. Phys.: Condens. Matter ,203101 (2014). P. Loche, C. Ayaz, A. Schlaich, D. J. Bonthuis, and R. R. Netz,J. Phys. Chem. Lett. , 6463 (2018). P. Loche, C. Ayaz, A. Wolde-Kidan, A. Schlaich, and R. R. Netz,J. Phys. Chem. B , 4365 (2020). A. Schlaich, E. W. Knapp, and R. R. Netz, Phys. Rev. Lett. , 048001 (2016). L. Fumagalli, A. Esfandiar, R. Fabregas, S. Hu, P. Ares, A. Ja-nardanan, Q. Yang, B. Radha, T. Taniguchi, K. Watanabe,G. Gomila, N. K. S., and G. A. K., Science , 1339 (2018). L. Bocquet, Nature Materials , 254 (2020). W. M. Latimer, K. S. Pitzer, and C. M. Slansky, J. Chem. Phys. , 108 (1939). T. T. Duignan, M. D. Baer, G. K. Schenter, and C. J. Mundy,Chem. Sci. , 6131 (2017). T. T. Duignan, M. D. Baer, G. K. Schenter, and C. J. Mundy,J. Chem. Phys. , 161716 (2017). R. C. Remsing and J. D. Weeks, J. Phys. Chem. B , 6238(2016). J. P. Bardhan, P. Jungwirth, and L. Makowski, J. Chem. Phys. , 124101 (2012). G. Hummer, L. R. Pratt, and A. E. Garc´ıa, J. Phys. Chem. ,1206 (1996). S. Rajamani, T. Ghosh, and S. Garde, J. Chem. Phys. , 4457(2004). R. Lynden-Bell and J. Rasaiah, J. Chem. Phys. , 1981 (1997). S. J. Cox and P. L. Geissler, J. Chem. Phys. , 222823 (2018). H. S. Ashbaugh and D. Asthagiri, J. Chem. Phys. , 204501(2008). R. C. Remsing and J. D. Weeks, J. Stat. Phys. , 743 (2019). H. J. C. Berendsen, J. R. Grigera, and T. P. Straatsma, J. Phys.Chem. , 6269 (1987). R. C. Remsing, M. D. Baer, G. K. Schenter, C. J. Mundy, andJ. D. Weeks, J. Phys. Chem. Lett. , 2767 (2014). H. S. Ashbaugh, J. Phys. Chem. B , 7235 (2000). J. ˚Aqvist and T. Hansson, J. Phys. Chem. B , 3837 (1998). M. A. Kastenholz and P. H. H¨unenberger, J. Chem. Phys ,124106 (2006). P. H¨unenberger and M. Reif,
Single-Ion Solvation , Theoreti-cal and Computational Chemistry Series (The Royal Society of Chemistry, 2011) pp. P001–664. C. C. Doyle, Y. Shi, and T. L. Beck, J. Phys. Chem. B ,3348 (2019). L. Horv´ath, T. Beu, M. Manghi, and J. Palmeri, J. Chem. Phys. , 154702 (2013). We define a molecule’s center of charge according to the chargedsites that specify a particular SPC model. In the case of SPC/Ewater, this center is displaced from the oxygen atom by approx-imately 0.029 nm along the molecular dipole. E. Harder and B. Roux, J. Chem. Phys. , 234706 (2008). D. Chandler, Nature , 640 (2005). While the vapor phase is very dilute at ambient temperature,its nonzero density does yield an average potential different fromthe vacuum environment of a volume-excluding cavity. Here weneglect this small distinction. Y. Shi and T. L. Beck, J. Chem. Phys. , 044504 (2013). T. P. Pollard and T. L. Beck, Curr. Op. Colloid Interf. Sci. ,110 (2016). D. Asthagiri, L. R. Pratt, and H. Ashbaugh, J. Chem. Phys. , 2702 (2003). L. R. Pratt, J. Phys. Chem. , 25 (1992). R. Kubo, J. Phys. Soc. Jpn. , 1100 (1962). F. Hirata, P. Redfern, and R. M. Levy, Int. J. Quantum Chem. , 179 (1988). A. Grossfield, J. Chem. Phys. , 024506 (2005). A. P. Willard and D. Chandler, J. Phys. Chem. B , 1954(2010). S. Vaikuntanathan and P. L. Geissler, Phys. Rev. Lett. ,020603 (2014). S. Vaikuntanathan, G. Rotskoff, A. Hudson, and P. L. Geissler,Proc. Natl. Acad. Sci. USA , E2224 (2016). B. Widom, J. Phys. Chem. , 869 (1982). T. T. Duignan, M. D. Baer, and C. J. Mundy, J. Chem. Phys. , 222819 (2018). T. P. Pollard and T. L. Beck, J. Chem. Phys. , 222830 (2018). R. Scheu, B. M. Rankin, Y. Chen, K. C. Jena, D. Ben-Amotz,and S. Roke, Angew. Chem. Int. Ed. , 9560 (2014). M. Maroncelli and G. R. Fleming, J. Chem. Phys , 5044 (1988). M. S. Skaf and B. M. Ladanyi, J. Phys. Chem. , 18258 (1996). P. L. Geissler and D. Chandler, J. Chem. Phys. , 9759 (2000). C. Zhang and M. Sprik, Phys. Chem. Chem. Phys. , 10676(2020). J.-P. Ryckaert, G. Ciccotti, and H. J. Berendsen, J. Comput.Phys. , 327 (1977). R. W. Hockney and J. W. Eastwood,
Computer simulation usingparticles (CRC Press, 1988). J. Kolafa and J. W. Perram, Mol. Sim. , 351 (1992). B. D¨unweg and W. Paul, Int. J. Mod. Phys. C , 817 (1991). T. Schneider and E. Stoll, Phy. Rev. B , 1302 (1978). S. Plimpton, J. Comput. Phys. , 1 (1995). F. Figueirido, G. S. Del Buono, and R. M. Levy, J. Phys. Chem.B , 5622 (1997). G. Hummer, L. R. Pratt, and A. E. Garc´ıa, J. Chem. Phys. ,9275 (1997). P. H. H¨unenberger and J. A. McCammon, J. Chem. Phys. ,1856 (1999). G. Hummer, L. R. Pratt, and A. E. Garc´ıa, J. Phys. Chem. A , 7885 (1998). M. R. Shirts, D. L. Mobley, J. D. Chodera, and V. S. Pande, J.Phys. Chem. B , 13052 (2007). For a rigid molecule Ω (0) is defined simply by the molecule’sorientation, e.g., a set of Euler angles. S. M. Kathmann, I.-F. W. Kuo, C. J. Mundy, and G. K. Schen-ter, J. Phys. Chem. B , 4369 (2011). M. A. Wilson, A. Pohorille, and L. R. Pratt, J. Chem. Phys. ,5211 (1989). Supporting Information
This document contains a detailed account of the issues faced when trying to isolate contributions to φ neut from localand distant sources. Also given are solvent density profiles ρ ( z ) in the presence of the neutral solute for the differentsystems studied, and the position of the solute at the interface is indicated in each instance. Solute-solvent radialdistribution functions g ( r ) are shown for q = − e , 0 and + e with the solute in the center of the slab. Details underlyingthe piecewise linear response model are also presented. A brief description of how P ( φ solv ) is obtained is given. S1. ELECTROSTATIC CONTRIBUTIONS FROM NEAR AND FAR
The challenge of identifying and interpreting a potential drop across the liquid-vapor interface can be viewed as anissue of partitioning molecules between distinct regions of space.Consider a macroscopic droplet of liquid bounded by an interface S with the vapor phase (as illustrated in Fig. S1).The origin of our coordinate system lies deep within the bulk liquid phase. We will aim to calculate the averageelectric potential (cid:104) φ (cid:105) at the origin, distinguishing contributions of molecules that are far from the probe (includingthose at the phase boundary) from those that lie nearer the origin. Specifically, we will divide the two populations atan imaginary surface B that is also deep within the bulk liquid. We will take B to be distant enough from the originthat liquid structure on this surface is bulk in character, even if the microscopic vicinity of the origin is complicatedby a solute’s excluded volume. probe atorigin liquid-vaporinterfaceimaginary dividingsurface S B L BS FIG. S1. Sketch of a macroscopic droplet of liquid (shaded region at right) surrounded by dilute vapor. The two phases contactat a macroscopically smooth interface S . The surface B within the droplet is a mathematical device to isolate the electrostaticcontribution of molecules residing near the phase boundary S . The droplet’s overall scale L is a macroscopic distance. Amagnified view of a microscopic region straddling B is shown at left. Molecules intersected by B (dashed white line) couldreasonably be assigned to either the near (inside B ) or far (outside B ) domains. A. Partitioning schemes
The vast majority of molecules in the droplet are unambiguously located either outside B (“far”) or inside B (“near”). A tiny fraction straddle the surface B . In the case of water this could involve a molecule’s oxygen atomlying on one side of B , while its hydrogen atoms lie on the other. One division scheme (an M-scheme) would judgethe molecule’s location based on the O atom; another M-scheme might base the classification on the molecule’s centerof charge. A still different scheme (a P-scheme) could divide the molecule in two, with some pieces “near” and otherpieces “far”. (The M-scheme and the P-scheme are well known in the literature. See e.g. Ref. 44.) The total potential φ at the probe site is not sensitive to which of these schemes is chosen. But its contributions φ near and φ far fromatoms/molecules in the near and far regions are sensitive, in an offsetting way.2Let’s first treat the M scheme, with the molecule’s near/far classification based on the position r (0) of some sitewithin the molecule (say, its O atom). The average far-field potential in this case is (cid:104) φ Mfar (cid:105) = N (cid:90) outside B d r (cid:90) d Ω p ( r , Ω ) (cid:88) α q α | r + ∆ r α ( Ω ) | , (S1)where N is the total number of molecules in the droplet and α indexes charged sites within each molecule. Here, p ( r , Ω ) = (cid:104) δ ( r − r (0) ) δ ( Ω − Ω (0) ) (cid:105) is the joint probability distribution of a molecule’s position (i.e., r (0) ) and in-tramolecular configuration Ω (0) (specified relative to the reference position r (0) , as indicated by the superscript). By ∆ r α = r α − r (0) we denote the displacement of charge q α from the reference point r (0) . This intramoleculardisplacement is entirely determined by Ω (0) .For the P-scheme, each charge α contributes to φ Pfar if r α lies outside B . The corresponding far-field potential is (cid:104) φ Pfar (cid:105) = N (cid:88) α (cid:90) outside B d r (cid:90) d Ω p α ( r , Ω ) q α | r | (S2)= N (cid:88) α (cid:90) outside B d r (cid:90) d Ω p ( r − ∆ r α , Ω ) q α | r | (S3)where p α ( r , Ω ) is the joint probability distribution for site position r α and intramolecular configuration of a solventmolecule. In Eq. S3 we have made use of the connection p α ( r , Ω ) = (cid:104) δ ( r − r α ) δ ( Ω − Ω (0) ) (cid:105) = (cid:104) δ ( r − ∆ r α − r (0) ) δ ( Ω − Ω (0) ) (cid:105) (S4)= p ( r − ∆ r α ( Ω ) , Ω ) (S5)between the distributions p and p α . B. Multipole expansion
Since the entire “far” region is macroscopically distant from the origin, small-∆ r α expansions of | r + ∆ r α | − and p ( r − ∆ r α , Ω ) are well justified. These yield (cid:88) α q α | r + ∆ r α | = (cid:32)(cid:88) α q α ∆ r α (cid:33) · ∇ r + 12 (cid:32)(cid:88) α q α ∆ r α ∆ r α (cid:33) : ∇∇ r + . . . (S6)and (cid:88) α q α p ( r − ∆ r α , Ω ) = −∇ · (cid:88) α q α ∆ r α p ( r , Ω ) + 12 ∇∇ : (cid:88) α q α ∆ r α ∆ r α p ( r , Ω ) + . . . (S7)where we have omitted leading terms proportional to (cid:80) α q α , which vanish by molecular charge neutrality. Whencarried through subsequent calculations, terms beyond quadrupole order in these expansions would vanish due eitherto symmetry or to the macroscopic scale of the droplet.Defining dipole and quadrupole densities as m ( r ) = N (cid:90) d Ω p ( r , Ω ) (cid:88) α q α ∆ r α (S8)and Q ( r ) = N (cid:90) d Ω p ( r , Ω ) (cid:88) α q α ∆ r α ∆ r α (S9)we can write (cid:104) φ Mfar (cid:105) = (cid:90) outside B d r (cid:18) m ( r ) · ∇ r + Q ( r ) : ∇∇ r + . . . (cid:19) (S10)3and (cid:104) φ Pfar (cid:105) = (cid:90) outside B d r r ( −∇ · m ( r ) + ∇∇ : Q ( r ) + . . . ) (S11)Integrating by parts, and noting that m ( r ) and ∇ : Q ( r ) vanish both on B and at infinity, (cid:104) φ Pfar (cid:105) = (cid:90) outside B d r (cid:18) m ( r ) · ∇ r − (cid:18) ∇ r (cid:19) · (cid:0) ∇ · Q ( r ) (cid:1)(cid:19) (S12)= (cid:90) outside B d r (cid:18) m ( r ) · ∇ r − ∇ · (cid:18) ∇ r · Q ( r ) (cid:19) + Q ( r ) : ∇∇ r (cid:19) (S13)Using the divergence theorem, (cid:104) φ Pfar (cid:105) = (cid:104) φ Mfar (cid:105) − (cid:90) B d R ˆ n ( R ) · (cid:18) ∇ R · Q ( R ) (cid:19) (S14)where R is a point on B and ˆ n ( R ) is the corresponding local inward -pointing normal vector. Since B lies within thebulk liquid, where the average quadrupole density Q liq is isotropic, Q ( r ) = I (Tr Q liq /
3) everywhere on this surface.As a result, (cid:104) φ Pfar (cid:105) = (cid:104) φ Mfar (cid:105) + Tr Q liq (cid:90) inside B d r ∇ r (S15)= (cid:104) φ Mfar (cid:105) − π Q liq (S16)These two measures of the far-field potential are thus different. Moreover, the quadrupole trace that determines thisdifference depends on the choice of r (0) . This ambiguity is a well-known feature of the so-called Bethe potential − (4 π/ Q liq . C. Dipole surface potential
To simplify the result for (cid:104) φ Mfar (cid:105) , note that Q ( r ) is isotropic everywhere outside B , except in the microscopic vicinityof S . In the bulk regions of the far domain, we then have Q ( r ) : ∇∇ r − ∝ δ ( r ) = 0. The final term in Eq. S10 thereforehas nonzero contributions only from a thin shell whose volume is proportional to L , where L is the macroscopic scaleof the droplet. Since ∇∇ r − ∼ L − in this shell, the quadrupolar contribution to (cid:104) φ Mfar (cid:105) has a negligible magnitude, L − . As a result, (cid:104) φ Mfar (cid:105) = (cid:90) outside B d r m ( r ) · ∇ r (S17)This integral similarly has nonzero contributions only from a microscopically thin shell of broken symmetry, centeredon the phase boundary S . Since the macroscopic surface is very smooth on this scale, and because the average dipoledensity points normal to the locally planar interface, the far-field potential may be written (cid:104) φ Mfar (cid:105) = (cid:90) S d R (cid:90) dz m ⊥ ( z ) ˆ n ( R ) · ∇ r , (S18)where R is the point on S nearest to r , the coordinate z = ( r − R ) · ˆ n ( R ) is the perpendicular displacement fromthe liquid-vapor interface, ˆ n ( R ) is the outward-pointing normal of S , and m ⊥ ( z )ˆ n ( R ) is the average dipole field at r .Neglecting contributions of O ( z/L ), we may replace r − by R − , and easily evaluate the surface integral, yielding (cid:104) φ Mfar (cid:105) = − π (cid:90) z vap z liq dz m ⊥ ( z ) , (S19)where the integral is performed in the direction from liquid ( z liq <
0) to vapor ( z vap > As they note, its value depends on the reference position r (0) defining the molecular referenceframe. In our calculation this dependence arises from the way molecules are classified relative to the dividing surface B .4 D. Near-field potential
In evaluating (cid:104) φ far (cid:105) , we have made no assumptions about the liquid’s structure near the probe. If the origin liesinside a solute’s excluded volume, then the near-field potential is complicated by the microscopically heterogeneousarrangement of solvent molecules in its vicinity. If, however, the probe is simply a point within the isotropic bulkliquid, then (cid:104) φ near (cid:105) can be easily determined.For a probe that resides in uniform bulk liquid, m ( r ) = 0 and Q ( r ) = Q liq everywhere inside B . In the P-schemewe can conclude immediately from the analogue of Eq. S11 that (cid:104) φ near (cid:105) = 0. In the M-scheme we have (cid:104) φ Mnear (cid:105) = (cid:90) inside B d r Tr Q liq ∇ r = − π Q liq (S20)In either case the total potential sums to (cid:104) φ (cid:105) = (cid:104) φ Mnear (cid:105) + (cid:104) φ Mfar (cid:105) (S21)= (cid:104) φ Pnear (cid:105) + (cid:104) φ Pfar (cid:105) (S22)= − π (cid:90) z vap z liq dz m ⊥ ( z ) − π Q liq (S23)These calculations of local and nonlocal contributions to the mean electrostatic potential resemble previous de-velopments of surface potential in many ways. Ours are somewhat more general than standardcalculations, in that we do not require a specific shape of the liquid domain. (The standard development presumesan idealized geometry of the liquid phase e.g. planar interface or spherical droplet, and integrates the resulting 1-dimensional Poisson equation.) More interestingly, it places the ambiguities surrounding surface potential in an easilyconceived context: The electrostatic bias of an interface is not well defined because there is no unique way to assignmolecules to that interface. Any attempt to do so carries an arbitrariness that (in the case of water) is comparable inmagnitude to the apparent surface potential itself.5
S2. SOLVENT DENSITY PROFILES AND SOLUTE-SOLVENT RADIAL DISTRIBUTION FUNCTIONS (a)(b)
FIG. S2. (a) Average solvent density ρ ( z ), plotted as a function of the coordinate z perpendicular to the liquid-vapor interface,with the solute ( R = 0 .
240 nm) located in the bulk ( z liq = 0 nm, solid blue line), and at the interface ( z int = 1 nm, dashedorange line). The dotted green line is drawn at z = z int . Only half ( z > g ( r ), plotted as a function of the distance r between the solute’s center and the oxygen atom of a watermolecule, with the solute at z = z liq with q = − e , 0, and + e . The vertical dot-dashed gray line is drawn at r = R . (a)(b) FIG. S3. (a) ρ ( z ) with the solute ( R = 0 .
317 nm) located in the bulk ( z liq = 0 nm, solid blue line), and at the interface( z int = 1 nm, dashed orange line). The dotted green line is drawn at z = z int . Only half ( z > g ( r ) with the solute at z = z liq with q = − e , 0, and + e . The vertical dot-dashed gray line is drawn at r = R . (a)(b) FIG. S4. (a) ρ ( z ) with the solute ( R = 0 .
415 nm) located in the bulk ( z liq = 0 nm, solid blue line), and at the interface( z int = 1 nm, dashed orange line). The dotted green line is drawn at z = z int . Only half ( z > g ( r ) with the solute at z = z liq with q = − e , 0, and + e . The vertical dot-dashed gray line is drawn at r = R . (a)(b) FIG. S5. (a) ρ ( z ) with the solute ( R = 0 .
75 nm) located in the bulk ( z liq = 0 nm, solid blue line), and at the interface( z int = 1 .
75 nm, dashed orange line). The dotted green line is drawn at z = z int . Only half ( z > g ( r ) with the solute at z = z liq with q = − e , 0, and + e . The vertical dot-dashed gray line is drawn at r = R . (a)(b) FIG. S6. (a) ρ ( z ) with the solute ( R = 1 . z liq = 0 nm, solid blue line), and at the interface( z int = 1 .
75 nm, dashed orange line). The dotted green line is drawn at z = z int . Only half ( z > g ( r ) with the solute at z = z liq with q = − e , 0, and + e . The vertical dot-dashed gray line is drawn at r = R . S3. EVALUATING PIECEWISE LINEAR RESPONSEA. Outline
Here we present details of the piecewise linear response (PLR) model discussed in the main article. The PLR modelis based on the observation that solvent response to charging a solute is linear for both anions and cations, but differsbetween the two cases.
In such a model, the average electrostatic potential due to the solvent at the center ofa charged cavity can be written as (cid:104) φ solv (cid:105) q = (cid:40) φ neut − βq (cid:104) ( δφ solv ) (cid:105) + ( q ≥ q c ) φ neut − βq (cid:104) ( δφ solv ) (cid:105) − − βq c (cid:2) (cid:104) ( δφ solv ) (cid:105) + − (cid:104) ( δφ solv ) (cid:105) − (cid:3) ( q < q c ) , (S24)where q c is the value of the ‘crossover charge’ between the two linear regimes, (cid:104) ( δφ solv ) (cid:105) + is the variance of φ solv for q ≥ q c , and (cid:104) ( δφ solv ) (cid:105) − is the variance of φ solv for q < q c . (As written, it is implicitly assumed that q c ≤
0, assuggested by simulations.) Let us define J = (cid:2) (cid:104) ( δφ solv ) (cid:105) + − (cid:104) ( δφ solv ) (cid:105) − (cid:3) . F chg is then, F chg ( q ) = (cid:40) qφ neut − βq (cid:104) ( δφ solv ) (cid:105) + ( q ≥ q c ) qφ neut − βq (cid:104) ( δφ solv ) (cid:105) − − βJ (cid:16) qq c − q c (cid:17) ( q < q c ) , (S25)and ψ is, ψ ( q ) = (cid:40) φ neut ( q ≤ | q c | ) φ neut − βJ q ( q − | q c | ) ( q > | q c | ) . (S26)In general, φ neut , q c and J will depend upon solute size, and whether or not the solute is located in bulk or at theinterface. B. Results
Figures S7, S8 and S9 show (cid:104) φ solv (cid:105) q vs q for R = 0 .
240 nm , .
317 nm and 0 .
415 nm, respectively, both for the solutein bulk and at the interface. Note that these results have not been corrected for the finite size of the simulation cell:we will correct φ neut for finite size effects when computing ∆ ads ψ (PLR) , where other finite size effects largely cancel. For R = 0 .
317 nm and R = 0 .
415 nm we can see that PLR is broadly reasonable for the solute in bulk, but some smalldeviations are seen. These deviations are more pronounced when the solute is at the interface. For R = 0 .
240 nm, theabove PLR model breaks down at large negative q , but it remains reasonable for smaller values of the absolute charge.By fitting straight lines to the anion and cation response, we can obtain values for q c , (cid:104) ( δφ solv ) (cid:105) + and (cid:104) ( δφ solv ) (cid:105) − .The results from using these in Eq. S26 to compute ∆ ads ψ (PLR) are presented in Fig. 4b in the main article. Resultsfor R = 0 .
75 nm and R = 1 nm are not shown because, while anion and cation response do still differ, the degreeof nonlinearity is much less on an absolute scale than for the smaller solutes. This makes it challenging to reliablyobtain q c .9 (a)(b) FIG. S7. (cid:104) φ solv (cid:105) q vs q for R = 0 .
240 nm with the solute located (a) in bulk and (b) at the interface. The dashed and dottedlines show linear fits to the left and right shaded regions, respectively. (a)(b)
FIG. S8. (cid:104) φ solv (cid:105) q vs q for R = 0 .
317 nm with the solute located (a) in bulk and (b) at the interface. The dashed and dottedlines show linear fits to the left and right shaded regions, respectively. (a)(b) FIG. S9. (cid:104) φ solv (cid:105) q vs q for R = 0 .
415 nm with the solute located (a) in bulk and (b) at the interface. The dashed and dottedlines show linear fits to the left and right shaded regions, respectively. S4. CONSTRUCTING P ( φ solv ) In order to compute F chg ( q ) from Eq. 7, we require P ( φ solv ), the probability distribution of φ solv . For the range of q of interest, i.e. − ≤ q/e ≤
1, sampling P directly (in the absence of solute charge) would yield grossly insufficientdata in the extreme wings of the distribution. Instead, we obtain P by histogram reweighting using MBAR. As anillustration, Fig. S10 (a) shows probability distributions P q ( φ solv ) of φ solv at the center of the solute ( R = 0 .
240 nm)with different values of q . Using data from simulations across the full range of q , we then construct P ( φ solv ), asshown in Fig. S10 (b). (a)(b) FIG. S10. (a) P q ( φ solv ) for q/e = − . , − . , . . . , . , . . . , . , . R = 0 .
240 nm. Solid lines indicate normalized Gaussiandistributions and are included as a guide to the eye. (b) P ( φ solv ) reconstructed from the set of P q using MBAR (solid line).The dashed line indicates a normalized Gaussian distribution with mean and variance obtained from the simulation at q/eq/e