Molecular Fingerprints of Hydrophobicity at Aqueous Interfaces from Theory and Vibrational Spectroscopies.

Hydrophobicity/hydrophilicity of aqueous interfaces at the molecular level results from a subtle balance in the water-water and water-surface interactions. This is characterized here via density functional theory-molecular dynamics (DFT-MD) coupled with vibrational sum frequency generation (SFG) and THz-IR absorption spectroscopies. We show that water at the interface with a series of weakly interacting materials is organized into a two-dimensional hydrogen-bonded network (2D-HB-network), which is also found above some macroscopically hydrophilic silica and alumina surfaces. These results are rationalized through a descriptor that measures the number of "vertical" and "horizontal" hydrogen bonds formed by interfacial water, quantifying the competition between water-surface and water-water interactions. The 2D-HB-network is directly revealed by THz-IR absorption spectroscopy, while the competition of water-water and water-surface interactions is quantified from SFG markers. The combination of SFG and THz-IR spectroscopies is thus found to be a compelling tool to characterize the finest details of molecular hydrophobicity at aqueous interfaces.

The specific molecular-level organization of water at aqueous interfaces is at the origin of many natural phenomena, ranging from biology and catalysis,[1−8] to pollutant transport in groundwater and mineral dissolution,[9,10] to atmospheric chemistry.[11−14] It is also of crucial importance in electrochemistry, phase-separation processes,[15] and many other technological applications.[16,17]

The balance between hydrophobic and hydrophilic interactions in particular dictates the microscopic arrangement at aqueous interfaces.[18−23] At an aqueous interface, water rearranges in response to the abrupt termination of the bulk water H-bond (HB) network in ways that depend on the local strength of the water–surface interactions. Water maximizes water–surface HBs above a locally hydrophilic surface, while dangling “free” OH groups, pointing out of the liquid phase, are observed in hydrophobic environments. Such dangling OH groups have been experimentally detected via surface-specific vibrational sum frequency generation (SFG) spectroscopy[24] for water in contact with air, oil, organic monolayers, silica, alumina, graphene, and boron nitride.[25−32] The spectroscopic signature of dangling OH groups has been historically interpreted as a molecular/local marker for hydrophobicity.[33] However, recent experimental and theoretical studies have shown that SFG-active dangling OH groups can be detected also at macroscopically hydrophilic surfaces, such as heat-treated silica[29,31,34−36] and 0001-α-alumina.[32,37] These studies pointed out how subtle the concept of hydrophobicity becomes when a molecular-level perspective is adopted. This calls for going beyond the knowledge of dangling OH groups and achieving a deeper rationalization of how the water HB-network rearranges once it is exposed to hydrophobic surfaces.

Progress in this direction has been possible thanks to molecular dynamics (MD) simulations, which showed how the final arrangement of water at complex inhomogeneous interfaces depends on not only the amount of hydrophilic and hydrophobic sites exposed to water but also their spatial distribution over the surface.[18−23] Giovambattista et al.[23] for instance demonstrated that the water surface density is considerably higher in contact with a hydrophobic “patch” surrounded by hydrophilic borders than it is at a purely hydrophobic surface. For water at the boundary with model surfaces with mixed hydrophilic–hydrophobic sites, Erte et al.[19] further found that the exposure of interfacial water to the hydrophobic areas is maximized when the hydrophobic and hydrophilic sites form separate patches over the surface, while it is minimized for homogeneous distributions. In a recent work, combining MD simulations with SFG experiments, some of us have shown that hydrophobic patches are formed on heat-treated silica surfaces. For these systems, the SFG marker-band for dangling OH groups can be observed, although these surfaces are classified as macroscopically hydrophilic from contact angle measurements.[31] When the pretreatment of the silica surface was changed, the attenuation of the free-OH band intensity in the SFG spectra was further shown to correlate with the decrease in the measured contact angles.[31] The correlation between structural, dynamical, and thermodynamic properties of interfacial water was also highlighted by Monroe et al.,[20] concluding that the surface patterns control the dynamics of hydration water and that the hydration water orientational entropy, diffusivity, and H-bonding properties are intrinsically connected.

All these works showed how the spatial variation in surface chemical and geometric topology can tune local and global surface water properties in a surprisingly complex way. The fundamental reason for this complexity is that water molecules have directional interactions and can form a wide variety of networks with neighboring surface groups and water molecules.[21]

In the present study we provide a detailed characterization of the HB-network which is formed by water at a set of distinct “hydrophobic” environments, by combining density functional theory-based MD (DFT-MD) simulations with theoretical SFG and experimental/theoretical THz-IR vibrational spectroscopies. While SFG is a natural probe of the vibrational properties of buried interfaces, we complement it with THz-IR absorption spectroscopy, which has been proved to be an extremely sensitive technique to reveal the intermolecular dynamics of water.[38,39] The strength of the THz low-frequency spectroscopic fingerprints is that they can be directly related to the corresponding hydrogen bond network motifs close to hydrophilic and hydrophobic moieties.[38−43]

We performed DFT-MD simulations (see Methods for all details) for liquid water in contact with (1) air, the prototype hydrophobe, (2) graphene, and (3) hexagonal boron nitride (BN), chosen as examples of non H-bonding surfaces; (4) heat-treated silica (i.e., amorphous silica with low ∼3.5 SiOH/nm2 silanol density[31]) and (5) (0001)-α-alumina, chosen as examples of H-bonding, macroscopically hydrophilic surfaces. The heat-treated silica, with its low degree of surface hydroxylation (∼3.5 SiOH/nm2),[29,31,35] is in particular considered because it is known from our previous work[31] to expose extended hydrophobic patches (with a local silanol density ≤1.5 SiOH/nm2) at the surface despite its macroscopic hydrophilicity (measured by contact angles). Conversely, the (0001)-α-alumina–water interface has a very high density of AlOH termination (15.4 AlOH/nm2), supposedly forming a very hydrophilic surface. However, it also exposes “hydrophobic” sites because less than 1/3 of the aluminols is found H-bonded to water (4.7 HBs/nm2, similar to the value of 4.3 SiOH–water HBs/nm2 formed by the heat treated silica), because of their very basic pKa.[44] The remaining 2/3 of aluminols are dangling −OH groups pointing toward the liquid.[37,44] All five interfaces have been considered at isoelectric conditions, where the BIL (binding interfacial layer[45,46]), in direct contact with the surface/air, is expected to be the only interfacial water layer, directly followed by bulk water. No diffuse layer is present at such isoelectric conditions, according to what we and others have shown in refs (36 and 45−51). These conditions allow us to maximize the SFG sensitivity to the topmost interfacial monolayer (BIL), which is ideal for the purposes of the present investigation.

DFT-MD to Characterize the 2D-HB-Network

We use our previously developed deconvolution scheme to identify the BIL region at the five aqueous interfaces,[45] based on three water–structure descriptors, i.e., the water density profile with respect to the vertical distance from the instantaneous water surface, the average number of HBs/water molecule, and the average orientation within the water HB-network. We find that all five interfaces are composed by a BIL-water monolayer of ∼3.5 Å thickness, directly followed by bulk liquid water.

As shown in Figure , water in the BIL has systematically a higher density than bulk water. There are however differences in the density profiles of the five systems, with the more planar BN–water, graphene–water, and alumina–water interfaces having a more intense and sharper first density peak than the more corrugated air–water and silica–water interfaces (Figure ). As revealed by the simulations, this correlates with the fluctuations of the instantaneous water surface, defined using the Willard and Chandler formalism.[52] The instantaneous surface is found with more oscillations (in time and space) at the more corrugated interfaces, while it is more “flat” at the more planar interfaces. The water spatial ordering thus follows the morphology of the surface. It is worth noting that the oscillations beyond the second peak in the density profiles observed in the bulk region for all interfaces were shown in our recent investigations[53,54] to be due to the limited DFT-MD box dimensions. They disappear for larger simulation boxes, with typically a liquid phase composed by 500 water molecules or more (less than 256 waters are used to model the liquid water in the present simulations).[53] This however has no impact on the structural and vibrational properties of interest, as shown in refs (28, 31, 53, and 54).

Figure 1

Water density profiles calculated from DFT-MD simulations for the five aqueous interfaces. The vertical black dashed line defines the separation between BIL, where the 2D-HB-network is formed, and bulk liquid. For all systems, the x-axis reports the r distance of the water molecules from the instantaneous water surface,[52] which is on average positioned at 3.0–3.5 Å from the four solid surfaces. By construction, r values are positive on the liquid side, while negative values identify waters that protrude out of the liquid phase toward the vacuum/solid.

The high density in the BIL-water interfacial layer, observed for all five systems, is somehow counterintuitive for hydrophobic interfaces: only hydrophilic interfaces would be expected to have such behavior as a result of water strongly interacting with the surface, and thus accumulating in the BIL. As demonstrated in refs (28, 48, 53, and 55), the high BIL-water density observed at the air–water interface is due to the specific HB-structure which is formed in order to maximize the water coordination at the boundary with the vapor. In particular, water at the interface with the air (vacuum) rearranges by forming an extended HB-structure connecting ∼90% of water molecules in the interfacial layer through oriented HBs, preferentially formed parallel to the instantaneous water surface.[28,53] This two-dimensional-H-bonded-network (2D-HB-network) is further made of adjacent 2D-HB-polygons.[54] As long as the water interactions with the other medium are weak enough, water–water HBs remain the dominant driving force for the interfacial water organization, and the formation of a 2D-HB-network in the BIL can be expected. We now demonstrate that the high BIL-water density has the same microscopic origin at all five interfaces, because the 2D-HB-network is systematically formed in the BIL.

Three descriptors are sufficient to reveal the presence of the 2D-HB-network at aqueous interfaces.[53] The first descriptor is the number of in-plane intra-BIL water–water HBs (i.e., HBs formed within the BIL only), which is greater than 1.6 HBs/molecule at interfaces where the 2D-HB-network is formed,[53] while it is lower than 1.1 HBs/molecule for hydrophilic interfaces where no 2D-HB-network is formed.[31,56]

As shown in Figure , more than 1.6 intra-BIL HBs/molecule (denoted “horizontal” because their orientation is parallel to the water surface) are formed in all the investigated systems. The number of these “horizontal HBs” is much greater than the number (≤0.8) of “vertical” HBs. These latter are mostly formed between BIL-water and bulk-water molecules (>85% for all the simulations), while only a fraction (<15%) is formed with H-bonding surface sites (at silica–water and alumina–water only). This behavior is opposite to the one observed for hydrophilic interfaces. For example, when water is in contact with highly hydroxylated (non-heat-treated) silica or quartz, more than 1.7 “vertical” HBs/molecule are formed (≥50% with surface OH-terminations, as detailed in refs (31, 44, and 56)).

Figure 2

Average number of HBs formed per water molecule in the 2D-HB-network, either considering intra-BIL (“horizontal”) HBs, i.e. forming the 2D-HB-network, or considering the remaining “vertical” HBs, which are formed by BIL-water molecules engaged in the 2D-HB-network either with bulk-water molecules (>85%) or with solid surface O–H terminations (<15%). A schematic illustration of “horizontal” (red) and “vertical” (blue) HBs formed by water molecules in the BIL is provided.

The second descriptor used to reveal the 2D-HB-network is the time-evolution of the most extended HB-structure made of intra-BIL (horizontal) HBs. At each MD step, all the possible horizontal HB-structures (by “structure” we mean a noninterrupted HBonded motif that connects several water molecules) made by BIL water molecules are identified, together with the number of water molecules composing each one of these motifs. We are interested in only the largest H-bonded-structure, composed of nmax water molecules (see Figure a that reports the evolution with time of the normalized value nmax/⟨NBIL⟩, where ⟨NBIL⟩ is the average number of water molecules located in the BIL). If a 2D-HB-network is formed in the BIL, nmax/⟨NBIL⟩ is expected to fluctuate around an average value ≥0.85.[53] For all interfaces, we find that nmax/⟨NBIL⟩ oscillates around an average value of 0.9, thus very close to the total number of BIL-water molecules (horizontal blue line), along the 50 ps simulation time. A 2D-HB-network composed by ∼90% of interfacial water molecules not only is observed for all the investigated interfaces but also is stably maintained in time. Illustrative MD-snapshots of the 2D-HB-network formed at the air–water and alumina–water interfaces are presented in Figure b, where the 2D-HB-polygons composing it are highlighted.

Figure 3

(A) Evolution with time (ps) of the number of water molecules (nmax) that are interconnected by intra-BIL HBs into one single 2D-HB-network, normalized by the average number of water molecules in the BIL, ⟨NBIL⟩. (B) MD-snapshots illustrating the 2D-HB-network (orange connections) formed by BIL-water molecules (red oxygens and white hydrogens) at the air–water interface (side and top views) and at the (0001)-α-alumina–water interface (top-view). The instantaneous water surface is also shown in the side view.

The third descriptor is the dynamic behavior of the horizontal intra-BIL HBs, which have been shown to be shorter-lived at interfaces where a 2D-HB-network is formed than any HB in bulk water. The hydrogen bond lifetime (τHB) is computed following the same formalism as used in ref (53) and found to be equal to 0.44 ± 0.03 ps for all the five interfaces. This is 0.8 times smaller than in bulk water (τHB = 0.56 ps obtained with the same computational setup, see ref (53)). The faster HBs dynamics within the 2D-HB-network can be rationalized by the density of HBs (donors and acceptors): this density is indeed found higher within the 2D-HB-network thickness than in an equivalent plane randomly cut in the bulk. This in turn promotes fast HB switching. This rationalization was discussed in ref (53). A similar result has been previously reported for graphene–water and BN–water interfaces by means of ab initio and classical MD simulations.[57,58] The preferential orientation and shorter lifetime of interfacial HBs imposed by the formation of the 2D-HB-network has been shown by recent works to strongly affect the proton hopping mechanism at the air–water interface,[48,59] preferentially occurring via water wires (i.e., HB-chains) running parallel to the surface, as well as the proton conductivity, which is greater at the air–water interface than in bulk liquid water.[60] The same rationalization can be applied also for the high conductivity and lateral hydroxide diffusivity at graphene–water and BN–water interfaces, as demonstrated by Grosjean et al.[61] Furthermore, for all these interfaces sharing the same water–water HB-structure, similar intermolecular pathways for ultrafast energy transfers and vibrational relaxation processes should be expected. As a consequence, we speculate that the presence of a similar 2D-HB-network at air–water and (0001)-α-alumina–water interfaces could be one major reason for the similar (∼250–300 fs) vibrational relaxation dynamics measured via time-resolved TR-SFG and 2D-SFG spectroscopies (at isoelectric conditions, where the BIL is probed by SFG[45,46]). Our hypothesis that the fast BIL-water relaxation measured by TR-SFG is characteristic of the formation of a 2D-HB-network is further supported by a recent study,[62] where we find that addition of ions at (non-heat-treated) silica–water interfaces induces the formation of a 2D-HB-network as well as acceleration of vibrational relaxation dynamics from ∼650–700 to ∼250–300 fs.

In summary, our simulations reveal that the same 2D-HB-network is formed at all five interfaces considered here, thus allowing the maximization of “horizontal” water–water HBs and increasing the connectivity between BIL-water molecules. Water engaged within the 2D-HB-network forms on average almost two intra-BIL HBs with other molecules in the same BIL-layer and less than one HB with water molecules located in the subsequent bulk, resulting in a total 3-fold coordination (sum of blue and red histograms in Figure ). As shown in the figure, such “horizontal ordering” is more marked for water at the interface with non H-bonding surfaces, like BN and graphene (1.9 horizontal HBs/molecule), than at the interface with H-bonding silica (1.7 horizontal HBs/molecule) and alumina (1.6 horizontal HBs/molecule) surfaces. The 2D-HB-network is thus weakened (i.e., less interconnected because of fewer horizontal HBs) by the increase in the number of water–surface HBs formed, with BIL-water making 0.0 HBs/nm2 with graphene and BN surfaces, 4.3 HBs/nm2 with silica, and 4.7 HBs/nm2 with the alumina surface.

HD-SFG Spectroscopy to Quantify Water–Surface Interactions

We now seek the spectroscopic signatures of the 2D-HB-network. We start by focusing on SFG spectroscopy, which is the natural spectroscopic tool to investigate buried interfaces. Theoretical heterodyne detected HD-SFG spectra, where the knowledge on the interfacial water orientation is gained from the sign of the imaginary component of χ(2)(ω),[63] have been calculated for all systems following our previous derivation.[28,56,64] Im(χ(2)(ω)) spectra are plotted in Figure .

Figure 4

(A) Theoretical BIL-SFG Im(χ(2)(ω)) spectra calculated for BIL-water at the interface with boron-nitride (BN), graphene, heat-treated amorphous silica, (0001)-α-alumina, and air. (B) Theoretical BIL-SFG Im(χ(2)(ω)) contribution of BIL-water molecules having one OH group H-bonded to bulk water. (C) Scheme illustrating the correlation between the position of the positive peak in the BIL-SFG Im(χ(2)(ω)) spectra and the strength of water–surface interactions.

The spectra have been calculated considering the contribution of water molecules in the BIL, because our target is to reveal these specific signatures arising from the 2D-HB-network. For the alumina and silica surfaces, the contributions of the O–H terminations should be included as well for a complete modeling of the SFG signal.[56]

Figure A shows very different Im(χ(2)(ω)) spectra for the five interfaces. In order to rationalize these differences we should consider the contributions of the BIL-water OH groups oriented toward the surface, toward the bulk, and parallel to the surface separately. The OH groups oriented parallel to the surface plane are unfortunately not SFG-active (in the common ssp/ppp polarizations) because of their orientation. All the OH groups involved in “horizontal” intra-BIL HBs and forming the 2D-HB-network thus provide a negligible Im(χ(2)(ω)) intensity, although they are the main components of the interfacial structure. By contrast, SFG is sensitive to the interfacial OH groups not engaged in the 2D-HB-network and oriented perpendicular to the surface. From the previous discussions, there is on average one OH group per BIL-water molecule, oriented either toward the air/surface (1/3 of the OHs) or toward the liquid bulk (2/3). In the latter case, these OH groups are systematically HB-donors to bulk water molecules; thus, their signature is expected to be the same in all the SFG spectra. This is confirmed in Figure B, where the SFG contribution of these OH groups is deconvolved and shown to provide a similar negative band for all systems. Thus, the differences observed in the total spectra in Figure A arise solely from the OH groups of interfacial water molecules pointing toward the air/surface. These OH groups modulate the shape of the negative band in the final SFG spectra because of compensation of positive/negative contributions from the two OH populations. More importantly, they systematically provide a positive Im(χ(2)(ω)) peak (the orientation toward the normal to the surface by convention is defined from liquid to solid/air). Such a positive peak is positioned at ∼3700 cm–1 for the air/water interface, where water is in contact with vacuum, i.e., the ultimate hydrophobic medium. This band frequency is then taken as the reference for the hydrophobic character of any system. The strength of the water–surface interactions formed at all other interfaces can be hence ranked by quantifying the extent of red-shift in the positive Im(χ(2)(ω)) peaks with respect to the ∼3700 cm–1 reference. We can hence deduce that BIL-water forms slightly stronger interactions with the nitrogen atoms of the BN than with the carbons of graphene, resulting in a ∼30 cm–1 shift between the relative positive peaks (∼3644 cm–1 for graphene vs ∼3616 cm–1 for BN), in agreement with a previous study.[30] We can also infer that BIL-water makes two kinds of interactions with silica surface terminations, providing two distinct Im(χ(2)(ω)) signatures at ∼3650 and ∼3450 cm–1. As demonstrated in ref (31), the ∼3450 cm–1 band is due to BIL-water OH groups donating HBs to SiOH silanol terminations, while the ∼3650 cm–1 band is due to BIL-water OH groups more weakly interacting with Si–O–Si siloxane bridges. Finally, the in-plane O–H terminations at the alumina surface[64] are able to receive much stronger HBs from BIL-water than silica and all the other surfaces considered, thus leading to the positive band located below 3200 cm–1. One should note here that part of the water–surface interactions, e.g., the HBs possibly received by BIL-water molecules from surface OH terminations, cannot be inferred from the present spectra, because their signatures would be obtained from the SFG activity of surface OH groups that are not taken into account in the present calculations (see e.g. refs (31 and 56) for such signatures). These have already been quantified at silica–water interfaces in ref (56), while they are negligible at the (0001)-α-alumina–water interface because of the pKa values and orientation of the AlOH terminations.[44]

The final ranking on the strength of water–surface interactions formed at the five interfaces obtained from the position of the positive band in Im(χ(2)(ω)) spectra results to alumina > silica > BN > graphene > air, as illustrated by the scheme in Figure C. This is opposite to the ranking obtained for the number of intra-BIL HBs/molecule in Figure , i.e., for the strength of the 2D-HB-network. The water arrangement at the interface is hence triggered by the competition between water–water and water–surface interactions; that is, the weaker the water–surface interactions, the higher the number of intra-BIL HBs in the interfacial layer and the stronger the 2D-HB-network. The formation of the 2D-HB-network in the BIL-water layer is thus a direct indicator of the hydrophobic character of the surface at the molecular level. However, we want to stress that HD-SFG spectra in the OH-stretching frequency domain (2700–4000 cm–1) do not provide any direct signature of the 2D-HB-network.

THz-IR Absorption Spectroscopy to Detect the 2D-HB-Network

The THz difference absorption spectra of aqueous dispersions of thin BN nanoplatelets have been recorded in the 80–240 cm–1 range (see the Supporting Information for further details). In this range, the intermolecular stretching of H-bonded water molecules contribute to the THz spectra discussed hereafter.[65] Three concentrations have been measured, corresponding to 20, 50, and 100 mg/mL, identified as high-, medium-, and low-water-content regimes. In Figure we plot the THz absorption difference spectra (Δα), which were obtained by Δα = , where αsample and αdryBN are the absorption coefficient of a BN aqueous dispersion and that of the corresponding dry BN in the considered dispersion, respectively, and xw is the molar fraction of water in each dispersion. The absorption coefficient α was deduced from the transmitted intensity, using the Lambert–Beer law (see the Supporting Information).

Figure 5

THz difference absorption spectra show Δα as obtained by subtracting the dry BN-nanoplatelets spectrum from the spectrum of each of three aqueous dispersions of BN nanoplatelets, in the high- (top), medium- (middle), and low- (bottom) water-content regimes. Each experimental spectrum (black line) was modeled as a sum of damped harmonic oscillator (red line), as described in the text and the Supporting Information. The damped harmonic oscillators were assigned to BIL-water forming the 2D-HB-network (orange) and to bulk-like water (blue), respectively. The theoretical DFT-MD THz signature of the 2D-HB-network at the BN–water interface, calculated considering the contribution of water molecules forming the 2D-HB-network, is also shown (green), rescaled to the experimental intensity for the sake of comparison. See the Supporting Information for more details on the normalization of the intensity in the experimental spectra.

As described in the Supporting Information and routinely done in our previous THz studies of water systems,[38,39] the Δα spectra are fitted with a sum of damped harmonic oscillators (DHOs). This analysis reveals that the experimental data are well described by one or two DHO(s) . As shown in the top panel of Figure , a DHO alone, centered at ∼180 cm–1 (blue line), dominates the spectrum, when the BN nanoplatelet dispersions are in the high-water-content regime. The spectrum of BN nanoplatelets dispersions in the medium-water-content regime results from two components at 162 ± 2 cm–1 (orange) and 197 ± 4 cm–1 (blue), while a single DHO at 160 ± 2 cm–1 is present in the spectrum of BN nanoplatelet dispersion in the low-water-content regime (orange line, bottom panel). We want to point out that the resonance at 180 cm–1 in the spectrum of the sample in the high-water-content regime is about at the same center frequency as that of the intermolecular stretching band of bulk water.[65,66] We can thus infer that the DHO centered at about 160 cm–1 (orange) component is mostly due to the water in the first interfacial layer, which is the dominant component in the low-water-content spectrum, while the DHO centered at about 197 cm–1 (blue) is due to bulk-like water. The partial contribution of the 197 cm–1 component is increased when additional water layers are introduced in the medium- and high-water-content regimes. We note that the center frequency of this band is similar to that found at ∼193 cm–1 in bulk liquid water at 273.2 K.[38]

Interestingly, the lower-frequency mode at 160 cm–1 has been already observed for water hydrating hydrophobic solutes, including alcohols,[67] clathrate, and semiclathrate hydrates.[38] These latter are cages formed by polygonal water structures around an apolar guest and thus serve as a model system for hydrophobic hydration.[38] We hence tentatively assign the 160 cm–1 mode to the 2D-HB-network at the BN-surface which mostly forms intra-BIL HBs. When more “bulk-like” hydration layers are added (medium-water-content regime), water molecules can now form both intra-BIL and BIL-bulk (or bulk–bulk) HBs, and both the intermolecular modes at ∼160 and ∼197 cm–1 (in orange and blue in Figure , respectively) contribute to the spectrum, analogously to alcohol–water solutions.[67] In the high-water-content regime, the HBs formed between bulk water molecules naturally become the major contributors to the THz spectrum, which is accordingly dominated by one broad band at ∼180 cm–1 (blue), similar to bulk water (see refs (65 and 66) and the Supporting Information).

To support the assignment, we also present in Figure the theoretical THz signal calculated for the 2D-HB-network at the BN–water interface (green line). The spectrum has been obtained by computing the theoretical far-IR spectrum of the water molecules forming the 2D-HB-network , so that the intermolecular stretching of intra-BIL HBs solely contributes in the 80–240 cm–1 frequency range of interest. As shown in the figure, an excellent agreement is obtained between theoretical and experimental spectra (green versus orange line), with the 2D-HB-network signal calculated from the DFT-MD simulation of the BN–water interface giving a maximum centered at ∼160 cm–1, as in the experiments. The ∼160 cm–1 band is hence the direct THz-marker of the 2D-HB-network.

Finally, it is worth noting that water in BN nanoplatelet dispersions shows a narrower line width with respect to that of bulk water at 293 K (see the Supporting Information). Such narrowing is slightly increasing upon decreasing the water content in the dispersion. Any decrease in line width has been ascribable to a reduced number of degrees of freedom, i.e., an entropic signature of a more restricted set of molecular configurations that are available to water molecules in confined environments.[66,68] Here, we suggest that this results from the interaction of water with the BN nanoplatelet surface. Interestingly, the same reduced values for the line widths were observed for water molecules in the hydration shells of alcohols, where the 2D-HB-network was also found.[67]

When DFT-MD simulations are combined with theoretical SFG and experimental/theoretical THz-IR spectroscopies, the molecular-level understanding of the water arrangement at aqueous interfaces and its spectroscopic markers can be obtained. We have demonstrated that, for water at the interface with weakly interacting surfaces, the air–water-like 2D-HB-network[28,53] or HB-wrap[28] is systematically formed, hence maximizing the number of water–water HBs oriented “horizontally” (i.e., parallel to the surface) within the topmost BIL water layer. These HBs have low SFG activity because of their orientation; hence, their spectroscopic signatures are not directly probed in static SFG spectra. However, we have shown that a direct marker band for the 2D-HB-network centered at ∼160 cm–1 is provided by THz-IR absorption spectroscopy experiments. A striking result is that water forms the 2D-HB-network not only at the interface with non H-bonding surfaces, like hexagonal boron nitride, graphene, and air, but also at the interface with H-bonding surfaces, like heat-treated amorphous silica and (0001)-α-alumina, where hydrophobicity arises at the nanoscale level from surface patches which have a local low density of H-bonding sites. Microscopically hydrophobic water, arranged in the 2D-HB-network, is thus present also at macroscopically hydrophilic surfaces. The 2D-HB-network is further found to be weakened by the increase in the strength of water–surface interactions, which is directly measured by HD-SFG spectroscopy.

Our work thus shows that beyond the well-known SFG spectroscopic signatures of the water free O–H groups, the interfacial molecular hydrophobicity is encoded into a planar 2D-H-bond-network formed by the water molecules in the BIL, which is directly revealed by one specific THz spectroscopic fingerprint. Generally speaking, the water HB-networks formed at any aqueous interface result from the subtle balance between horizontal water–water (intra-BIL) and vertical water–surface interactions, which reflects (i) the local density of surface sites that can interact with water molecules, (ii) the strength of the vertical water–surface H-bond interactions (i.e., pKa’s), and (iii) how much the structural pattern(s)/patch(es) formed by the surface sites can be commensurate to the H-bonded structures that can be formed horizontally by interfacial water above the surface. The 2D-H-bond-network revealed by the present MD simulations and THz spectroscopy at various solid–water interfaces (and air–water interface) is the ultimate horizontal ordering that can be obtained. The subtle balance among points i–iii can thus lead to the appearance of the 2D-H-bond-network even at macroscopically hydrophilic interfaces, as shown here.

We further infer that the high connectivity and fast HB-dynamics within the 2D-HB-network revealed here could be used to rationalize the ultrafast vibrational relaxation processes observed at interfaces where a 2D-HB-network is expected.[37,62] In addition, the lateral diffusivity of H3O+/OH– ions should be enhanced within the 2D-HB-network interconnected plane. We argue that such enhanced diffusivity contributes to the high conductivity of interfaces such as graphene–water, BN–water, and air–water, as discussed in refs (60 and 61) as well as to the acidity/basicity of, for example, silica–water and air–water interfaces, as discussed in refs (48, 59, and 69).

As a final note, the complementarity between THz and SFG spectroscopies is proved to be a compelling tool to shed light on the horizontal/vertical HB balance and reveal molecular hydrophobicity, providing direct marker bands for both the formation of a 2D-HB-network (THz) and the strength of water–surface interactions (SFG). This tool opens perspectives in the rationalization and optimization of local hydrophobic effects and therefore in the design of aqueous interfaces.

Methods

For all systems, Born–Oppenheimer DFT-MD simulations have been carried out using the CP2K package.[70,71] The BLYP[72,73] functional plus Grimme D2 correction[74] for van der Waals interactions are adopted, with a combination of Gaussian (DZVP-MOLOPT-SR) plus plane waves basis sets (400 Ry) and GTH pseudopotentials.[75] The nuclei displacements have been predicted through classical Newton’s equations of motions integrated through the velocity-Verlet algorithm, with a time-step of 0.4 fs. Three-dimensional periodic boundary conditions have been applied.

All the simulations have been performed for 50 ps in the NVE ensemble, after equilibration dynamics in the NVE ensemble, with possible rescaling of velocities. All the simulation details are reported in Table . Note that for the silica/water interface, the model for the silica surface is taken from ref (76).

Table 1

Computational Details of the Simulations Performed

simulation	Nmol(H2O)	Natoms(solid)	size (Å3)
air/water[28]	256		19.76 × 19.76 × 35.0
graphene/water	256	180	21.49 × 22.22 × 35.0
boron nitride/water	256	180	21.85 × 22.71 × 35.0
silica/water[31]	116	198	12.67 × 13.27 × 37.0
alumina/water	168	408	14.20 × 16.40 × 35.0

For all structural analyses, the H-bond definition proposed by White and co-workers[77] has been adopted, with O(−H)···O ≤ 3.2 Å and the O–H···O angle in the range of [140-220]°. For the analysis reported in Figure , the HBs formed by BIL water molecules are classified as “horizontal” if they are oriented parallel to the instantaneous water surface within ±30° fluctuation and as “vertical” otherwise.

The SFG signal, coming from the imaginary component of the total resonant electric dipole nonlinear susceptibility χ(2)(ω), has been calculated following the time-dependent method introduced by Morita et al.[78,79] The model from ref (64) is used, largely validated on a variety of water–vapor and water–solid interfaces.[28,31,45,47,48] The THz signal is calculated using the atomic polar tensors (APT)-based method, as developed, applied, and validated in refs (67 and 80). The THz calculated spectrum includes self-correlation terms and all cross-correlation terms between the water molecules (BIL-BIL and BIL-Bulk) (see ref (67)). The 2D-HB-network marker-band at ∼160 cm–1 is due to the sum of the self-correlation and cross-correlation terms between the water molecules that belong to the 2D-HB-network. In ref (67) we have shown that the remaining cross-correlation terms between the 2D-HB-network and the adjacent liquid bulk provide a distinct band centered at ∼190 cm–1, which has a more bulk-like character and therefore cannot be used as a marker for the formation of the 2D-HB-network. The theoretical assignment of such a band is thus not discussed in the present Letter, while details can be found in ref (67).

THz absorption spectra of BN nanoplatelet aqueous dispersions were recorded in the frequency range from 50 to 240 cm–1 by FTIR absorption spectroscopy at room temperature. THz-FTIR measurements were performed using a Bruker Vertex 80v spectrometer equipped with a mercury arc lamp as source and a liquid helium-cooled bolometer from Infrared Laboratories as detector. The samples were placed in a temperature-controlled liquid transmission cell from Harrick with z-cut quartz windows and a 25 μm thick Kapton spacer, and 64 scans with a resolution of 2 cm–1 were averaged for each spectrum, which was smoothed with a 3 point-wide moving average. The apparent center frequencies of the damped resonances have been corrected for the damping factors of the modes to obtain the corresponding undamped peak frequencies, which are reported in the text. The error on the fit parameters is less than 5%.

BN nanoplatelet aqueous dispersions were purchased from Sigma-Aldrich (at a concentration of 20 mg/mL), and all the other concentrations investigated here were obtained by evaporation of the dispersion directly in the measurement cell. The size distribution of the thin nanoplatelets in diluted dispersions was measured by dynamic light scattering (DynaPro Nanostar, Wyatt Technology), resulting in a 90% size fraction with a 140 nm equivalent radius and a 10% size fraction with a 1500 nm equivalent radius for the nanoplatelets. The thickness of the thin BN nanoplatelets is about 2.4 nm, as measured in ref (81). For further details, see the Supporting Information.