The nonlocal dielectric response of water in nanoconfinement
Geoffrey Monet, Fernando Bresme, Alexei Kornyshev, Hélène Berthoumieux
TThe nonlocal dielectric response of water in nanoconfinement
G. Monet and H. Berthoumieux
Sorbonne Universit´e, CNRS, Laboratoire de Physique Th´eorique dela Mati`ere Condens´ee (LPTMC, UMR 7600), F-75005 Paris, France
F. Bresme and A. Kornyshev
Department of Chemistry, Molecular Sciences Research Hub,Imperial College London, W12 0BZ 2AZ London, United Kingdom
Recent experiments reporting a very low dielectric permittivity for nanoconfined water have re-newed the interest to the structure and dielectric properties of water in narrow gaps. Here, wedescribe such systems with a minimal Landau-Ginzburg field-theory composed of a nonlocal bulk-determined term and a local water-surface interaction term. We show how the interplay between theboundary conditions and intrinsic bulk correlations encodes dielectric properties of confined water.Our theoretical analysis is supported by molecular dynamics simulations and comparison with theexperimental data.
Introduction -
Interest in the dielectric propertiesof confined water has been boosted by the remarkedmeasurement of the dielectric permittivity of nanomet-ric water layer confined between hydrophobic surfaces[1]. Fumagali et al. reported an anomalously low di-electric constant in the direction perpendicular to thesurface. [2] Water permittivity in the vicinity of a sur-face is inhomogeneous[3, 4] leading to a significant in-crease of the electrostatic interactions, as postulated inthe 1950’s by Schellman,[5] and observed experimentallyand in simulations [6–8]. The stability of emulsions andcolloidal solutions [9, 10], ion transport and reactivityin channels of proteins,[11], in subsystems of geologi-cal interest [12] or in nanotechnologic devices [13] arestrongly influenced by electrostatic properties of confinedwater. However, a fundamental analytic theory connect-ing the dielectric response to the properties of the con-fining surfaces, namely chemical composition, degree andgeometry of confinement, is still outstanding.[14] At themolecular scale, the relative dielectric permittivity ten-sor (cid:15) α,β ( ~r − ~r ) of bulk water is non local.[15–17] Thestructuration in the fluid at an interface induced by thisnonlocality has been widely studied at the atomic scaleusing molecular dynamics (MD) simulations [18–22]. Ata coarse-grained scale, continuum nonlocal electrostat-ics provide a useful framework to quantify the dielectricproperties of confined correlated fluids. [19] This can bebased on phenomenological energy functionals that arewritten in terms of the polarization field ~m . They arethe sum of the electrostatic energy depending on the dis-placement field ~D and of a correlation term [23–26]. Itreads U bulk [ ~m, ~D ] = 12 (cid:15) Z d~r (cid:16) ~D − ~m ( ~r ) (cid:17) + 12 (cid:15) Z d~rd~r m α ( ~r ) K α,β ( ~r, ~r ) m β ( ~r ) , (1)where (cid:15) is the vacuum dielectric permittivity. We specify the kernel K α,β ( ~r, ~r ) to mimic the simu-lated nonlocal dielectric properties of bulk water. We fur-ther introduce a phenomenological interaction energy be-tween the surface and the fluid as a sum of harmonic po-tentials. We show that this framework reproduces bothMD simulations for two hydrophobic surfaces, grapheneand hexagonal boron nitride (hBN), and an experimentaldata.[1] In addition, it formalizes the effect of the con-fining material on the dielectric properties of ’interfacialwater’. Bulk water -
The dielectric properties of bulk water areencoded in the two-points susceptibility tensor χ α,β ( ~r − ~r ) = δ α,β ( ~r − ~r ) − (cid:15) − α,β ( ~r − ~r ). This nonlocal kernel canbe expressed through the correlations of the polarization ~m using the classical approximation for the fluctuation-dissipation theorem[27], χ α,β ( ~r − ~r ) = h m α ( ~r ) m β ( ~r ) i (cid:15) k B T . (2)The correlations h m α ( ~r ) m β ( ~r ) i can be written in termsof the experimentally measured partial HH, OH, OOstructure factors of water[28] under the assumptionof simple point charges localized at the atoms ofmolecules.[17] The q − dependence of longitudinal partof the susceptibility in the Fourier space ˆ χ k ( q ) illus-trates the nonlocal nature of dielectric properties wa-ter (see Fig. 1). The main peak of ˆ χ k ( q ) (centeredat q = 30 nm −
1) exceeds 1, corresponding to a rangeof wavelengths associated with a negative permittivity (cid:15) k ( q ) = 1 / (1 − χ k ( q )). This overscreening zone is a con-sequence of the H-bonding network in water[17].To model these properties, we follow a Landau-Ginzburg (LG) approach which proved its value in thestudy of critical surface phenomena.[29] We choose thefollowing form of the second item in Eq.(1), U m [ ~m ] = Z d~r (cid:15) (cid:18) K ~m + K l (cid:16) ~ ∇ ~m (cid:17) + β (cid:16) ~ ∇ (cid:16) ~ ∇ ~m (cid:17)(cid:17) (cid:19) , (3) a r X i v : . [ phy s i c s . c h e m - ph ] F e b
10 20 30 40 50 600 . . . . . . ˆ χ k ( q ) / ˆ χ k , m a x q [nm − ]20 40 6002040 ˆ χ k ( q ) ModelFrom exp. structure factors
FIG. 1. Dielectric susceptibility of bulk water. Blackdots are recovered from inelastic neutron scattering datafor oxygen-oxygen, hydrogen-hydrogen, and oxygen-hydrogenstructure factors.[17, 28] Red curve has been computed fromEq. (4) with K = 1 / K l = − . × − nm − , β =1 . × − nm − . The inset shows the susceptibility whichhas not been normalized to 1. which includes terms up to second spacial derivative ofthe field and leads to the longitudinal susceptibility,ˆ χ k ( q ) = 11 + K + K l q + βq . (4)For derivation and discussion see the supporting material(SM). The model parameters ( K , K l and β ) are chosento capture: (i) the permittivity of bulk water at q = 0,(ii) the position of the first peak and (iii) its width athalf height of the simulated or experimentally recovered χ k ( q ). The theoretical susceptibility is plotted in red inFig. 1. Its poles define a decay length λ d and a period λ o , λ d = 2 √ β q p ( β (1 + K ) + K l , λ o = 4 π √ β q p β (1 + K ) − K l , (5)characterizing the polarization correlations in bulk. Theyare equal to λ d = 2 . (cid:6) A and λ o = 2 . (cid:6) A for the chosenparametrization.
Theoretical model for interfacial water -
We considerwater delimited by a planar interface infinite in the xy plane and located at z =0 (See Fig. 2a). A static ho-mogeneous external field ~D = D ~u z is applied in thez-direction. According to the symmetry of the problem,this field excites exclusively the longitudinal polarizationthat depends on z : ~m ( ~r ) = m ( z ) ~u z . We write the energyof the system per unit area U [ m, D ] = U bulk + U s , thesum of the bulk-determined term, U bulk , derived from(Eqs. 1,3), and a surface term U s as U bulk = Z ∞ z =0 dz (cid:15) h ( D − m ) + Km + K l ˙ m + β ¨ m i U s = k m m (0) − m )) + k ρ ρ (0) − ρ )) (6)where the upper dot stands for the spatial derivationalong z . In the spirit of the LG development used to ex-press the kernel K (Eq. (1)), U s is written as an expansionof elastic energies[29, 30] depending on the polarizationfield and its derivative ˙ m ( z ), equal to minus the boundcharge, ρ ( z ).[31] The major contribution promotes a sur-face polarization m and the corrective second term fa-vors a water charge density ρ at the interface. The stiff-nesses k m and k ρ quantify the strength of the boundaryconditions. In the strong interaction limit ( k m , k ρ ) → ∞ ,the surface fixes both polarization and charge density atinterface.The partition function of the system, Z [ D ] = R D [ m z ] exp [ − ( U bulk [ m, D ] + U s ) /k B T ] , can be split inthe form Z [ D ] = Z d ¯ md ¯ ρ exp (cid:20) k B T (cid:18) k m m − m ) + k ρ ρ − ρ ) (cid:19)(cid:21) Z m ( z →∞ )=0˙ m ( z →∞ )=0 m (0)= ¯ m ˙ m (0)= − ¯ ρ D [ m ] exp (cid:20) − k B T U bulk [ m, D ] (cid:21) . (7)This includes a partition of the fields m ( z ) satisfying theboundary conditions (right integral), then a sampling ofthe z = 0 boundary conditions ( ¯ m, ¯ ρ ) (left integral). Wefind the mean field solution, m ( z ), by first minimizing U bulk [ m, D ] with respect to m ( z ) leading to(1 + K ) m ( z ) − K l m (2) ( z ) + βm (4) ( z ) = D , (8)with m (0) = ¯ m , ˙ m (0) = − ¯ ρ , m ( z → ∞ ) = 0, ˙ m ( z → ∞ ) = 0. The solution of which is m ( z ) = D K (cid:18) − e − zλd (cid:18) cos( q o z ) + q d q o sin( q o z ) (cid:19)(cid:19) + e − zλd (cid:18) ¯ m (cid:18) cos( q o z ) + q d q o sin( q o z ) (cid:19) − ¯ ρq sin( q o z ) (cid:19) . (9)with q o = 2 π/λ o and q d = 1 /λ d , the wavenumbers of thebulk correlations. Second, we extremalize the total en-ergy of the system, U = U bulk + U s , with respect to ( ¯ m, ¯ ρ )obtained by injecting m ( z ) in Eqs. (6) and performingthe integral over z (see SM). The nature of the extremumdepends on the dimensionless stiffness constants (˜ k m , ˜ k ρ ),given in the SM. U ( ¯ m, ¯ ρ ) admits a minimum for ˜ k m and˜ k ρ belonging to the pointed zone represented in Fig. 2b,to which we restrict our study in the following. Themean field polarization m ( z ) is given by Eq. (9) for¯ m = m s , ¯ ρ = ρ s , the boundary conditions minimizing U ( ¯ m, ¯ ρ ). Their expressions are given in SM.To study the dielectric properties of interfacial wa-ter, we introduce the real space susceptibility, χ ( z ) = dm ( z ) /dD derived from Eq. (9). It quantifies the re-sponse to a homogeneous external field D and is con-stant and equal to χ b = ˆ χ (0) for bulk water.Fig. 2c, d show typical mean field polarization m ( z )and susceptibility χ ( z ) in the interfacial water. We ob-serve a nonvanishing polarization and a nonconstant χ ( z )that are oscillating functions of period λ o in an exponen-tially decaying envelope of range λ d . The surface inducesa layering of the fluid that extends over about 1 nm, alengthscale consistent with many previous simulations ofinterfacial water[3, 32]. The susceptibility shows alterna-tion of underresponding ( χ ( z ) (cid:28) χ b ) and overrespond-ing layers ( χ ( z ) (cid:29) χ b ), typical for overscreening effect.Whereas the amplitude of m ( z ) is a non-trivial functionof the bulk properties, the four parameters of the surfaceinteraction and D , the amplitude of χ ( z ) does not de-pend on ( m , ρ , D ). The interface affects the dielectricproperties of water only through the stiffnesses (˜ k m , ˜ k ρ ).We first study the case of vanishing ˜ k ρ for which χ ( z ) ˜ k ρ =0 can be written as χ ( z ) ˜ k ρ =0 χ b = 1 + ˜ k m e − q d z − ˜ k m (cid:18) cos( q o z ) + q d − q o q d q o sin( q o z ) (cid:19) . (10)The amplitude of χ ( z ) ˜ k ρ =0 decreases with ˜ k m and tendsto a finite value for ˜ k m (cid:29)
1. This case is represented inFig. 2e (blue curve).[33] Then we consider the correctiveeffect of ˜ k ρ in the limit of a large ˜ k m by studying χ ( z ) ˜ k m (cid:29) χ b = χ ( z ) ˜ k m (cid:29) , ˜ k ρ =0 χ b − q d + q o q d q o ˜ k ρ e − q d z k ρ sin( q o z ) . (11)An increasing ˜ k ρ induces a dephasing and an amplitudedecrease up to a factor 2 of χ ( z ) (See 2e). The behaviorof χ ( z ) as a function of (˜ k m , ˜ k ρ ) illustrates that differ-ent surfaces, having stronger or weaker influence on po-larization and partial charge, induce different dielectricproperties of interfacial water. Comparison with MD simulations -
We performed MDsimulations of pure water confined in a slab geometry us-ing the GROMACS MD simulation package.[34] Watermolecules are described with the SPC/E model and thewalls are made up of atoms of frozen positions. We con-sidered graphene and hBN surfaces ( details in the SM). Water0 z D ~e z a ˜ k m ˜ k ρ b − m ( z ) / m c ˜ k m = 9˜ k ρ = 1 χ ( z ) / χ b u l k z [nm] d χ ( z ) / χ b u l k z [nm]˜ k m (cid:29) e (cid:28) − (cid:29) ˜ k ρ FIG. 2. Dielectric properties of water in the vicinity of asurface. a. Sketch of the system. b. Diagram presenting thezone of finite minimum (dotted zone) as a function of ˜ k m and˜ k ρ . Profile of the polarization m ( z ) (c) and the normalizedsusceptibility χ ( z ) (d) computed for (˜ k m = 9 , ˜ k ρ = 1) and( m = −
10 V / nm , ρ /q o = −
10 V / nm). e. Susceptibilitynormalized to the bulk susceptibility computed from Eq. (11)with different values of ˜ k ρ . We analyze the polarization, m MD ( z ) = − R z dzρ MD ( z ) dz , with ρ MD the charge den-sity of water, and the susceptibility χ MD ( z ) =( m MD ( z, D + δD ) − m MD ( z, D )) /δD [3] with δD = 0 . V /nm , in the vicinity of the surfaces. Theprofiles are similar for both surfaces (Fig. 3): first, avacuum layer ( m MD ( z ) = 0, χ MD ( z ) = 0) between thesurface and the liquid, due to the repulsive part ofthe surface-fluid Lennard-Jones (LJ) interaction, thendecaying oscillations over about 1 nm before reaching thebulk value. The theoretical decay λ d and the period λ o are in very good agreement with the simulated ones (seeSM). This validates the derivation of the characteristiclengths of interfacial water from the bulk dielectricsusceptibility, ˆ χ ( q ).In MD simulations, the position of the interfaces is notas clear-cut as in theory due to thermal capillary fluctua-tions and the non-infinitly sharp repulsion of the surface-fluid LJ interaction.[35] This is taken into account byapplying a smearing to the theoretical predictions,˜ f ( z ) = ( G ∗ ( θf )) ( z + z ) , G ( z ) = e − z / η η √ π , (12)with θ being the Heaviside function and f standing for m or χ . The position z and the standard deviation η are determined for each surface by fitting the first oxy-0 . . . − − m e p s i l o n [ V / n m ] a Graphene
Sim. z < z Sim. z ≥ z Model . . . − − b hBN . . . . . χ ( z ) c z [nm] 0 0.50 . . . . . d z [nm] FIG. 3. Comparison between model (in red) and MD simula-tions (in black) for a graphene layer (left panels) and a hBNlayer (right panels). Top (respectively bottom) panels showthe polarization (respectively the susceptibility). Simulationcurves for z ≤ z are represented with dotted lines. gen density peak with a Gaussian G ( z ) which positionand width define z and η (see SM for details). ThehBN surface is characterized by a deeper LJ potential andconsequently a smaller dispersion η than the graphene.Correspondly, m MD ( z ) amplitude is smaller in interfacialwater for graphene than for hBN Figs. 3a-b.We validate the theoretical model in three steps. First,we adjust the simulated susceptibilities with ˜ χ ( z ) definedin Eq. (12). If we choose (˜ k m (cid:29)
1, ˜ k ρ = 0) for grapheneand (˜ k m (cid:29)
1, ˜ k ρ = 0 .
2) for hBN, we obtain a good agree-ment between the calculated and the simulated value ofthe susceptibilities as shown in figures 3c-d. Next, wefit the simulated polarization for graphene surface with˜ m ( z ) by fixing m , the single left unknown parameter forgraphene as ˜ k ρ = 0. Finally, we fit the simulated polar-ization for a hBN surface. Taking the surface polarization m previously determined in the case of graphene, we fix ρ . The comparison between theoretical and simulatedpolarization are presented in Figs. 3a-b. The dotted partof the simulated curves correspond to the vacuum gapand the contribution of hydrogen located in z < z . Thetheoretical model describes this zone as a vacuum gap(see Eq. (12)).Graphene and hBN surfaces are parametrized by ˜ k m (cid:29)
1, thus both surfaces freeze the interfacial polarization to m ( z ) = m which does respond to D . At the micro-scopic scale, this result can be interpreted as the effect ofthe vacuum gap on the organization in the first layer ofwater which optimizes the number of H-bonds.[36] Mostlikely, ˜ k m is very large for a wide variety of surfaces, bothhydrophobic and hydrophilic, as they impose a layout in 10 L [nm]10 (cid:15) e ff ( L ) -L/2 0 L/2 z [nm]02 (cid:15) − ( z ) . . FIG. 4. Effective dielectric permittivity (cid:15) eff of water nanocon-fined in a channel of width L . Comparison between exper-imental measurements reproduced from [1] and theoreticalmodel. the first hydration layer.[3, 37] For a non-vanishing cor-rective term ˜ k ρ , the surface has an effect on the interfacialcharge, ρ s ( z ), and its variation under D . We investi-gate the microscopic origin of this effect by performingMD simulations for artificial surfaces associated with hy-brid properties between graphene and hBN surfaces (seeSM). We find out that it is induced by a large mean depthof the LJ minimum. A non-vanishing ˜ k ρ is related withimportant variations of the interaction energy betweenthe surface and a water molecule in the ( xy ) plane for z = z that constrains the position of water molecules inthis plane. Nanoconfined water -
We use now this theoreticalmodel to derive the dielectric properties of a confinedwater layer. The experimental measurements report aneffective dielectric permittivity up to (cid:15) eff = 2 for a chan-nel of about 1 nm.[1] (reported on Fig. 4). The authorssuggest the existence in the channel of three water lay-ers of homogeneous dielectric properties: two interfaciallayers ( (cid:15) = 2 .
1, thickness: 0 . (cid:15) = 78). We compute the effective permittiv-ity (cid:15) eff = L/ R L (1 − χ ( z )) dz as a function of L for twographene surfaces. Our model can be seen as two vacuumgaps and an inhomogeneous water layer. This inhomo-geneity is not implemented ad hoc but is the signatureof the nonlocal dielectric properties of water, revealedby the boundary conditions. The results are presented inFig. 4. The model reproduces the experimental measuresand catches in particular a non-homogenous behavior ofthe permittivity as a function of L for small L as shownin the insert that cannot be described by a three homo-geneous layer model.[1, 4, 38] Conclusion -
Nanoconfined water is a non-homogeneous dielectric material which properties differdramatically from the bulk. Water-surface interactionand the bulk properties of the fluid combine to inducespecific dielectric profiles. The complexity of this systemis captured by a minimal phenomenological Hamiltoniandepending on the polarization field and composed of (i) aLG forth order development for a bulk-determined termand (ii) a harmonic surface-water term. We show thatthe dielectric susceptibility of interfacial water may bestrongly affected by the coarse-grained parameters (˜ k m ,˜ k ρ , η ) characterizing local surface-water interaction. Itgives a framework to compare graphene and hBN thatcould be predictive for other surfaces and also to derivethe dielectric properties confined water in other geome-tries, such as nanotubes.[39]This work was supported by Sorbonne Sciences undergrant Emergences-193256. HB and GM thank B. Delam-otte for fruitful discussions. [1] L. Fumagalli, A. Esfandiar, R. Fabregas, S. Hu, P. Ares,A. Janardanan, Q. Yang, B. Radha, T. Taniguchi,K. Watanabe, G. Gomila, K. S. Novoselov, and A. K.Geim, Science , 1339 (2018).[2] S. V. Kalinin, Science , 1302 (2018).[3] D. J. Bonthuis, S. Gekle, and R. R. Netz, Langmuir ,7679 (2012).[4] C. Zhang, The Journal of Chemical Physics , 156101(2018).[5] J. A. Schellman, Journal of Physical Chemistry , 472(1953).[6] V. Ballenegger and J.-P. Hansen, The Journal of Chem-ical Physics , 114711 (2005).[7] S. Chen, Y. Itoh, T. Masuda, S. Shimizu, J. Zhao, J. Ma,S. Nakamura, K. Okuro, H. Noguchi, K. Uosaki, andT. Aida, Science , 555 (2015).[8] T. Sato, T. Sasaki, J. Ohnuki, K. Umezawa, andM. Takano, Physical Review Letters , 206002 (2018).[9] V. Bergeron, Journal of Physics: Condensed Matter ,R215 (1999).[10] N. E. Levinger, Science , 1722 (2002).[11] E. Gouaux, Science , 1461 (2005).[12] P. Fenter, S. Kerisit, P. Raiteri, and J. D. Gale, TheJournal of Physical Chemistry C , 5028 (2013).[13] A. Siria, M.-L. Bocquet, and L. Bocquet, Nature ReviewsChemistry , 91 (2017).[14] D. Mu˜noz-Santiburcio and D. Marx, Physical ReviewLetters , 056002 (2017).[15] A. A. Kornyshev, Journal of Electroanalytical Chemistryand Interfacial Electrochemistry , 79 (1986).[16] P. A. Bopp, A. A. Kornyshev, and G. Sutmann, The Journal of Chemical Physics , 1939 (1998).[17] P. A. Bopp, A. A. Kornyshev, and G. Sutmann, PhysicalReview Letters , 1280 (1996).[18] M. A. Vorotyntsev and A. A. Kornyshev, Soviet Electro-chemistry , 560 (1979).[19] A. A. Kornyshev, Electrochimica Acta , 1 (1981).[20] A. A. Kornyshev, E. Spohr, and M. A. Vorotynt-sev, Encyclopedia of Electrochemistry: Online10.1002/9783527610426.bard010201 (2007).[21] A. A. Kornyshev, The Chemical Physics of Solvation ,355 (1988).[22] C. Schaaf and S. Gekle, The Journal of Chemical Physics , 084901 (2016).[23] A. Hildebrandt, R. Blossey, S. Rjasanow, O. Kohlbacher,and H.-P. Lenhof, Physical Review Letters ,10.1103/physrevlett.93.108104 (2004).[24] A. C. Maggs and R. Everaers, Physical Review Letters , 230603 (2006).[25] H. Berthoumieux and A. C. Maggs, The Journal ofChemical Physics , 104501 (2015).[26] M. Vatin, A. Porro, N. Sator, J.-F. Dufrˆeche, andH. Berthoumieux, Molecular Physics , e1825849 (2020).[27] There is a more general formulation taking into accountquantum corrections[16] that is not considered here forsimplicity.[28] A. K. Soper, The Journal of Chemical Physics , 6888(1994).[29] R. Lipowsky and W. Speth, Physical Review B , 3983(1983).[30] A. Ajdari, B. Duplantier, D. Hone, L. Peliti, and J. Prost,Journal de Physique II , 487 (1992).[31] J. D. Jackson, Classical electrodynamics (Wiley, NewYork, 1975).[32] A. Schlaich, E. W. Knapp, and R. R. Netz, Physical Re-view Letters , 048001 (2016).[33] Note that the diverging value ˜ k m = 1 corresponds tothe boundary of the stability region of the phase spaceparameter (see Fig. 2b).[34] Lindahl, A. , Hess, and V. D. Spoel, Gromacs 2019 sourcecode (2018).[35] P. Yang, Z. Wang, Z. Liang, H. Liang, and Y. Yang,Chinese Journal of Chemistry , 1251 (2019).[36] S. Varghese, S. K. Kannam, J. S. Hansen, and S. P.Sathian, Langmuir , 8159 (2019).[37] Q. A. Besford, A. J. Christofferson, J. Kalayan, J.-U.Sommer, and R. H. Henchman, The Journal of PhysicalChemistry B , 6369 (2020).[38] P. Loche, C. Ayaz, A. Wolde-Kidan, A. Schlaich, andR. R. Netz, The Journal of Physical Chemistry B ,4365 (2020).[39] P. Loche, C. Ayaz, A. Schlaich, Y. Uematsu, and R. R.Netz, The Journal of Physical Chemistry B , 10850(2019). upplementary Materials forThe nonlocal dielectric response of water in nanoconfinement G. Monet and H. Berthoumieux
Sorbonne Universit´e, CNRS, Laboratoire de Physique Th´eorique dela Mati`ere Condens´ee (LPTMC, UMR 7600), F-75005 Paris, France
F. Bresme and A. Kornyshev
Department of Chemistry, Molecular Sciences Research Hub antThomas Young Center for Theory and Simulations of Materials,Imperial College London, W12 0BZ 2AZ London, United Kingdom
CONTENTS
The theoretical model 1Bulk water 1Polarization for one-wall geometry system 2Water confined in between two surfaces : a slab geometry 4The details of molecular dynamics simulations 5Periodic water box 6Water confined in slab geometry 7Validation of the fitting strategy of the bulk susceptibility ˆ χ k ( q ) 8Electrostatic field for large slab geometry 9Smearing the surface 9Molecular dynamics simulations as function of the slab width L k ρ THE THEORETICAL MODELBulk water
The functional energy U bulk of the dielectric medium modeling bulk water is written as U bulk [ m , D ] = 12 (cid:15) Z d r ( D ( r ) − m ( r )) + 12 (cid:15) Z d r d r m α ( r ) K α,β ( r , r ) m β ( r ) (S1)= 12 (cid:15) Z d r (cid:0) D ( r ) − D ( r ) · m ( r ) (cid:1) + 12 (cid:15) Z d r d r m α ( r ) ( δ α,β ( r − r ) + K α,β ( r , r )) m β ( r ) , (S2)where we use the Einstein summation convention over repeated indices.The dielectric material is homogeneous and isotropic so K α,β ( r , r ) = K α,β ( | r − r | ) and we can write the non-local contribution of the functional energy (Eq. (S1)) as function of the Fourier transform of the polarization fieldˆ m ( q ) = π ) / R d r m ( r ) e − i qr and introduce the dielectric susceptibility tensor ˆ χ − α,β ( q ) : Z d r d r m α ( r ) ( δ α,β ( r − r ) + K α,β ( r , r )) m β ( r ) = 1(2 π ) / Z d q ˆ m α ( q ) ˆ χ − α,β ( q ) ˆ m β ( − q ) . (S3) a r X i v : . [ phy s i c s . c h e m - ph ] F e b We then write the functional of bulk energy U bulk in Fourier space: U bulk [ ˆm , ˆD ] = 12 (cid:15) π ) / Z d q h ˆD ( q ) ˆD ( − q ) − ˆm ( − q ) ˆD ( q ) + ˆ m α ( q ) ˆ χ − α,β ( q ) ˆ m β ( − q ) i . (S4)The mean field polarization, ˆm MF , is defined as the minimum of the functional energy: δ U bulk δ ˆm [ ˆm MF ] = 0 . (S5)By performing the functional derivative of the equation S4, we show:ˆ D ( q ) = ˆ χ − ( q ) ˆ m MF ( q ) ⇔ ˆ m MF ( q ) = ˆ χ ( q ) ˆ D ( q ) (S6)Since the dielectric material is homogeneous and isotropic, the dielectric susceptibility tensor ˆ χ − has two compo-nents : ˆ χ − α,β ( q ) = ˆ χ − k ( q ) q α q β q + ˆ χ − ⊥ ( q ) (cid:18) δ α,β − q α q β q (cid:19) , (S7)where ˆ χ − k is the longitudinal component and ˆ χ − ⊥ , the transverse one. It can be seen that developing the longitudinalcomponent to the fourth order, ˆ χ − k ( q ) = 1 + K + K l q + βq (S8)ˆ χ − ⊥ ( q ) = K, (S9)is equivalent to the Landau Ginzburg approach used in the main article: U bulk [ m , D ] = 12 (cid:15) Z d r ( D − m ( r )) + 12 (cid:15) Z d r (cid:16) K m + K l ( ∇ · m ) + β ( ∇ ( ∇ · m )) (cid:17) . (S10) Polarization for one-wall geometry system . . . . z [nm] m ( z ) Water I n t e r f a ce FIG. S1. System with water close to a wall. The polarization field is schematically shown with the blue curve (computed fromEq. S11).
As stated in the main text, there is a polarization field m ( z ) that minimizes the bulk part of the functional of theenergy U bulk defined in equation 8 : m ( z ) = D K (cid:18) − e − zλd (cid:18) cos( q o z ) + q d q o sin( q o z ) (cid:19)(cid:19) + ¯ me − zλd (cid:18) cos( q o z ) + q d q o sin( q o z ) (cid:19) − λ o ¯ ρe − zλd sin( q o z ) . (S11)with m (0) = ¯ m and ˙ m (0) = − ¯ ρ . We introduce this expression of polarization into the total functional of energy U = U bulk + U s = 12 (cid:15) Z + ∞ z =0 dz ( D − m ( z )) + 12 (cid:15) Z + ∞ z =0 dzKm ( z ) + K l ( ˙ m ( z )) + β ( ¨ m ( z )) + k m m − m ) + k ρ ρ − ρ ) , (S12)where dots stand for d/dz . Performing the integral over z in U bulk , we express U as a function of the boundaryconstants ( ¯ m, ¯ ρ ): U ( ¯ m, ¯ ρ ) = k m m − m ) + k ρ ρ − ρ ) + q d D ( q d + q o )(1 + K ) (cid:15) + D ( q d + q o ) (cid:15) ¯ ρ + q d (1 + K )( q d + q o ) (cid:15) ¯ ρ − q d D ( q d + q o ) (cid:15) ¯ m + q d (1 + K )( q d + q o ) (cid:15) ¯ m − K ( q d + q o ) (cid:15) ¯ ρ ¯ m. (S13)We find the stationary point of the energy U ( ¯ m, ¯ ρ ) by cancelling its gradient: ∂U∂ ¯ m ( m s , ρ s ) = 0 ∂U∂ ¯ ρ ( m s , ρ s ) = 0 ⇒ − (cid:0) q d + q o (cid:1) q d (cid:18) m s − D K (cid:19) + ρ s + ˜ k ρ ( ρ s − ρ ) = 02 q d (cid:18) m s − D K (cid:19) − ρ s − q d − q o q d ˜ k m ( m s − m ) = 0 (S14)where ˜ k m = k m /k m, bulk and ˜ k ρ = k ρ /k ρ, bulk are dimensionless stiffness constants with k m, bulk = − (cid:0) q d − q o (cid:1) (1 + K )2 q d ( q d + q o ) (cid:15) (S15) k ρ, bulk = 2 q d (1 + K )( q d + q o ) (cid:15) . (S16)The expression of the stationary point ( m s , ρ s ) is linear function of the external field D and be written as ρ s ( D ) = a ρ D K + b ρ (S17) m s ( D ) = a m D K + b m , (S18)with a ρ = − (cid:0) q d − q o (cid:1) (cid:0) q d + q o (cid:1) ˜ k m q d (cid:16) (3 q d − q o ) (cid:16) ˜ k m (cid:16) k ρ (cid:17) − (cid:17) − q d ˜ k ρ (cid:17) (S19) a m = − (cid:0) q d − q o (cid:1) + 4 q d ˜ k ρ (3 q d − q o ) (cid:16) ˜ k m (cid:16) k ρ (cid:17) − (cid:17) − q d ˜ k ρ (S20)and b ρ = (cid:0) q d − q o (cid:1) ˜ k m (cid:16)(cid:0) q d + q o (cid:1) m + 2 q d ρ ˜ k ρ (cid:17) − q d ρ ˜ k ρ q d (cid:16) (3 q d − q o ) (cid:16) ˜ k m (cid:16) k ρ (cid:17) − (cid:17) − q d ˜ k ρ (cid:17) (S21) b m = (cid:0) q d − q o (cid:1) m ˜ k m (cid:16) k ρ (cid:17) − q d ρ ˜ k ρ (3 q d − q o ) (cid:16) ˜ k m (cid:16) k ρ (cid:17) − (cid:17) − q d ˜ k ρ . (S22)The stationary point ( m s , ρ s ) is not always a minimum as evidenced by the determinant Det( H U w ) of the Hessianmatrix, H = ∂ U bulk ∂ ¯ m∂ ¯ m ∂ U bulk ∂ ¯ ρ∂ ¯ m∂ U bulk ∂ ¯ ρ∂ ¯ m ∂ U bulk ∂ ¯ ρ∂ ¯ ρ ! . (S23)We find that Det( H )(˜ k m , ˜ k ρ ) = ( q o − q d )(1 + K ) (cid:15) ( q d + q o ) (cid:18) ˜ k m (cid:16) k ρ (cid:17) − q d q d − q o ˜ k ρ − (cid:19) (S24)can be negative for given values of (˜ k m , ˜ k ρ ). Under these conditions, the functional energy U ( ¯ m, ¯ ρ ) extremum is asaddle point. This saddle point is actually shown in figure S2a for a non-interaction surface (˜ k m = 0 , ˜ k ρ = 0) and avanishing value of D . The polarization and the layering of the fluid by a non-interacting interface, i.e. by a puregeometrical constraint, is a consequence of the pronounced overscreening associated with specific wavelengths. Thisnon-physical limit could be avoided by adding a non-linear term in m in Eq. (S10) that would ensure a saturationeffect in the fluid. In Figure 2a of the main article and Figure S2b, the gray curve separating the region for which theextremum of the energy functional U ( ¯ m, ¯ ρ ) is a minimum (gray dotted area) and that for which this extremum is asaddle point (white area) is computed by cancelling equation S24. − − m [V / nm] − − ¯ ρ / q o [ V / n m ] a ˜ k m = ˜ k ρ = 0 − − U / k B T [ n m − ] × . . . . k m ˜ k ρ MinimumSaddlepoint b FIG. S2. a . Function U ( ¯ m, ¯ ρ ) for (˜ k m = 0 , ˜ k ρ = 0) and a vanishing value of D . b . The dotted area represents a zone ofparameters ˜ k m and ˜ k ρ for which the function U ( ¯ m, ¯ ρ ) possesses a minimum for finite value of ¯ m and ¯ ρ . Outside this area thefunctional energy shows a saddle point as depicted in the left panel. Water confined in between two surfaces : a slab geometry -L/2 0 L/2 m ( z ) Water I n t e r f a c e I n t e r f a c e z [nm] FIG. S3. System with water confined between two walls separated by a length L . The blue curve is a schematic representationof the polarization highlighting the oddity of this function for vanishing D . We now consider a system composed of two identical surfaces, located in z = − L/ z = L/ U = U bulk + U s , (S25)the sum of the electrostatic energy of water U bulk and the interaction energy U s with U bulk [ m ] = 12 (cid:15) Z L/ − L/ ( m ( z ) − D ) dz + 12 (cid:15) Z L/ − L/ (cid:0) Km ( z ) + K l ˙ m ( z ) + β ¨ m ( z ) (cid:1) dz (S26) U s = k ρ ρ − − ρ ) + k ρ ρ + − ρ ) + k m m − − m ) + k m m + + m ) (S27)and ¯ ρ − = ρ ( z = − L/
2) ¯ ρ + = ρ ( z = + L/ m − = m ( z = − L/
2) ¯ m + = m ( z = + L/ . (S28)As one sees, the confinement induces an interaction term U s composed of two contributions, one in z = − L/ L/
2. Minimizing the total energy with respect to the polarization field still leads to a fourth order differentialequation : (1 + K ) m ( z ) − K l m (2) ( z ) + βm (4) ( z ) = D . (S29)The solutions of this differential equation are written in the following form: m ( z ) = D K + ¯ C cos ( q o z ) cosh ( q d z ) + ¯ C sin ( q o z ) cosh ( q d z )+ ¯ C cos ( q o z ) sinh ( q d z ) + ¯ C sin ( q o z ) sinh ( q d z ) . (S30)Written in this manner, we underline the fact that the polarization m ( z ) is odd (and the bound charge ρ ( z ) = − ˙ m ( z )is even) for a vanishing D .Following the approach presented in the previous section, the solution of the differential equation (Eq. S30)is introduced into the functional energy U . One then gets an energy U ( ¯ C , ¯ C , ¯ C , ¯ C ) which stationary point( C s , C s , C s , C s ) can be determined by cancelling its gradient. Results are given in the appendix. Note thatin the case where the external field D = 0, the problem is much simpler. Indeed, the system is then symmetricalaccording to z : m ( z ) = − m ( − z ) which leads to having ¯ C = ¯ C = 0. The energy U ( ¯ C , ¯ C ) is then a function of twovariables which nature of the extremum can be assessed by studying the sign of the determinant of the associatedHessian matrix. Figure S4 shows that as soon as the distance between the two walls is greater than ∼ L < k ρ ˜ k m − L [ n m ] FIG. S4. The curves are plotted for a null-Hessian matrix determinant of the functional energy and for D = 0. The color ofthe curves is related to the distance L between the two walls given on the colorbar. The black curve is drawn for a system withonly one wall. THE DETAILS OF MOLECULAR DYNAMICS SIMULATIONS
We performed molecular dynamics simulations using GROMACS 2019.6 package.[1] Dispersion interactions aremodelled using effective Lennard-Jones (LJ) potentials truncated at 1 . r ij is V LJ ( r ij ) = (cid:15) ij (cid:18)(cid:16) σ ij r ij (cid:17) − (cid:16) σ ij r ij (cid:17) (cid:19) for r ij ≤ . r ij > . , (S31)where (cid:15) ij is the depth of the potential well (or dispersive energy) and σ ij is the distance at which the particle-particlepotential energy is zero (or size of the particle). Only atom pairs involving an oxygen atom interact through theLennard-Jones potential. Associated parameters, (cid:15) X − O and σ X − O , are given in the following Table I. X Partial charge [e] σ X − O [nm] (cid:15) X − O [kJ / mol] SourceH 0.4238 0 0O -0.8476 0.3166 0.650C 0 0.319 0.392 Werder et al.[2]B 0.37 0.331 0.508 Won and Aluru[3]N -0.37 0.326 0.628 Won and Aluru[3]TABLE I. Partial charges and Lennard-Jones parameters for the interaction between the atoms of the wall X and the oxygenof the water molecules. The Coulomb force is treated using a real-space cutoff at 1 . τ = 0 . Periodic water box n m n m n m FIG. S5. The simulated system after the thermalization in NVT ensemble: a 5 nm side cube, periodic along the 3 directions ofspace and filled with 4055 SPC/E water molecules.
The simulated system is a 5 × × cube containing 4055 SPC/E water molecules and with periodicity followingthe 3 directions of space. A first simulation carried out in an NVT canonical ensemble and of short duration allowsto ensure the thermalization of the system. Then we let the system adjust the volume of the box at constant pressureequal to 1 bar with a Berendsen barostat [6] ( τ = 2 . ∼
998 kg / m )thanks to a second simulation performed in a NPT ensemble with a long range dispersion corrections for energy andpressure. At the end of this relaxation, the side of the box is slightly different and equal to 4 .
963 nm. Finally, wecarry out long-term simulations (2 ns) in a canonical ensemble from which we compute H-H, O-H and O-O structurefactors and deduce the susceptibility ˆ χ k ( q ) with the classical approximation for the fluctuation dissipation theorem(see black curves on Figure S8).[7]Figure S8 shows the two strategies to fit the dielectric susceptibility of bulk water with a 4 order model:ˆ χ k ( q ) = 11 + K + K l q + βq . (S32)The first one aims at fitting the maximum of the susceptibility while the second one focuses on fitting the width ofthe susceptibility normalized to 1 (the aspect ratio) bearing in mind that in both cases ˆ χ k ( q = 0) is set by the valueof the relative dielectric permittivity of the water SPC/E (cid:15) r = (1 − ˆ χ k ( q = 0)) − ∼ q [nm − ]01020304050 ˆ χ k ( q ) a SPCE waterMax fittingAspect ratio fitting
10 20 30 40 50 60 70 q [nm − ]0 . . . ˆ χ k ( q ) / ˆ χ k , m a x b FIG. S6. Susceptibility χ ( q ) ( a ) and normalized susceptibility χ ( q ) /χ max ( b ) computed from pure SPC/E water moleculardynamic simulation (black curve). The blue curve arises from the fit of the maximum of the first component of the susceptibilitywhile the orange curve shows the result of the fit of the aspect ratio. Water confined in slab geometry L . n m . n m FIG. S7. The simulated system after the thermalization in NVT ensemble: SPC/E water molecules confined between twographene sheets.
We performed molecular dynamics simulations of SPC/E water confined in a slab geometry. The walls, perpendic-ular to the z direction, are made up of atoms which positions are fixed and arranged in a hexagonal lattice. L is thedistance between the two walls. The simulation box is extended in the z-direction until it reaches a length equal to3 L . Thus, even if we considered a periodic system in the 3 directions of space, the periodized slab system is separatedby a 2 L -thick void layer along the z axis.In our study we considered two types of hydrophobic walls: graphene and a hBN wall. The dispersion interactionsbetween the wall and the water molecules are modelled using the truncated Lennard-Jones potential between theatoms of the wall and the oxygen atoms of the water molecules. We recall that the parameters of this potential aregiven in the Table I. The walls are electrically neutral overall but, in the case of hBN, the atoms carry partial chargesleading to an electrostatic interaction with the water molecules.As for the previous system, we conducted a first simulation in a NVT canonical ensemble and of short duration thatallows to ensure the thermalization of the system. Then, we check that the number N of water molecules introducedbetween the two walls, separated by a distance L , leads to a density close to bulk water ( ∼ / m ) with a secondsimulation performed in a NPT ensemble with a semi-isotropic Parrinello-Rahman barostat [8] ( τ = 2 . / nm. Besides, we have computed simulations with two separated graphene walls with a length L ranging from1 to 10 nm with E ext = 0 and 0 . / nm.From the molecular trajectories, we compute the charge density ρ ( z ) as a histogram of the point charge distributionwith a step of 0 .
01 nm along z . Then, the polarization field m ( z ) calculate by performing a numerical derivativeof the charge density, m ( z ) = − dρ ( z ) /dz . Finally, the susceptibility is computed from the polarization field of twosimulations in which the amplitude of the external field is different by a value δD = 0 . / nm : χ ( z, D ) = m ( z, D + δD ) − m ( z, D ) δD (S33) Validation of the fitting strategy of the bulk susceptibility ˆ χ k ( q ) . . . . . z [nm]10 − − − − l og | m ( z ) | SimulationMax fittingAspect ratio fitting
FIG. S8. The black curve is the polarisation field computed from MD simulation of water in the vicinity of graphene layer atvanishing D . Color curves have the equation m ( z ) = Ae − q d z where q d is computed from the fitting of the bulk susceptibilityˆ χ ( q ). The blue one arises from the fit of the maximum of the first component of the susceptibility while the orange curve showsthe result of the fit of the aspect ratio. The parameter A is chosen to match the back curve. Electrostatic field for large slab geometry
The simulations for which the results are shown in figures S9 were performed with a 5 nm distance between the twographene walls. This value is much larger than the characteristic decay length λ d . Therefore, the water structurationin the vicinity of one wall is not affected by the presence of the second wall. We then choose to display the results forthe left wall only.Figure S9 shows the polarization and susceptibility in such a system for increasing values of the external field D up to 2 V / nm. Considering the variations in amplitude of the susceptibility for interfacial water, one can estimatethat the regime is linear up to D ∼ / nm. This value is however lower than the one reported in previous work.[9](HB: Cette valeur a aussi ´et´e estim´ee dans les Si de ... ???)0 . . . − − − m [ V / n m ] a Graphene . . . − − − b hBN . . . z + L [nm]0 . . . . . χ ( z ) c . . . z + L [nm]0 . . . . . d D [ V / n m ] FIG. S9. Polarization m ( z ) (top panel) and susceptibility χ ( z ) (bottom panel) from molecular dynamics simulations of waternear a graphene (left panel) or hBN layer (right panel) for different external field strengths (see color code). The part of thecurves for z < z are shown as dotted lines. Smearing the surface
The Figure S10 represents the density of Oxygen ρ O and Hydrogen ρ H in the vicinity of graphene (green curve) andhBN (red curve) layers. The surface induces a layering of the fluid that is more pronounced for hBN as for graphenewhich is indicated by a higher and narrower first hydration peak for this surface. In order to take into account thatthe numerical interface is not 2D but is associated with a small width, we convolute the theoretical expressions ofpolarization and susceptibility with a Gaussian G ( z − z ) = e − z / η η √ π . (S34)The position z and the width η of the Gaussian are determined for each surface by fitting the first peak of the Oxygendensity with G ( z − z ). The inset of Fig S11 shows the Gaussian G (full red zone) and the simulated density foroxygen in the vicinity of graphene. The parameters η and z do not vary with the application of an excitation field D as shown in Fig S11.00 . . . ρ O [ n m − ] a GraphenehBN0 . . . z [nm]050100150 ρ H [ n m − ] b FIG. S10. Oxygen ρ O ( z ) (top panel) and Hydrogen ρ H ( z ) density near a graphene or a hBN layer. . . . . . . . . . . z [ n m ] a . .
5z [nm]050 ρ O [ n m − ] . . . . . D [V/nm]0 . . . . . . . η [ n m ] b GraphenehBN
FIG. S11. Smearing parameters z (a) and the associated standard deviation η (b) as a function of the external field and fordifferent types of walls: graphene (blue), hBN (orange) and hBN without partial charges (green). The inset illustrates theGaussian fit of the first oxygen layer from which the parameters z and η are derived. Molecular dynamics simulations as function of the slab width L The theoretical model that we have developed can be used to determine the profile of the polarization and thesusceptibility of a nanometric water layer of varying thickness between graphene surfaces. The Figure S12 shows m ( z )(a) and χ ( z ) for water thickness going from 5 nm to 1 nm. The theoretical model allows to retrieve the organizationin successive layers of over ( χ ( z ) (cid:29) χ b ) and underresponding ( χ ( z ) (cid:28) χ b ) region in interfacial water.10 . . . . z + L [nm] − m ( z ) [ V / n m ] a Model 0 . . . . z + L [nm]02468101214 χ b FIG. S12. Polarization m ( z ) (a) and susceptibility χ ( z ) (b) from molecular dynamics simulations of water in a slab geometrywith graphene on both sides of the system and without external field. The width L of the geometry is given on the right sideof the figure. The colored curves, coming from MD simulations, are shown as dashed line for z < − L/ z and z > L/ − z .Curves have been shifted for sake of clarity. Physical origin of a non-zero ˜ k ρ In this section, we probe the influence of partial charge q and dispersion energy (cid:15) of the surface’s atoms on thedielectric properties of interfacial water. We consider artificial surfaces for which one of these parameters has beenmodified with respect to graphene or hBN. Impact of partial charges
First, we study the influence of partial charge of surface’s atoms by considering a ’charged’graphene as well as ’non-charged’ hBN. The parameters of these two surfaces are given in tables S13a and S13b. Fromthose custom MD simulation, we show that the dielectric susceptibility of interfacial water is the same for ’charged’graphene and graphene. It is also the same for ’non-charged’ hBN and hBN. We conclude that it is not the partialcharge of the surface’s atoms that trigger the differences of the dielectric susceptibility of graphene and hBN interfacialwater.2ABAB σ (cid:15) q G r a ph e n e a ABAB σ (cid:15) q h B N b . . . z [nm]0 . . . . . χ ( z ) . . . z [nm]0 . . . . . χ ( z ) FIG. S13. Dielectric susceptibility of water in the vicinity of graphene (a) or hBN (b) layer without and with partial charges q . σ is given in nm, (cid:15) in kJ / mol and the partial charge q in eV. The curves are made transparent for z < .
33 nm.
Impact of heterogeneity
Secondly, we evaluate the impact of the bi-atomic nature of hBN on the susceptibility. Todo so, we consider a ’bi-atomic’ like graphene surface associated with LJ parameters given in Table S14a. The twoartificial atoms have the same size σ as carbon’s in graphene and do not carry partial charges. However, they areassociated with two different (cid:15) which average is equal to the carbon’s dispersion energy in graphene. In contrast, wesimulate a ’mono-atomic’ like hBN surface (see Table S14b). This surface is composed of two atoms carrying oppositepartial charges equal to ± .
37 e. These two atoms are associated with the same LJ parameter set. The size and thedispersion energy are equal to the average of those for Boron and Nitride. From those custom MD simulations, weshow that the dielectric susceptibility of interfacial water is the same for ’bi-atomic’ like graphene and graphene aswell as ’mono-atomic’ like hBN and hBN. This aspect has no influence on the dielectric susceptibility of interfacialwater. ABAB σ (cid:15) q G r a ph e n e a ABAB σ (cid:15) q h B N b . . . z [nm]0 . . . . . χ ( z ) . . . z [nm]0 . . . . . χ ( z ) FIG. S14. Focus on the impact of a bi-atomic wall on the dielectric susceptibility of water. The orange curves are derived froma mono-atomic wall while the black dashed curves result from a bi-atomic wall. σ is given in nm, (cid:15) in kJ / mol and the partialcharge q in eV. The curves are made transparent for z < .
33 nm.
Average dispersion energy influence
Finally, we increase the dispersion energy of carbon atoms in the graphenesurface (see Table S15a). In contrast, we decrease the average dispersion energy for hBN by keeping the ratiobetween the two atoms (see Table S15b). We can see that the greater the average dispersion energy, the smaller thesusceptibility amplitude.3ABAB σ (cid:15) q G r a ph e n e a ABAB σ (cid:15) q h B N b . . . z [nm]0 . . . . . χ ( z ) . . . z [nm]0 . . . . . χ ( z ) FIG. S15. Dielectric susceptibility of water in the vicinity of a surface when the Lennard-Jones interaction between the waterand this surface is deeper (from orange curves to black dashed curves). σ is given in nm, (cid:15) in kJ / mol and the partial charge q in eV. The curves are made transparent for z < .
33 nm.
The susceptibility of interfacial water for graphene surface shows a higher amplitude than the one for hBN surface.Thanks to this investigation, we conclude that this difference is triggered by the average dispersion energy which issmaller for graphene.
Appendix : Stationary point in the case of a slab geometry C s = (cid:18) m − D K (cid:19) a b + c (S35) C s = ρ a b − c (S36) C s = ρ a b − c (S37) C s = (cid:18) m − D K (cid:19) a b + c (S38)4with a = 2 k m (cid:15) (cid:0) q d + q o (cid:1) (cid:18) sinh (cid:18) Lq d (cid:19) (cid:18) (cid:15) k ρ q o (cid:0) q d + q o (cid:1) cos (cid:18) Lq o (cid:19) + ( K + 1) ( q d − q o ) ( q d + q o ) sin (cid:18) Lq o (cid:19)(cid:19) + cosh (cid:18) Lq d (cid:19) (cid:18) (cid:15) k ρ q d (cid:0) q d + q o (cid:1) sin (cid:18) Lq o (cid:19) + 2( K + 1) q o q d cos (cid:18) Lq o (cid:19)(cid:19)(cid:19) (S39) a = 2 (cid:15) k ρ (cid:0) q d + q o (cid:1) (cid:18) sinh (cid:18) Lq d (cid:19) (cid:18) (cid:15) k m (cid:0) q d + q o (cid:1) cos (cid:18) Lq o (cid:19) + ( K + 1) q o sin (cid:18) Lq o (cid:19)(cid:19) + ( K + 1) q d cosh (cid:18) Lq d (cid:19) cos (cid:18) Lq o (cid:19)(cid:19) (S40) a = − (cid:15) k ρ (cid:0) q d + q o (cid:1) (cid:18) cosh (cid:18) Lq d (cid:19) (cid:18) (cid:15) k m (cid:0) q d + q o (cid:1) sin (cid:18) Lq o (cid:19) − ( K + 1) q o cos (cid:18) Lq o (cid:19)(cid:19) + ( K + 1) q d sinh (cid:18) Lq d (cid:19) sin (cid:18) Lq o (cid:19)(cid:19) (S41) a = − k m (cid:15) (cid:0) q d + q o (cid:1) (cid:18) sinh (cid:18) Lq d (cid:19) (cid:18) (cid:15) k ρ q d (cid:0) q d + q o (cid:1) cos (cid:18) Lq o (cid:19) − K + 1) q o q d sin (cid:18) Lq o (cid:19)(cid:19) + cosh (cid:18) Lq d (cid:19) (cid:18) − (cid:15) k ρ q o (cid:0) q d + q o (cid:1) sin (cid:18) Lq o (cid:19) + ( K + 1) ( q d − q o ) ( q d + q o ) cos (cid:18) Lq o (cid:19)(cid:19)(cid:19) (S42)and b = − K + 1) (cid:15) q d q o (cid:0) q d + q o (cid:1) cos ( Lq o ) (cid:0) k ρ (cid:0) q d + q o (cid:1) − k m (cid:1) + q d sin ( Lq o ) (cid:0) (cid:15) k m k ρ (cid:0) q d + q o (cid:1) − ( K + 1) (cid:0) q d − q o (cid:1)(cid:1) (S43) c = 2( K + 1) (cid:15) q d q o (cid:0) q d + q o (cid:1) cosh ( Lq d ) (cid:0) k ρ (cid:0) q d + q o (cid:1) + k m (cid:1) q o sinh ( Lq d ) (cid:0) (cid:15) k m k ρ (cid:0) q d + q o (cid:1) + ( K + 1) (cid:0) q d − q o (cid:1)(cid:1) (S44) [1] Mark James Abraham, Teemu Murtola, Roland Schulz, Szil´ard P´all, Jeremy C. Smith, Berk Hess, and Erik Lindahl. GRO-MACS: High performance molecular simulations through multi-level parallelism from laptops to supercomputers. SoftwareX ,1-2:19–25, September 2015.[2] T. Werder, J. H. Walther, R. L. Jaffe, T. Halicioglu, and P. Koumoutsakos. On the water-carbon interaction for use inmolecular dynamics simulations of graphite and carbon nanotubes.
The Journal of Physical Chemistry B , 107(6):1345–1352,February 2003.[3] Chang Y. Won and N. R. Aluru. Water permeation through a subnanometer boron nitride nanotube.
Journal of theAmerican Chemical Society , 129(10):2748–2749, March 2007.[4] In-Chul Yeh and Max L. Berkowitz. Ewald summation for systems with slab geometry.
The Journal of Chemical Physics ,111(7):3155–3162, August 1999.[5] Giovanni Bussi, Davide Donadio, and Michele Parrinello. Canonical sampling through velocity rescaling.
The Journal ofChemical Physics , 126(1):014101, January 2007.[6] H. J. C. Berendsen, J. P. M. Postma, W. F. van Gunsteren, A. DiNola, and J. R. Haak. Molecular dynamics with couplingto an external bath.
The Journal of Chemical Physics , 81(8):3684–3690, October 1984.[7] Philippe A. Bopp, Alexei A. Kornyshev, and Godehard Sutmann. Static nonlocal dielectric function of liquid water.
PhysicalReview Letters , 76(8):1280–1283, February 1996.[8] Giovanni Bussi, Tatyana Zykova-Timan, and Michele Parrinello. Isothermal-isobaric molecular dynamics using stochasticvelocity rescaling.
The Journal of Chemical Physics , 130(7):074101, February 2009.[9] Douwe Jan Bonthuis, Stephan Gekle, and Roland R. Netz. Profile of the static permittivity tensor of water at interfaces:Consequences for capacitance, hydration interaction and ion adsorption.