Molecular Structure and Modeling of Water–Air and Ice–Air Interfaces Monitored by Sum-Frequency Generation

From a glass of water to glaciers in Antarctica, water–air and ice–air interfaces are abundant on Earth. Molecular-level structure and dynamics at these interfaces are key for understanding many chemical/physical/atmospheric processes including the slipperiness of ice surfaces, the surface tension of water, and evaporation/sublimation of water. Sum-frequency generation (SFG) spectroscopy is a powerful tool to probe the molecular-level structure of these interfaces because SFG can specifically probe the topmost interfacial water molecules separately from the bulk and is sensitive to molecular conformation. Nevertheless, experimental SFG has several limitations. For example, SFG cannot provide information on the depth of the interface and how the orientation of the molecules varies with distance from the surface. By combining the SFG spectroscopy with simulation techniques, one can directly compare the experimental data with the simulated SFG spectra, allowing us to unveil the molecular-level structure of water–air and ice–air interfaces. Here, we present an overview of the different simulation protocols available for SFG spectra calculations. We systematically compare the SFG spectra computed with different approaches, revealing the advantages and disadvantages of the different methods. Furthermore, we account for the findings through combined SFG experiments and simulations and provide future challenges for SFG experiments and simulations at different aqueous interfaces.


INTRODUCTION
Liquid water−air and solid-state water (ice)−air interfaces play a critical role in many biological, chemical, physical, and atmospheric processes as well as environmental science. Water and ice interfaces provide unique platforms for evaporation and sublimation of water, respectively, controlling the phase diagram of water. 1,2 These interfaces interact with small molecules, governing important atmospheric and environmental processes. For example, the water interface interacts with CO 2 molecules in the air emitted from fossil fuels and captures them, increasing the acidity of sea, 3 and organic molecules are absorbed onto ice interfaces, leading to the growth of aerosols in the air. 4 From Faraday's measurement, it is known that the ice−air surface has a thin premelted water layer, 5 causing ice blocks to stick together. Moreover, ice can also adhere to biological tissue (e.g., fingers or tongues) when placed in contact with the ice surface. Although such processes critically affect our daily life, health, and climate change, molecular-level understanding of the structure of water and ice interfaces is far from complete. To probe water and ice interfaces, surface sensitive techniques such as atomic force spectroscopy, 6,7 X-ray spectroscopy, 8,9 and sum-frequency generation (SFG) spectroscopy 10−13 have been used. Among these techniques, SFG is a powerful technique, as it can probe the vibrational mode of water. 10,14−17 An SFG signal, which is described by the second-order nonlinear susceptibility (χ (2) ), generated by infrared (IR) and visible pulses, is enhanced when the molecular vibration at interfaces is resonant with the IR frequency. As such, vibrational spectra of the ice or water interface are acquired. The frequency of a vibrational mode of water provides information on the local environment of the interfacial water molecules. For example, the O−H stretch frequency of water is sensitive to the hydrogen-bond (H-bond) strength; a low (high) O−H stretch frequency indicates that a water molecule forms a strong (weak) H-bond. 11,18 The SFG intensity spectra (|χ (2) | 2 ) for O−H (O−D) stretch mode at the H 2 O(D 2 O)−air interface was first measured by Shen and coworkers, using the homodyne-detection technique, 19,20 and has been repeatedly reported by other researchers. 21− 23 The O−H stretch spectrum contains two features; a broad 3000−3500 cm −1 band and a sharp 3700 cm −1 peak (Figure 1a). The broad 3000−3500 cm −1 feature evidences the H-bonding of the interfacial water, while the sharp 3700 cm −1 peak arises from the dangling O−H groups sticking out of the water interface. 17 One of the challenges associated with interpreting these spectral features originates from the fact that different contributions to the real and imaginary components of the complex χ (2) can interfere in the |χ (2) | 2 SFG intensity spectra.
In 2005, this realization led to a new technique called phasesensitive SFG measurements, 24−26 which was soon thereafter applied to the water−air interface. 26−28 The first reports of phase-sensitive SFG spectroscopy date back to the early 1990s. 29 Phase-sensitive SFG has, in its implementation using femtosecond pulses, also been called heterodyne detected SFG. In this technique, Im(χ (2) ) can be determined from the phase of the signal, relative to a signal of known phase. Im(χ (2) ) can be either positive or negative, reflecting the absolute orientation of the transition dipole moment and thereby the molecules (i.e., pointing up or down, relative to the surface plane). The early experimental Im(χ (2) ) data at the water−air interface showed positive 3700 cm −1 , negative 3400 cm −1 , and positive 3200 cm −1 bands for the O−H stretch mode (Figure 1b). 24,26,27,30−37 However, in 2015, Tahara and co-workers corrected such a positive−negative−positive feature of the O−H stretch mode to a negative−positive signal by recognizing the contamination in the reference signal ( Figure 1c). 38−40 Whether an SFG spectrum of water at the water−air interface should show a positive− negative−positive feature or negative−positive feature is still under debate. 40−44 This discussion highlights the importance of a critical check between experimentally measured and simulated SFG spectra, as simulations can provide an ideal condition.
Indeed, SFG spectra simulations have been shown to be a very powerful means to obtain molecular-level information about the interfacial water structure, 45−49 interpret different spectral features, 50−52 disentangle the layer-by-layer contributions to the SFG spectra, 53−55 and examine the effects of coupling between the water molecules 56 and between the layers 57 on the SFG spectra. The strength arises from the ability of simulations to separate the SFG spectra into the contributions from each water molecule or intra/intermolecular coupling. Accordingly, there has been increased interest in the simulation of SFG spectra of interfacial water.
It has turned out that the SFG spectra of water are sensitive to the models used in the molecular dynamics (MD) simulations such as the force field models 58−60 and calculation methods in the density functional theory (DFT). 58,61,62 Here, we provide an overview of the recent progress on molecular modeling of SFG spectra of water, with particular focus on the progress after 2015. Reviews on modeling SFG spectra before 2015 can be found elsewhere. 63−65 This review article will discuss not only the O− H stretch mode but also the other modes such as H−O−H bending mode and librational mode of water at the water−air interface.
On the basis of the successfully modeling of the SFG spectra at the water−air interface, the microscopic structure of the ice−air interface has been investigated using combined SFG experiments and simulations. This surface has been measured by Shen and co-workers 66,67 as well as Shultz and co-workers, 68−71 but the physical origin of the SFG spectra of the O−H stretch mode and their bulk contribution were not clear. 72 Recently, the SFG spectra of water at the ice−air interface have been computed and compared with experimental data. The direct comparison of the SFG spectra between experiments and simulations provided a rich understanding of the underlying physics of the surface premelting of the ice−air interface. 73,74 These will be also discussed.
In this review, we provide an overview for the simulation protocols for the SFG spectra by introducing the SFG principles (sections 2.1.1, 2.1.2, 2.1.5, and 2.1.6), an efficient algorithm for response function formalism (section 2.1.3), simulation settings (section 2.1.4), various force field models (section 2.2), ab initio MD simulation (section 2.3) and the comparison of SFG experiments and simulations (section 2.4). We also introduce useful simulation techniques to characterize the microscopic interfacial water structure (section 2.5). Subsequently, we will explain the O−H stretch mode of water at the liquid water−air interface (sections 3.1.1 and 3.1.3). The effects of the intra/ intermolecular vibrational couplings (section 3.1.2), the differential SFG spectra due to the temperature changes (section 3.1.4), and the orientation of free O−H groups (section 3.1.5) will be discussed. We then move to the discussion of the water bending (section 3.2) and librational modes (section 3.3). We further focus on the ice−air interface in section 4. Unresolved questions and the future challenges of SFG at these interfaces, including simulation of time-resolved and two-dimensional (2D) SFG and X-ray SFG, will be given in section 5. section 6 provides a conclusion.

SIMULATION FRAMEWORK
2.1. Simulation Protocols for Computing SFG Spectra 2.1.1. Principles of SFG. An SFG signal is generated by overlapping IR and visible pulses, an IR pulse generates a coherent state between the ground state and vibrationally excited state, and a sequential visible pulse up-coverts this vibrationally excited state to the virtual state, producing a signal with a frequency equal to the sum of the IR and visible frequencies (see Figure 2a). Such operations constitute a second-order nonlinear optical process.
The second-order susceptibility (χ cba (2) ) vanishes in a centrosymmetric medium. 75 This can be shown as follows. The second-order nonlinear polarization in the c-direction is given by where E a is the electric field applied in the a-direction, and a, b, and c = x, y, and z. In the centrosymmetric medium, the flipping of the pulse directions such as E → −E should give rise to the flipping of the polarization P → −P, namely; To satisfy eqs 1 and 2 simultaneously, P c (2) should be zero. As such, a SFG signal excludes the contribution from the region where the medium is centrosymmetric; in other words, the molecular response is nonzero in the region where the centrosymmetry will be broken. We note that this is true in the electric dipole approximation. Bulk quadrupole contributions to the SFG intensity remain possible in centrosymmetric media.
How can the sensitivity to the breaking of the centrosymmetry be connected to the interface specificity? When a molecule A is present at the interface of mediums B and C, the molecular interaction between molecule A and medium B differs from the molecular interaction between molecule A and medium C. In most cases, the probabilities of up-and down-orientation of the molecule A are no longer the same, as the A−B and A−C interactions are different (Figure 3a), leading to noncentrosymmetry at the interface. In contrast, when a molecule A is present Figure 2. (a) Double-sided Feynman diagram representing SFG optical process. |g⟩, |v⟩, and |vir⟩ denote the key state of the ground state, vibrationally excited state, and electronically excited virtual state, respectively. (b) Time-independent vs (c) time-dependent scheme for SFG optical process. In the time-independent picture, the transition energy (ℏω gν ) is known, while it cannot be determined with much accuracy in the time-dependent picture because of the Heisenberg uncertainty principle. in the bulk of medium B, there is no difference in the up-oriented and down-oriented probabilities (Figure 3b). In this case, the centrosymmetry holds. As a result, SFG can detect the molecular response near the interfaces. In fact, an SFG signal is dominated by the response of interfacial water molecules within the 1−2 topmost layers at the water−air interface. This will be discussed in more detail in sections 2.1 and 3.
2.1.2. Time Correlation Function Approach. The total SFG signal, χ (2) (ω), is the sum of the resonant contribution (χ (2),R (ω)) and nonresonant contribution (χ (2),NR ); Because the origin of the χ (2),NR signal has not been fully understood, the χ (2),NR signal cannot be computed directly from the simulation by now. Nevertheless, from the SFG measurement, it is empirically known that the χ (2),NR signal can be assumed to be independent of the frequency ω in a limited frequency region. Because the χ (2),R (ω) signal arises from the vibrational mode of the molecule, it can be accessed by MD simulation. Calculation of the χ (2),R (ω) spectrum can be classified into two categories: time-independent approach and time-dependent approach (see parts b and c of Figure 2, respectively). 45 In the time-independent, transition energy (ℏω gν )-specific approach, MD simulation is often carried out using the rigid-body force field models of water such as SPC/E 76 and E3B models, 77,78 and vibrational frequencies, ω gν , are often computed through the frequency mapping technique. 79−81 Here, a rigid-body model denotes the model with the two fixed O−H bond lengths and the one fixed H−O−H angle. In contrast, the time-dependent approach requires the flexiblebody model of water.
In the time-dependent approach, which is currently being used in most simulation groups, an SFG spectrum χ cba (2),R (ω) at the polarization directions of a, b, and c for IR, visible, and SFG pulses, respectively, is given by the Fourier transformation of the time-correlation function R cba (2) (t): where Q(ω) is the quantum correction factor 82−85 and A bc (t) (M a (t)) denotes the bc-element of the polarizability (acomponent of the dipole moment) of the system at time t. ⟨···⟩ represents the thermal average. The dipole moment and polarizability of the system are the sum of the molecular dipole moment (μ) and molecular polarizability (α): The first and second terms of the right side in eq 6 represent the autocorrelation term for ith molecule and cross-correlation term between ith and jth molecules, respectively. The localized vibrational feature and the effects of the intramolecular vibrational couplings are included in the first (autocorrelation) term, while the contribution of the intermolecular vibrational couplings leading to the delocalization of the vibrational chromophores 86,87 is included in the second (cross-correlation) term. 54,88 This means that one can realize the contribution of the intermolecular vibrational couplings to the SFG spectra by evaluating the second term of the right side in eq 6. On the other hand, to examine the effects of the intramolecular vibrational couplings on the spectra, one needs to compare simulations with HOD and H 2 O (or D 2 O) molecules. These are detailed in section 3.1.2. Finally, we note that eq 6 indicates a vanishing SFG signal in a centrosymmetric medium. Here, let us think about the rotation of the molecular dipole moment and polarizability. The rotation of the molecular dipole moment by 180°changes the sign (μ → −μ), whereas the molecular polarizability is invariant to the 180°r otation (α → α). As such, in a centrosymmetric medium, ⟨αμ⟩ is equal to zero.
2.1.3. Truncated Time Correlation Formalism. In eq 6, we decomposed the time correlation function of the systemdipole moment (M) and system-polarizability (A) into the autocorrelation terms and cross-correlation terms of the molecular dipole moment (μ) and molecular polarizability (α). The numbers of the autocorrelation terms ) and in the cross-correlation terms (α i, bc (t)μ j, a (0)) are N and N(N-1), respectively, where N is the number of the molecules. This indicates that the number of the crosscorrelation terms is much larger than that of the autocorrelation terms. In contrast, the contribution of the autocorrelation terms to the spectral amplitude is larger than, or at least comparable to, the contribution of the cross-correlation terms. 54,59,89,90 The cross-correlation terms are largely canceled out due to the phase and frequency mismatch between the different vibrational chromophores in the large number sampling limit. This is intuitive, as chromophores that are physically far apart, should not contribute, on average, to the cross-correlation term. In a finite sampling scheme, however, the cross-correlation terms can provide substantial noise. As a result, obtaining the converged cross-correlation terms requires a much longer trajectory than obtaining the converged autocorrelation terms.
To evaluate the cross-correlation terms accurately with a limited length of an MD trajectory, we have introduced the truncated time correlation function formalism. 54 In this approach, by introducing the cross-correlation cutoff of r t , we take into account the cross-correlation of the ith and jth chromophores in the time correlation calculations, when the i−j distance is shorter than the cutoff radius r t . Inversely, we exclude the i−j cross-correlation when the i−j distance is larger than r t . Such a time-correlation function can be given by , ,  (8) where g t (r ij (t);r t ) is the function to control the cross-correlation terms, and r ij (t) is the distance between chromophores i and j at time t. Note that eq 7 is the same as eq 6 when r t is infinite. For water, it would be suitable to set r ij as the intermolecular distance of the oxygen atoms of water (see the schematics in the left panel of Figure 4). The choice of r t is in principles arbitrary. However, it is convenient to determine r t based on the minima of the radial distribution function (RDF) of water, as it allows us to connect the hydration shell of water with an SFG signal. For example, when the SFG signal computed from the truncated correlation function formalism (eq 7) is converged with r t = 3.3 Å (5.6 Å), the SFG signal contribution can be unambiguously attributed to the contribution of the water molecules within the first (second) hydration shell of water (see the right panel of Figure 4). A similar analysis has been used in ref 86. 2.1.4. Slab Model. To model water−air (water−vapor) or ice−air (ice−vapor) interfaces in the simulation cell, we need to insert a vacuum region into the cell and isolate the condensed water molecules in the simulation box. As a result, both the topmost and bottommost layers of water are exposed to the vacuum. Typically, the dimension of the cell perpendicular to the interface (the length of the box along z-axis) is larger than those along x-and y-axes. This is called the slab model. In the slab model, the water−air interface is repeated due to the threedimensional (3D) periodic boundary condition.
Because the two interfaces in the simulation cell are oppositely oriented, one needs to calculate the SFG responses at the two interfaces separately. To do so, we introduce a screening function g sc (z i ) in eq 6; (10) z i denotes the z-coordinate of the water molecule i. The region of |z i | > z c2 determines the regions where the dipole moment and polarizability are included in the time correlation function (eq 9), and the region of z c1 < |z i | ≤ z c2 defines the region where the dipole moment and polarizability are switched off smoothly. The origin point is set to the center of mass of the system. For the SFG calculation in the slab model, one needs to pay attention to the following two issues. First, because of the periodic boundary condition, the one surface artificially interacts with the other surface in the image cell through the air region. In particular, when the macroscopic surface dipole is generated at the interface, it may affect the interfacial structure, as the electrostatic interaction is long-range. To avoid such an artifact, we need to have an air region that is thick enough to make this artifact negligibly small in the 3D periodic boundary condition. The in-house data indicate that a thickness of the vapor region in excess of 1.5 nm would be sufficient for removing the artificial interactions between the original and image cells at the water− air interface. The thickness of the air region may need to be enlarged when one simulates the water−air interface with the presence of ions or water-charged interfaces. An alternative approach for removing this artifact is to employ the twodimensional Ewald summation method for computing the electrostatic energy in a two-dimensional periodic cell 91,92 instead of the standard Ewald summation for the threedimensional periodic cell.
Second, the water layer in the slab model should be thick enough to have a region where the net orientation of water is zero (centrosymmetric). This centrosymmetric water region separates the two interfaces, thus being essential for SFG calculations. The thickness of the centrosymmetric water region can be examined by calculating the axial profile of dipole moments of water. The net dipole moment should be zero in the middle of the layer. In this region, physical properties such as the density and the H-bond network should be very similar to those in bulk water. The density profile of water in the slab model simulated with the POLI2VS model is displayed in Figure 5. This indicates that a 2 nm thick water slab contains a 1 nm thick region of the bulk water. Note, however, that the situation is much more complicated when the surface is charged; it is known that charged surfaces generate two noncentrosymmetric layers, the Stern layer and the diffuse layer. The thickness of the diffuse layer easily reaches ∼10 nm, much longer than that of the Stern layer. Thus, the SFG spectra calculation at a charged interface requires, in principle, a very long side of the cell along. 93,94 Although the electrical double layer model is well-established for describing the electrostatic potential near a charged interface, the effects of the Stern and diffuse layers on the SFG signal have not been fully tested in a cell with its longest side of ∼10 nm. Only an approximate approach to estimate χ (3) -contributions is available. 95 where ssp (ppp) denotes s-(p-), s-(p-) and p-polarized SFG output, visible input, and IR input, respectively. Here, e a denotes the component along a direction of the unit polarization vector of the optical field, L aa (ω) represents the frequency-dependent Fresnel factors. ω SFG , ω Vis , and ω are the frequencies of SFG, visible, and IR, respectively, and β SFG , β Vis , and β are their incident/reflection angles formed by the surface normal and SFG, visible, and IR beams, respectively, as shown in Figure 6a. Note that ω SFG = ω Vis + ω.
Generally, the SFG and visible beams are not in resonance with a vibrational transition, while the IR beam is resonant with a vibrational transition of the molecule. Thus, the Fresnel factors of visible and SFG are featureless, while the Fresnel factors of IR shows significant frequency dependence (see Figure 6b). That means, to obtain the χ (2),R (ω) spectrum, one needs to account for the local field effects of the IR beam as described by the Fresnel factors.
The Fresnel factor is defined as a macroscopic or mesoscopic property, which describes the relationship of the electric field in two different media. However, in the few surface layers which give rise to the SFG signal, the dielectric constant is not welldefined. To obtain the local field at the surface, a microscopic local field correction might be needed. This factor is included in a factor, ε′, and is incorporated in the equations of the Fresnel factor. Therefore, there are two methods to calculate Fresnel factors. One is ignoring the microscopic local field correction, in this case, ε′ adopts the dielectric constant value of one of the two media. The other is considering the microscopic local field correction. To obtain an estimation of ε′, a slab model is presented in ref 98. They obtain ε′ = ε(ε + 5)/(4ε + 2) by assuming the dielectric constant of the upper medium is 1. Figure 6b shows the spectra of the IR Fresnel factor, e z L zz , calculated with ε′ = ε and ε′ = ε(ε + 5)/(4ε + 2). The difference of the SFG spectra of the experimental data using these two corrections are shown in section 3.1.1.   with the information on the angle Ω formed by the free O−H group and the surface normal. The amplitudes of the free O−H stretch mode A xxz and A zzz in the χ xxz (2) and χ zzz (2) spectra, respectively, are related to ⟨cos 3 (13) in the slow motion limit, 99 and is linked to in the fast motion limit. 100 For the slow (fast) motion limit, the decay of the orientational memory of the free O−H group is much slower (faster) than its dephasing time. So far, it is recognized that the slow-motion limits can reproduce the data in a self-consistent manner. r is given by the ratio of the transition polarizability α̇ζ ζ /α̇ξ ξ , where ξ and ζ denote the directions parallel to and perpendicular to the O−H bond. This r can be easily calculated using ab initio calculations. In the case of the water's free O−H group, r can be calculated from the water trimer, giving 0.15 at the B3LYP/aug-cc-pVTZ level of theory.
Together with the ratio of A xxz and A zzz of 0.43, one can get ⟨ cos Ω⟩/⟨cos 3 Ω⟩ = 1.55 in the slow motion limit (1/⟨cos 2 Ω⟩ = 1.55 in the fast motion limit).
To compute the average angle ⟨Ω⟩ from ⟨cos 3 Ω⟩/⟨ cos Ω⟩ or ⟨cos 2 Ω⟩, one needs to assume an orientational distribution. This will be discussed in section 3.1.5.

Force Field Model of Water: Empirical vs ab Initio-Based
Above, we explained the basic formalism to compute the SFG spectra. Currently, most research groups use the timecorrelation function approach with the truncated timecorrelation formalism. Thus, there is no difference in the calculation methodology. Nevertheless, one can see a variety of the SFG spectra of water below. These differences arise solely from the model chosen to describe water. When one uses an accurate model of water, one can reproduce the correct experimentally measured spectrum of water. Thus, the accuracy of the water model is essential for accurate prediction of the water spectrum. The model of water for computing vibrational spectra can be categorized into the force field model and ab initio model. In this section, we present an overview of different force field models of water.
So far, many force field models of water have been proposed such as simple point charge (SPC) model, 101 extended simple point charge (SPC/E) model, 76 three-point transferable intermolecular potential (TIP3P) model, 102 and four-point transferable intermolecular potential (TIP4P) model. 102 However, these SPC or TIP water models assume a fixed geometry of the water molecules (rigid-body model), which do not have O− H stretch vibrations or H−O−H bend vibrations. Therefore, the rigid-body models can be used only for computing the spectra for librational modes and intermolecular vibrational modes such as the H-bond stretching mode in the time-dependent approach, as is mentioned in section 2.1.2. 103−106 For computing the intramolecular vibrational mode spectra, one needs to use flexible body models of water. The intramolecular potential in the flexible body model was designed for an isolated water molecule and provides access to the explicit molecular vibrations of water such as the O−H stretch mode and H−O−H bending mode. The modulation of the intramolecular potential due to the surrounding molecules, e.g., H-bonds, is accounted for through intermolecular interactions of molecules. However, because the O−H stretch mode is primarily affected by the anharmonicity of the potential, it is known that the flexible-body water model with a harmonic potential term for the O−H stretch vibration (the lowest-order term for vibrations) such as flexible simple point charge water (SPC/Fw) model 107 is insufficient for reproducing the line shape of the O−H stretch spectrum. Modeling with higherorder (anharmonic) terms is required for accurately reproducing vibrational spectra.
The water models with an explicit intramolecular potential can be classified into "empirical" and "ab initio-based" models. The empirical force field models have been parametrized based on fitting to experimental data, while the ab initio-based force field models have been parametrized based on fitting to the intramolecular potential energy surface obtained from ab initio calculation. 108 Rahman and co-workers proposed the first "empirical" flexible-body water model, including the anharmonic terms for the intramolecular degrees of freedom, in a way that the water model can reproduce the experimental observables. 109 Such empirical flexible-body force field models of water have further been improved 110−112 and have been used not only for IR spectra calculation 113,114 but also for SFG spectra calculation. 52,53 The typical criterion for an "empirical force field model" for water's O−H stretch mode is the accurate prediction of the O−H stretch frequency, without nuclear quantum effects. 37,[52][53][54]113,115 Because the O−H covalent bond is stiffer without nuclear quantum effects than with nuclear quantum effects, 116,117 the O−H stretch frequency without nuclear quantum effects should exceed the experimentally observed value if the potential energy surface for the O−H stretch mode is properly described. In other words, the typical signature for "ab initio-based" force field model is the overestimation of the O−H stretch mode by 150−200 cm −1 when nuclear quantum effects are not included. [58][59][60]118 This explains why the inclusion of the nuclear quantum effects or corresponding frequency shift provides a frequency of 3700 cm −1 free O−H peak for several ab initio-based force field models such as the polarizable water model for intramolecular and intermolecular vibrational spectroscopies (POLI2VS) 119 and the many-body potential for water (MB-pol), 60,121 while the empirical force field model underestimates the free O−H frequency with the nuclear quantum effects. 43 Above, we explained the force field models of water, based on the classification of "rigid-body" vs "flexible-body". In addition to this classification, a commonly used classification is the "polarizable" vs "nonpolarizable" force field models. In the polarizable force field model, induced effects are explicitly taken into account through charge or dipole moments. We have thus four categories for the force field models of water: rigid-body nonpolarizable model, rigid-body polarizable model, flexiblebody nonpolarizable model, and flexible-body polarizable model. These are depicted in Figure 7a. The ab initio-based force field models constitute a specific class of the flexible-body, polarizable model.
The ab initio-based force field models are constructed such that the predicted conformational energy of water molecules with the force field model matches that obtained from the highlevel ab initio calculation results. Thus, one cannot adjust parameters intentionally to reproduce the experimental data with the ab initio-based force field model, unlike the empirical force field model.  126 Luber, 127 Cho, 44 Klein/Bourget, 62,128 Saalfrank, 129,130 Chandra, 131 and Car/Selloni 132 groups as well as Nagata/Bonn group 58,133−136 have used AIMD simulation for computing the SFG spectra. Because the forces acting on the atoms are calculated within the electronic structure theory in the AIMD simulation, one does not need force field models. However, the quality of the data is governed by the different levels of the exchange-correlation (XC) functionals. In fact, within the framework of the DFT, the many-body interaction between electrons is mapped to a single-particle problem with a practice form of the XC functional approximation. By including more and more constraints and corrections, the accuracy of the XC functional will be improved by climbing the so-called Perdew's metaphorical Jacob's ladder of XC functionals where each rung introduces more physics. 137 The Jacob's ladder is shown in Figure 7b. The XC functionals in the higher rung are especially valuable for water because the structural, electronic, and dynamical properties of water are improved dramatically. 138 Furthermore, because the description of the van der Waals (vdW) interaction is insufficient in the DFT approach, the vdW correction schemes have been often used. 139 We have included the similar Jacob's ladder for the different level of vdW corrections in Figure 7b as well.
Because force field modeling is not required in AIMD simulations, AIMD simulations are very suitable to model complicated interfaces such as liquid−solid interfaces. 97,125,135 On the other hand, the computational cost of AIMD simulation is at least 100 times larger than that of MD simulation using force field models, while a very high level of AIMD could be even by 100000 times more expensive than classical MD simulation. AIMD has been well investigated for predicting the structure and dynamics properties of bulk water, 138,140−148 while the accuracy of the interfacial water properties predicted by AIMD simulations has not been systematically examined not only at the generalized gradient approximation (GGA) level of theory with the empirical vdW corrections 61 but also at the meta-GGA 149−151 and hybrid-GGA level of theory 152−160 with the nonlocal vdW corrections. 161,162 For example, we have reported that the values of the surface tension of water predicted using different AIMD simulations vary substantially due to the vdW interactions, while most of the force field models predict reasonable values for the surface tension of water. 163,164 Likewise, the predicted details of the free O−H group at the water−air interface changes drastically when we use different XC functionals and vdW correction schemes. 61,161 The SFG spectra of water at the water−air interface are also different when we use different XC functionals. In fact, the choice of the XC functional has a substantial influence on the frequency of the O−H stretch mode. The Perdew−Burke− Ernzerhof (PBE) GGA functional 157 is known to show a huge red-shift of the O−H stretch vibrational modes. 58, 165 Becke-Lee−Yang−Parr (BLYP) GGA functional 153,166 provides a reasonable frequency of the O−H stretch mode, but without nuclear quantum effects, indicating that the reasonable frequency arises from the cancellation of the missing nuclear quantum effects and self-interaction correction error. As is mentioned section 2.2, the accurate ab initio technique should predict a higher O−H frequency than the experimental data when the nuclear quantum effects are not included. Because a frequency shift of the different DFT method may cause a misleading peak assignment of water, one needs to be careful about the choice of the XC functional and calculation method. This will be discussed in section 4.2. Furthermore, we note here that Car−Parrinello-type AIMD using a fictional electron mass may cause the systematic red-shift, unlike Born−Oppenheimer type AIMD simulations. 167 A small electron mass should be chosen, in particular for computing the vibrational spectra.

Comparison of Simulation with Experiment
Simulated SFG data can be compared with the experimental data of either homodyne-detected SFG intensity signal (|χ (2) (ω)| 2 ) or heterodyne-detected data (Im(χ (2) (ω))). However, the comparison of the simulation and experimental data is not straightforward due to the nonresonant contribution χ (2),NR . As is mentioned in section 2.1.2., the origin of the nonresonant contribution is unclear and thus cannot be directly addressed by MD simulation. This nonresonant contribution complicates the comparison of simulation and experiment in the homodynedetected measurement. On the other hand, homodyne detection is much easier to perform experiments than heterodyne detection, and the signal-to-noise ratio is much better: heterodyne detection requires interference between the local oscillator and the SFG from the sample, unlike homodyne measurement. As such, although the heterodyne detection is much better suited for comparing the simulated and experimental SFG spectra, the heterodyne measurement requires much more efforts in detecting an SFG signal when the signal is weak. Below, we overview how the simulation and Chemical Reviews pubs.acs.org/CR Review experimental spectra can be compared in the homodyne detection and heterodyne detection. For the heterodyne detection, the total SFG signal, Im(χ (2) (ω)), is the sum of the resonant contribution (Im(χ (2), R (ω))) and nonresonant contribution (Im(χ (2), NR )), as is mentioned in eq 3; The approximation of the frequency-independent χ (2),NR is valid when considering a limited frequency range. The Im(χ (2),R (ω)) contribution can be addressed through MD simulation, whereas the Im(χ (2),NR ) contribution cannot be addressed directly from MD simulations. Thus, comparison of the experimentally obtained Im(χ (2) (ω)) spectrum and the simulated Im(χ (2),R (ω)) spectrum can be achieved by shifting the simulated response of Im(χ (2),R (ω)).
On the other hand, the homodyne-detected signal can be written as , The term of |χ (2),NR | 2 only shifts the whole spectra in a vertical direction because it is frequency independent, while the crossterm of (χ (2),R (ω))*χ (2),NR + c.c. is frequency-dependent, modulating the spectral lineshapes significantly. As such, the homodyne-detected data gives more uncertainty on the agreement of the simulation data and experimental data than the heterodyne-detected data. 168 For a comparison of the homodyne-detected experimental data with the simulation data, χ (2),NR is often treated as a parameter for the fitting 169 because the simulation does not provide a value for χ (2),NR . This gives rise to an uncertainty in the consistency between experiment and simulation. Furthermore, the information on the absolute orientation of the transition dipole moment is automatically included in the sign of Im(χ (2),R (ω)), whereas it is missing when one measures |χ (2) (ω)| 2 .
2.5. Tools for Analysis of the Water Interface 2.5.1. Averaged vs Instantaneous Liquid Interface. One important step in analyzing the simulation results is defining the interface. In one approach, the molecular information is computed for the normal of the macroscopic surface and then averaged over the ensemble of the conformations of water. This is named averaged interface description. Because the surface normal of the interface is defined in the lab-frame, the averaged interface description allows us to correlate the microscopic structure of interfacial water with macroscopic observables like the SFG response. The depth profile of the density of water at the water−air interface with the averaged interface description is plotted in red in Figure  8a.
The axial profile of the water density at the water−air interface is known to be well-described with the sigmoid function: where ρ represents the density of water in the bulk in a zeropressure condition, z G the z-coordinate of the Gibbs dividing surface, and δ the thickness parameter for the interface. Note that the thickness parameter δ reflects the magnitude of the nanoroughness caused by capillary waves. Because capillary waves are highly dependent on the long-range vdW interaction and the size of the simulation cell (surface area), 163,170,171 δ depends not only on the water force field models but also on the simulation cell size (see Figure 8b). 172 As such, one needs to pay special attention to the comparison of δ for different simulations.
If one wants to minimize the effect of the surface roughness on the density profile and extract the intrinsic density fluctuation of water at the interface, the instantaneous interface description has to be used. 173 In this approach, the molecular information is computed for the axis, which is defined as the normal to a specific position of the temporal interface. The instantaneous interface description employs the coarse-grained view of water; a water molecule is treated as a single sphere. The collection of the spheres generates the gradient of the instantaneous density field of the nuclear positions. The location of the interface is then defined as the isodensity surface corresponding to half the bulk density of water. 173 As such, one can construct the nonuniform instantaneous interface. Because this instantaneous interface has already included the surface nanoroughness, the computed depth profile along the normal of the instantaneous interface can exclude the effects of the surface roughness to some extent. The depth profile of the interfacial water obtained with this instantaneous description plotted in blue in Figure 8a. This figure demonstrates that the instantaneous interface description visualizes the density oscillation spread over at least two water layers at the water−air interface, which is smeared out in the averaged interface description due to the surface roughness. Furthermore, this instantaneous interface description has the potential to provide rich information on the water interface. 57,124,174,175 For example, the Kuḧne 174 and Gaigeot groups 124,175 have pointed out that in the topmost layer defined by the instantaneous liquid interface description, the H-bond network is expanded in the plane of the water−air interface. To find this ensemble, one needs to define the geometrical criteria. 178 The O−H group defined by the optimized free O−H group criterion should include the whole 3700 cm −1 positive SFG feature and exclude the negative SFG feature in the frequency region of ω ≤ 3600 cm −1 . Furthermore, the free O−H groups captured by this criterion should reproduce (1) the experimentally obtained orientational parameter, d ≡ ⟨cos Ω⟩/ ⟨cos 3 Ω⟩ of 1.55 ± 0.02, 179 and (2) the free O−H group lifetime of 1.1 ps. 180 An O−H group can be categorized into free O−H when it cannot find any H-bond partners. The H-bond partner has been defined as the geometrical relation of the water dimer, using an intermolecular O···O distance R or intermolecular O···H distance r, and an angle of H−O···O or O−H···O, which are denoted as β and θ, respectively (see Figure 9a). 181 An H-bond is formed when R < R c and β < β c for the R-β definition, when R < R c and θ > θ c for the R-θ definition, and when r < r c and θ > θ c for the r-θ definition. Our goal is to obtain the optimal values of (R c , β c ), (R c , θ c ), and (r c , θ c ) for a free O−H group subensemble, which reproduces the positive 3700 cm −1 SFG feature.
We optimized the criteria for the geometrical criteria (R c , β c ) in a manner that the selected O−H groups can reproduce the positive 3700 cm −1 peak in the SFG spectra of O−H in D 2 O, and obtained the optimal criterion of (R c , β c ) = (3.5 Å, 50°). 182 The obtained criterion is different from the commonly used criterion (R c , β c ) = (3.5 Å, 30°) for defining a H-bond. 183 Figure 9b displays the free O−H group contributions to the SFG spectra with both criterion. The criterion of (R c , β c ) = (3.5 Å, 30°) contains not only the free O−H peak at 3700 cm −1 but also a lower frequency tail around 3600 cm −1 , while the criteria of (R c ,  Figure 9c. For the criteria using other pairs of the angle and distance, we obtained the optimal value of (R c , θ c ) = (3.5 Å, 110°) and (r c , θ c ) = (3.0 Å, 90°). 182 These can be compared with literature suggesting the criteria for defining the free O−H groups. Campen and co-workers proposed the criteria of (R c , θ c ) = (3.5 Å, 140°), based on the comparison of time-resolved SFG data with the SPC/E and TIP4P/2005 water model simulation. 184 This criterion is slightly tighter for H-bonding and thus slightly looser than our criteria for defining a free O−H group of (R c , θ c ) = (3.5 Å, 110°). 182 Richmond and co-workers also use their own definition of (r c , θ c ) = (2.5 Å, 150°), 185 which is apparently much tighter for H-bonding and thus much looser for defining free O−H group than our criteria of (3.0 Å, 90°).
The H-bonded O−H/free O−H group definition is frequently used for categorizing interfacial water molecules based on specific hydrogen bonding criteria. Such categorization provides a powerful means to interpret the SFG spectra. 28,50,51,90,121 However, it depends critically on the precise (geometric) definition of H-bonded O−H/free O−H groups, the different force field models, and different functionals used in AIMD simulations. Checking the consistency of the interpretation, which is partially done in the body of this review, is essential for an accurate evaluation of simulated data.
2.5.3. Sign of Heterodyne-Detected Spectra. A water molecule has three intramolecular vibrational modes: two O−H stretching modes and one H−O−H bending mode. The relation between the sign of their Im(χ (2),R (ω)) feature and the orientation of the transition dipole moment is mode-dependent. To understand these mode-dependent relations, one needs to think about the phase of the transition dipole moment and the transition polarizability. Here, we discuss the relations by taking into account two typical conformations of the water molecules at the water−air interface.
One typical conformation of water at the water−air interface is DAA (or DA)-type, where one of the two H atoms in a water (2) )), 89 as well as the simulated Im(χ ssp,fOH ) with the optimal criteria of (R c , β c ) = (3.5 Å, 50°) 182 and (R c , β c ) = (3.5 Å, 30°). 184 Table 1. The other typical conformation is a DDA-type, where both two H atoms have H-bond partners. Here, for a representative case, we assume that these two O−H stretch modes have the same frequency and thus the transition dipole moment for the symmetric (antisymmetric) mode is parallel (perpendicular) to the H−O−H angular bisector. When H−O−H angular bisector is parallel to the z-axis, the symmetric stretch is SFG active and the antisymmetric stretch mode is SFG inactive. However, the value of ∂μ z /∂Q · ∂α xx /∂Q for the antisymmetric mode becomes positive once an H−O−H bisector of a water molecule is no longer parallel to the surface normal and the O−H group pointing more down to the bulk form a stronger H-bond than the other O−H group. This antisymmetric mode will be revisited at the water−air interface (section 3.1) as well as the ice−air interface (section 4.2).
This table also indicates that the positive (negative) sign of the molecular susceptibility for the bending mode arises from the transition dipole moment pointing up (down). However, one should be very careful about the orientation of the transition dipole moment for the bending mode, as the permanent dipole moment and induced dipole moment have different phases in the condensed phase. 186 In fact, the SFG simulation indicates that the positive sign of the SFG peak corresponds to the transition dipole moment pointing down to the bulk, at odds with this table. This will be detailed in section 3.2. (ω)) features can be determined using the heterodyne-detected SFG technique. Until 2015, many experimentally reported SFG spectra show a positive−negative-positive feature: a positive 3700 cm −1 free O−H peak, a negative 3500 cm −1 H-bonded O−H peak, and a positive 3200 cm −1 peak. Because the physical origin of the positive 3200 cm −1 peak was unknown, MD simulations have been conducted to reproduce these SFG spectral features and address the molecular origin of the 3200 cm −1 peak. The reproduction of the 3200 cm −1 positive SFG band has been reported by the Morita group, 52,115,187 the Skinner group, 51,188 and the Sulpizi/Sprik/Gaigeot group. 123 In contrast, MD simulation with the POLI2VS model as well as AIMD simulation in Nagata/Bonn group did not reproduce such a positive 3200 cm −1 peak. 58 In 2015, Yamaguchi 38 and Tahara group 39 reported that the presence of the 3200 cm −1 peak arises from the artifacts of the measurement; the SFG spectra of water at the water−air interface should not show a positive 3200 cm −1 peak. These results were in agreement with the data of POLI2VS and AIMD SFG simulations. 58 Since 2016, many papers have been published confirming the absence of the 3200 cm −1 peak. Experimentally, Wen and co-workers showed a similar SFG spectrum, 189 but Tian and co-workers noted that the careful assessment of the phase for D 2 O may be needed. 40 Theoretically, Paesani and co-workers calculated the SFG spectra with the MB-pol model 60  Sulpizi and co-workers reported that the 3200 cm −1 positive band might arise artificially from insufficient convergence of the SFG time correlation function due to the shortage of the trajectory. 122,123,191 These lessons illustrate that proper setting of the SFG calculation and sufficient sampling of the trajectory Table 1. Signs of Molecular Susceptibilities of DAA and DDA-Type Water Molecules a a α xx represents the xx-component of molecular polarizability of a target water molecule, μ z the dipole moment along z-axis, and Q the normal mode. Calculation of the transition dipole moment and transition polarizability for the highlighted water molecule was done in the gas phase.

THE LIQUID WATER−AIR INTERFACE
Chemical Reviews pubs.acs.org/CR Review are essential to access the correct SFG spectra, in addition to the accurate modeling of the force fields.
Most of the simulated spectra reported after 2015 including those from the Gaigeot 124 and Kuḧne 57 groups do not contain a 3200 cm −1 band. However, the spectral shapes of these simulated SFG spectra of water differ significantly. The selected simulated spectra together with the experimentally measured spectra 39,40 (with and without the Fresnel factor corrections using ε′ = ε(ε + 5)/(4ε + 2) 98 and ε′ = ε 40 ) are shown in Figure  10. First, we focus on the ratio of the positive 3700 cm −1 peak amplitude and negative 3400 cm −1 peak amplitude. Simulation data show the scattered ratio of the amplitudes of the 3700 and 3400 cm −1 peaks. Next, the presence of the 3600 cm −1 shoulder peak does not agree. For example, this shoulder peak is missing in the data of Skinner and co-worker 190 as well as the data of Cho and co-workers. 44 We shall revisit this shoulder peak in the next section.
We then focus on the spectra simulated with the POLI2VS and MB-pol models. The simulated spectrum with the POLI2VS model by Nagata and co-workers shows a double-peak feature for the negative 3300−3400 cm −1 signal, while the spectrum with the MB-pol model by Paesani and co-workers shows a single negative peak at 3400 cm −1 . A comparison of these data with the experimental spectra after correction for the Fresnel factor indicates that the MB-pol data by Paesani and co-workers underestimates the 3300 cm −1 negative amplitude, while the POLI2VS data by Nagata and co-workers enhances the 3300 cm −1 negative amplitude excessively. However, the agreement of the simulation and experimental data does highly depend on the dielectric constant used in the Fresnel factor correction. Further examination of the accuracy of the simulated spectra will be made in section 3.1.4.
The quantum simulations using centroid MD technique with the MB-pol model 192 and the partially adiabatic centroid MD technique with the quantum version of the flexible TIP4P model 117 show very broad free O−H peak unlike the experimental data, 57,60 which is associated with the intrinsic problem of the centroid MD simulations. 193 The improved quantum description 116,194−196 of the SFG spectra remains challenging. Furthermore, studying other polarization combinations such as sps, is also essential for checking the agreement between the simulation and experiment. The simulated spectra presented here have been assumed to arise from the interface, while a possible contribution of, for example, the quadrupole terms to the SFG spectra 197 cannot be strictly excluded.
3.1.2. Signatures of Intra-/Intermolecular Vibrational Coupling in O−H SFG Spectra. It is well-recognized that bulk IR and Raman spectra of neat H 2 O differ significantly from those of isotopically diluted water, in particular at ∼3300 cm −1 . 11,56 This evidences that the O−H stretch chromophores of water are intramolecularly and/or intermolecularly coupled with other O−H stretch chromophores in the bulk. 11,56,87,198−202 The intramolecular vibrational coupling includes the Fermi resonance between the overtone of the H−O−H bending mode and O−H stretch mode and the vibrational energy splitting of the nearly identical two O−H stretch chromophores. In contrast, the intermolecular vibrational coupling allows the O−H stretch chromophores to be delocalized over several vibrational chromophores (vibrational exciton). Here, a question is how such intra-/intermolecular vibrational couplings affect the SFG spectra of water.
We quantify the effects of the intra-/intermolecular vibrational couplings on the O−H stretch SFG spectra of water. Here, we present the following computed three SFG spectra of water: the spectra with both the inter-and intramolecular vibrational couplings of water (H 2 O spectrum), with intramolecular vibrational couplings but without intermolecular vibrational couplings (H 2 O spectrum without intermolecular vibrational couplings), and without any couplings (HOD in D 2 O spectrum). The intramolecular vibrational couplings can be removed by replacing H 2 O by HOD molecule, while the intermolecular vibrational coupling of the O−H vibrational chromophores can be controlled by neglecting the crosscorrelation terms in eq 9, so that  58,190,60,44,122,57,40, and 39, respectively. For the spectra of "Nagata (Sim.)" and "Paesani (Sim.)", the frequency is red-shifted to compensate for the missing nuclear quantum effects, while it is not done in the spectrum of "Morita (Sim.)", as the nuclear quantum effects are already included. "Paesani (CMD Sim.)" represents the centroid MD simulation data without shifting. For the data of "Skinner (Sim.)", the frequency mapping technique is used, which includes such effects implicitly. For the data of "Kuḧne (Sim.)", the partially adiabatic centroid MD simulation was used. For the experimental data (Tian and Tahara), the dotted, dashed, and dashed−dotteded lines denote the measured spectrum, the spectrum with Fresnel factor correction of ε′ = ε(ε + 5)/(4ε + 2), and the spectrum with the correction of ε′ = ε, respectively. 40 The arrows indicate the maximum peak position of the positive free O−H peak and center-of-mass peak position of the negative H-bonded O−H peak, while the filled circle indicates the crossing point. Note that we normalized all the spectra at the peak of the H-bonded OH at around 3400 cm −1 .

Chemical Reviews
The simulated spectra are displayed in Figure 11a. 89 The corresponding experimental data are presented in Figure  11b,c. 203 First, we shall focus on the impact of the intermolecular vibrational coupling on the SFG spectra of water by comparing the H 2 O spectrum (in red in Figure 11a) and H 2 O spectrum without intermolecular vibrational coupling (in green). The striking difference of these spectra can be seen in the H-bonded O−H stretch region; the intermolecular vibrational coupling enhances the amplitude of the negative peak at ∼3300 cm −1 . A similar trend can be seen in the IR and Raman spectra of bulk water. 56,86,88 The delocalization of the vibrational chromophores contributes to the spectral signature at ∼3300 cm −1 both in the bulk and at the interface.
Subsequently, we focus on the impact of the intramolecular vibrational coupling on the SFG spectra. This can be seen from the comparison of the H 2 O spectrum without intermolecular vibrational couplings (in green in Figure 11a) and the HOD spectrum (in blue). Figure 11a shows that the 3600 cm −1 shoulder is missing in the "HOD in D 2 O" spectrum, while it is present in the H 2 O spectrum without intermolecular vibrational couplings. This signature is not clear in the experimental data of Figure 11b but is clear in Figure 11c. The computational analysis using the H 2 O spectrum without intermolecular vibrational couplings shows that the shoulder peak arises from the intramolecular vibrational coupling not from the intermolecular vibrational coupling. 89 Shen and co-workers have pointed out this shoulder arises from the intramolecular coupling, and subsequently 31 Skinner and co-workers also assigned this shoulder peak to the antisymmetric mode induced by intramolecular vibrational coupling of the DDA water molecules. 28  Note that a different assignment has also been presented, 203 indicating that the discussion of the molecular origin of the 3600 cm −1 shoulder peak is not completed.  58,135 Here, we tested the quality of the spectra with the AIMD trajectories at the different levels of theory. For computing the O−D in H 2 O SFG spectra, we used the surfacespecific velocity−velocity correlation function formalism 58 and the frequency axis was shifted by using a factor 0.96 50,59,119 to take the nuclear quantum effects into consideration. The results with various DFT levels are displayed in Figure 12. 204 All the simulated O−D stretch mode spectra show a sharp positive peak at 2550−2750 cm −1 and a broad negative peak centered at 2300−2500 cm −1 . Compared with POLI2VS data 89 and experimental data, 37 we conclude that the hybrid-GGA revPBE0 functional 156 combined with Grimme's D3(0) vdW correction 205 is the best for reproducing the SFG features at the water−air interface. For the GGA XC functionals (BLYP, PBE, revPBE), the simulated spectra commonly show the red-shift of the spectra. This means that within the GGA, the O−D covalent is artificially weakened owing to the self-interaction error. 206 This can be corrected by including the Hatree−Fock exchange term, resulting in the hybrid GGA functionals. Except the overall red-shift, the revPBE together with Grimmes's D3(0) vdW corrections reproduce the experimentally obtained line shape well.
Next, we focus on the details of the spectral shapes. The meta-GGA (SCAN, 151 M06-L, 149 and B97M-rV 150,207 ) commonly lacks a plateau region in 2620−2650 cm −1 . These are well reproduced in the BLYP-family and PBE-family. In particular, the SFG feature simulated with the M06-L functional is largely distorted despite the excellent reproduction of the molecular geometry. Overall, the revPBE0+D3(0) reproduces the trend of the experimental SFG spectrum. These indicate that the SFG spectra at the water−air interface can be a critical check for reproducing the description of the XC functional. In fact, the water−air SFG includes the O−H groups with an H-bond donor as well as the free O−H group, unlike in the bulk water, where a free O−H group is missing. The collections of the SFG spectra can be a good guide for simulating the SFG spectra of interfacial water. 62 3.1.4. Temperature Dependence. Modulation of the topmost layer of water at the water−air interface due to the variation of the physical conditions such as temperature and pressure is linked to the stability of the interface. For instance, the surface stability changes with temperature, resulting in the variation of the surface tension of water. 208 It has been, however, unclear how the molecular structure of interfacial water is modulated by temperature. The variation of the SFG spectra with temperature may provide information on the variation of the microscopic structure of interfacial water. Furthermore, understanding the differential SFG spectra is useful to estimate the heat contribution in the time-resolved SFG data; 36,180,209,210 a pump pulse excites the O−H stretch modes of water, which results in the heating of the local environment after ∼700 fs. 200 To disentangle the vibrational energy relaxation from the heat contribution in the pump−probe data, one needs to know the temperature-dependent SFG spectra. 36,180,209−211 Parts a and b of Figure 13 display the temperature-dependent SFG spectra calculated with the POLI2VS model 59 and the MBpol model, 121 respectively. Both data indicate that the SFG amplitude at ∼3700 cm −1 is rather insensitive to the temperature Chemical Reviews pubs.acs.org/CR Review change, while the SFG amplitude at ∼3300 cm −1 decreases significantly with increasing temperature. In addition, the MD simulation results indicate that the positive shoulder at ∼3600 cm −1 becomes featureless with increasing temperature. Let us focus more on individual SFG features. We start with the 3300 cm −1 SFG signature. It is known from the bulk Raman spectra that the 3300 cm −1 O−H stretch mode is sensitive to the temperature. 11 Thus, we compare the SFG spectra with the IR (I IR ), Raman (I Raman ), and constructed IR/Raman spectra ( I I IR Raman ), 59,89,212,213 which are plotted in parts c, d, and e Figure 13, respectively. The 3300 cm −1 O−H stretch mode is very sensitive to the temperature, in particular in the VV-Raman spectra. As a result, the constructed IR/Raman spectra also vary drastically at ∼3300 cm −1 . The variation of the constructed spectra shows a similar trend to that of the SFG spectra, indicating that the variation of the SFG spectra of water is not characteristic for the interfacial ordering of water.
Next, we discuss the insensitivity of the 3700 cm −1 peak to temperature variations. This seems to indicate that the free O− H group is unchanged even when the temperature changes. However, this is not the case. From the analysis using the free O−H definition, we found that the surface density of the free O−H group increases with increasing temperature, enhancing the SFG amplitude. This seems reasonable, given that the increase in temperature enhances the breaking of the H-bonds. At the same time, the average angle between the free O−H group and the surface normal increases, resulting in the reduction of the SFG amplitude. As such, the effects of the enhanced average angle on the free O−H peak amplitude and the effects of the enhanced surface density canceled out, making the free O−H SFG peak amplitude insensitive to the temperature.
Are these predicted trends in agreement with the experimental data? To examine this, we collect the differential SFG spectra from Δχ ssp (2),R (ω) ≡ χ ssp (2),R (ω; T = ΔT + T 0 ) − χ ssp (2),R (ω; T = T 0 ) with the POLI2VS model, 59 the MB-pol model, 121 and the CRK model including the short-range damping function, 115,214 and then compare them with the experimentally measured differential spectra by Bonn et al. 36,59 and Tahara et al. 215 The data is shown in Figure 14. Note that the experimental data of refs 36 and 215 arise from pump−probe SFG measurements with a 1.5 ps delay time. These differential spectra ΔIm(Δχ ssp (2),R (ω)) can be reasonably assumed to be governed by the heated spectra, as it is estimated to take ∼700 fs for converting from the O−H stretch vibrational energy to the heat. 200 The POLI2VS simulation, the MB-pol simulation, and the experimental data show a similar trend: a broad positive band at 3300 cm −1 and a large negative peak around 3550 cm −1 in addition to a small positive peak at 3670 cm −1 . However, one can see a small difference between the POLI2VS simulation, MB-pol simulation, and the experiment; the 3300 cm −1 ΔIm(Δχ ssp (2),R (ω)) feature is more pronounced in the simulation than in the experiment. In contrast, a striking difference exists between the spectrum simulated with the CRK model by Morita and co-workers 115,214 and the experimental data. This demonstrates that the differential spectra are more sensitive to the quality of the molecular modeling than the static spectra. This is consistent with the claim that the higher-order vibrational spectra are more sensitive to the details of modeling. 216 3.1.5. Orientational Distribution of the Free O−H Group. Above, we show how sensitive the SFG spectrum of water at the water−air interface is to the choice of the force field model. Does this mean that the predicted interfacial structure of water is similarly sensitive to the choice of the model? To examine this, we recently computed the average angle of the free O−H groups with respect to the surface normal (⟨Ω⟩), with the help of the free O−H group definition (see section 2.5.2). 182 These are summarized in Figure 15a.
The average angles of the free O−H groups ⟨Ω⟩ are very similar between the empirical nonpolarizable force field models (SPC/E 76 and TIP4P/2005 217 ) and ab initio-based polarizable force field models (POLI2VS 118 and MB-pol 218 ). This indicates that the simulated structures of the interfacial free O−H groups are very similar among these models. Likewise, it has also been reported that the H-bonded structure of the interfacial water is similar. 175 These are surprising, as the molecular dipole moments of water differ at the water−air interface and in the bulk and thus it is expected that the interfacial water structure predicted with the polarizable force field models differs from that predicted with the water models designed for the bulk water (nonpolarizable force field models). The insensitivity of the structure of the interfacial water molecules indicates that the Hbond strengths of water molecules at the water−air interface and in the bulk are similar despite the apparent difference in the magnitude of the dipole moment of water. This is consistent with the SFG probe of the H-bond strength of the interfacial water molecules. 28 These suggest that, at least for understanding the fundamental structure of the interfacial water at the hydrophobic interfaces, a nonpolarizable force field model can be used.
Subsequently, we turn our focus on the distribution of angle Ω. This is plotted in Figure 15b. An important issue understood from the simulation results is that the distribution function of the  228 It is also essential when discussing the possible overlap of the bending overtone and O−H stretch modes because the overlap dictates vibrational energy relaxation of H 2 O. 180,[209][210][211]215,229,230 So far, the H−O−H bending mode of the interfacial water molecules has been experimentally probed by three groups: the Benderskii and Bonn groups performed the homodyne-detected SFG measurements, 90,231 while the Tahara group did the heterodyne-detected measurement. 232 For the SFG simulations, the Nagata/Bonn group, 90 Skinner group, 169 Sulpizi group, 122 and Paesani group 121 computed the bending SFG features, with different techniques. The experimental SFG intensity spectra are summarized in Figure 16a. All the data commonly show a maximum intensity around 1680 cm −1 , a drastic intensity change in 1630 cm −1 < ω < 1660 cm −1 and a moderate decrease in the intensity with increasing frequency above 1700 cm −1 . In the frequency region below 1630 cm −1 , the intensity spectra are almost flat. These data can be reproduced with the simulation data using a parameter of χ (2),NR (see section 2.4.) However, when we turn our focus on Im(χ ssp (2),R (ω)) spectra, one can see a non-negligible difference between the experimental and simulation data. Figure 16b displays the simulated spectra of Nagata/Bonn group 90 and Skinner   90 and Tahara's group (heterodyne detection). 232 All the spectra are offset by increments of 0.4 for clarity. (b) Experimental and simulated Im(χ ssp (2),R (ω)) spectra of the water bend mode with the POLI2VS model 90  (ω)) peak to the H−O−H bisector pointing down arises from the opposite sign of the transition dipole moment and transition polarizability with respect to the bending mode of water. 90 The simulation data indicates that interfacial water molecule can have both up-and down-orientations for the H−O−H bisector. This mixed positive−negative peak for the bending mode is consistent with the homodyne-detected bending mode data at the negatively charged surfactant−water interface. 233 On the other hand, Tahara and co-workers attributed the positive peak in the experimental data to the quadrupole contributions in the bulk; according to their claim, the SFG signal arises from the higher-order contribution than the interfacial dipole contribution, and thus the bulk SFG does not contain the information on the interfacial water. 232 Probably, the best way to check whether this arises from the bulk or interface, measuring the SFG spectra at the watercharged lipid interfaces, and examining the variation of the frequency and sign of χ (2) . If an SFG signal is dominated by the bulk, the sign and frequency of χ (2) will not be changed due to the different lipids with positive/negative charges. Such a trial has been partially done by Benderskii group, 233,234 whereas disentangling χ (2) contribution from χ (3) contribution at the charged interfaces 94,95,235,236 has not been performed yet, except in one very recent study. 237 These should be a future challenge.
Probing the H−O−H bending mode provides complementary information on the structure and dynamics of interfacial water. For example, in section 2.5.1, we mention the H-bond network expanded in the direction parallel to the water−air interface. 124,174,175 These H-bonded O−H groups are, however, SFG inactive, when we probe the O−H stretch mode. On the other hand, these in-plane O−H groups can be probed through the bending mode, as the direction of the bending mode transition dipole moment is not parallel to that of the O−H stretch. As such, analyzing multiple vibrational modes proves to be particularly valuable.

The Librational Mode
The librational mode of water is composed of three different reciprocating rotational motions (see Figure 17a). The frequency of the librational mode is sensitive to the H-bond connectivity, thus to the molecular conformation of water; the librational mode frequencies of bulk liquid water, microporous amorphous ice, annealed ice, and crystalline ice are 670, 770, 810, 840 cm −1 , respectively. 238,239 As such, measuring the librational modes of the interfacial water can clarify how the Hbond network is formed at the aqueous interfaces.
The first and only measurement so far for the librational mode of the interfacial water has been conducted by Campen and coworkers. 240 The SFG intensity spectrum at the water−air interface is shown in Figure 17b. This indicates that the peak frequency is ∼850 cm −1 , much higher than the peak frequency of 670 cm −1 in the bulk liquid water and close to the frequency of 840 cm −1 in crystalline ice. Because the peak frequency may be shifted by the nonresonant contribution in the SFG intensity spectra (see section 2.4.), a question here is whether this peak frequency can be considered as the characteristic frequency of the librational mode of the interfacial water.
To examine this, we have performed the MD simulation with several force field models (POLI2VS, TIP4P, and SPC/E). 120 Note that one can use the rigid-body force field models such as TIP4P and SPC/E model for calculating the vibrational spectra of librational mode (see the discussion in section 2.2), as the librational mode does not need any intramolecular vibrational degrees of freedom, like the O−H stretch and H−O−H bend modes. 103,104,106 The simulated Im(χ ssp (2),R ) in Figure 17c indicates that the peak is located at ∼760 cm −1 , not 850 cm −1 . However, once we fit the simulated χ ssp (2),R (ω) data to the homodyne-detected experimental data of SFG intensity |χ (2) (ω)| 2 via eq 16 using χ ssp (2),NR as a parameter for the fitting, one can see excellent agreement between simulation and experiment, having a peak frequency of 850 cm −1 in the |χ (2) (ω)| 2 spectra (see Figure 17b). This clearly demonstrates that the nonresonant contribution χ ssp (2),NR may shift the peak frequency by ∼100 cm −1 . 120 This demonstrates that the SFG heterodyne detection of the water librational mode is required to confirm the peak frequency of the librational mode. Furthermore, although the discussion so far is limited to the peak frequency of the librational motion, an analysis using MD simulations would be also needed to understand the properties of the librational motion of the interfacial water in the near future.
Here, we note a pioneering work on SFG simulations by the Space group. 47,48,241 They have already predicted the whole spectral range including bending and librational mode in 2005, 47 while the spectral shape of the well-known O−H stretch mode differs from the reported literature; they reported the same sign One can expect that the reorganization of the topmost water layer can be captured by probing the free O−H stretch mode. In fact, in the first SFG measurement of the ice−air interface, Shen and co-workers discussed the variation of the free O−H stretch intensity in the ssp and ppp SFG signals. 66,67 In these works, the basal plane of the crystal ice was measured. The change in the ssp/ppp signals with varying temperature was attributed to the variation of the free O−H angle (see section 3.4.2) and it was concluded that disordering of the topmost layers starts above 200 K. 66,67 In contrast, it was not clear how the topmost layer reorganizes their structures by varying temperature (see Figure  18a). SFG spectra contain not only the information on the ssp/ ppp ratio but also the frequency shift of the free O−H group. Below, we outline how the frequency shift of the free O−H peak can be used to unveil the structural change of the topmost water layer of the ice surface, by combining the POLI2VS MD simulation with different temperature as well the experimental data. Parts b and c of Figure 18 display the Im(χ ssp (2),R ) spectra of free O−H stretch mode at ice−air interface obtained from the heterodyne-detected SFG measurement and the POLI2VS MD simulation, respectively. 119 When focusing on the temperature variation of the peak frequencies for the free O−H stretch mode, we can see that both experiment and simulation exhibit a blueshift of the peak frequencies from 3694 to 3700 cm −1 when the temperature increases from 185 to 245 K (see Figure 18d). 119 Here, a question is how such a frequency shift occurs for the free O−H groups, although the free O−H group has no H-bond partner. There are two possible mechanisms that may alter the free O−H group frequency; one is the intramolecular vibrational coupling between the free O−H stretch mode and the other O− H stretch mode, 28,242 while the other mechanism is the interconversion of the DAA water species to DA water species. 243,244 Let us focus on the first mechanism. The intramolecular vibrational coupling alters the free O−H frequency by 10 cm −1 , when the other H-bonded O−H stretch frequency is changed from ∼3400 cm −1 to the O−D frequency of ∼2500 cm −1 (∼900 cm −1 shift). 28 However, the temperature rise from 150 to 245 K changes of the H-bonded O−H stretch mode frequency only 150 cm −1 , which is then expected to redshift the free O−H peak frequency only by 1.2 cm −1 . 119 Apparently, this 1.2 cm −1 shift caused by intramolecular vibrational coupling is insufficient to account for the 6 cm −1 shift observed in the SFG spectra.
Subsequently, we examine the second mechanism of the interconversion of the DAA species to DA species.  Chemical Reviews pubs.acs.org/CR Review This interconversion of DAA to DA is rather complicated; below 200 K, the fraction of DAA water species decreases with increasing temperature, while the fraction of DA water species is almost unchanged (see Figure 18e). As a result, the total number of the free O−H groups, which can be counted as a sum of the DA and DAA species, decreases with increasing temperature up to around 200 K. On the other hand, above 230 K, the DAA fraction is almost constant, whereas the DA fraction increases with temperature. As a result, the number of the free O−H groups is minimal, and thus the number of the excess H-bonds is maximized at around 200 K. This is counterintuitive, as the entropic contribution increases with increasing temperature, which is expected to decrease the H-bond number of the interfacial water. This anomaly of the interfacial water can be rationalized as follows. At temperatures below 170 K, the topmost ice layer forms hexagonal ice. When the temperature is increased up to 200 K, the entropic contribution becomes increasingly important and the structure of the topmost water layer becomes disordered. With this disordered structure, the topmost ice surface generates nonhexagonal ice structures such as pentagonal, heptagonal, and octagonal ice structures, 245 allowing the interfacial water molecules to increase the H-bond numbers. When the temperature is increased above 200 K, the weak H-bonds are broken, thus producing more DA water species. These processes are schematically depicted in Figure  18f.

Continuous Change between Premelting Water and Liquid Water: 3530 cm −1 O−H Stretch Mode
A feature in the SFG spectra at the ice−air interface, which is distinct from the water−air interface, is the presence of a 3530 cm −1 peak. Figure 19 shows the experimentally measured and simulated SFG spectra of this 3530 cm −1 peak at the ice−air interface, together with the individual contributions of the DDAA, DDA, and sum of DA and DAA contributions (DAA + DA). One can see that this 3530 cm −1 peak arises from the competing contributions of the DDAA molecules and the DDA molecules with opposite signs; the DDAA water molecules show a positive contribution at around 3475−3530 cm −1 , while the DDA water molecules have a negative contribution at around 3440−3480 cm −1 . Because the amplitude for the positive contribution from the DDAA water molecules is larger than that of the negative contribution from the DDA water molecules, the competing contributions result in the 3530 cm −1 positive peak. The contributions of the DDAA water molecules originate from their asymmetric O−H stretch mode, while the contribution of the DDA water molecules results from their symmetric O−H stretch mode (see the right panel of Table 1). As such, the negative sign of the DDA water molecules contribution reflects the net orientation of the DDA water molecules with both O−H groups pointing down to the bulk.
Because the DDA water molecules are located exclusively in the topmost layer, it is interesting to explore the temperature dependence of the spectral contribution of the DDA water molecules, in addition to the contribution of the water molecules with the free O−H group (DA + DAA, see section 4.1). Simulation data shows that the amplitude of the negative ∼3450 cm −1 peak (DDA symmetric mode) decreases when the temperature increases. At the same time, a positive 3590 cm −1 peak for the DDA antisymmetric mode emerges; this positive peak is absent at 150 K because the antisymmetric DDA mode is oriented parallel to the surface for the perfect ice surface, resulting in being SFG inactive. With increasing temperature, the positive peak emerges above 200 K, demonstrating that the disordering of the DDA water molecules starts from ∼200 K. 50 This is consistent with the previous SFG study of ice. 66,67 This variation of the contribution of the DDA water species provides a logical connection between ice−air and liquid water− air SFG signals. The transition dipole moment of the DDA symmetric O−H stretch mode points down to the bulk, providing a 3450 cm −1 negative peak (see the right panel of Table 1). This negative peak is weakened with increasing temperature, which can be attributed to disordered DDA conformations; because the transition dipole moment is more tilted, the 3450 cm −1 negative peak amplitude is reduced. By  Table 1). This 3600 cm −1 DDA antisymmetric peak eventually appears as a 3600 cm −1 shoulder peak at the liquid water−air interface, as discussed in section 3.1.2. The 3530 cm −1 SFG feature can be captured not only by POLI2VS MD simulations but also by AIMD simulations. Galli and co-workers performed AIMD simulation at the ice−air interface. 126 The SFG spectrum simulated by AIMD simulation at the PBE level of theory (without any red-shift) 126 is shown in the broken blue line of Figure 20a. This spectrum seems to differ significantly from the spectrum simulated with the POLI2VS simulation (red line). However, as is discussed in section 3.1.3, the frequency can be shifted with different XC functional. To obtain the relation of the O−H stretch frequency experimentally measured (ω exp ) and that predicted at the PBE level of theory (ω PBE ), it is convenient to use the water−air SFG data, as it is well-established. Figure 20b displays a simulated SFG spectrum at the PBE level of theory (broken blue line) 58 along with an experimental spectrum (red line) 38 obtained at the liquid water−air interface. Comparing with the experimental SFG spectrum, the PBE functional tends to provide blue-shifted free O−H stretch frequency and significantly red-shifted H-bonded O−H stretch frequency. To correct the frequency mismatch, we used the relation of: This relation leads to a reasonable agreement between the simulated (solid blue line) and experimental (red line) SFG spectra (Figure 20b). By utilizing eq 19, we replot the SFG spectrum (solid blue line) at the ice−air interface in Figure 20a. The spectra shifted with eq 19 and the POLI2VS data (red line) are in good agreement, indicating that the POLI2VS and AIMD can provide qualitatively consistent results at the ice−air interface. At the same time, this may raise another question about the molecular origin of a negative ω PBE = 3200 cm −1 peak, which is out of scope in the current review.

Challenge for Static SFG Spectra
Above, we have overviewed the recent SFG measurement and modeling at the water−air and water−ice interfaces, mainly by focusing on MD simulation. At the water−air interface, recent technical developments in experimental SFG allow us to cover the frequency range from 700 cm −1 (librational motion of water) to 4000 cm −1 (free O−H stretch region). Further measurement reached the combination mode of bending mode and stretch mode of interfacial water at 5000−5300 cm −1 . 246 In contrast, vibrational modes below 700 cm −1 have not been investigated. For example, the H-bond stretch mode and H-bond bending mode have characteristic frequencies of 180 and 50 cm −1 , respectively. These are neither experimentally probed nor theoretically analyzed. Intense THz light source may pave a path to explore such a low-frequency vibrational mode.
At the ice−air interface, following the pioneering SFG works by the Shen 66,67 and Shultz 68−71,247,248 groups, new approaches have appeared such as heterodyne detection 50,119,249−251 and combined SFG measurements and simulation. 50,55,119,250 These attempts started just 1−2 years ago. Although we have learned much about the reorganization of the topmost bilayer of the ice surface from the high-frequency O−H stretch modes (ω > 3400 cm −1 ), substantial controversy remains concerning the O−H stretch mode in the 3000 cm −1 < ω < 3400 cm −1 frequency range. The heterodyne-detected measurement in ref 249 shows a large positive 3100 cm −1 peak and a small negative 3250 cm −1 at ∼130 K, whereas ref 50 shows a small positive peak at 3100 cm −1 and a large negative peak at 3150 cm −1 (see Figure 21). The simulated SFG data with ab initio-based force field models is not available yet, while the MD simulation of the CRK model of water by Morita and co-workers predicts a small positive 3200 cm −1 peak and a large negative 3250 cm −1 peak (see Figure  21). 252 Note that because the CRK model tends to predict higher frequency for the ice spectra, 225,252 the displayed spectrum was red-shifted. The simulation data seems to resemble the experimental data of Bakker and co-workers but disagrees with the experimental data of Yamaguchi and coworkers. 50 This controversy should be solved in the near future. 72 Further ice spectra have been presented using isotopically diluted water; Matsumoto and co-workers found the positive 3250 cm −1 band in addition to the negative 3300 Together with these, a simulated spectrum shifted with eq 19 is also plotted. Note that the original data in ref 126 was red-shifted by 100 cm −1 from the spectrum in broken blue line. (b) Experimentally measured SFG spectrum 37 as well as the spectrum simulated at the PBE level of theory for the isotopically diluted liquid water−air interface. 58 The spectrum with the frequency corrected with eq 18 is also displayed. Note that the experimental measurement overestimates the signal around ∼3200 cm −1 , as mentioned in section 3.1.1. Because there is no available corrected data for HOD−air interface, we adopted it from ref 37. Note that the experimental data are not corrected for the Fresnel factor. Chemical Reviews pubs.acs.org/CR Review cm −1 band at 120 K, whereas the heterodyne-detected data of Bakker and co-workers have only a negative 3300 cm −1 peak and do not have a positive 3250 cm −1 peak at 245 K. 251 For obtaining unified view of the ice surface, one needs to confirm that the emergence of the positive 3250 cm −1 band arises from the temperature change. We note that these are also nicely summarized in ref 72, with a particular focus on the experimental data.
Although modeling of the SFG spectra at the ice interface has been started for a simple ice basal model which has random proton ordering and is composed of Ice Ih structure, 50,55,119,252 much more complicated behavior of the ice interface has been reported. Theoretical studies have suggested that the ice surface has a proton ordering pattern, a so-called Fletcher pattern, rather than random proton ordering. 248,253,254 Furthermore, the ice surface may be composed of mixed hexagonal ice (Ih) and cubic ice (Ic) structures through the stacking disorder. 255 These questions can be addressed only through combined SFG measurement and simulations; simulations are potentially able to calculate the SFG spectra of ice for various situations and testify to the agreement of the spectral shape between the simulation and experiments. Such a pioneering attempt has been conducted by Buch and co-workers. 248 To make these approaches successful, again the accurate force field models are essential. Furthermore, the computed spectra assume that the contribution arises solely from the interface. This assumption has been checked carefully.  119,209,230,256 and two-dimensional SFG measurements 210,257−261 have been developed, and an increasing number of studies has been undertaken to understand the vibrational dynamics of water 36,215,229,230,262−269 and protein structure characterization performed 270−272 at the interfaces. However, the simulation of the 2D SFG has not been reported so much due to the complex modeling of the two-dimensional spectra together with the surface-specific modeling. 273−277 For computing the 2D spectra, one can have two routes: the route based on frequency mapping technique and the route of dipole/polarizability time-correlation function, as accounted for in section 2.1.2. The mapping technique does not require the explicit intramolecular vibrations, but it requires the additional treatment for including the coupling effects in the spectra. The time-correlation function approach is more straightforward than the mapping technique because coupling effects are automatically included in the time-correlation function. On the other hand, the computational cost for the time-correlation function approach is huge not only for the 2D-SFG spectra calculation but also for general 2D spectra calculation. Here, we outline the time-correlation approach.
The quantum expression of the fourth-order response function for 2D-SFG is given by  20) where [A, B] is the quantum commutator and t 1 , t 2 , and t 3 are the time intervals between the first IR and the second IR, between the second IR and the third IR, and between the third IR and visible pulses, respectively. The 2D-SFG susceptibility is then given by As is seen, the time intervals of t 1 and t 3 were Fourier transformed, converting the time domain into the frequency domain (ω 1 and ω 3 , respectively). t 2 is the waiting time. By monitoring the correlation of ω 1 and ω 3 as a function of delay time t 2 , one can monitor the vibrational dynamics. Equation 20 is the quantum expression for the fourth-order response function, and thus one cannot connect eq 20 with the classical MD simulation. This equation can be recast in the classical limit as  is equal to ⟨{M b (t 1 ), Ṁa(0)} P.B. ⟩ requires many sampling of the trajectories. Rather than evaluating eq 22 directly, one can use one more efficient technique based on the hybrid equilibrium nonequilibrium technique, which avoids the second Poisson bracket calculation, by using the external electric field in the nonequilibrium MD simulation. 279,280 In fact, this is the major route for computing the 2D-IR spectra with explicit vibrations in the MD simulation (without frequency mapping technique). 281−283 5.2.2. Excitation with Specific Frequency. With the above formalism, the molecules are excited through a small perturbation of position and momentum in the Poisson bracket. In this scheme, all the vibrational modes are excited by using a delta-function pulse. However, often, we only need to know the vibrational dynamics for excitation at a specific frequency. For this purpose, it is straightforward to excite a specific vibrational mode and monitor the vibrational dynamics for an excited (2),R ) spectra of the ice−air interface, from Bakker's group at 150 K, 50 Yamaguchi's group at 130 K, 249 and Morita's group at 130 K. 252 Note that we shifted Morita's spectrum to the low frequency side by 120 cm −1 . Note that the Fresnel factor is not corrected in experimental data. Chemical Reviews pubs.acs.org/CR Review vibrational mode. This can be done by applying the alternating current electric field with a targeted frequency. 106 The interaction of the electric field and molecules excites the vibrational modes in the nonequilibrium MD simulation. By computing the time-correlation functions as a function of the time passing after the excitation (delay time), one can monitor the variation of the spectra. Describing the interactions between the time-dependent electric field and molecules in the simulations is sometimes not straightforward. In particular, when we use the AIMD simulations, the interaction of the electric field is treated in the Berry phase approach. 284,285 Thus, one may avoid the explicit modeling of interactions between the electric field and molecules. The vibrational mode excitation without an explicit electric field can be achieved, for example, by using the velocity swapping technique 88 or the specialized thermostat. 286 In this velocity swapping technique, the atom velocities are swapped with a target frequency so that the specific vibrational mode acquires excess kinetic energy. In these approaches, it is essential that the perturbation should be small enough to ensure that the linear response approximation works and the input energy for modulating the system is much smaller than the sum of the total kinetic energy of the whole system. Experimentally, the temperature increase due to the vibrational excitation of the system is estimated to be 3 K. 200 All approaches using the explicit electric field, velocity swapping, and thermostat methods should increase the kinetic energy very moderately, following the spirit of the Poisson bracket where the small perturbation is required. This is linked with the fact that the energy levels in classical mechanics are continuous, unlike in quantum mechanics.
On the other hand, when we consider the quantum excitation of the O−H stretch mode, the excitation energy will be ∼16 kT, which is orthogonal to the small perturbation explained above. When one wants to have the high excitation energy in the classical MD simulation, the above-mentioned approaches using velocity swapping, thermostatting, or explicit electric field does not work; an approach is required to add a large excess of energy to the one vibrational mode. This excess energy can be added by simply pulling the O−H bond or adding the additional velocity to the O−H group. Subsequently, by monitoring the kinetic energy of the vibrational mode, one can monitor the release and flow of the excess vibrational energy. This technique has been initiated by Hynes and co-workers, 287,288 and has been extended to the interfacial water. 289 Here, we note that the excitation (excess) energy mimics the experiment, while it breaks the quantum-classical correspondence of iℏ[A, B] → {A, B} P.B. /kT. Comparing these two techniques is a future challenge.

5.
3. X-ray SFG Spectroscopy 5.3.1. X-ray Techniques on Water and Ice. So far, the review has focused on sum-frequency processes involving IR and visible pulses to probe vibrations. We now discuss SFG processes involving X-ray pulses. While IR and optical frequencies efficiently probe vibrational and electronic transitions, respectively, they lack the element selectivity that X-rays can offer. Specific elements can be excited by tuning the X-ray frequency to a desired core electronic transition. Different local environments can also lift the degeneracy of the core transition (chemical shift) to address atoms located in different local environments. Adding the spatial resolution obtained from offresonant scattering, X-rays are naturally an excellent probe of molecular, crystalline structure and can advantageously be used to investigate the water interfaces.
Nonlinear X-ray spectroscopies have experienced dramatic progress due to the advent of X-ray free electron laser (XFEL), third-generation synchrotron light and high harmonic generation (HHG) sources. 290 In particular, XFEL pulses typically have many improved parameters such as fluence, spatial and temporal coherence, and time resolution 291 that make nonlinear coherent X-ray spectroscopies a reality. For example, the X-ray pump−probe (XPP) end station of the Linear Coherent Light Source (LCLS, USA) delivers 80 fs pulses between 5 and 25 keV, 2 mJ pulse energy, and 3 × 3 μm size. XFELs have a brilliance about 25 orders of magnitude higher than that of X-ray tube lamps and ∼8 orders higher than at synchrotrons. Moreover, the spatial and temporal coherences of XFELs allow envisioning the implementation of multidimensional coherent spectroscopies, 292 involving sequences of several pulses. XFEL can deliver femtosecond pulses, down to the attosecond regime, 293 providing an excellent resolution for time-resolved spectroscopies. There are currently seven large-scale XFELs facilities: , and European X-ray FEL (EXFEL, Germany). These sources extend coherent ultrafast spectroscopy to the X-ray regime. X-ray photons have enough energy to photoexcite core−electrons either to bound states (i.e., core excited states) or to the continuum (photoionization). The absorption cross-section increases dramatically near the bound core electronic states and leads to absorption edges. In Figure 22a, we display the X-ray absorption spectrum at the oxygen K-edge in ice and liquid water measured by Nilsson et al. 294 The K-edge corresponds to excitation of the 1s orbital, while L-edge refers to the 2s, 2p orbitals. Being a light element, the oxygen K-edge occurs in the soft X-ray regime, the main edge features for water and ice are around 538 and 541 eV, respectively (see Figure 22a). Core orbitals of different atoms have specific binding energies, and thus the incoming photon frequencies required to generate core-excited states are very different. X-ray absorption and in general resonant X-ray processes thus have high element selectivity. 294−297 For example, the sudden increase in absorption observed at ∼540 eV is due to the excitation of 1s orbitals in oxygen while the carbon 1s is at 284 eV. This permits addressing of local sites in a material and to follow local dynamics in pump−probe schemes. Below and near the absorption edge, one can observe pre-edge transitions that carry useful information regarding the electronic structure. This is known as X-ray absorption near-edge spectroscopy (XANES, also called NEXAFS for near edge X-ray absorption fine structure). At higher energies, the cross-section is dominated by photoionization and backscattering of the ionized electrons by the nearest neighbors. This region, named extended X-ray absorption fine structure (EXAFS), displays oscillations due to the interference of different backscattered electrons and carries local structure information. 298 As shown in Figure 22b, precise information on the O−O distance can be obtained from such measurements. 299 The experimental and theoretical study of Xray absorption at the oxygen K-edge has led to a debate on the local structure of liquid water assuming either H-bonded rings and chains or the standard local tetrahedral structure. 294,300 Moreover, the NEXAFS signal can be made surface-sensitive by measuring it in total ion yield mode. 301 303 Molecules with an uncoordinated O−H group display strong pre-edge X-ray absorption features, and one can deduce from these measurements the structure of the ice surface. 302 Another unique property of X-ray light is its extremely short wavelength. In the hard X-ray (∼10 keV), the wavelength is roughly 1 Å, and the light can be diffracted from charge densities at this scale. In ordered samples, e.g., crystals, this leads to the well-known Bragg diffraction that has been used for a long time in crystallography, while disordered samples show diffraction rings. 267,268 Figure 23 displays diffraction patterns from an ice crystal 304 and liquid water. 305 X-ray diffraction can be made surface-sensitive by using grazing incidence geometries in techniques such as grazing incidence X-ray diffraction (GIXD) or diffuse X-ray scattering (DXS). 306,307 The low momentum transfer measured with these techniques imply that the structural information is obtained for large scales (>10 Å).
The first step for nonlinear X-ray techniques was implementing pump−probe experiments, either with an absorption or a diffraction event as an observable after an IR or visible excitation. 308 Just like in the visible regime, being a third-order technique, 309 the pump−probe process is not surface-sensitive. Frequency domain nonlinear experiments have also been carried out including two-photon absorption, three-and four-wave mixing. 310,311 Pump−probe XANES is highly sensitive to changes in oxidation states and electronic configuration, 312 while pump− probe EXAFS or XRD measure accurately triggered nuclear dynamics. 313,314 These techniques, and the nonlinear spectroscopies involving X-rays currently under development, are usually Chemical Reviews pubs.acs.org/CR Review hybrid techniques that merge the sensitivity of different photon probes. These techniques mixing a low energy probe in the THz, IR, visible, or UV regimes launch a low-energy dynamic, subsequently probed by XANES, EXAFS, and XRD. At the oxygen K-edge, X-rays are highly sensitive to the local electronic structure in the vicinity of the oxygen atoms but probe bulk properties due to the large penetration depth of X-rays (530 nm in water and 570 nm in ice at 545 eV). Surface sensitivity is routinely achieved by grazing incidence experiments, 315 allowing the X-ray to penetrate only tens of nanometers in the material. Below we provide an expression for these various observables for X-ray SFG (XSFG) techniques. The X-ray SFG signals have the potential to offer a new window into surfacesensitive techniques.

Minimal Coupling
Hamiltonian for X-ray Interactions. The signals expressions of the SFG signals in the former sections rely on the nonlinear optics formalism in the dipole approximation that has been developed over the past decades. 316,317 This approximation is not always suitable in the X-ray regime due to the short wavelength of X-ray that allows for structure determination.
For nonlinear X-ray spectroscopies, and because we aim at discussing both off-resonant (diffraction-like) and resonant (absorption-like) processes observables, we describe lightmatter interaction in the minimal coupling Hamiltonian: 318 where j(r) and σ(r) are the current and charge density operators and A is the electromagnetic vector potential. Equation 23 leads to two observables where the final process involves the j · A term or with the σA 2 one. Moreover, signals can be either homodyne or heterodyne detected as discussed earlier. The minimal coupling Hamiltonian includes naturally diffraction detected signals through the charge density operator, while the current density operator includes implicitly all multipoles involved in resonant techniques. 319 5.3.3. Off-Resonant X-ray SFG. We first focus on offresonant homodyne-detected SFG. In this case, the first interaction brings the molecule into a vibrational or electronic excited state that is then probed by an off-resonant X-ray pulse; see diagrams (d) and (e) in Figure 24 for the heterodyne and homodyne cases, respectively. Thus, the off-resonant SFG is linear in the j · A and in the σA 2 interaction. Off-resonant SFG signals are dominated by σ(q) where q is the momentum transfer between the incoming X-ray beam and the diffracted photon. The diffraction-detected SFG then reads: This represents a diffraction event with the field A X after a linear interaction with an infrared, visible, or even X-ray field noted A. F 2 (q) is the structure factor, q is the momentum transfer, ϵ s is the detected field polarization, and ε ω = ℏ Ω A /2 s s 0 . The material quantity shows up as two-point ⟨σj⟩ matter correlation function (dropping the space and time arguments for brevity). It is obtained by a second-order perturbative expansion in the incoming field of the observable and is indeed a sum-frequency process. This is clearly seen by expressing the homodyne signal in the frequency domain:  (25) which contains the sum over the pumped manifold of states. The surface sensitivity of this SFG technique is complicated by the fact that the charge density is a scalar field that does not belong to an irreducible representation of the inversion group. This means that depending on the transition matrix element σ ge probed, the scalar field can have odd or even parity, or a combination of the two. The even parity term leads to an odd response tensor ⟨σj⟩ and is thus surface-sensitive following the standard argument on a χ (2) tensor. The technique can be made surface-sensitive under various conditions. First, even σ gc can be Chemical Reviews pubs.acs.org/CR Review preferably considered. Because σ gc (r) = ψ g *(r)ψ c (r), where ψ g and ψ c are the differing orbitals between states g and c, an even σ gc (r) is obtained when the states c and g are either both even or both odd with respect to the local frame of coordinate.
Otherwise, X-ray experiments are routinely done in grazing incidence, allowing to probe only surface properties of materials despite the large penetration depth of X-rays. Being diffraction detected, this type of XSFG signals should be able to follow in time-domain the time evolution of the charge density of the surface molecules at an interface following an IR excitation triggering a vibrational motion for example. The time-domain XSFG in the impulsive limit reads: Because the matter correlation function involves three vector operators, the resonant signal is surface sensitive. In this technique, the interactions with the various pulses usually involve hybrid infrared/X-ray scheme, offering element selectivity and broadband detection. A resonant all X-ray technique is hard to envision because it would involve doublecore holes that are extremely short-lived, screened by various other incoherent recombination processes and carry a too localized information to be surface sensitive.
5.3.5. Homodyne vs Heterodyne Detection for Resonant X-ray SFG. Summarizing the possible signals in the minimal coupling framework, we are left with four possible detection modes: Homodyne-detected signals are the most commonly implemented in the X-ray regime up to now. This is due to the difficulty of generating interferences with an X-ray local oscillator because of the short coherence length of X-ray pulses.
However, heterodyne detection schemes should be feasible in the near future. In Figure 24, we present a general level diagram for the X-ray hybrid SFG. A low energy pump excites the system in a |e⟩⟨g| coherence where e is a vibrational or an electronic excited state. Then, an X-ray beam probes the sample. At resonance, a core-excited state coherence |c⟩⟨g| is created and can be probed in a heterodyne or homodyne detection scheme, as shown in parts b and c of Figure 24, respectively. Light-matter interaction is given by electric transition dipole moments μ in the dipolar coupling or by transition current density j in the minimal coupling. In the off-resonant case, no core excitation is created, and the transition is described by the transition charge density matrix element between the e and g states (transition polarizability), as discussed above. Both heterodyne and homodyne detection are possible diagrams in parts d and e of Figure 24, respectively. A third approach to XSFG prepares a low-energy wavepacket through two initial interactions with an X-ray pump, similarly to a Raman process, which is subsequently probed by a visible or infrared probe. On resonance, this is again described by a ⟨jjj⟩ multipoint correlation function, while off-resonance, the signal is given by a sum of ⟨σj⟩ and ⟨jjj⟩ response functions. The resonant excitation can make use of the X-ray element specificity, while the offresonant one can use the large broadband offered by X-ray pulses.
Chemical Reviews pubs.acs.org/CR Review 5.3.6. Outlook. We have surveyed the characteristics of Xray spectroscopies (element sensitivity, atomic resolution) and the recent progress in implementing ultrafast and nonlinear Xray spectroscopies. We then provided a theoretical framework that expresses these signals using a sum-overstates expression involving transition current and charge densities. A truncated multipolar expansion of the current densities can recover the infrared−visible expressions discussed earlier in the electric dipole approximation.
Time-domain XSFG signals have not been demonstrated yet, and their pursuit has just started recently with the development of XFEL and new generation synchrotron sources. Secondharmonic generation has been demonstrated at the carbon Kedge. 310 With the same symmetry argument used in lower energy regime, resonant XSFG is expected to be surface specific. At the oxygen K-edge, it is clear that the ability to probe at the ice or water surface the oxidation state of the oxygen atoms and the O−O distance is appealing for the observation of intermediate states triggered by a low energy pump beforehand.
Off-resonant XSFG may still lead to a nonvanishing signal from the bulk due to the fact that transition charge densities do not belong necessarily to an irreducible representation of the inversion group.
We have focused on hybrid XSFG techniques. Recent experimental developments 320 show that doubly core-excited state may be accessible and observable. However, the transition from there to all X-ray XSFG seems challenging in the near future and the obtained information is expected to be highly localized on the oxygen atom and may not be sensitive to the surface as a whole.

CONCLUSIONS
Through the comparison of theoretical SFG spectra with experimental data, we obtained critical feedbacks both for the simulation and the experimental communities. The message to the simulation community is that calculating SFG spectra offers a unique avenue for accessing the accuracy of the current simulation's methodologies, which are often inaccessible with the standard benchmarks. For example, among the force field models/DFT functionals which can closely reproduce the experimental bulk IR spectra of water, they show a large variation in the SFG spectra at the water−air interface because the SFG spectroscopy probes heterogeneous interactions at the interface. Newly developed models of water, based on machine learning or traditional parametrization, can be well-testified by calculating the SFG spectra. To the experimental community, we highlight that the state-of-the-art simulation methodologies can be used not only for connecting the microscopic structure of interfacial water with SFG spectra but also for complementing the experimental SFG spectra. Currently available force field models and DFT simulations for water are accurate, thanks to the advances in modeling methodologies and increased computational resources. By combining the algorithms outlined in this review with these advanced models, the SFG spectra can be accurately predicted.
Finally, we would like to emphasize that combined SFG experiment and simulation have a significant potential to access the chemistry/physics/biology at the interface. Some of the examples are presented in this review (water orientation and surface nanoroughness of water-air interface), while other examples can be found elsewhere (for example, ice friction, 74,328 ). Tight collaboration between experiment and simulation will expand the potential application field of SFG spectroscopy.