Deconvolution of BIL-SFG and DL-SFG spectroscopic signals reveals order/disorder of water at the elusive aqueous silica interface

Simone Pezzotti *ab, Daria Ruth Galimberti ab and Marie-Pierre Gaigeot *ab
aLAMBE CNRS UMR8587, Laboratoire Analyse et Modélisation pour la Biologie et l'Environnement, Université d'Evry val d'Essonne, Blvd F. Mitterrand, Bat Maupertuis, 91025 Evry, France. E-mail:;
bUniversité Paris-Saclay, France

Received 15th May 2019 , Accepted 8th August 2019

First published on 8th August 2019

Through the prism of the rather controversial and elusive silica/water interface, ab initio DFT-based molecular dynamics simulations of the structure and non-linear SFG spectroscopy of the interface are analysed. Following our recent work [Phys. Chem. Chem. Phys., 2018, 20, 5190–5199], we show that once the interfacial water is decomposed into BIL (Binding Interfacial Layer) and DL (Diffuse Layer) interfacial regions, the SFG signals can be deconvolved and unambiguously interpreted, and a global microscopic understanding on silica/water interfaces can be obtained. By comparing crystalline quartz/water and amorphous (fused) silica/water interfaces, the dependence of interfacial structural and spectroscopic properties on the degree of surface crystallinity is established, while by adding KCl electrolytes at the quartz/water interface, the chaotropic effect of ions on the interfacial molecular arrangement is unveiled. The evolution of structure and SFG spectra of silica/water interfaces with respect to increasing surface deprotonation, i.e., with respect to pH conditions, is also evaluated. Spectroscopic BIL-SFG markers that experimentally allow one detect the water order/disorder in the BIL as a function of surface hydroxylation and ion concentration are revealed, while the pH-induced modulations in the experimentally recorded SFG spectra are rationalized in terms of changes in both BIL and DL SFG signatures.

1 Introduction

Because silica is one of the most abundant minerals on Earth, its interaction with water is of pivotal importance in many aspects of fundamental and applied sciences. Geochemistry and environmental chemistry are certainly two domains in high demand of fundamental knowledge for understanding mineral/water structural organisation and dynamics, in particular because of the underground water navigating the Earth crust, interacting with the minerals and transporting organic molecules and pollutants. Silica dissolution is also of utmost importance for the environment, e.g., in the context of nuclear waste storage. Improving the knowledge of oxides at the interface with aqueous electrolytes is of high relevance for designing new materials, for instance in the context of heterogeneous catalysis, photovoltaics, batteries, while nano-/meso-porous silica nanoparticles are relevant for molecular drug encapsulation and more generally for biotechnological applications.

The structure of the silica/water interface is probably one of the most controversial and elusive of the oxide/water family. Vibrational Sum Frequency Generation (SFG) non-linear spectroscopy has certainly been the experimental technique the most in use within the past two decades in order to unravel such structure.1–11 Because SFG is a probe of buried interfaces, probing the non-centrosymmetric part of these media, it is one of the most powerful techniques to interrogate mineral/water structures within a limited number of water layers over a few Å thickness. While the experimental SFG spectra recorded for silica/water interfaces are similar between the groups, there has always been disagreement in their microscopic interpretation. This can possibly be due to experiments performed on crystal vs. fused (amorphous) silica, with non-equivalent water interfacial structures (see for instance ref. 12–14), or due to the use of conventional SFG spectroscopy without phase knowledge which is more complex to interpret than phase-resolved SFG. Coming recently into light, interferences between the emitted photons in the SFG experiment15,16 also have to be taken into account in the final experimental spectra, possibly shifting band positions and reshaping the bands, which complicates further comparison and interpretation of spectra published in the literature.

In the H-bonded spectral region, four main scenario have been presented for the microscopic assignment of the 3200 and 3400 cm−1 SFG bands at the silica/water interface (conventional |χ2(ω)|2 and phase-resolved Imχ2(ω)). The earlier interpretation6 of liquid-like (3400 cm−1) and ice-like (3200 cm−1) bands has been replaced by ‘two water populations at the interface’, being either described as disordered/ordered water structures6 or as water donors/acceptors of H-bonds with surface silanols12 (going in line with the two-pKa silanol surface sites initially inferred from Eisenthal's pioneering SHG experiments12,17). Fermi resonance, and more generally, collective coupled vibrational modes between water molecules have also been used for the interpretation of the 3200 cm−1 band.7 A scenario based on water molecules close to the surface/located further away from the surface responsible for the 3400/3200 cm−1 bands has also been proposed.18 Richmond further suggested that there exits an additional interfacial water population also contributing to the two bands.19 Note also that the origin of the free O–H peak at ∼3700 cm−1 normally present at hydrophobic silica/water interfaces and suppressed at hydrophilic interfaces (e.g., quartz) is still a source of controversy.8,9,20,21 On the other hand, molecular dynamics (MD) simulations12,22–24 of silica/water interfaces have not yet provided a sufficiently clear picture of how water is organized at silica interfaces, and they have not provided an unambiguous structure-SFG fingerprints relationship to date.

The latter is not straightforward to establish from experiments alone while the MD results have always suffered from ambiguous definitions of the water layers at the interface. Most of the interpretations present in the literature (experiments and simulations,25–28 including ours for quartz/water for instance12) have mostly focused on the water molecules directly in contact with the surface. In a pioneering work by Tian et al.,29 it became clear that such interpretation is insufficient and possibly misleading, and that two water layers are however universally responsible for the SFG signals, i.e., the BIL (Binding Interfacial Layer) and the DL (Diffuse Layer). Beyond these two layers, whose total thickness depends on the surface chemical nature and pH/electrolytic conditions, bulk liquid water is recovered without SFG activity. In our recent work,30 we have demonstrated how to define these water layers from MD simulations, solely based on the structural properties of water (using three structural descriptors: density, H-bond network, orientation), and we have furthermore demonstrated that only these two layers are responsible for SFG activity. In ref. 30 and in a subsequent paper,31 we have shown that the DL-SFG signal is universal and systematically provides a 3200–3400 cm−1 two-band structure, whose amplitude and sign are surface charge dependent, and that the BIL-SFG signal is the only one that carries relevant structural information on the organisation of water at charged interfaces. Although the DL-SFG signal does not carry any direct information on the water structure at the direct interface, it does carry a wealth of complementary information such as the surface charge, the isoelectric point, the EDL (electric double layer) formation, and the direct relationship between aqueous electrolyte solution pH and surface hydroxylation state.31 The DL-SFG knowledge has also been shown to be crucial to reveal the proton location at the acidic air–water interface.32 Other authors have also attempted to distinguish the double layer conformation at aqueous interfaces, see e.g., ref. 33 and 34.

In the present work, we apply the same strategies for defining the BIL and DL SFG active water layers at the interface with silica surfaces, for calculating SFG nonlinear spectroscopic signatures of these layers, and for extracting spectroscopic markers specific to the BIL and DL, and hence provide an unambiguous interpretation of SFG features.

DFT-based molecular dynamics simulations (DFT-MD/AIMD for Ab Initio MD) are performed on three silica/water interfaces, chosen because they represent three limiting cases: the crystalline neat fully hydroxylated (0001) α-quartz/water interface under PZC conditions (denoted QW), the fully hydroxylated amorphous silica/water interface (denoted AW) under PZC conditions, and the fully hydroxylated (0001) α-quartz surface interfaced with 2 M KCl aqueous solution (denoted QW + KCl). The latter two cases represent well defined perturbations introduced upon the water interfacial ordering observed at the neat QW interface (see ref. 12 and 14), either by decreasing the solid surface regularity (when we deal with the AW interface) or by the adsorption of electrolytes which induces a surface electric field35 (when we look at the QW + KCl interface). We furthermore made the choice to consider a highly hydroxylated amorphous silica surface with 7.6 SiOH per nm2, whose hydroxylation state is close to the value of 9.6 SiOH per nm2 at the crystalline (0001) α-quartz surface. This reflects our goal to investigate the effect of one parameter at a time, i.e., the surface morphology (crystalline vs. amorphous silica), while the surface hydroxylation state is kept more or less constant. As regards the QW + KCl interface, K+ and Cl ions are known to be enriched at the interface with cations like K+ and Na+ being inner-sphere adsorbed,24,36–38 hence directly altering the water HB-network that was initially existing in the BIL (in the absence of the electrolytes). The rather high 2 M ionic concentration considered in our work allows an easier detection of the effect of electrolytes on interfacial structural and spectroscopic properties when they are located in the BIL.

The interfacial structures of these three systems have already been presented in some of our previous works;12,13,35 they are hereby discussed in terms of BIL and DL interfacial water layers and associated SFG spectroscopy.

The details of the DFT-MD simulations and the structural definitions of BIL-water and DL-water are given in the next two sections, and their application to the three chosen silica/water interfaces is presented next. The SFG calculations including the instrumental role of χ(3)bulk(ω) third order nonlinear susceptibility of bulk water (see also ref. 30) as well as the photon interference effects are discussed. Unambiguous analyses of SFG spectra and of spectroscopic marker bands are revealed for the three chosen silica/water interfaces. The separation of the interfacial water into SFG-active BIL and DL layers is instrumental in the spectroscopic interpretation, as one can calculate, deconvolve and finally interpret the SFG features arising separately from each of these structurally distinct waters.

2 Methods and computational details

Generalities on the DFT-MD simulations

Three silica–water interfaces are investigated in the present work: the neat (0001) α-quartz/water interface (denoted QW), the amorphous/water interface (denoted AW), and 2 M KCl electrolyte at the (0001) α-quartz/water interface (QW + KCl). See below for the details of the simulation boxes and surface charge conditions. State-of-the-art first principles molecular dynamics simulations in the framework of the electronic DFT representation (DFT-MD) are carried out with the CP2K package,39,40 adopting the BLYP41,42 functional, augmented by D2 Grimme corrections,43,44 in combination with the mixed Gaussian/Plane wave (PW) basis sets and GTH pseudopotentials.45 A cutoff of 280 Ry for the PW combined with the TZV2P Gaussian basis set for O, H, Si, and Cl, and DZVP-MOLOPT-SR for K+, is adopted as a good compromise between accuracy and computational cost in reproducing structural and vibrational properties (see ref. 12, 35 and 46).

Born–Oppenheimer molecular dynamics simulations are performed with a time-step of 0.4 fs, and the average temperature is in the range of 315–320 K. An equilibration dynamics is firstly performed for a duration of 5 ps (NVE ensemble, with possible rescaling of velocities) followed by a 15–100 ps NVE production run on which analyses are done. We investigate the dependence of the SFG spectroscopy signal on the length of the DFT-MD trajectory, therefore the 100 ps trajectories mentioned above. Note that (for reasons that will become clear when we analyse these 100 ps trajectories at the end of this section) the 100 ps trajectories have been accumulated for the QW and QW + KCl interfaces. Anticipating our results for the time-dependence, 15 ps trajectories are shown to be long enough to extract SFG signals.

The (0001) α-quartz and amorphous surfaces are considered as fully hydroxylated in our simulations. According to ref. 47–49, more than 75% of quartz surface sites are in the SiOH neutral form throughout the pH range 0–10, and slightly more than 15% of the surface sites are in the SiO deprotonated form at pH 7. Our present simulations are mostly concerned with PZC conditions (pH ∼2–4 for silica), where the surface does not have a net charge as a result of the surface terminations being mostly in the SiOH protonation state, while the less probable SiO and SiOH2+ sites balance each other. The simulated fully hydroxylated silica surface (where all silanols are in the SiOH protonation state) is hence more than a reasonable representation. As demonstrated in ref. 36, the neat (0001) α-quartz/water interface (denoted QW) is however not at the isoelectric point: at the interface with liquid water, the SiOH silanols acquire negative charges36 (surface charge slightly less than −0.1 C m−2),31 which creates an electric field at the surface, as will be shown again in Section 3. The amorphous/water interface (denoted AW, the amorphous model is from Bernasconi et al.50) is at the isoelectric point, as will be shown in Section 3. The quartz/water/2 M KCl interface (denoted QW + KCl) is a positively charged interface because of the inner-sphere adsorption of the cation (see Section 4 for the structural details and also ref. 35 and 36). The simulated box dimensions are 9.820 × 8.504 × 32.165 Å3 for QW and QW + KCl, and 9.117 × 16.342 × 32.000 Å3 for AW. There are 8 Si–OH at the surface of the (0001) α-quartz/water in the simulation box, and 11 at the surface of the amorphous/water, which corresponds to silanol densities of respectively 9.6 and 7.6 SiOH per nm2.

Instrumental in our work is the calculation of the instantaneous surface following the definitions of Willard and Chandler,51 at each time step of the dynamics, hence providing the (average) evolution of the water density with respect to the vertical distance from the (instantaneous) surface. Once the water density profile is obtained, water layers above the surface can be defined, and in each water layer the structural properties adopted by the water molecules (H-bonds, orientations, average density) can be unraveled. These individual structural properties are used as descriptors in order to define BIL (Binding Interfacial Layer), DL (Diffuse Layer) and bulk liquid water. The methodology is reviewed in the next section, and it has been presented in ref. 30, and has been successfully applied in ref. 31 and 32. For the characterization of hydrogen bonds, the definitions by Galli et al.52 have been adopted, with O(–H)⋯O ≤ 3.2 Å and the O–H⋯O angle in the range of 140–220°. 3D-plots of the joint probability for a given H-bond orientation and H-bond distance are calculated, hence providing the full characterization of H-bond networks in between water molecules. The 3D-plot for pure bulk liquid water has a homogeneous distribution of orientations for H-bonds 2.65–2.80 Å in length. In these 3D-plots, the orientation of the O–H group donating the H-bond is given by the cosine of the angle between the vector connecting O to H and the normal to the surface (θ, see Fig. 3 for the direction of [n with combining right harpoon above (vector)]). For spectroscopic calculations, the vector normal [n with combining right harpoon above (vector)] to the solid/liquid surface is taken as oriented towards the solid, and defines the z direction (lateral x, y directions are within the oxide surface).

Theoretical SFG spectroscopy

The SFG signal, coming from the total resonant electric dipole non-linear susceptibility χ(2)(ω) (real and imaginary components), is calculated following the time-dependent approach of Morita et al.53,54 For the water contribution to χ(2)(ω), individual molecular dipole moments and polarisability tensors are calculated with the model proposed by Khatib et al.,55 which supposes that in the high frequency region (>3000 cm−1) only the OH stretching motions are contributing to the spectrum. By neglecting the intermolecular terms, at it has been shown reasonable in ref. 56 through the truncation scheme, and as it is done in the current literature for theoretical SFG, χ(2)(ω) becomes:
image file: c9cp02766a-t1.tif(1)
where (P, Q, R) are any x,y,z direction in the laboratory frame, and kb and T are the Boltzmann constant and temperature of the simulated system. 〈⋯〉 is a time-correlation function, αPQ(t) and μR(t) are respectively the individual O–H bond contribution to the total polarization and dipole moment of the system and [small alpha, Greek, dot above]PQ(t) and [small mu, Greek, dot above]R(t) are their time derivatives. M is the number of water molecules and n1 and n2 are the two indices that identify each of the two O–H oscillators per molecule. Using the direction cosine matrix (D) projecting the molecular frame (x, y, z) onto the fixed laboratory frame (P, Q, R), and assuming that the O–H stretching is much faster than the modes involving a bond reorientation, one can write
image file: c9cp02766a-t2.tif(2)
image file: c9cp02766a-t3.tif(3)
The D matrix and the projection of the velocities on the O–H bond axis (vz) can be readily obtained from the DFT-MD trajectory, while image file: c9cp02766a-t4.tif and image file: c9cp02766a-t5.tif are parameterized.55,57 In this work, ssp SFG signals are calculated, i.e., P = x, Q = x and R = z in eqn (1)–(3). Eqn (2) and (3), now including velocities along the bond axis, follow a similar strategy as in ref. 58.

As demonstrated in ref. 30–32, 46 and 59, the interfacial region at charged interfaces is composed of BIL (Binding Interfacial Layer) and DL (Diffuse Layer) layers, both being SFG-active, due to the breaking of centrosymmetry induced by specific water–surface interaction and field-induced water reorientation respectively, so that

χ(2)(ω) = χ(2)BIL(ω) + χ(2)DL(ω)(4)

While the molecular arrangement in the BIL is specific to each interface, and therefore the χ(2)BIL(ω) SFG signal is specific to each interface, the DL is always composed of bulk water reoriented by the field generated at the charged surface. Only the magnitude and sign of the field (equivalently of ΔϕDL) differ from one interface to another while there is a common structural origin for χ(2)DL(ω) at any aqueous interface. See our two most recent works dedicated to this issue.31,32 The χ(2)DL(ω) SFG signal therefore always has the same shape and position arising from the two-band structure inherited from χ(3)bulk(ω) (3200 and 3400 cm−1, see Tian et al.29 and theoretical papers from Joutsuka et al.60 and from us30) modulated in sign and amplitude by ΔϕDL (specific to each interface). More discussions and demonstrations can be found on this subject in ref. 31 and 32. The χ(2)DL(ω) signal can thus be rewritten as

image file: c9cp02766a-t6.tif(5)
where z is the distance from the surface, zBIL/DL is the lower limit of the DL water region (i.e., plane boundary separating BIL and DL water layers, roughly one water monolayer of ∼3 Å thickness above the surface in the systems investigated in the present work as well as in our previous work30), and z formally stands for the Debye-length limit over which EDC(z) is nonvanishing. image file: c9cp02766a-t7.tif is the potential difference across the DL, ΔϕDL, while Δkz is the phase mismatch between the emitted photons (the optical beams in the experiment vary in phase as they propagate into the aqueous phase), giving rise to interferences. Following Roke's detailed investigation,16 interferences become non-negligible for large Debye lengths, typically for ionic concentrations lower than ∼10−2 M; consequently they might alter both SFG intensities16 and band-shapes.15 See ref. 16 for Δkz expressions in relation with experimental set-ups. At the high ionic concentrations of most of the MD simulations in the literature, including the present ones at ≫10−2 M, the Debye length is only a few Å (∼3 Å Debye length at 2 M electrolyte concentration); interference effects are thus absent over such a short distance close to the surface and z in eqn (5) reduces to the DL layer extracted from our DFT-MD simulations (see Section 3 and Fig. 2). χ(2)BIL(ω) and χ(2)DL(ω) are calculated from the DFT-MD simulations through eqn (1)–(3), where the sum over m (eqn (1)) is reduced to the water belonging to BIL or DL regions only. Despite the fact that χ(2)DL(ω) could alternatively and equivalently be calculated through χ(2)DL(ω) = χ(3)bulk(ωϕDL, with ΔϕDL extracted from the DFT-MD trajectories, it arises from a pure second order process. At low ionic concentrations and at neat interfaces, interference effects have to be taken into account, either by including them in the theoretical SFG spectrum or by removing them from the experimental spectrum (as suggested by Roke, Geiger and collaborators15,16). One of these two actions has to be accounted for in order for theoretical and experimental SFG signals to be comparable. To the best of our knowledge, this has not been done in the literature yet, with the exception of model calculations done by Geiger's group.15

The rewriting of the SFG signal through eqn (2) and (3) for the molecular dipole/polarizability derivatives is based on the general strategy of APT (Atomic Polar Tensor, made of image file: c9cp02766a-t8.tif elements, with l one of the 3-cartesian coordinates of space) and Raman tensors (made of image file: c9cp02766a-t9.tif elements), exploited in these equations only along the O–H direction of each O–H oscillator (thus only well-designed for O–H stretchings). The image file: c9cp02766a-t10.tif and image file: c9cp02766a-t11.tif terms included in eqn (2) and (3) have been parameterized based on conformations generated for liquid water (see details in ref. 55). The derivatives were numerically calculated fully ab initio taking into account the variation of the total dipole/polarisability of the system following the elongation/contraction of one individual O–H bond-length (of one chosen molecule), making use of numerous decorrelated snapshots. This parameterization strategy takes into account not only the N-body polarisation effects acting on all water dipoles and polarisabilities due to the water molecules undergoing O–H motions (i.e., inter-molecular couplings), as the wavefunction of the modified system is systematically recalculated, but also charge fluxes in between water molecules. Within the dipole derivatives used for the parameterization, one has indeed (in the following equation, i represents the elongated/contracted water molecule and μ the total dipole moment vector of the system (image file: c9cp02766a-t12.tif with μi the dipole of molecule i)) image file: c9cp02766a-t13.tif, where image file: c9cp02766a-t14.tif are charge fluxes in between molecules i and j. The same is true for the polarizability derivatives.

With charge fluxes hence included in the parameterized derivatives, intermolecular couplings and cross-correlation terms are also included in a ‘mean field approach’, even though only self-correlation terms for the SFG signal (eqn (1)) are formally calculated. Also of importance in the model employed here are the velocities (eqn (2) and (3)) that take into account intra- and inter-molecular couplings from the dynamics, hence again introducing cross-correlation terms in a ‘mean field approach’ into the SFG self-correlation function of eqn (1).

We have also included solid surface O–H (silanols, Si–OH groups) contributions to the SFG theoretical spectra, as surface silanols vibrate in the 2800–3700 cm−1 stretching region,61 gaining SFG activity as soon as they have an out-of-plane (OP) orientation from the silica (planar) surface. The SFG signals arising from the O–H oscillators of the silica surface and of the water are added in order for the theoretical spectra to be comparable to experiments. The silanols’ contributions are calculated via the formalism in eqn (2) and (3) above. The parameterization of image file: c9cp02766a-t15.tif and image file: c9cp02766a-t16.tif (respectively of units Å2 and D/Å) can be found in ref. 59: image file: c9cp02766a-t17.tif, image file: c9cp02766a-t18.tif, image file: c9cp02766a-t19.tif, image file: c9cp02766a-t20.tif, image file: c9cp02766a-t21.tif, image file: c9cp02766a-t22.tif, image file: c9cp02766a-t23.tif, image file: c9cp02766a-t24.tif, image file: c9cp02766a-t25.tif. Smoothing procedures are applied to the calculated SFG spectra (Gaussian smoothing with σ = 35 cm−1).

Convergence of the SFG spectra over time-length

Time-convergence of vibrational spectroscopic signals is always a critical issue in theoretical spectroscopy, all the more when spectroscopy is extracted from rather short time-scale DFT-based MD simulations. In order to evaluate the simulation time needed to converge SFG spectra for the interfaces investigated in this work, Fig. 1 reports the final SFG spectra calculated at the quartz–water interface (QW, left side) and at the 2 M KCl electrolytic quartz–water interface (QW + KCl, right side), comparing the signals extracted over 15 ps (blue lines) and 100 ps (red lines) DFT-MD trajectories. Also shown is the comparison of the deconvolved SFG signals for water-BIL, water-DL, and water-bulk. As expected, the water-bulk SFG signal is systematically equal to zero. At the two interfaces shown here, the BIL-SFG signal is unchanged in going from the 15 ps to the 100 ps trajectory, which is due to the fact that these water molecules are interacting rather strongly with the oxide surface sites (especially at the QW interface), hence providing a rather static water-BIL (see Sections 3 and 4 for a detailed description of the water molecules that belong to the water-BIL). The final SFG signal at the neat QW interface is however changed from the 15 to the 100 ps trajectory, because of the change in the DL-SFG signal over the longer time-scale. The DL-SFG is entirely due to the non-linear third order susceptibility χ(3)bulk(ω) signal, which is a universal two-band signal with peaks at 3400 and 3200 cm−1.29,30,60 While the two-band structure is already correctly obtained from the 15 ps trajectory at the QW + KCl charged interface (see blue and red lines for QW + KCl in Fig. 1), it is not the case from the 15 ps trajectory of the neat QW interface where the 3400 cm−1 band is almost non-existent. The DL-SFG signal at charged interfaces stems from the surface field that induces reorientation of the water HB network within the DL (see Section 3 for all details): this signal is therefore dependent on the water orientational dynamics, which requires longer time-scales than 15 ps for a statistical representative sampling, as shown here at the neat QW interface. Once 100 ps trajectory is achieved, the correct two-band structure is recovered. The situation is different at the QW + KCl interface where the surface charge is larger than at the neat QW interface (which is revealed in the amplitude of the DL-SFG that will be discussed in more detail in Section 4), once the K+ cation is adsorbed at the surface: it thus reorients more strongly the water molecules in the DL, hence reducing the influence of the orientational dynamics effects on the convergence of the DL-SFG signal at this interface.
image file: c9cp02766a-f1.tif
Fig. 1 Convergence of the Imχ(2)(ω) SFG signal over trajectory length scale. Imχ(2)(ω) SFG signal (‘TOT’) and its decomposition into BIL-SFG (‘BIL’), DL-SFG (‘DL’) and bulk-SFG (‘bulk’) calculated for the neat quartz–water interface (neat QW, left side) and the 2 M KCl electrolyte at the neat quartz–water interface (QW + KCl, right side), comparing the results for a 15 ps trajectory (blue lines) and a 100 ps trajectory (red lines).

Summarizing, Fig. 1 shows that the 15 ps time-scale of DFT-MD trajectories is enough for BIL-SFG spectra to be converged, while longer time-scales improve the DL-SFG final spectrum for interfaces where there is a rather low surface (reorienting) field as is the case here at the neat quartz–water (QW) interface. Despite this improvement, the interpretation of the SFG features over the 15 ps trajectories is fully consistent with the one from 100 ps trajectories. This will be re-emphasized in Section 4 where the SFG spectra are analysed in detail for all three silica–water interfaces investigated here. Note also that the interfaces at the isoelectric point, as is the case for the amorphous–water interface (AW) investigated in this work, have no DL and therefore 15 ps trajectories are enough for full convergence of the interfacial (BIL-)SFG signal. Such short time-scale for SFG spectra calculation has already been emphasized in our previous works where experiment–theory comparisons were successfully presented,30,32,46 and is due to the rewriting of the dipole–polarisability correlation function into a velocity based correlation function together with the approach employed for the parameterization of the APT and Raman tensors that introduce charge fluxes in a ‘mean-field’ approach. This is not the strategy adopted in other works where longer time-scales for SFG calculations are proposed.62–65

3 Unraveling BIL and DL layers at various silica/water interfaces

We now apply our methodology to three silica/water interfaces in order to reveal the structural organisation of water into layers, namely BIL (Binding Interfacial Layer), DL (Diffuse Layer), and bulk, starting from the case of the quartz/water interface (hereafter denoted QW). Once the water layers are identified it is straightforward to calculate their associated spectroscopic SFG signatures, and hence deconvolve the SFG spectrum of the whole interface into contributions arising separately from each of these layers. The methodology is presented and validated in ref. 30.

The structural definitions employed for revealing BIL/DL/bulk layers are solely taken from the point of view of the water molecules, as they are the ones predominantly giving rise to the SFG signals at aqueous interfaces. Consequently, the structural definitions of the BIL and DL possibly differ from those of the Stern and diffuse layers defined in electrochemistry of electrolytic interfaces, as these ones take the viewpoint of the electrolytes.

Three descriptors of the water structure are used to reveal BIL/DL/bulk water layers at interfaces: water density evolution with respect to the distance from the surface, water coordination, and strength and orientation of water–water HBs. We use the Willard–Chandler instantaneous surface51 for the definition of the surface and its oscillations in space and time.

Definition and characterization of interfacial L0–L3 water layers

In Fig. 2a, the evolution of the water density with respect to the distance from the solid–liquid boundary (instantaneous surface) is plotted. Four different water layers can hence be identified in the case of quartz/water (QW) and quartz/water + 2 M KCl (QW + KCl) interfaces and three in the case of the amorphous silica/water (AW) interface, corresponding to the four/three peaks in the density profile. They are labelled L0, L1, L2 and L3 (see Fig. 2 for QW). In all three interfaces, water in layer L0 is characterized by a maximum in the water density profile (Fig. 2), with the density being found to be almost twice the bulk water density at the QW interface, 1.5 times the bulk water density at the QW + KCl interface and ∼1.2 times the bulk water density at the AW interface.
image file: c9cp02766a-f2.tif
Fig. 2 Water organisation at three silica/water interfaces: (a) quartz/water (QW, upper panel), (b) quartz/water + 2 M KCl electrolyte (QW + KCl, middle panel), (c) amorphous silica/water (AW, lower panel). Left: Schemes of the composition of the liquid phase in terms of BIL, DL and bulk (definitions detailed in the text). Right: Water density profiles evaluated with respect to the distance from the surface (based on the instantaneous surface definition51). In the density profiles, the water density is shown in black and in the case of the QW + KCl interface, the density profiles for Cl (red) and K+ (orange) are added to the plot, multiplied in intensity by a factor of five for the sake of clarity. In plot (a) for the QW interface, we also identify/illustrate the different peaks in the density profile as the different water layers, labelled L0, L1, L2, and L3. The correspondence with BIL, DL and bulk is made in the text using these density plots and the orientational 3D-plots in Fig. 3 and 4. The instantaneous surface (blue surface) and the simulation box (red oxygen, white hydrogen, yellow silicon) for the AW interface are illustrated in panel (d) (taken from one snapshot along the DFT-MD).

The number of water–water H-bonds formed in these layers (Table 1) shows that only layer L0 at these three interfaces has a lower value than the one in the bulk. The number of water–water H-bonds formed in L0 is specific to each oxide surface and differs significantly from the liquid reference value (i.e., less than the 3.4 reference value in bulk liquid water from DFT-MD35,46,52). In the case of the QW interface, the water molecules in L0 make on average 2.8 HBs with other water molecules, plus 0.5 HBs with solid surface O–H groups. This number decreases to 2.2 and 2.3 in the case of the QW + KCl interface and AW interface respectively. In L1, L2 and L3 layers, the water molecules instead recover the average 3.4 HBs/water obtained in bulk liquid water.

Table 1 Number of water–water H-bonds per water molecule reported for each layer. Layers are denoted here as L0–L3; see text for definitions. This table of H-bonds per molecule serves for the final classification of interfacial layers into BIL (Binding Interfacial Layer), DL (Diffuse Layer) and bulk regions. The reference value for bulk liquid water at the BLYP-D2 level of theory is 3.4. QW is the (0001) α-quartz/water interface, QW + KCl is the (0001) α-quartz/water interface accommodating a nominal aqueous 2 M KCl concentration, and AW is the amorphous silica/water interface
L0 2.8 2.2 2.3
L1 3.4 3.4 3.4
L2 3.4 3.4 3.4
L3 3.4 3.4

The joint probability for the water molecules to donate H-bonds with a given O–O distance and a given H-bond orientation with respect to the normal to the surface is then calculated layer per layer for each interface, and compared to the reference for bulk water (Fig. 3b; see also the Methods and computational details section). In these plots, one focuses on the maxima that appear in red-yellow colors. Fig. 3a reports the 3D-plots obtained for the L0–L3 layers at the QW interface as a test case. Starting with layer L3, the plot displays a uniform distribution in cosine values for HB distances in the range of ∼2.5–2.9 Å (bluish background), as for bulk liquid water (Fig. 3b), i.e., all possible orientations of the water–water H-bonds with respect to the normal to the surface are equally probable. It shows that the isotropy in the HBs typical of bulk water is recovered in layer L3 (e.g., the average of the cosine values is 0.01, within our numerical error for the zero value expected for centrosymmetric bulk water). Note that the centrosymmetry of this layer will be confirmed by its SFG spectrum in Fig. 6a. Conversely, the 3D-plots for all the other water layers show that the uniform distribution in cosine values is lost, as the bluish distribution is not uniform for all possible cosine values or even red-spots are also present.

image file: c9cp02766a-f3.tif
Fig. 3 Characterization of the water H-bond networks at the neat quartz/water (QW) interface. (a) 3D-plots of joint probability to donate one H-bond with a given O–O distance and given orientation relative to the normal to the water surface (cos[thin space (1/6-em)]θ) for the O–H groups of the water molecules in each layer, i.e., L0, L1, L2, and L3 as defined in Fig. 2a. The maxima in the joint probability appear in yellow-red colors in the plots. (b) Reference 3D-plot for bulk water (left), and scheme for the definition of the θ angle and the normal to the surface (oriented out of the aqueous phase and towards the solid).

Definition of the interfacial and bulk water regions

The first conclusion to be drawn is hence that the bulk isotropic structure is recovered on average for distances beyond 8.4 Å from the quartz surface at the neat interface. This distance is reduced to 6.0 Å in the case of the QW + KCl interface, where layers L2 + L3 are already centrosymmetric (average cosine values of −0.01 in our simulations), and to 3.0 Å for the AW interface, where even L1 is already bulk for our descriptors (average cosine values of 0.01 in our simulations).

Definition of the DL interfacial water layer

Coming back now to the QW interface, in layers L1 and L2 the average number and strength of HBs are found to be the same as in the bulk (see Table 1), and the water density (Fig. 2a) is on average close to bulk liquid water, but a net orientation in the HB-network is clearly observed (Fig. 3a), with a preferential orientation of the water–water H-bonds pointing towards the surface (red spot in the 3D plots at cosine values of ∼+0.6/+1.0, with an average cosine value of +0.13). These structural properties correspond to water in the DL layer (as defined in the works from ref. 29 and 30), i.e., bulk liquid water reoriented by an electric field. There is indeed an electric field at the neat QW interface, despite the fact that the surface of quartz is fully hydroxylated, as shown in our previous work.35,36 This is due to Si–OH surface silanols carrying a partial negative charge when in contact with liquid water.36 At the QW interface, the DL is thus composed of L1 and L2 water layers, located between 2.6 and 8.4 Å from the instantaneous surface (Fig. 2a), with a net orientation of H-bonds towards the solid surface (red spot at positive cosine values in the 3D-plot of the DL, see Fig. 3 and 4). The opposite is observed when the 2 M KCl electrolyte is present at the interface, as shown by the 3D-plot of the DL for the QW + KCl interface (Fig. 4), which displays a red spot at negative cosine values (the average cosine value is −0.18). In the latter case, only L1 is reoriented by the field while layer L2 is centrosymmetric (as in the bulk). There is on the other hand no DL at the AW interface, and layers L1 and L2 have all properties of bulk liquid water. This is because the amorphous surface is at the isoelectric point. It is remarkable that only the QW and QW + KCl interfaces have such a DL. It extends between 2.6–8.4 Å and 2.6–6.0 Å above the surface respectively for QW and QW + KCl. By integrating the electric field across the DL boundaries one obtains the difference in the Hartree electric potential across the DL (ΔΦDL): values of 0.01 V for QW and of −0.02 V for QW + KCl are obtained. Beware that these DL potentials cannot be directly compared to Galvani potentials.
image file: c9cp02766a-f4.tif
Fig. 4 Characterization of the water H-bond networks at QW and QW + KCl interfaces. Joint probability to donate one H-bond with a given O–O distance and given orientation for the O–H groups of the water molecules in the DL (left) and in the bulk (right), calculated for the QW (up) and QW + KCl (down) interfaces. The orientation of the O–H group donating the H-bond is given by the cosine of the θ angle between the vector connecting O to H and the normal to the surface. The maxima in the joint probability appear in yellow color in the plots.

Definition of the BIL interfacial water layer

Finally, the water structure in layer L0 differs significantly from the one in the bulk and in the DL for the three structural descriptors used in this work. In the case of the QW interface the water density is ∼1.8 times larger than in the bulk, the number of water–water HBs per molecule is 2.8, and a net orientation in the water–water HB-network is observed. The 3D-plot in Fig. 3 shows that there are two separate maxima in the joint-probabilities for water in L0, one at negative cosine values (−0.6/−1.0) due to water molecules donating H-bonds to water located in L1, and the other at cosine values between −0.2 and +0.2 due to water–water H-bonds parallel to the surface. A similar maximum of probability to form HBs parallel to the surface is not observed in any of the other layers. Together with these water–water H-bonds, there is an additional 0.5 water–surface silanol H-bond per water molecule (not seen in this water–water 3D-plot), already described in ref. 12, 35 and 36. There is indeed an alternate motif of in-plane/out-of-plane surface silanols respectively being a H-bond acceptor/donor with water molecules located in L0. See the illustration in Fig. 7b. We have shown in ref. 12 and 36 that the quartz surface at the interface with water is made of a 2D motif of in-plane (IP) and out-of-plane (OP) silanols. These two populations are respectively acting as a H-bond acceptor and donor with the water molecules located just above the surface in the BIL. Consequently, the structural organisation of BIL-water follows the order induced by the IP/OP silanols, with a regular alternate motif of water molecules being a H-bond donor (denoted ‘DS’ for Donor to Solid) and acceptor (denoted ‘AS’ for Acceptor from Solid) to the O–H terminations (geminal-Q2 type). On top of these oxide–water H-bonds, a third water population (WB) bridges the DS/AS water molecules making H-bonds oriented parallel to the surface. See Fig. 7b and the joint probabilities of angles in the 3D-plot (Fig. 3) of layer L0 (=BIL) at cosine values −0.2/+0.2. There are on average 9.3 oxide–water H-bonds per nm2 formed at the QW interface.

In layer L0 the structural organisation of water is, not surprisingly, dictated by the direct interactions with the surface sites. This is the layer that has been called ‘Binding Interfacial Layer’ (BIL) in ref. 29 and 30. The BIL is thus composed of one water monolayer (L0), vertically located between 0.0 and 2.6 Å from the surface at the neat QW interface (Fig. 2a).

Such regular organization is disordered not only at the AW (amorphous/water) but also at the QW + KCl (quartz/water + 2 M KCl) interface, as shown in ref. 13, 35 and 36. Due to the corrugated amorphous surface and less acidic surface silanols (Q2 and Q3 types, geminals, vicinals, isolated),13,14,66 the water molecules can no longer interact so strongly with the silanols and form the ‘ideal’ quartz–water repeated network of Fig. 7b. Less and weaker H-bonds are hence found at the AW interface. There are 7.0 oxide–water H-bonds per nm2 formed at the AW interface (vs. 9.3 at the QW interface). Furthermore, once the 2 M KCl electrolyte is added at the QW interface, K+ and Cl ions form solvent separated ion-pairs,35 the cation being adsorbed at the surface and using surface silanols to complete its solvation shell (inner-sphere adsorption), while the anion lies farther away from the surface,35 fully solvated as in pure liquid water (see Fig. 2b). The inner-sphere cation affects the structural organization of both surface silanols (i.e., inducing more IP silanols) and water (i.e., reorienting water molecules and reducing surface–water and water–water numbers of H-bonds), all detailed in ref. 14, 35 and 36. These findings are in agreement with previous theoretical studies using classical MD simulations from Dewan et al.24 and Machescky et al.37,38 concluding that small alkali cations undergo inner-sphere adsorption at the silica surface, altering the water organization both in the topmost interfacial layer (i.e., the BIL) and in the DL, and from Bouhadja and Skelton,67 concluding that cations like K+ and Na+ are attracted by the quartz-surface with consequences on interfacial water dynamics. There are hence 7.2 oxide–water H-bonds per nm2 formed at the QW + KCl interface (vs. 7.0 at the AW and 9.3 at the neat QW). The K+ cation has therefore a chaotropic effect, disordering the water of the (initially ordered) BIL (at the neat interface), by breaking both solid–water and water–water H-bonds. Note that the solvent separated KCl ion-pair, organized with the cation localized in the BIL and the anion in the DL, is maintained over the 100 ps DFT-MD trajectory performed in this work. Although ref. 36 has shown that taken separately, the K+ BIL-inner-sphere adsorption is energetically favored at the QW interface while it is not the case for the Cl anion (hence pushed farther away into the DL), other ion-pair conformations should also be considered (e.g., a contact ion-pair located in the BIL, etc.) if one were interested in the full structural characterization of the electrolytic quartz–water interface in order to make a 1-to-1 assignment of the QW + KCl theoretical SFG spectrum to the SFG experiments (which would require spectral weighting over all possible ion-pair configurations at the interface). As our goal is to characterize changes in the SFG spectroscopy of the QW interface once an aqueous molar concentration of KCl electrolyte is introduced, and especially to assess the changes in the BIL and DL spectra, without making a full comparison to the experiments, one single interfacial KCl structure is enough for our purpose.

4 Rationalization and comprehensive view of the silica/water SFG spectra

In the following, we do not make a 1-to-1 theory–experiment comparison but instead we qualitatively use our theoretical SFG spectra, deconvolved into BIL-SFG and DL-SFG, in order to provide a rational understanding of the spectral features that have been previously recorded and presented in the literature for silica–water interfaces. The presentation of results in this section uses acronyms for the molecular systems that are investigated, which we recall once again here: QW is the crystalline neat fully hydroxylated (0001) α-quartz/water interface under PZC conditions, AW is the fully hydroxylated amorphous silica/water interface, and QW + KCl is the fully hydroxylated (0001) α-quartz surface interfaced with 2 M KCl aqueous solution.

We make the choice to discuss in this section the theoretical spectra calculated over 15 ps trajectories (common time-length for all three investigated interfaces of this work). We have demonstrated in Section 2 that the only modification from the 15 ps to 100 ps trajectory has been observed for the neat QW interface DL-SFG signal, because of the longer time-scale for the orientation dynamics of the DL-water molecules under the rather low surface field that arises at the neat interface. Despite this difference, the interpretation done using the 15 ps trajectory is fully consistent with the one obtained from the 100 ps trajectory.

Let us first comment the SFG experimental data on silica/water interfaces reported in Fig. 5. One main set of experimental data is from Shen et al.6 with phase-resolved SFG at the neat quartz/water interface measured at pH values of 1.5 and 6.5 (Fig. 5a). The second set of experiments (Fig. 5b) is for the amorphous (fused) silica/water interface, i.e., conventional |χ(2)(ω)|2 SFG by Tyrode et al.8 (two curves reported because of two surface treatments applied in the experiment, pH ∼ 6), conventional (black) and phase-resolved (red) SFG by Tahara et al.7 (pH ∼ 2), and conventional SFG by Shen et al.1 (pH ∼ 1–3). All measurements are for the ssp polarization. Briefly, the quartz–water experimental spectra from Shen's group (Fig. 5a, pH = 1.5 and pH = 6.5) have been chosen since they are the closest experimental data in the literature to our simulations, i.e., alpha-(0001) quartz–water interface and PZC conditions of pH = 2–4 (thus the fully hydroxylated surface in the simulation box of that system). The fused/amorphous silica–water experimental spectra (Fig. 5b) have been selected from the literature following two criteria: (1) experiments performed under pH conditions as close as possible to PZC (i.e., again corresponding to fully hydroxylated amorphous surfaces in our simulation boxes), and (2) experiments where the 3650 cm−1 band can be clearly observed. As pointed out in ref. 8 and 59, such a band is a characteristic feature of water located above hydrophobic patches of the amorphous silica surface that are indeed present at the nanoscale level. As shown in ref. 8 the intensity of such a band depends on the preparation of the silica surface (for instance, a more intense 3650 cm−1 band is obtained by heat-treating the surface). We hence selected experiments where such a band is more clearly visible, since it has a central role in relation to the water ordering/disordering in the BIL at silica–water interfaces (as discussed in this manuscript). These selected experiments are thus from Tyrode's group (Fig. 5b, top), Tahara's group (Fig. 5b, middle) and Shen's group (Fig. 5b, bottom).

image file: c9cp02766a-f5.tif
Fig. 5 Silica–water interfaces. (a) Comparison of our theoretical Imχ(2)ssp(ω) SFG spectrum under PZC conditions (DFT-MD PZC), calculated by including both water and solid surface contributions, with the experiments from Shen's group6 at pH values of 1.5 and 6.5 (below and above PZC), for the (0001) α-quartz–water interface. (b) Experimental evidence for the quasi-free OH (3600–3650 cm−1) band at an amorphous (fused) silica/water interface. Upper panel: Tyrode's group experiments,8i.e., |χ(2)ssp(ω)|2, heat surface treatment (green) and plasma surface treatment (black). Middle panel: Tahara's group experiments,7i.e., Imχ(2)ssp(ω) in red and |χ(2)ssp(ω)|2 in black. Lower panel: Shen's group experiments,1i.e., |χ(2)ssp(ω)|2 in black.

image file: c9cp02766a-f6.tif
Fig. 6 Phase-resolved SFG spectra calculated from DFT-MD simulations (O–H stretching region), including only the water contribution. (a) Imχ(2)(ω) spectra (calculated from 15 ps trajectories) for neat quartz–water (QW, black), amorphous silica–water (AW, cyan) and quartz–water + 2 M KCl (QW + KCl, red) interfaces at PZC (i.e., pH = 2–4, corresponding to fully hydroxylated silica surfaces in the simulations), calculated including all the water molecules from all layers (TOT), including only the water molecules located in the BIL (BIL), including only the water molecules located in the DL (DL), and including only the bulk water molecules (bulk). (b) Imχ(2)(ω) spectra calculated for QW interfaces with 0% (black), 3% (cyan) and 12% (red) of surface SiOH sites being deprotonated. The simulations for 3% and 12% deprotonation are taken from ref. 31. The total spectra (TOT) as well as the deconvolved BIL/DL/bulk spectra are shown.

image file: c9cp02766a-f7.tif
Fig. 7 Silica surface and BIL-water contributions to SFG. (a) Imχ(2)(ω) spectra calculated for the neat quartz/water interface (QW, 100 ps trajectory) and the amorphous-silica/water interface (AW, 15 ps trajectory) at PZC (i.e., pH = 2–4, corresponding to fully hydroxylated silica surfaces in the simulations), either including only the contribution from water molecules (upper panels) or adding up the contribution from surface out-of-plane silanols (lower panels). (b) Left: Schematic representation of the repeated structural unit observed in the BIL, with ‘DS’ being a water molecule donating an H-bond to an in-plane (IP) silanol, ‘AS’ a water molecule accepting an H-bond from one out-of-plane (OP) silanol, and ‘WB’ a water molecule non H-bonded to any surface silanol but bridging AS/DS water molecules. Right: Individual contributions to Imχ(2)(ω) for AS (blue), DS (red) and WB (orange) water populations.

As mentioned in the previous sections, when interference effects are present and not compensated for, as is the case of all experimental SFG spectra recorded to date for silica–water interfaces, a direct comparison between SFG spectra obtained from different experimental set-ups cannot be strictly done. In the same way, also a direct comparison between theory and experiments is not possible, since in theory χ(2)BIL(ω) and χ(2)DL(ω) are not taking into account the interference terms. See Section 2 for a full discussion.

Despite the challenges in comparing SFG experimental spectra of silica/water interfaces because of the interference term in eqn (5) that should be accounted for, either experimentally or theoretically, it remains that a two-band structure with maxima located at ∼3200 and ∼3400 cm−1 is systematically recorded for the SFG spectra of neat silica–water interfaces, with the addition of a supplementary quasi-free O–H band around 3600–3650 cm−1 for fused silica/water (Fig. 5b). As shown by Shen et al.,6 the intensity ratio between the ∼3200 and ∼3400 cm−1 bands is pH dependent (Imχ(2)(ω), red lines in Fig. 5a), with the ∼3400 cm−1 band being more intense at low pH, while it is the ∼3200 cm−1 band that is more intense at high pH. This ratio could however be modified once the SFG intensity is corrected for the photon interference term.

Revealing the microscopic origin of BIL-SFG and DL-SFG signatures at the three interfaces

In Fig. 6a we present our theoretical phase-resolved SFG spectra for the three silica–water systems investigated here, including only the contribution coming from water, together with the decomposition into individual contributions from BIL/DL/bulk layers (as unraveled in the previous sections). The total SFG spectra including both water and surface silanol contributions to SFG are presented in Fig. 7a for both neat QW and AW interfaces, as surface out-of-plane silanols are SFG active and contribute to the experimentally recorded spectra (see ref. 61 and 7 and 68, where silanols’ contribution to the final SFG is inferred). The calculated SFG spectra are given with their absolute intensities (m2 V−1 unit).

The same two bands structure as in the experiments is present also in our theoretical SFG spectrum of the neat QW interface at PZC (∼3400 and 3150 cm−1 positions, Fig. 6 and 7a), with the ∼3400 cm−1 band being the most intense (this is better appreciated in the 100 ps trajectory SFG-signal in Fig. 1, because of the better convergence of the DL-SFG signal at this interface over the longer trajectory, see Section 2). This shows that the two-bands structure of the SFG spectra of silica–water interfaces is maintained also when the phase mismatch affecting the experimental data is not taken into account, i.e., it is due to the intrinsic structural properties of the interface.

In order to provide a full rationalization of the molecular origin of these two bands, we have deconvolved our theoretical spectrum of the QW interface in terms of the contributions arising from the BIL and DL. To this aim, we present in Fig. 6aχ(2)(ω) spectra separately for the water molecules belonging to the BIL and DL layers (the spectra calculated for the QW interface at PZC are colored in black). With the decomposition into contributions from different layers, we show that only the band located at 3400 cm−1 arises from the BIL-water, with a much lower contribution arising from the DL at the same frequency (see also Fig. 1 with the 100 ps trajectory), while the band at 3200 cm−1 arises solely from the DL-water.

As discussed in Section 3, three different water populations are present in the BIL at the QW interface. Two of these water populations are respectively ‘AS’, i.e., accepting H-bonds from the surface OP silanols, and ‘DS’, i.e., donating H-bonds to IP silanols, with DS/AS water molecules being connected by water bridge molecules (‘WB’) forming in-plane HBs. See Fig. 7b for an illustrative scheme. Strikingly, despite three different water populations being present in the BIL, only one dominant vibrational band is observed in the Imχ(2)(ω) calculated for this interfacial layer. On evaluating the SFG contribution of AS, DS and WB water molecule populations separately (Fig. 7c), we clearly see that the BIL-SFG signal is dominated by the O–H groups of the DS molecules giving H-bonds to in-plane surface silanols (SiOH). As these DS water molecules are strongly oriented along the normal to the surface, they lead to the intense SFG positive signal (Fig. 7c, red line) with the maximum at 3400 cm−1 in the total SFG spectrum. The WB population forms H-bonds oriented parallel to the surface, thus not contributing to the ssp SFG signal (orange line), while the AS population provides a low intensity negative band below 3300 cm−1 (blue line), mostly compensated by the tail of the positive band arising from the DS population. It is important to stress here that the negative band is observed in the Imχ(2)(ω) experimental spectra from Shen and Tahara in Fig. 5.

As already concluded in recent investigations7,61,68 on aqueous fused-silica surfaces, also out-of-plane surface silanols contribute to the negative Imχ(2)(ω) in this region. This assignment is fully confirmed in our simulations, as our theoretical spectra have the same intensity for the negative band as in experiments, for both QW and AW interfaces, once the SFG contribution from the solid surface SiOH groups is included in the theory (Fig. 7a, bottom). The negative band at ∼3200 cm−1 in experiments is thus arising from both AS–water O–H and surface out-of-plane SiOH groups.

If we now compare the Imχ(2)(ω) calculated for the water-BIL at the neat QW interface (Fig. 6a, black line) with the ones of the QW + KCl interface (red line) and of the AW interface (cyan line), we observe that the quasi-free O–H band (∼3650 cm−1) appears at these two interfaces. Moreover the H-bonded O–H positive band (3400 cm−1) is reduced in intensity and blue-shifted to ∼3500 cm−1, thus now overlapping with the tail of the high frequency band.

As already pointed out in our structural analyses, the interfacial regular structure in the BIL observed for the QW interface is partially lost for both QW + KCl and AW interfaces, and the average number of H-bonds per water molecule in the BIL is reduced from 3.3 (counting water–water + water–oxide H-bonds) for the QW interface to 2.7 and 2.6 respectively for the AW and QW + KCl interfaces. A new BIL-water population is found at these two interfaces, i.e., having one O–H group not H-bonded (or weakly H-bonded) although pointing towards the solid oxide. These dangling O–H groups do not have the geometrical requirement(s) for H-bonds with surface silanols, but they still weakly interact with the solid surface oxygen sites. Consequently, the SFG positive peak associated to these dangling OH stretching motions is red-shifted with respect to a “pure” free OH stretching motion (3700 cm−1), leading to the signal observed at ∼3650 cm−1. In a recent investigation,59 we have shown that the 3650 cm−1 band at amorphous/fused silica–water interfaces is in particular due to water O–H groups pointing toward surface Si–O–Si (siloxane bridge) sites, and its intensity correlates with macroscopic contact angle measurements: the less the silica surface is hydrophilic (i.e., the larger the contact angle), the more the Si–O–Si sites that are exposed at the surface and the higher is the 3650 cm−1 intensity in the Imχ(2)(ω) spectra. The Si–O–Si sites have been further shown to be inhomogeneously distributed on the amorphous silica surface, leading to the formation of water hydrophobic patches above the areas with higher Si–O–Si density.59 Interestingly, we find here that the same disordered water structure in the BIL, as observed above the hydrophobic patches of the amorphous silica, is recovered also when electrolytes are adsorbed at the crystalline quartz–water interface: it is indeed measured by the similar decrease in BIL-water coordination at the QW + KCl and AW interfaces with respect to QW (discussed in Section 3), and experimentally probed by the similar 3650 cm−1 intensity for both QW + KCl and AW interfaces (Fig. 6a).

We can therefore conclude that the SFG peak located at 3650 cm−1 is a marker band for the degree of water ordering/disordering in the BIL: no such band is observed for crystalline silica/water interfaces where a regular pattern of water molecules at the interface with silanol sites is built, but the marker band is systematically present once a perturbation to that water ordering is introduced (irregular surface and/or adsorbed electrolytes as shown here).

Chaotropic effect of inner-sphere adsorbed KCl electrolytes at the quartz–water interface

To further investigate the relationship between the 3650 cm−1 dangling OH band and the ion adsorption at the QW interface, Fig. 8 reports the integrated intensities of the 3400 cm−1 band (“order band”, from DS water) and of the 3650 cm−1 band (“disorder band”, from water with dangling OH) as a function of KCl concentration, from our DFT-MD simulations. The points at 0 M and 2 M have been obtained from the QW and QW + KCl interfaces investigated in this work, while the points at 0.2 M and 1 M have been obtained from DFT-MD simulations (with identical computational set-up) of the same QW interface with 0.2 and 1 M KCl concentration taken from our recent investigation in ref. 31. The SFG-spectra for QW interfaces with 0.2 M and 1 M KCl electrolytes, from which the integrated intensities are calculated, have been already reported by our group in ref. 31. As shown in the figure, the intensities of the order/disorder bands follow a clear trend with ion concentration: the higher the KCl concentration, the more cations are inner-sphere adsorbed in the BIL (see ref. 35), the more extended are the disordered water areas within the BIL, the higher is the disorder-band intensity and the lower is the order-band intensity.
image file: c9cp02766a-f8.tif
Fig. 8 Order/disorder marker bands vs. KCl concentration. Integrated Imχ(2)(ω) intensity of the 3400 cm−1 band (“order band”) and of the 3650 cm−1 band (“disorder band”) as a function of KCl concentration. The points at 0 M and 2 M have been obtained from the QW and QW + KCl interfaces investigated in this work, while the points at 0.2 M and 1 M have been obtained from DFT-MD simulations (with identical computational set-up) from ref. 31.

The order/disorder marker bands rationalize the SFG signals from Tahara et al.7 (phase-resolved SFG) and from Tyrode et al.8 (conventional SFG) on fused silica interfaces (Fig. 5b). In Tyrode's work, the effect of two surface treatments on the amorphous silica surface has been exploited, leading to quantitatively different SFG signals, but both spectra display the same band around 3650 cm−1. Similarly, the HD-SFG spectrum from Tahara et al. at pH ∼ 2 has the positive 3650 cm−1 band. Interestingly, while such a quasi-free O–H marker band is present in the Imχ(2)(ω) spectrum, it is hard to recognize it in the |χ(2)(ω)|2 spectrum (black line in the figure, Tahara's data), in which it is hidden into the convolution with the Reχ(2)(ω) part. The same remark can be applied to the experiments reported by Shen et al.6 at pH ∼ 1–3 (lower panel of Fig. 5b), where the |χ(2)(ω)|2 spectra of the fused silica/water interface at low pH values do not show a prominent peak at 3650 cm−1.

Rationalization of the SFG spectra of quartz–water interfaces as a function of pH/surface protonation state

Coming back to the decomposition of our theoretical SFG spectra into BIL/DL contributions (Fig. 6) at the quartz–water interface, now that we have shown that the 3400 cm−1 Imχ(2)(ω) band arises predominantly from the BIL-water, while the one at 3200 cm−1 arises from the DL-water, we further investigate the evolution of the two bands with pH (as in Shen's experiments in Fig. 5a). To this end, the theoretical Imχ(2)(ω) spectra (as well as deconvolved BIL/DL/bulk signals) calculated for the quartz–water interface under increasing pH conditions are reported in Fig. 6b. The increasing pH is mimicked in the simulations by introducing an increasing number of deprotonated (SiO) silanol sites at the surface. The systems considered are QW interfaces with 0% (neat QW interface), 3% and 12% of surface SiOH sites being deprotonated. The simulations with 3% and 12% deprotonation are taken from our recent investigation in ref. 31. As can be seen from BIL-SFG and DL-SFG intensities reported in the second and third panels of Fig. 6b, at increasing surface deprotonation (i.e., increasing pH conditions), the 3200 cm−1 DL-SFG intensity increases, while the 3400 cm−1 BIL-SFG intensity decreases (the integrated area of DL-SFG becomes larger than that of BIL-SFG). The increase of DL-SFG intensity is simply due to the increase in the negative surface charge resulting from the deprotonated sites: the more the quartz surface is deprotonated, the higher is the negative surface charge, the higher is the potential difference across the DL, and the higher is the DL-SFG intensity (see eqn (5)). The decrease of the 3400 cm−1 BIL-SFG band intensity with increasing pH is due to the appearance of a new SFG-active water-OH-population, which is composed of BIL-water molecules donating one H-bond to a SiO site instead of to an in-plane SiOH site. Once some of the surface SiOH groups are deprotonated, BIL-water molecules interact with the resulting SiO groups by donating stronger HBs compared to the BIL-water molecules donating to in-plane SiOH groups (DS water population). As a consequence, the associated positive signature in BIL-SFG Imχ(2)(ω) spectra shifts in intensity from 3400 cm−1 to 3200 cm−1. The 3200 cm−1 BIL-SFG band has been already experimentally detected in ref. 68 and it has been assigned to water donating HBs to SiO groups at the silica surface, in agreement with our findings.

In the light of these results, we are finally in a position to rationalize the evolution of the ∼3200 and ∼3400 cm−1 bands with pH for the QW experimental spectra in Fig. 5a. At low pH, when the surface charge and thus the potential across the DL are at their minimum values, the BIL contribution (the ∼3400 cm−1 band) is dominant in the SFG spectrum, while on increasing the pH the quartz surface becomes more and more deprotonated so that the potential across the DL increases. The DL contribution to the SFG thus becomes more prominent in the total SFG signal, while the ∼3400 cm−1 band loses intensity due to the appearance of BIL-water OH-groups donating strong HBs to SiO surface sites, shifting their vibration down to ∼3200 cm−1. At high pH values (higher than 6.5), the ∼3200 cm−1 band hence becomes the most intense in the whole SFG spectrum. The here revealed distinct pH dependence of the SFG signatures from BIL and DL is in nice agreement with a recent investigation from Darlington et al.,69 concluding that two interfacial water populations coexist at the silica–water interface, i.e., water molecules located in the Stern and Diffuse layer, and their SFG signatures have distinct pH dependence.

The present assignment of the positive BIL-SFG intensity to water molecules having one OH-group pointing to the silica surface (and vibrating at frequencies that depend on the interacting surface Si–OH (3450 cm−1)/Si–O (3200 cm−1) sites) furthermore allows us to make a microscopic rationalization of two recent SFG observations made by Tahara's group based on isotopic dilution and pH/electrolyte conditions. (1) In ref. 7 related to pH ∼ 6 silica–water interfaces, the 3200 cm−1 band has been shown to almost disappear upon isotopic dilution, while the 3450 cm−1 band remains. At this pH, the 3450 cm−1 band is indeed not affected by isotopic dilution since, as shown here, it mostly arises from ‘DS’ BIL-water molecules (Fig. 7) having one OH-group H-bonded to an IP-silanol, while the other O–H group is H-bonded either to a BIL-water and hence being in-plane and not SFG active or to a DL-/bulk-water molecule, hence contributing to the negative SFG-band at frequencies ≤3300 cm−1. The two OH groups of these ‘DS’ BIL-water molecules are therefore vibrationally uncoupled, so that their signatures do not change upon isotopic dilution. The 3200 cm−1 band in the SFG spectrum, which dramatically changes upon isotopic dilution, is in contrast solely due to DL-water. The isotopic dilution effect seen experimentally on such a band is thus expected, since the DL-SFG activity is due to the reoriented bulk water, which is well known to exhibit strong intra- and inter-molecular couplings, as also pointed out in ref. 29. (2) In ref. 68, where pH > 10 conditions are applied, the DL-SFG contribution is now screened by the 2 M NaCl electrolyte. The SFG positive band observed at 3200 cm−1 is known from our current work to be due to BIL-water molecules donating H-bonds to surface SiO, instead of the 3450 cm−1 SFG signature obtained at lower pH values from the BIL-water HB-donor to surface IP/Si–OH.

Before concluding, it is important to mention that no “disorder band” (i.e., a positive BIL-SFG band centered at 3650 cm−1) is observed when increasing the pH from 2 to 10 (corresponding to surface deprotonation from 0% to 12% in our DFT-MD simulations), contrary to what is observed for the increase in the concentration of the KCl electrolyte (see Fig. 6 and 8). As already shown, inner-sphere adsorbed ions locally reorient the SiOH surface sites in-plane, and hence break the crystal-induced ordering initially observed in the BIL at the neat interface. It seems then reasonable to hypothesize that a similar effect on the BIL structure can also be obtained because of a large deprotonation of SiOH sites at very high pH values. Following this hypothesis, pH induced BIL enhanced disordering for higher pH values than the ones already explored in this work should be observed. Such order/disorder transition could be tracked by BIL-SFG via the relative intensities of the positive 3450/3650 cm−1 bands (i.e., order/disorder bands).

This is indeed what a DFT-MD simulation of the QW interface with 24% of SiOH silanol sites being deprotonated (pH ∼ 13) provides. Fig. 9 reports the BIL-SFG signal calculated for all deprotonated QW interfaces investigated in this work (0%, 3%, 12%, 24%); the BIL-SFG spectra including only water contributions have already been presented in Fig. 6 for the series 0–12% deprotonation states. As shown in Fig. 9, two major changes are observed in the BIL-SFG spectra as a function of increasing bulk pH/increasing deprotonated surface silanols: (1) the negative band at frequencies <3200 cm−1 present in the spectrum of the neat QW interface (i.e., 0% deprotonation) loses intensity with increasing pH, and (2) the positive band at 3450 cm−1 present in the spectrum of the neat QW interface blue-shifts for the most basic conditions (i.e., 24% deprotonated silanols). Point (1) is mostly due to the out-of-plane silanols, that contribute to the negative band at low pH values (as already discussed in this manuscript), while out-of-plane silanols are progressively reoriented in-plane when pH increases, up to the point where virtually no more out-of-plane silanols are statistically present at the QW interface in the most basic pH conditions. Point (2) is due to the breaking of the regular alternate motif of ‘DS’ and ‘AS’ water molecules initially observed in the BIL under PZC conditions, with the consequent appearance of the “disorder band” and decrease in the intensity of the “order band”, as already discussed in this text for QW interfaces with increasing KCl concentrations. This confirms that the BIL-disordering induced by inner-sphere adsorbed ions is also promoted by the most basic pH conditions/higher number of surface deprotonated silanols. The microscopic mechanism leading to BIL disordering is therefore found as follows: by increasing pH conditions, the percentage of deprotonated SiOH increases, each SiO inducing local in-plane reorientation of adjacent silanols (similarly to the in-plane reorientation induced by the inner-sphere adsorbed KCl electrolyte), locally breaking the regular pattern of in-plane (IP)/out-of-plane (OP) SiOH surface sites. This in turn prevents the templating effect of the surface IP/OP motif to act as the ordering motif of the BIL-water structure, up to the point that an order/disorder transition is observed for the higher surface deprotonation state (here 24%), which leads to the appearance of the “disorder” SFG-band and overall blue-shift of the positive BIL-SFG Imχ(2)(ω) band.

image file: c9cp02766a-f9.tif
Fig. 9 Imχ(2)(ω) BIL-SFG spectra calculated from DFT-MD simulations (O–H stretching region) of QW interfaces with 0% (black), 3% (green), 12% (magenta) and 24% (orange) of surface SiOH sites being deprotonated. The BIL-SFG spectra are calculated including both BIL-water and solid contributions (BIL), only BIL-water contributions (BIL-water) and only solid surface contributions (BIL-solid).

Importantly, such blue-shift at increasing pH values has never been reported in the long literature of SFG spectroscopy at QW interfaces. Confirmation of our DFT-MD results could be obtained through BIL-resolved SFG spectra, by applying the recent developments from ref. 15 and 29. As a final note, same pH-induced order/disorder transition is expected to occur at fused silica/water interfaces, however starting at lower surface deprotonation states than for crystalline QW. As shown in this work, far less ordered BIL-motifs are seen at amorphous surfaces, hence the lower pH conditions under which disorder can be triggered.

Summarizing, the SFG spectra of silica–water interfaces are thus found to be modulated by changes in the degree of crystallinity of the silica surface, by changes in the surface protonation state (pH dependence) and by the presence of electrolytes adsorbed in the BIL. The SFG marker bands identified for each of the conditions investigated in the present work are summarized in the scheme presented in Fig. 10.

image file: c9cp02766a-f10.tif
Fig. 10 Schematic representation of the interfacial water and surface silanols OH-populations responsible for the SFG-bands in the experimental and theoretical Imχ(2)(ω) spectra of silica–water interfaces. The O–H groups providing the SFG-signatures are highlighted with blue circles, and the position (frequency, cm−1) and sign (+,−) of their marker-bands in Imχ(2) spectra are also reported.

5 Conclusions

With the present work, by combining microscopic characterization from DFT-MD simulations with theoretical SFG spectroscopy, we provide a rationalization on the molecular organization at silica–water interfaces together with the structure–spectroscopy relationships. By comparing a crystalline quartz/water (QW) interface with an amorphous (fused) silica/water (AW) interface, the dependence of interfacial structural and spectroscopic properties on the degree of surface crystallinity is established, while by adding KCl electrolytes at the QW interface, the effect of ions on the interfacial molecular arrangement is investigated. The evolution of the structure and SFG spectra of QW interfaces with respect to increasing surface deprotonation, i.e., with respect to pH conditions, is also evaluated. The ability to disentangle SFG Imχ(2)(ω) contributions coming from the two distinct interfacial water regions, namely BIL (Binding Interfacial Layer) and DL (Diffuse Layer), is instrumental for an unambiguous assignment of all spectroscopic signatures to the water and/or surface OH-populations.

The water arrangement in contact with silica is thus shown to strongly depend on the degree of crystallinity of the surface, as well as on pH and ionic strength of the aqueous solution. Firstly, the BIL-SFG spectra of silica–water interfaces are found to carry precious marker-bands for the ordered/disordered organisation of water at the interface, which allows one distinguish between crystalline and amorphous silica surfaces: the positive band at 3400 cm−1 is the signature of the crystal-induced ordering of water in the BIL, while the positive band at 3650 cm−1 reflects the loss of such crystal-induced ordering. As already shown in ref. 59, the 3650 cm−1 band at amorphous silica–water interfaces specifically arises from water above hydrophobic patches of the silica surface, rich in Si–O–Si sites. We here find that similar disordered areas within the BIL, where the regular HB-structure of the QW interface is locally broken leading to the appearance of the same 3650 cm−1 band, are also formed at the QW interface either when electrolytes are inner-sphere adsorbed at the quartz surface or when extreme basic pH conditions result in a massive deprotonation of the surface SiOH terminations.

Secondly, both BIL and DL contributions are found to determine the SFG spectral changes with respect to variations in pH conditions. At low pH, when the surface charge and thus the potential across the DL are at their minimum values, the BIL-SFG contribution dominates over the DL-SFG one, and the ∼3400 cm−1 positive band coming from BIL-water molecules donating one HB to in-plane SiOH sites is more intense than the ∼3200 cm−1 band arising from DL-water. On increasing the pH, more and more SiOH surface sites are deprotonated, and the BIL-water molecules donating strong HBs to the appearing SiO sites (instead of to SiOH sites) now provide a positive SFG band at ∼3200 cm−1, reducing the intensity of the ∼3400 cm−1 band initially observed at the neat interface. Simultaneously, the increased negative surface charge leads to a more intense DL-SFG intensity. As a result of both these effects, the positive ∼3200 cm−1 band becomes the most intense in the experimental SFG spectrum at high pH values (above 6.5).

Finally, the negative band below ∼3200 cm−1, appearing in the experimental Imχ(2)(ω) spectra of silica–water interfaces at pH ≤ 7, is shown to arise both from surface out-of-plane SiOH (silanol) groups donating HBs to BIL-water and from O–H groups of BIL-water molecules donating one HB to water in the subsequent bulk layer.

Further studies are however needed to complete the picture presented in this paper, for instance including scenarios where the degree of hydroxylation of the silica surface is varied, as it strongly alters the balance between order and disorder marker-bands. In the same way, the ions’ effect on the interfacial organization has been exploited in the present work only under low pH conditions and for crystalline quartz–water interfaces. An understanding of how the ion adsorption process is modulated by surface hydroxylation, crystallinity and pH conditions is still missing. Moreover, since different interfacial HB-structures are formed at each given combination of these three parameters, the disordering “chaotropic” effect of KCl electrolytes on the BIL-water arrangement, revealed here for crystalline quartz in contact with water in acidic conditions, could possibly be not observed for amorphous silica under other hydroxylation and pH conditions.

Conflicts of interest

There are no conflicts to declare.


This work was performed under Grant ANR DYNAWIN ANR-14-CE35-0011-01 and using HPC resources from GENCI-France Grant 072484 (CINES/IDRIS/TGCC). All the authors acknowledge a France-Berkeley fund for all discussions with Prof Y. R. Shen's group at the University of California at Berkeley, USA. Flavio Siro Brigiano and Louis Potier are gratefully acknowledged for fruitful discussions and sharing code developments in the group.

Notes and references

  1. Y. R. Shen and V. Ostroverkhov, Chem. Rev., 2006, 106, 1140–1154 CrossRef CAS PubMed.
  2. C. Tian and Y. R. Shen, Surf. Sci. Rep., 2014, 69, 105–131 CrossRef CAS.
  3. P. Covert and D. Hore, Annu. Rev. Phys. Chem., 2016, 67, 233 CrossRef CAS PubMed.
  4. P. E. Ohno, S. A. Saslow, H. F. Wang, F. M. Geiger and K. B. Eisenthal, Nat. Commun., 2016, 7, 13587–13591 CrossRef CAS PubMed.
  5. F. Geiger, Annu. Rev. Phys. Chem., 2009, 60, 61 CrossRef CAS PubMed.
  6. V. Ostroverkhov, G. A. Waychunas and Y. R. Shen, Phys. Rev. Lett., 2005, 94, 046102 CrossRef PubMed.
  7. A. Myalitsin, S.-h. Urashima, S. Nihonyanagi, S. Yamaguchi and T. Tahara, J. Phys. Chem. C, 2016, 120, 9357–9363 CrossRef CAS.
  8. E. P. Laetitia Dalstein and E. Tyrode, Phys. Chem. Chem. Phys., 2017, 19, 10343–10349 RSC.
  9. O. Isaienko and E. Borguet, Langmuir, 2013, 29, 7885–7895 CrossRef CAS PubMed.
  10. Z. Yang, Q. Li and K. C. Chou, J. Phys. Chem. C, 2009, 113, 8201–8205 CrossRef CAS.
  11. K. C. Jena and D. K. Hore, J. Phys. Chem. C, 2009, 113, 15364 CrossRef CAS.
  12. M. Sulpizi, M.-P. Gaigeot and M. Sprik, J. Chem. Theory Comput., 2012, 8, 1037–1047 CrossRef CAS PubMed.
  13. A. Cimas, F. Tielens, M. Sulpizi, M.-P. Gaigeot and D. Costa, J. Phys.: Condens. Matter, 2014, 26, 244106 CrossRef PubMed.
  14. M. Pfeiffer-Laplaud, D. Costa, F. Tielens, M.-P. Gaigeot and M. Sulpizi, J. Phys. Chem. C, 2015, 119, 27354–27362 CrossRef CAS.
  15. P. E. Ohno, H. F. Wang and F. M. Geiger, Nat. Commun., 2017, 8, 1032–1040 CrossRef PubMed.
  16. G. Gonella, C. Lutgebaucks, A. G. F. de Beer and S. Roke, J. Phys. Chem. C, 2016, 120, 9165–9173 CrossRef CAS.
  17. S. Ong, X. Zhao and K. B. Eisenthal, Chem. Phys. Lett., 1992, 191, 327 CrossRef CAS.
  18. M. C. Gurau, G. Kim, S. M. Lim, F. Albertorio, H. C. Fleisher and P. S. Cremer, Chem. Phys. Chem., 2003, 4, 1231 CrossRef CAS PubMed.
  19. G. L. Richmond, Chem. Rev., 2002, 102, 2693 CrossRef CAS PubMed.
  20. Z. Yang, Q. Li, M. R. Gray and K. C. Chou, Langmuir, 2010, 26, 16397–16400 CrossRef CAS PubMed.
  21. E. Tyrode and J. F. D. Liljeblad, J. Phys. Chem. C, 2013, 117, 1780–1790 CrossRef CAS.
  22. A. A. Skelton, P. Fenter, J. D. Kubicki, D. J. Wesolowski and P. T. Cummings, J. Phys. Chem. C, 2011, 115, 2076 CrossRef CAS.
  23. A. V. Bandura, J. D. Kubicki and J. O. Sofo, J. Phys. Chem. C, 2011, 115, 5756–5766 CrossRef CAS.
  24. S. Dewan, V. Carnevale, A. Bankura, A. Eftekhari-Bafrooei, G. Fiorin, M. L. Klein and E. Borguet, Langmuir, 2014, 30, 8056–8065 CrossRef CAS PubMed.
  25. T. Ohto, H. Tada and Y. Nagata, Phys. Chem. Chem. Phys., 2018, 20, 12979–12985 RSC.
  26. G. Melani, Y. Nagata, J. Wirth and P. Saalfrank, J. Chem. Phys., 2018, 149, 014707 CrossRef PubMed.
  27. M. F. Calegari Andrade, H.-Y. Ko, R. Car and A. Selloni, J. Phys. Chem. Lett., 2018, 9, 6716–6721 CrossRef CAS PubMed.
  28. M. J. DelloStritto, S. M. Piontek, M. L. Klein and E. Borguet, J. Phys. Chem. Lett., 2019, 10, 2031–2036 CrossRef CAS PubMed.
  29. Y.-C. Wen, S. Zha, X. Liu, S. Yang, P. Guo, G. Shi, H. Fang, Y. R. Shen and C. Tian, Phys. Rev. Lett., 2016, 116, 016101 CrossRef PubMed.
  30. S. Pezzotti, D. R. Galimberti, Y. R. Shen and M.-P. Gaigeot, Phys. Chem. Chem. Phys., 2018, 20, 5190–5199 RSC.
  31. S. Pezzotti, D. R. Galimberti, Y. R. Shen and M.-P. Gaigeot, Minerals, 2018, 8, 305–321 CrossRef.
  32. S. Pezzotti and M.-P. Gaigeot, Atmosphere, 2018, 9, 396–411 CrossRef.
  33. L. B. Dreier, Y. Nagata, H. Lutz, G. Gonella, J. Hunger, E. H. G. Backus and M. Bonn, Sci. Adv., 2018, 4, eaap7415 CrossRef PubMed.
  34. S. N. Wren, B. P. Gordon, N. A. Valley, L. E. McWilliams and G. L. Richmond, J. Phys. Chem. A, 2015, 119, 6391–6403 CrossRef CAS PubMed.
  35. M. Pfeiffer-Laplaud and M.-P. Gaigeot, J. Phys. Chem. C, 2016, 120, 14034–14047 CrossRef CAS.
  36. M. Pfeiffer-Laplaud and M. Gaigeot, J. Phys. Chem. C, 2016, 4866–4880 CrossRef CAS.
  37. N. Allen, M. L. Machesky, D. J. Wesolowski and N. Kabengi, J. Colloid Interface Sci., 2017, 504, 538–548 CrossRef CAS PubMed.
  38. M. Predota, M. L. Machesky, D. J. Wesolowski and P. T. Cummings, J. Phys. Chem. C, 2013, 117, 22852–22866 CrossRef CAS.
  39. J. Hutter, M. Iannuzzi, F. Schiffmann and J. VandeVondele, WIREs Comput. Mol. Sci., 2014, 4, 15–25 CrossRef CAS.
  40. J. VandeVondele, M. Krack, F. Mohamed, M. Parrinello, T. Chassaing and J. Hutter, Comput. Phys. Commun., 2005, 167, 103–128 CrossRef CAS.
  41. A. D. Becke, Phys. Rev. A: At., Mol., Opt. Phys., 1988, 38, 3098–3100 CrossRef CAS PubMed.
  42. C. Lee, W. Yang and R. G. Parr, Phys. Rev. B: Condens. Matter Mater. Phys., 1988, 37, 785–789 CrossRef CAS PubMed.
  43. S. Grimme, J. Comput. Chem., 2006, 27, 1787–1799 CrossRef CAS PubMed.
  44. S. Grimme, J. Comput. Chem., 2004, 25, 1463–1473 CrossRef CAS PubMed.
  45. S. Goedecker, M. Teter and J. Hutter, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54, 1703–1710 CrossRef CAS PubMed.
  46. S. Pezzotti, D. R. Galimberti and M.-P. Gaigeot, J. Phys. Chem. Lett., 2017, 8, 3133–3141 CrossRef CAS PubMed.
  47. J. N. Malin, J. G. Holland and F. Geiger, J. Phys. Chem. C, 2009, 113, 17795 CrossRef CAS.
  48. Y. Duval, J. A. Mielczarski, O. S. Pokrovsky, E. Mielczarski and J. J. Ehrhardt, J. Phys. Chem. B, 2002, 106, 2937 CrossRef CAS.
  49. P. M. Dove, Am. J. Sci., 1994, 294, 665–712 CrossRef CAS.
  50. S. Iarlori, D. Ceresoli, M. Bernasconi, D. Donadio and M. Parrinello, J. Phys. Chem. B, 2001, 105, 8007–8013 CrossRef CAS.
  51. A. Willard and D. Chandler, J. Phys. Chem. B, 2010, 114, 1954–1958 CrossRef CAS PubMed.
  52. J. A. White, E. Schwegler, G. Galli and F. Gygi, J. Chem. Phys., 2000, 113, 4668–4673 CrossRef CAS.
  53. A. Morita and J. T. Hynes, J. Phys. Chem. B, 2002, 106, 673–685 CrossRef CAS.
  54. A. Morita and T. Ishiyama, Phys. Chem. Chem. Phys., 2008, 10, 5801–5816 RSC.
  55. R. Khatib, E. H. G. Backus, M. Bonn, M.-J. Perez-Haro, M.-P. Gaigeot and M. Sulpizi, Sci. Rep., 2016, 6, 24287 CrossRef CAS PubMed.
  56. Y. Nagata and S. Mukamel, J. Am. Chem. Soc., 2010, 132, 6434 CrossRef CAS PubMed.
  57. S. A. C. J. L. Skinner, J. Phys. Chem. A, 2005, 109, 6154–6165 CrossRef PubMed.
  58. T. Ohto, K. Usui, T. Hasegawa, M. Bonn and Y. Nagata, J. Chem. Phys., 2015, 143, 124702 CrossRef PubMed.
  59. J. D. Cyran, M. A. Donovan, D. Vollmer, F. Siro Brigiano, S. Pezzotti, D. R. Galimberti, M.-P. Gaigeot, M. Bonn and E. H. G. Backus, Proc. Natl. Acad. Sci. U. S. A., 2019, 116, 1520–1525 CrossRef CAS PubMed.
  60. T. Joutsuka, T. Hirano, M. Sprik and A. Morita, Phys. Chem. Chem. Phys., 2017, 20, 3040–3053 RSC.
  61. M.-P. Gaigeot, M. Sprik and M. Sulpizi, J. Phys.: Condens. Matter, 2012, 24, 124106 CrossRef PubMed.
  62. G. R. Medders and F. Paesani, J. Am. Chem. Soc., 2016, 138, 3912–3919 CrossRef CAS PubMed.
  63. Q. Wan and G. Galli, Phys. Rev. Lett., 2015, 115, 246404 CrossRef PubMed.
  64. Y. Ni and J. L. Skinner, J. Chem. Phys., 2016, 145, 031103 CrossRef PubMed.
  65. T. Ishiyama and A. Morita, Annu. Rev. Phys. Chem., 2017, 68, 355–377 CrossRef CAS PubMed.
  66. A. Rimola, D. Costa, M. Sodupe, J.-F. Lambert and P. Ugliengo, Chem. Rev., 2013, 113, 4216–4313 CrossRef CAS PubMed.
  67. M. Bouhadja and A. A. Skelton, J. Phys. Chem. C, 2018, 122, 1535–1546 CrossRef CAS.
  68. S.-h. Urashima, A. Myalitsin, S. Nihonyanagi and T. Tahara, J. Phys. Chem. Lett., 2018, 9, 4109–4114 CrossRef CAS PubMed.
  69. A. M. Darlington, T. A. Jarisz, E. L. DeWalt-Kerian, S. Roy, S. Kim, S. Azam, D. K. Hore and J. M. Gibbs, J. Phys. Chem. C, 2017, 121, 20229–20241 CrossRef CAS.

This journal is © the Owner Societies 2019