Effective dielectric constant of water at the interface with charged C 60 fullerenes
Setare Mostajabi Sarhangi, Morteza M. Waskasi, Seyed Majid Hashemianzadeh, Dmitry V. Matyushov
aa r X i v : . [ phy s i c s . c h e m - ph ] J a n Effective dielectric constant of water at the interface with charged C fullerenes Setare Mostajabi Sarhangi, Morteza M. Waskasi, Seyed Majid Hashemianzadeh, and Dmitry V. Matyushov Molecular Simulation Research Laboratory, Department of Chemistry, Iran University of Science and Technology,Tehran 16846-13114, Iran School of Molecular Sciences, Arizona State University, PO Box 871504, Tempe,AZ 85287-1504 Department of Physics and School of Molecular Sciences, Arizona State University, PO Box 871504, Tempe,AZ 85287-1504 a) Dipolar susceptibility of interfacial water and the corresponding interface dielectric constant were calculatedfrom numerical molecular dynamics simulations for neutral and charged states of buckminsterfullerene C .Dielectric constants in the range 10–22, depending on temperature and solute charge, were found. Thehydration water undergoes a structural crossover as a function of the solute charge. Its main signaturesinclude the release of dangling O-H bonds pointing toward the solute and the change in the preferentialorientations of hydration water from those characterizing hydrophobic to charged substrates. The interfacedielectric constant marks the structural transition with a spike. The computational formalism adopted hereprovides direct access to interface susceptibility from configurations produced by computer simulations. Therequired property is the cross-correlation between the radial projection of the dipole moment of the solvationshell and the electrostatic potential of the solvent inside the solute. I. INTRODUCTION
The dipolar susceptibility of a bulk material χ is mea-sured by the dielectric experiment in terms of the elec-trostatic free energy stored in a plane capacitor. Thesusceptibility defines the bulk dielectric constant ǫ =1 + 4 πχ . Whether this material property can be appliedto interfaces of molecular or mesoscopic dimensions haslong been a subject of contention. It has long beensuggested that an effective dielectric constant of a mi-croscopic interface has to be introduced, and most re-searchers have agreed that this effective dielectric con-stant has to be reduced from the bulk value.
Theextent of reduction has mostly remained unknown. Fol-lowing earlier indications, direct measurements havebeen recently reported for the dielectric constant of wa-ter in contact with a graphite substrate as a functionof the water film thickness. The dielectric constant inthe direction perpendicular to the graphite plane wasfound to be as low as ∼ ∼ ǫ ⊥ , and parallel, ǫ k , to the sub-strate plane. Ref. 8, therefore, reports ǫ ⊥ ∼ ǫ ⊥ ∼ . while ǫ ⊥ ∼ In this study, we have addressed the problem of the in-terfacial dipolar response for a somewhat related type ofthe interface formed between buckminsterfullerenes C z a) Electronic mail: [email protected] and SPC/E water. Here, z denotes the total chargeof the fullerene, which, in practical applications, canbe altered by electrochemistry. We apply moleculardynamics (MD) simulations to study the effects of so-lute’s charge ( − ≤ z ≤
1) and water’s temperature(240 ≤ T ≤
360 K) on the interface susceptibility. Thelinear dipolar susceptibility of the interface is used todefine the interface dielectric constant ǫ int , which is aproperty characterizing the interface and distinct fromthe bulk dielectric constant ǫ ( ǫ ∼
71 for SPC/E waterat 300 K ). We find that ǫ int ∼ −
22 is not sig-nificantly affected by temperature, but is much strongeraffected by the solute charge. The interface susceptibilityas a function of charge z passes through a spike markinga structural crossover of the hydration shell. II. INTERFACE SUSCEPTIBILITY
The dielectric constant of a bulk dielectric is a ma-terial property, independent of sample’s shape, becauseof the locality of the Maxwell field E . The Maxwellfield is defined as an ensemble average of the microscopicfield E m , followed by a coarse-graining protocol averagingout molecular-scale oscillations of the microscopic field. The postulated locality of the Maxwell field allows oneto relate it to the local polarization density P throughthe scalar susceptibility, P = χ E .Despite this widely adopted reasoning, E itself is neveraccessible experimentally and only the line integral ∆ φ = R E · d l connects the Maxwell field to experimentally ac-cessible voltage difference ∆ φ between the end points. Dielectric experiments take advantage of the fact that E z is uniform in a plane capacitor and find the z -projectionof the Maxwell field E z = ∆ φ/d in terms of the sepa-ration d between the plates ( z -axis is perpendicular tocapacitor’s plates).The Maxwell field must be non-uniform in the inter-face and, by that fact, it becomes a property not accessi-ble to measurements. Also, locality of an inhomogeneousMaxwell field has never been established. The fundamen-tal difficulty of measuring local fields in non-uniformlypolarized dielectrics has long been recognized. The onlyknown resolution is to either measure the field producedby a polarized dielectric in vacuum or to measure fieldsinside small cavities carved within the dielectric. The second strategy is realized in experiments record-ing solvent-induced shifts of optical lines giving access tolocal cavity fields. The locality of polar response disappears for inter-faces, which cannot be characterized by a well-definedscalar susceptibility. In fact, the locality of the responseis gone even in the bulk at microscopic length-scales,when bulk dielectric susceptibility χ is converted to anonlocal 2-rank tensor response function χ ( r − r ′ ). Itreduces to longitudinal and transverse dielectric projec-tions χ L,T ( k ) depending on the scalar wavevector k in re-ciprocal space. For inhomogeneous polarization en-countered in solvation and interfaces, the polar suscep-tibility, χ ( k , k ), loses its isotropic symmetry and be-comes a function of two wavevectors. A complete microscopic solution for the interface po-larization is clearly complex and is only remotely relatedto dielectric properties of the bulk. One still wondersif a coarse-grained description in terms of an effectivesusceptibility of the interface can be formulated. It isclear that any such definition will not be unique and islikely to apply to a set of problems for which the responseof the interface is well defined. Our focus is on inter-face polarization in terms of the field it produces withina cavity in response to a probe charge. Since theCoulomb law applies also to microscopic fields, this goalcan be achieved by a proper reformulation of the bound-ary value problem. The interface susceptibility givingaccess to the field inside a void is not necessarily trans-ferrable to other electrostatic problems and, specifically,to another well-established problem where dielectric con-stant is prominent: the screening of ions in solution.It is easy to realize that a length-scale is involvedin dielectric screening: the potential of mean force be-tween ions is an oscillatory function of the distance atmolecular scale, obviously not reducible to a sin-gle screening parameter. Likewise, the average polariza-tion density in the interface induced by a spherical ion h P r i = h ˆr · P i ( ˆr = r /r is the radial unit vector) is anoscillatory function of the radial distance r as found insimulations by Ballenegger and Hansen and in a num-ber of follow-up simulation studies of spherical and pla-nar interfaces. There is obviously a length-scale,dictated by the microscopic structure, specific to the lo-cal polarization density of the interface. A definition ofa scalar parameter of interfacial polarity based solely on P ( r ) is not possible.The polarization of the interface P ( r ) is not directlyaccessible experimentally and is not even required if the focus is on the measurable local field inside a solute orcavity in the dielectric. From this perspective, the ques-tion at hand is what is the integrated response of the in-terface to a probe charge. Posed in this way, the problemof an effective susceptibility of the interface can poten-tially be formulated without a length-scale involved, if aproper coarse-graining formalism is formulated. Cast inthis way, the problem becomes somewhat similar to theproblem of surface tension, which is clearly a microscopicinterfacial property, but without a length-scale involved.One asks the question of what is the integrated responseof the interface to altering its surface area. Similarly, weare asking what is the integrated polarization of the in-terface creating a certain field inside a void. Presentinga formalism to address this question and its applicationto a realistic interface of charged buckminsterfullerenesin water is the goal of this study.The microscopic electric field E m = −∇ φ m is ex-pressed in terms of the microscopic electrostatic poten-tial φ m in the presence of the external charge density ρ ( r ) = qδ ( r ) and the instantaneous density of boundcharge ρ b = −∇ · P . The fluctuating electrostatic poten-tial satisfies the Laplace equation at each configurationof the liquid ∇ φ m = − π [ ρ + ρ b ] . (1)The density of bound charge, a scalar field, is in turnexpressed through the divergence of the fluctuating po-larization field P = P d − ∇ · Q + . . . , which includesthe dipolar field P d and spatial derivatives of densities ofhigher multipoles, starting from the quadrupolar 2-ranktensor Q . When a solute is immersed in a polar liquid, ρ b = 0 inside the solute and φ m is determined from thestandard Laplace equation, ∇ φ m = − πρ . A signif-icant result here is that one can find a solution for anensemble-averaged potential h φ m i provided the bound-ary conditions are additionally supplied. Therefore, find-ing the electric field inside a cavity in the dielectric isreduced to the question of formulating proper bound-ary conditions that preserve, in a coarse-grained man-ner, some information about the molecular structure ofthe interface. The boundary conditions account for the discontinuityof the electric field at some dividing surface between thesolute and the liquid E n − h E n i = 4 π [ σ + h P n i ] . (2)In this equation, E n = − b n · ∇ φ is the normal pro-jection of the electric (vacuum) field inside the soluteand h E n i = − b n · ∇h φ m i is the normal projection of theensemble-averaged electric field inside the solvent, bothare taken at a nonspecified dividing surface separatingthe solute from the solvent. The surface normal b n is di-rected outward from the solvent into the solute (Fig. 1).Equation (2) averages over the fluctuations of the solute-solvent system by taking ensemble averages h . . . i . Fur-ther, σ in Eq. (2) is the density of free charges at the di- ˆn a R m j ˆr j Figure 1. Schematics of the C -water interface. The spheri-cal dividing surface with the radius a (solid line) has the nor-mal unit vector b n pointing outward from the liquid toward thesolute. The spherical region in water with the radius R is usedto calculate the total dipole moment of the water moleculeswithin the R -sphere, Eq. (10). The scalar normal projectionof the dipole moment is calculated by projecting each indi-vidual dipole moment m j (blue arrow) of the water moleculewith r j < R on the radial unit vector b r j = − b n j , b r j = r j /r j (red arrow). viding surface. This component of the electrostatic prob-lem is important for our simulations of charged fullerenesC z since the solute charge is distributed over the surfaceatoms when z = 0. This charge configuration is differentfrom solvation of small ions for which only P n enters theboundary condition. The normal projection of the instantaneous polariza-tion density of water P n = b n · P is the main focus ofour formalism and of the simulations presented below.Its ensemble average σ b = h P n i on the right-hand side ofEq. (2) is the density of the bound (water) charge at thedividing surface. The surface charge density is a param-eter quantifying the preferential alignment of dipoles inthe interface. Its relation to the bulk properties of the di-electric material is different for liquid and solid dielectricsas can be highlighted by first looking at the results fol-lowing from the theories of continuum (macroscopic) di-electrics.The result for the surface charge density is particu-larly simple for the spherical symmetry of the dielec-tric interface and the polarizing electric field. When aprobe charge q is placed at the center of a spherical cav-ity with the radius a carved from a dielectric, the sur-face charge density becomes σ b = − (1 − ǫ − )( q/S ) =(4 π ) − (1 − ǫ − ) E n , where S = 4 πa is the surface area. The electrostatic potential of the bound charge inside thedielectric becomes φ b ( r ) = − (cid:0) − ǫ − (cid:1) ( q/r ) . (3)The surface charge is, therefore, opposite to the probecharge and screens it. When combined with the vac-uum potential φ ( r ) = q/r , Eq. (3) yields the standardCoulomb potential screened by the dielectric, h φ m ( r ) i = q/ ( ǫr ). !" % & $ ,-./0 Figure 2. Density of bound charge of water ρ b ( r ) [Eq. (4)] forC z with the solute charges indicated in the plot. Equation (3) assumes that a bulk material property,the dielectric constant here, can define a property of theinterface, the dipolar polarization of the interface in ourcase. This assumption strictly applies only to solid di-electrics, which can propagate bulk stress through theentire material by means of a uniform strain when theuniform stress/field is applied. One can view the polar-ization of the interface as preferential alignment of inter-facial dipoles uniformly propagating from the bulk. Thepolarization density in the interface P is then the sameas in the bulk when the dielectric is uniformly polarized(plane capacitor), which is the meaning of the notion ofa continuum dielectric described by the boundary condi-tions of Maxwell’s electrostatics.Liquids do not maintain bulk stress, and the polariza-tion in response to an external field must form in a sur-face layer of molecular dimension. As mentioned above,the issues involved are similar to the distinction betweenthe surface tension, a macroscopic property characteriz-ing interface only, and the cohesive energy of the bulk.Drawing from this analogy, bulk dielectric constant doesnot necessarily describe surface polarization. Two dif-ferent susceptibilities are, therefore, required: the inter-face susceptibility to describe polarizability of the inter-face and bulk dielectric constant. The former describesorientational preferences of the interfacial dipoles, whichare strongly affected by the local interfacial structure.The latter describes the buildup of dipolar correlationsby chains of mutually induced dipoles producing long-ranged, ∝ r − , correlations ultimately responsible for di-electric screening. There is no direct link between thesetwo susceptibilities since the alignment of dipoles in theinterface is a function of both the liquid and the sub-strate.In contrast to the dielectric constant of the bulk,the dipolar susceptibility of the interface has not beenuniquely defined. The vector field of the dipole momentdensity P ( r ) is highly oscillatory in the interface and doesnot provide a scalar susceptibility, even if a specific pro-jection is taken. This is illustrated in Fig. 2 where !" ! % & $$$ ’()*+ ,-,./0 ,1$, ,2$, Figure 3. Interface susceptibility χ n ( a ) for C z with the so-lute charges indicated in the plot. The dashed lines (nearlyindistinguishable on the scale of the plot) show coarse-grained χ n calculated from the slopes of χ I in Eq. (10). the density of bound charge of water, ρ b ( r ) = N X i =1 X a =1 , h q i,a δ ( r − r i,a ) i (4)in the interface of buckminsterfullerenes C z is shown(the sum here runs over all i = 1 , . . . , N water moleculeswith atomic charges q i,a at r i,a , where a = 1 − E n and h P n i provides such aclosure h P n i = χ n E n . (5)This route is applied here to calculate χ n from molec-ular dynamics trajectories and to evaluate the interfacedielectric constant ǫ int = ǫ based on the input from nu-merical simulations.Solutes studied here are neutral and charged fullerenesC z , with the charge z varied between z = − z =1. The entire charge z of a charged fullerene is spreadover its surface with the charge density 4 πσ = ze/a = − E n ( a + ) (Fig. 1), where E n ( a + ) is the electric field ofthe solute charges at the outer surface of the fullerene.In contrast, the electric field inside the charged fullerenein zero, E n ( a − ) = 0, if the surface charge is assumedto be uniformly spread. It is this field, E n = E n ( a − )that enters the boundary condition on the left-hand sidein Eq. (2).In the case of the neutral fullerene with z = 0, we as-sume that the solute field E n is produced by the probecharge q placed at the center of the solute. Both casesof charged and neutral fullerenes can be combined in Eq.(2) to obtain a closed-form equation for the dipolar sus-ceptibility χ n connecting h P n i to the ensemble-averaged(Maxwell) field h E n i : h P n i = χ n h E n i . The same ex-pression for the interface dielectric constant ǫ int − !" ! ! " $ % & ’ ( ) * +, - ( Figure 4. Correlation of the radial projection of the shelldipole moment M r ( R ) [Eq. (12)] with the electrostatic po-tential of the solvent φ s at the center of C z as a functionof the radius R of the spherical shell used to calculate thewater dipole moment [Eq. (11)]. Calculations are performedat different values of charge z (plots at different temperaturesare collected in Figs. S5–S8 in supplementary material). Theslope of χ I ( R ) ∝ h δM r ( R ) δφ s i with the radius R determinesthe interface susceptibility χ n [Eq. (10)]. The dashed linesare linear fits through the simulation points. πχ n follows in both cases ǫ int = [1 − πχ n ] − . (6)The perturbation theory gives h P n i as the statisticalcorrelation of the fluctuation δP n = P n − h P n i with thefluctuation of the solute-solvent Coulomb energy U C h P n i = − β h δP n δU C i , (7)where δU C = U C − h U C i . Further, U C = Qφ s is givenas the product of the fluctuating electrostatic potential φ s of the solvent and the solute charge Q : Q is either q (for z = 0) or ze (for z = 0). By combining Eqs. (5) and(7), the solute charge can be eliminated from the linearsusceptibility χ n = βa h δP n δφ s i , (8)where a is the radius of the spherical dividing surface and δφ s = φ s − h φ s i . Because of the spherical symmetry ofthe problem, φ s is taken at the center of C z for bothcharged and neutral solutes.Equation (8) is the integrated form of the equationgiven by Ballenegger and Hansen, in which the two-point correlation function between radial projections ofthe polarization density P r = ˆr · P needs to be integrated χ n = 4 πa β Z ∞ a h P r ( a ) P r ( r ′ ) i dr ′ . (9)The susceptibilities χ n in both equations rely on a spe-cific value of the radius a for the dividing surface. Thisdefinition is not computationally robust since the divid-ing surface separating dielectrics, and, even more so, theseparation between a liquid and a molecular-scale solute,is not well defined in dielectric theories. This problem isshared not only by solvation theories, where the “dielec-tric cavity” can only be empirically established, butalso by the local polar response discussed in the compu-tational literature. Specifically, χ n ( a ) oscillatesas a function of the radius a of the dividing sphere ontowhich both the polarization density and the external fieldare projected (Fig. 3 and Fig. S5 in supplementary ma-terial). Therefore, in order to arrive at a single robustparameter, averaging oscillations of the polarization den-sity out, a coarse-graining protocol was suggested in Ref.46.Instead of using χ n ( a ) calculated at a specific a in Eq.(8), an average χ n over a range of dividing surface radii a it is calculated. Coarse graining of molecular scale in-terfacial oscillations is performed by calculating the slopeof the integrated susceptibility χ I ( R ) vs the radius R ofthe spherical region chosen around the solute (Fig. 1) χ n = dχ I /dR. (10)In turn, the integrated susceptibility4 πχ I ( R ) = − β h δM r ( R ) δφ s i (11)is calculated by correlating the fluctuations of the electro-static potential produced by the solvent at the positionof the probe charge with the radial projection of the totaldipole moment of the solvent within the R -sphere (Fig.1) M r ( R ) = X r j The geometries and charge distributions of C z werecalculated with the density functional theory as de-scribed in supplementary material. These charges werecombined with CHARMM22/OPLS parameters for car-bon atoms and used as the force field for classicalMD simulations performed with NAMD 2.9 softwareprogram. The nonuniform distribution of atomic chargein charged fullerenes is caused by Jahn-Teller distor-tions of icosahedral symmetry characteristic of the neu-tral fullerene (see supplementary material for dis-cussion). However, we found that DFT charges andthe uniform distribution of atomic charge z/ 60 pro-duce indistinguishable results for the interfacial struc-ture of hydration water (Fig. S1 in supplementary ma-terial). Therefore, in addition to integer charges z =+1 , . . . , − z/ z = − . , − . , − . , − . , − . , − . T = 300 K. This set of simula-tions allowed us to obtain the dependence of the inter-facial dipolar susceptibility on the solute charge and toconnect it to the structural crossover of hydration wateroccurring with increasing | z | (see below).The C z solutes were hydrated with 2413 SPC/Ewater molecules. Initial equilibration with NPT wasfollowed by NVT simulations at different temperatures(240, 260, 280, 300, 320, 340 and 360 K). A typical sim-ulation length was 110 ns. A detailed list of simulationtimes is provided in Table S1 in supplementary material,along with other details of the simulation protocol.The main finding of our simulations and their analy-sis is a significant reduction of ǫ int compared to the bulkdielectric constant ( ǫ ∼ 71 for SPC/E ). This resultis relevant for the electrostatic boundary value problemfor which the bulk value ǫ is often used. In contrast,our formulation suggests that the interface dielectric con-stant ǫ int , carrying molecular properties of the interface,should be used in place of ǫ in the boundary conditions !" % &’ ( ) .’ ( .&/) .0 Figure 6. Radial distribution function for the hydrogens ofSPC/E water at different charges z of C z indicated in theplot ( T = 300 K). (Eq. (2)) for the Laplace equation. The resulting valuesof ǫ int depending on temperature and solute charge aresummarized in Fig. 5 and in Table S2 in supplementarymaterial. There are some noticeable changes with tem-perature, particularly for the neutral fullerene (a dashedblue line in Fig. 5). However, a much stronger variationof the interface dielectric constant is observed with thesolute charge z , reflecting a structural crossover in thefirst hydration shell for charged solutes.The structure of water interfacing charged fullerenes issignificantly altered compared to the bulk. Water’s den-sity in the first hydration layer increases with increasingcharge. More importantly, interfacial water undergoes astructural transition from preferential orientations spe-cific to hydrophobic interfaces to an orientational struc-ture characteristic of charged substrates. Signatures ofthe structural transition are seen already at z = − , − z = − , − ∼ 40 first-shell waters at z = − , − ∼ to a largepopulation of dangling O-H bonds was recorded by x-rayabsorption of water on gold substates under a negativebias. The release of dangling O-H bonds at the pointof crossover can be characterized by the order parametergiven by the fraction n OH s of dangling O-H in the first hy-dration shell. It is calculated from the relative area ofthe closest peak in the solute-hydrogen pair distributionfunction (Fig. 6). This order parameter is plotted in the !" ! $ % & ’ (" () ( + - . $ / (" () ( Figure 7. Upper panel: Interface dielectric constant at T =300 K vs the fullerene charge z . The dashed lines are fits toCurie-type functions a + b/ ( z − z z . The dashed line is a fit to a hyperbolic tangent functionoften appearing in mean-field theories of phase transitions. lower panel of Fig. 7 vs the solute charge z (also see TableS3). The dashed line fitting through the points is a hy-perbolic tangent function often appearing in mean-fieldtheories of phase transitions. The structural crossover,carrying some phenomenology of bulk phase transitions,is accompanied by a spike in the interface dielectric con-stant shown in the upper panel of Fig. 7. The dashedlines fitting the points are of Curie type: a + b/ ( z − z ).This functionality appears in the Landau theory of phasetransitions when the quadratic term in the free energyfunctional is taken in the form ∝ ( z − z ) (the standardCurie law follows from ∝ ( T − T )). One, however, shouldnot anticipate a Curie-type singularity found for suscep-tibilities characterizing bulk phase transitions because ofa small number of water molecules involved. Neverthe-less, the phenomenology of Curie’s law is approximatelyretained for the interface dielectric constant approachingthe crossover point.Our simulations employ a non-polarizable water modelaccounting only for the response of the nuclear coordi-nates. As a minimum, electronic susceptibility ǫ ∞ − ǫ int to account for electronic po-larizability of water molecules. Here, ǫ ∞ is the high-frequency dielectric constant . Its precise value for wa-ter is not known: values between squared refractive in-dex, ǫ ∞ ∼ . 8, and ǫ ∞ ∼ . The situation is potentially more complex,as is seen for the water-air interface. The average, mean-field dipole of water changes in the water-air interfacefrom a bulk value, enhanced relative to the gas phase, tothe gas-phase dipole. While the interface with a soluteor with a planar substrate is obviously distinct from thewater-air interface, a nonuniform interfacial electric fieldmight lead to an effective water dipole different from thebulk.The interface dielectric constant considered here is agauge of the integral ability of the interface to polarizein response to a probe charge placed inside a void in apolar liquid. In contrast, thermally-driven fluctuationsof the overall shell dipole M ( R ) can be gauged by ei-ther its variance h [ δ M ( R )] i or by the scalar product, h δ M ( R ) · δ M s i , with the total dipole moment of the sam-ple M s . By this measure, dipole fluctuations actuallyincrease, and not decrease, compared to the bulk for hy-drophilic solutes. The variance of the shell dipolesscales with the local density and, consistent with thislogic, a recent paper reports a drop of water’s dielectricconstant at model hydrophobic interfaces with reducedsurface density (dewetting). However, the Kirkwood for-mula was incorrectly applied to the calculations of thedielectric tensor in the interface. The Kirkwood for-mula is derived by tracing the dipolar susceptibility overits longitudinal and transverse components. For theslab geometry, those correspond to linear responses per-pendicular and parallel to the slab. The Kirkwoodequation, therefore, cannot be applied to components ofthe dielectric tensor, which can carry different symme-tries in respect to the longitudinal and transverse dipolarsusceptibilities. The variance of the shell dipole moment, which can beused to characterize dipolar fluctuations in the interface,does not directly enter the electrostatic boundary-valueproblem. It is only the normal projection of the dipolemoment that defines the interface susceptibility andis required for electrostatics. Only this susceptibility islower in the interface than in the bulk. This result im-plies suppression of the interfacial response in the normaldirection. Dipoles in the interface, frustrated by the localfields and geometric constraints, do not developthe complete dielectric screening of the bulk material.The deviation between ǫ int and ǫ is less pronounced forthe fluid of dipolar hard spheres interfacing a repulsivevoid. The origin of a large difference between ǫ int and ǫ for water is likely a signature of its specific interfacialorientational structure, which is difficult to characterizein more detail without expanding the interface suscepti-bility in basis functions sensitive to orientational dipolarorder. IV. CONCLUSIONS The computational formalism used here allows directaccess to interface susceptibility from configurations pro-duced by computer simulations. The required property isthe cross-correlation of the radial projection of the dipolemoment of the solvation shell with the electrostatic po-tential of the solvent inside the solute.This computational formalism has been applied to sim-ulations of a realistic interface chemically similar to thewater-graphite interface studied experimentally in Ref. 8, where dielectric constant ǫ ⊥ ∼ ∼ ǫ int ∼ − 22 are reported here for chargedfullerenes interfacing SPC/E water. Interface dielectricconstants in the same range, ∼ − 9, were previouslycalculated from simulations of model Lennard-Jones so-lutes of different sizes in TIP3P water. The interfacedielectric constant ∼ − 15 ( T = 240 − 360 K) for theneutral fullerene is consistent with ǫ int ∼ ∼ . ∼ . 75 ˚A for C ). Since the non-polarizable forcefields for water miss the response of the electronic polar-izability, our results likely constitute the lower bound forthe interface susceptibility. The interface dielectric con-stant is the property of an interfacial layer of water ofmolecular scale and is physically distinct from the bulkdielectric constant reflecting dipolar correlations in thebulk ( ǫ ∼ 71 for SPC/E water ).Increasing the charge | z | of hydrated fullerenes leadsto a structural transition of hydration water. It is char-acterized by breaking the interfacial network of hydrogenbonds and the release of dangling O-H bonds pointing to-ward the solute. The interface dielectric constant marksthis structural crossover with a spike. SUPPLEMENTARY MATERIAL See supplementary material for the simulation proto-cols, data analysis, and additional plots of the interfacesusceptibility at different temperatures. ACKNOWLEDGMENTS This research was supported by the National ScienceFoundation (CHE-1800243). This work used the Ex-treme Science and Engineering Discovery Environment(XSEDE) through allocation TG-MCB080071. REFERENCES J. D. Jackson, Classical Electrodynamics (Wiley, New York,1999). I.-C. Yeh and M. L. Berkowitz, J. Chem. Phys. , 7935 (1999). O. Teschke and E. F. de Souza, Phys. Chem. Chem. Phys. ,3856 (2005). J.-J. Velasco-Velez, C. H. Wu, T. A. Pascal, L. F. Wan, J. Guo,D. Prendergast, and M. Salmeron, Science , 831 (2014). B. Shi, M. V. Agnihotri, S.-H. Chen, R. Black, and S. J. Singer,J. Chem. Phys. , 164702 (2016). C. Zhang, J. Chem. Phys. , 031103 (2018). L. B. Dreier, Y. Nagata, H. Lutz, G. Gonella, J. Hunger, E. H. G.Backus, and M. Bonn, Sci. Adv. , eaap7415 (2018). L. Fumagalli, A. Esfandiar, R. Fabregas, S. Hu, P. Ares, A. Ja-nardanan, Q. Yang, B. Radha, T. Taniguchi, K. Watanabe,G. Gomila, K. S. Novoselov, and A. K. Geim, Science ,1339 (2018). P. Debye and L. Pauling, J. Am. Chem. Soc. , 2129 (1925). B. E. Conway, J. O. Bockris, and I. A. Ammar, Trans. Farad.Soc. , 756 (1951). L. S. Palmer, A. Cunliffe, and J. M. Hough, Nature , 796(1952). P. J. Stiles, J. B. Hubbard, and R. F. Kayser, J. Chem. Phys. , 6189 (1982). B. W. Ninham and V. Yaminsky, Langmuir , 2097 (1997). G. E. Brown, V. E. Henrich, W. H. Casey, D. L. Clark, C. Eggle-ston, A. Felmy, D. W. Goodman, M. Gr¨atzel, G. Maciel, M. I.McCarthy, K. H. Nealson, D. A. Sverjensky, M. F. Toney, andJ. M. Zachara, Chem. Rev. , 77 (1999). O. Teschke, G. Ceotto, and E. F. de Souza,Phys. Rev. E , 011605 (2001). P. J. Lenart, A. Jusufi, and A. Z. Panagiotopoulos, J. Chem.Phys. , 044509 (2007). M. C. F. Wander and A. E. Clark, J. Phys. Chem. C , 19986(2008). E. Lima, D. Horinek, R. Netz, E. C. Biscaia, F. W. Tavares,W. Kunz, and M. Bostr¨om, J. Phys. Chem. B , 1580 (2008). M. D. Boamah, P. E. Ohno, F. M. Geiger, and K. B. Eisenthal,J. Chem. Phys. , 222808 (2018). H. A. Stern and S. E. Feller, J. Chem. Phys. , 3401 (2003). V. Ballenegger and J.-P. Hansen, J. Chem. Phys. , 114711(2005). S. Gekle and R. R. Netz, J. Chem. Phys. , 104704 (2012). C. Zhang, F. Gygi, and G. Galli, J. Phys. Chem. Lett. , 2477(2013). S. De Luca, S. K. Kannam, B. D. Todd, F. Frascoli, J. S. Hansen,and P. J. Daivis, Langmuir , 4765 (2016). H. J. C. Berendsen, J. R. Grigera, and T. P. Straatsma, J. Phys.Chem. , 6269 (1987). C. A. Reed and R. D. Bolskar, Chem. Rev. , 1075 (2000). C. J. Fennell, L. Li, and K. A. Dill, J. Phys. Chem. B , 6936(2012). L. D. Landau and E. M. Lifshitz, Electrodynamics of ContinuousMedia (Pergamon, Oxford, 1984). L. Eyges, The Classical Electromagnetic Field (Dover Publica-tions, Ney York, 1972). W. Thompson Lord Kelvin, Reprint of Papers on Electrostaticsand Magnetism , 2nd ed. (MacMillan and Co., London, 1884, sec479). D. R. Martin, A. D. Friesen, and D. V. Matyushov, J. Chem.Phys. , 084514 (2011). S. D. Fried and S. G. Boxer, Acc. Chem. Res. , 998 (2015). O. V. Dolgov, D. A. Kirzhnits, and E. G. Maksimov, Rev. Mod.Phys. , 81 (1981). T. Fonseca and B. M. Ladanyi, J. Chem. Phys. , 8148 (1990). F. O. Raineri, H. Resat, and H. L. Friedman, J. Chem. Phys. , 3068 (1992). P. A. Bopp, A. A. Kornyshev, and G. Sutmann, Phys. Rev. Lett. , 1280 (1996). D. Chandler, Phys. Rev. E , 2898 (1993). D. V. Matyushov, J. Chem. Phys. , 1375 (2004). S. E. Huston and P. J. Rossky, J. Phys. Chem. , 7888 (1989). A. A. Rashin, J. Phys. Chem. , 4664 (1989). J. S. Bader and D. Chandler, J. Phys. Chem. , 6423 (1992). C. J. Fennell, A. Bizjak, V. Vlachy, and K. A. Dill, J. Phys.Chem. B , 6782 (2009). D. J. Bonthuis, S. Gekle, and R. R. Netz, Langmuir , 7679(2012). C. Schaaf and S. Gekle, Phys. Rev. E , 032718 (2015). M. Dinpajooh and D. V. Matyushov, J. Chem. Phys. , 014504(2016). D. V. Matyushov, J. Chem. Phys. , 224506 (2014). J. C. Maxwell, A Treatise on Electricity and Magnetism , Vol. 1(Dover Publications, New York, 1954, sec. 63). M. S. Wertheim, J. Chem. Phys. , 4291 (1971). J. S. Høye and G. Stell, J. Chem. Phys. , 562 (1974). L. Horv´ath, T. Beu, M. Manghi, and J. Palmeri, J. Chem. Phys. , 154702 (2013). C. Zhang and M. Sprik, Phys. Rev. B , 144201 (2016). B. Roux, H.-A. Yu, and M. Karplus, J. Phys. Chem. , 4683(1990). S. W. Rick and B. J. Berne, J. Am. Chem. Soc. , 3949 (1994). D. J. Bonthuis, S. Gekle, and R. R. Netz, Phys. Rev. Lett. ,166102 (2011). M. A. Wilson, A. Pohorille, and L. R. Pratt,J. Chem. Phys. , 5211 (1989). J. C. Phillips, R. Braun, W. Wang, J. Gumbart, E. Tajkhorshid,E. Villa, C. Chipot, R. D. Skeel, L. Kal´e, and K. Schulten,J. Comput. Chem. , 1781 (2005). J. Niklas, K. L. Mardis, and O. G. Poluektov, J. Phys. Chem.Lett. , 3915 (2018). S. M. Sarhangi, M. M. Waskasi, S. M. Hashemianzadeh, andD. Matyushov, Phys. Chem. Chem. Phys. , 27069 (2018). C. Y. Lee, J. A. McCammon, and P. J. Rossky, J. Chem. Phys. , 4448 (1984). Y. R. Shen and V. Ostroverkhov, Chem. Rev. , 1140 (2006). J. G. Davis, B. M. Rankin, K. P. Gierszal, and D. Ben-Amotz,Nat. Chem. , 796 (2013). H. E. Stanley, Introduction to phase transitions and critical phe-nomena (Oxford University Press, New York, 1987). C. J. F. B¨ottcher, Theory of Electric Polarization, Vol. 1: Di-electrics in Static Fields (Elsevier, Amsterdam, 1973). U. Kaatze, J. Soln. Chem. , 1049 (1997). F. S. Cipcigan, V. P. Sokhan, A. P. Jones, J. Crain, and G. J.Martyna, Phys. Chem. Chem. Phys. , 8660 (2015). J. Mart´ı, G. Nagy, E. Gu`ardia, and M. C. Gordillo, J. Phys.Chem. B , 23987 (2006). S. Mukherjee, S. Mondal, and B. Bagchi, J. Chem. Phys. ,024901 (2017). A. D. Friesen and D. V. Matyushov, Chem. Phys. Lett. , 256(2011). T. Sato, T. Sasaki, J. Ohnuki, K. Umezawa, and M. Takano,Phys. Rev. Lett. , 206002 (2018). G. Stell, G. N. Patey, and J. S. Høye, Adv. Chem. Phys. , 183(1981). P. Madden and D. Kivelson, Adv. Chem. Phys. , 467 (1984). D. V. Matyushov, J. Chem. Phys. , 087102 (2016). D. V. Matyushov, in Nonlinear Dielectric Spectroscopy , editedby R. Richert (Springer, Cham, 2018) pp. 1–34. Y.-C. Wen, S. Zha, X. Liu, S. Yang, P. Guo, G. Shi, H. Fang,Y. R. Shen, and C. Tian, Phys. Rev. Lett. , 016101 (2016). Q. A. Besford, M. Liu, and A. J. Christofferson, J. Phys. Chem.B122