# Microscopic view on the polarization-resolved S-SHG intensity of the vapor/liquid interface of pure water

Guillaume Le Breton, C. Loison, Kevin Vynck, E. Benichou, O. Bonhomme

► To cite this version:

Guillaume Le Breton, C. Loison, Kevin Vynck, E. Benichou, O. Bonhomme. Microscopic view on the polarization-resolved S-SHG intensity of the vapor/liquid interface of pure water. The Journal of Chemical Physics, 2024, 161 (15), pp.154712. 10.1063/5.0231240. hal-04756407

HAL Id: hal-04756407

https://hal.science/hal-04756407v1

Submitted on 4 Nov 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Microscopic view on the polarization-resolved S-SHG intensity of the vapor/liquid interface of pure water

G. Le Breton,¹, a) C. Loison,¹ K. Vynck,¹ E. Benichou,¹ and O. Bonhomme¹

Institut Lumière Matière, UMR5306 – UCBL – CNRS, 10 rue Ada Byron, 69622, Villeurbanne CEDEX, France

(*Electronic mail: oriane.bonhomme@univ-lyon1.fr)

Second harmonic generation (SHG) is a nonlinear optical phenomenon where two photons at the frequency $\omega$ combine to form a single photon at the second-harmonic frequency $2\omega$. Since that second-order process is very weak in bulk isotropic media, optical SHG responses of interfaces provide a powerful and versatile technique to probe the molecular structure and dynamics of liquid interfaces. Both local dipole contributions and non-local quadrupole contributions can be interesting to investigate different properties of the interface, such as the molecular orientation or the charge density. But a major difficulty is to comprehend the link between the S-SHG intensity and molecular details. This article reports a numerical approach to model the polarization-resolved SHG intensities of a model vapor/liquid interface of pure water. The influence of the interfacial local environment on the hyperpolarizability is taken into account using Quantum Mechanical / Molecular Mechanics (QM/MM) calculations. The numerical predictions are in very good agreement with experiments. We detail the hypotheses made during the modeling steps, and discuss the impact of various factors on the modeled SHG intensities, including the description of the exciting field in the interfacial layer, the effect of neighboring molecules on the second-harmonic polarization, and the presence of an additional static electric field at the interface.

## I. Introduction

Due to intrinsic symmetry constraints, the second-order non-linear optics effects are significantly stronger in non-centrosymmetric media. Exploiting this specificity, Y. Ron Shen was a pioneer in the development and application of second-order nonlinear optics nonlinear techniques to study interfaces and Surface-Second Harmonic Generation (S-SHG) has been developed in the last decades to monitor various phenomena at many kinds of interfaces¹⁻³. The principle of a S-SHG experiment on a liquid/vapor interface is schematized on Fig. 1(a). The interfaces of aqueous solutions have notably been extensively studied using S-SHG because of their importance in many physical, chemical or biological domains. For example, S-SHG spectroscopy can be used to monitor the orientation of adsorbed molecules, the adsorption of ions and the interaction of biological species²,⁴⁻⁶.

Liquid interfaces are complex, with specific molecular structures and a discontinuity in the refractive index⁷. The liquid composition, the hydrogen bond network and the molecular orientations are modified at the interface⁸,⁹. An intrinsic static electric field arises at the interface due to differences in refractive indexes, to anisotropic dipole orientations, or ionic layers near the interface¹⁰⁻¹². The S-SHG intensity depends on all these features: the refractive index variation, the atomic structure, the local static electric field and its gradient¹³,¹⁴. Experimental innovations alongside modelings may help to disentangle the different contributions to the intensity, and to interpret them at the microscopic scale¹⁰,¹⁵⁻¹⁷.

But even for pure water, a quantitative prediction of the po-

<figure><img src="image_2.png" /><figcaption>FIG. 1. (a) Schematic view of S-SHG experiments on a liquid/vapor interface. The altitude z corresponds to the distance to the interface, with negative values in the liquid. (b) Synopsis of the theoretical procedure used to describe numerically the S-SHG intensity collected. See details in Sec. II.</figcaption></figure>

larization resolved S-SHG intensities of the vapor/pure water interface remains complex¹⁸,¹⁹. Using a multipolar expansion, this response contains different contributions. Concerning the dipolar term, the molecular SHG responses quantified through the molecular hyperpolarizability $\beta^{dd}(-2\omega; \omega, \omega)$, strongly depend on the solvation, which in turn depends on the position of the molecules relative to the interface (noted

<footnote>a) Now at Institut de Chimie Moléculaire de l'Université de Bourgogne, Univ. Dijon, France</footnote>
as the altitude $z$, where the plane $z = 0$ defines the Gibbs dividing surface$^{19-23}$. In addition, the dipolar hyperpolarizability of a water molecule at the interface is not independent on its orientation$^{22}$. After some debate, it is now clear that the dipolar contributions alone are not enough to account for the experimental polarization resolved S-SHG measurements$^{13,18,19,24}$.

Shiratori et al.$^{18}$ used classical molecular dynamics (MD) simulations to obtain the structure of the interface, and exploited it to calculate altitude-dependent rotational matrices that describe the molecular orientation. They calculated a value for the dipolar and quadrupolar hyperpolarizability tensors using Density Functional Theory (DFT) calculations on a single water molecule, adding a static electric field to mimic solvation effects. They concluded that the quadrupolar contributions from the bulk are dominant.

Also using classical MD simulations followed by DFT calculations on large water assemblies representing a liquid/vapor interface, Foucaud et al. proposed a model where solvation and local-field effects are intrinsically included in the calculated response$^{24}$. Again, they demonstrated that the bulk quadrupolar contributions nearby the interface dominates the interface dipolar ones. More recently, they developed an alternative molecular-based model including both dipolar and quadrupolar terms$^{19}$, using the values of the polarizability and hyperpolarizability obtained by density functional theory (DFT) calculations on a water molecule surrounded by an ice-like solvation layer. They successfully reproduced the features of the experimental intensities, but some discrepancy remains. The authors suggested that a more realistic description of the solvation effect may increase the accuracy of their S-SHG prediction.

Along this line, this work presents a molecular-based model of the polarization-resolved S-SHG intensities of the air/water interface, taking into account (1) more realistic short-range solvation effects, (2) the correlations between the molecular orientation and the dipolar SHG response, and (3) the presence of a static electric field generated by the interface.

The article is structured as follows: In Sect. II, we develop the theoretical basis and hypotheses of the different steps, represented in Fig 1(b), leading to the S-SHG intensities. First, the electric field within the sample induced by the exciting laser is calculated analytically using Fresnel coefficients and local field effects (step A). Then the SHG polarization of the molecules is deduced from QM/MM calculations (step B), and summed up into a mesoscopic polarization (step C). Finally, the light emitted by this interfacial SHG polarization is estimated in the framework of the three-layer model (step D). Sect. III presents the experimental setup and the numerical methods. Section IV presents the experimental and numerical results, focusing first on the dipolar and quadrupolar hyperpolarizabilities obtained by QM/MM. Beyond the bare QM/MM results that include solvation effects up to a few nm only, the static electric field generated by the water dipoles at the interface is also taken into account in these contributions as its crucial effect on dipolar hyperpolarizability was already demonstrated$^{20,25}$. The resulting polarization resolved S-SHG intensity curves are then compared to the experimental ones.

Finally, we also discuss how some approximations involved in steps A to D may affect the results.

## II. Theoretical elements

In this first section, we present the different steps and corresponding formalisms used to describe the SHG intensity created by a vapor/liquid interface based on the responses of individual molecules, see Figure 1(a). The system studied reproduces the standard geometry of a polarization-resolved S-SHG experimental setup on a vapor-liquid interface in the $xy$ plane at an altitude $z=0$. The linearly-polarized fundamental incident electric field $E$ arrives at the interface with an angle of incidence $\theta$ relative to the surface normal. The harmonic field $E^{2\theta}$ is detected in the direction of the specular reflection. In a polarized-resolved system, the polarizations of the incident and harmonic fields are controlled: the S-polarized electric field is perpendicular to the incident plane whereas the P-polarized field is in the plane of incidence. As mentioned previously, four steps, labelled A to D, as illustrated in Fig 1(b) are mandatory. In the first step (A), the incident electric field in the vapor phase is transformed into a local electric field in the liquid that varies as a function of the altitude $z$. In the second step (B), the individual response of the molecule at the second harmonic frequency due to this electric perturbation is described, leading to an induced dipole $\mu$ and quadrupole $\nu$. In the third step (C), the mesoscopic polarization is introduced as the collective response of the molecules. This allows us to define the susceptibility tensor probed by S-SHG experiments. In the final step (D), the electric field generated at the harmonic frequency is computed from the mesoscopic polarization, and the resulting SHG intensity can be compared to experimental measurements.

### A. Perturbation: fundamental field across the liquid-gas interface

Experimentally, the fundamental field $E$ is applied with a grazing incident angle $\theta$ of about 70 degrees. This angle is fixed, while the incoming light polarization angle, noted $\alpha$, varies from 0 to 180 degrees. $\alpha = 0^{\circ}$ (resp. $\alpha = 90^{\circ}$) corresponds to a polarization in (resp. out) the plane of incidence, namely p-polarization (resp. s-polarization).

Due to the optical gap between the gas and the liquid phase, the field amplitude along the different $x$, $y$ and $z$ directions evolves through its propagation. The relation between the gas and bulk phase amplitudes are provided by the classical Fresnel coefficients, leading to 2 different field amplitudes between each part of the liquid-gas interface. The field felt by a molecule within the interface, often called “local field” and noted $e$, evolves from one limiting value to the other. In practice, the fundamental electric field that reaches the sample is expressed in the laboratory frame $\{x, y, z\}$ as a vector $E$ and the local field $e$ in the laboratory frame via a transfer matrix
$\mathcal{T}(z)$ as

$$e(z) = \mathcal{T}(z)E. \qquad (1)$$

where the $\mathcal{T}$ matrices depend on $\theta$ and $E$ on $\alpha$ ($E = E\{\cos\alpha, \sin\alpha, \cos\alpha\}$ with $E$ the field amplitude). Here, the field $e$ at the fundamental frequency depends solely on the altitude $z$, i.e. the surface respects a $C_{\infty\nu}$ symmetry.

Morita et al. developed a formalism for the matrix $\mathcal{T}(z)^{18,26}$ including two contributions: the Fresnel factors and the effect of neighboring molecules on the field exciting a molecule. Accordingly, the two values of $\mathcal{T}$ just above the surface and in the bulk liquid can be written:

$$\lim_{z \to 0^+} \mathcal{T} = \mathcal{L}_s, \qquad (2)$$

$$\lim_{z \to -\infty} \mathcal{T} = \frac{2+n_2^2}{3} \mathcal{L}_l,$$

The $(2+n_2^2)/3$ prefactor represents the local field effect in the liquid, where $n_2$ is the refractive index of the liquid phase ($n_2 = 1.33$ at 800nm). The continuity of the electric displacement is accounted for in the surface $\mathcal{L}_s$, and the liquid phase $\mathcal{L}_l$ matrices. The matrices $\mathcal{L}$ are:

$$\mathcal{L}_s = \begin{pmatrix} (1-r^p)\cos\theta & 0 & 0 \\ 0 & (1+r^s) & 0 \\ 0 & 0 & (1+r^p)\sin\theta \end{pmatrix}, \\ \mathcal{L}_l = \begin{pmatrix} t^p \cos\theta_2 & 0 & 0 \\ 0 & t^s & 0 \\ 0 & 0 & t^p \sin\theta_2 \end{pmatrix}. \qquad (3)$$

where the angle $\theta_2$ is defined between the direction of the refracted light in the liquid phase, and the normal to the surface ($z$ direction). The real numbers $r^p, r^s, t^p, t^s$ are the Fresnel reflection and transmission coefficients, describing the evolution of an electric field polarized in $p$ and $s$ directions respectively.

Between these two known values, we use a shift proportional to the molecular density $\rho(z)$, as suggested by Morita and co-workers$^{18}$:

$$\mathcal{T}(z) = \frac{\rho(z) - \rho_v}{\rho_l - \rho_v} \left( \frac{2+n_2^2}{3} \mathcal{L}_l - \mathcal{L}_s \right) + \mathcal{L}_s, \qquad (4)$$

where $\rho_l$ is the density in the liquid phase, and numerically, the density in the vapor $\rho_v$ was neglected. Fig. 2 reports the evolution of $\mathcal{T}(z)$ and its gradient along $z$. Here, the electric field gradients extend from about -4 to +4 Å. To summarize, the electric field experienced by the different molecules is dictated by the diagonal components of $\mathcal{T}$, which depends on the incidence angle $\theta$ and the electric field $E$ which depends on the polarization angle $\alpha$. We shall discuss the influence of $\theta$ uncertainties on the S-SHG polarization curves in Sect. IV.

## B. Molecular response: hyperpolarizability tensors

The second harmonic generated by the interface emerges from the electron's displacement upon the fundamental field

<figure><img src="image_3.png" /><figcaption>FIG. 2. (Top) Evolution of the transformation matrix $\mathcal{T}$ between the input electric field and the local electric field as a function of the altitude. (Bottom) Evolution of the first derivative along $z$ of $\mathcal{T}$ as a function of the altitude. The incident angle $\theta$ was fixed to 70°. Color curves correspond to the different components and symbols are guides to identify curves.</figcaption></figure>

perturbation$^{24}$. In a molecular description, the time-dependent Hamiltonian is noted $\hat{H}(t) = \hat{H}_0 + \hat{V}(t)$, where $\hat{H}_0$ is the Hamiltonian of the molecule without the light perturbation – but taking into account solvation effects – and $\hat{V}(t)$ is the potential describing the interaction with the electromagnetic field $e(z,t)$. For optical frequencies and small molecules, $\hat{V}(t)$ is often described within the dipole approximation: $\hat{V}(t) = -\hat{\mu} \cdot e(z,t)$, with $\hat{\mu}$ the dipolar operator. This expression is valid if the field $e$ is spatially homogeneous in the molecule vicinity.

At the interface, the local field evolves along with the density, and some strong electric field gradients appear (see Fig. 2). Hence, a quadrupolar term is added to the interaction potential:

$$\hat{V}(t) = -\hat{\mu} \cdot e(z,t) - \frac{1}{2} \hat{Q} : \nabla e(z,t), \qquad (5)$$

with $\hat{Q}$ the quadrupole operator in its canonical convention (i.e. not the traceless one). Within the response formalism, the induced dipole and quadrupole moments of a molecule $n$, noted $\mu_n$ and $\nu_n$, are obtained as tensorial products between the perturbation and response tensors:

$$\mu_n = \frac{1}{2} \beta_n^{dd} : ee + \frac{1}{2} \beta_n^{dq} : e\nabla e, \qquad (6)$$

$$\nu_n = \frac{1}{2} \beta_n^{qd} : ee. \qquad (7)$$

$\beta_n^{dd}$ is the molecular first hyperpolarizability describing the
induced dipole moment at the second harmonic upon a fundamental local field. $\beta_n^{dq}$ is the molecular hyperpolarizability expressing the induced dipole moment upon an electric field and electric field gradient and $\beta_n^{qd}$ the one expressing the induced quadrupole moment upon an electric field. As the notation and convention may differ in the literature, some properties for these 3x3x3x3 tensors are recalled in Supp. Mat. Sect. A. Note that $\beta^{dq}$ and $\beta^{qd}$ are of the same order of magnitude. We note $\beta$ the hyperpolarizabilities expressed in the molecular frame, and $\mathcal{B}$ the ones expressed in the laboratory frame.

These hyperpolarizability tensors are properties of the molecule that depend on the Hamiltonian $\hat{H}_0$ – and thus on the molecular environment. In this work, the hyperpolarizabilities of all molecules are computed taking into account their local electrostatic environment with the response framework of Density Functional Theory calculations (see Methods part) as a function of the altitude, and their values are described in part IV A. Note, however, that the induced moments as predicted from Equations (6) and (7) and used hereafter to obtain a mesoscopic polarization neglect certain aspects of intermolecular coupling. Specifically, the radiation from a molecule at the second harmonic frequency can polarize the neighboring molecules by radiative coupling; in other words, the dipole and quadrupole moments of a molecule $n$ should depend on the dipole and quadrupole moments of the other molecules $m \neq n$. This effect is less described in the literature notably because of the difficulty to properly handle the emitters’ heterogeneity (for instance, the induced dipole moment fluctuates a lot from one molecule to another at the interface$^{22}$). This aspect constitutes a possible future improvement of the work, as will be discussed in Section IV D.

### C. Collective response: susceptibility tensors

The next step is thus to define the mesoscopic SHG polarization arising from a collection of emitting molecules that constitute the polarization sheet. We define the surface density of polarization at the harmonic frequency as $P_s^{SHG} \equiv \int P^{SHG}(z)dz$ with $P^{SHG}(z)$ the volume polarization density at the altitude $z$.

#### a. Dipolar level:

Within the dipolar approximation, we have

$$P^{SHG}(z) \approx P^{dip}(z) = \rho(z)\mu(z). \qquad (8)$$

Numerically, $\mu(z)$ corresponds to the average of the molecular QM/MM observable $\mu_n$, over the molecules that are located between $z$ and $z + dz$. This is equivalent to approximating the surface polarization density as the sum of the induced dipole moment of each individual molecule per area, that is $P_s^{SHG} \approx \frac{1}{A} \sum \mu_n$, where $A$ is the area of the interface formed by the $N$ molecules.

Since the local electric field depends solely on the altitude $z$ in Eq. 1, the average molecular dipole in the laboratory frame is written as

$$\mu(z) = \frac{1}{2} \mathcal{B}^{dd}(z) : \mathcal{T}(z) \mathcal{T}(z) \mathcal{E} \mathcal{E}, \qquad (9)$$

where $\mathcal{B}^{dd}(z)$ is the molecular dipolar hyperpolarizability averaged over the ensemble of molecules located at the altitude $z$. Because it is expressed in the laboratory frame, it depends on both the individual molecular hyperpolarizability $\beta^{dd}$ and the molecular orientation.

Finally, the usual formalism in surface-SHG defines the volume susceptibility tensor $\chi$ as the dependence of the mesoscopic polarization on the fundamental light in the polarization sheet, such as $P_s^{SHG} \propto \chi : ee$. We generalize this definition to define the dipolar susceptibility tensor as a function of the altitude. Inserting the dd contribution of Eq. 6 into Eq. 8 leads to:

$$\chi^d(z) = \frac{1}{2} \rho(z) \mathcal{B}^{dd}(z) \qquad (10)$$

Integrating this susceptibility tensor along $z$ will allow computing the surface susceptibility tensor classically used to describe S-SHG experiments$^{24,27}$.

We also define an effective susceptibility tensor that links the mesoscopic polarization to the incident fundamental light such as $P^{dip}(z) \propto \chi^{\text{eff},d} : EE$. Using Eq. 10 leads to

$$\chi^{\text{eff},d}(z) = \chi^d(z) : \mathcal{T}(z) \mathcal{T}(z). \qquad (11)$$

Physically, this effective susceptibility contains both the susceptibility tensor and the local field correction of the fundamental electric field. This quantity will be used to compare the amplitude of the dipolar and quadrupolar contributions.

#### b. Quadrupolar level:

Including the quadrupolar terms impacts the molecular induced dipole moment, see Eq.(6), but also the definition of the mesoscopic polarization $P_s^{SHG}$. Eq. 8 is modified into

$$P^{SHG}(z) \approx \rho(z)\mu(z) - \frac{1}{2}\nabla \cdot \rho(z)\nu(z), \qquad (12)$$

with $\nu(z)$ the molecular quadrupole, averaged over the ensemble of molecules located at the altitude $z$.

Combining the microscopic equations with the mesoscopic ones, similarly to the procedure described for the dipolar term, one can show that the polarization in the quadrupolar approximation has two additional terms compared to the one in the dipolar approximation (see Supplementary Information Sect. B):

$$P^{SHG}(z) = P^{dip}(z) + P^{grad}(z) + P^{quad}(z). \qquad (13)$$

The first term is the one obtained in the dipolar approximation defined in Eq.(8).

The second involves all the contributions related to the gradient of the electric field:

$$P^{\text{grad.}}(z) = \frac{1}{2} \rho(z) \mathcal{B}^g(z) : e(z) \frac{\mathrm{d}e(z)}{\mathrm{d}z}, \qquad (14)$$

where $\mathcal{B}^g$ is a 3x3x3 tensor defined as $\mathcal{B}_{ijk}^g = \mathcal{B}_{ijkz}^{dq} - \mathcal{B}_{zijk}^{qd}$ where $\mathcal{B}^{dq}$ and $\mathcal{B}^{qd}$ are the molecular hyperpolarizabilities,
and the indices i, j and k can take the values of the three directions x, y and z. Accordingly, we approximate the effective gradient susceptibility tensor as:

$$ \chi^{\text{eff},g}(z) = \frac{1}{2} \rho(z) \mathcal{B}^g(z) : \mathcal{T}(z) \frac{d\mathcal{T}(z)}{dz}. \quad (15) $$

The last term is the quadrupolar one, i.e. the one due to the spatial evolution of the quadrupolar hyperpolarizability at the microscopic scale, and can be approximated as follows:

$$ P^{\text{quad.}}(z) = -\frac{1}{4} \frac{d[\rho(z)\mathcal{B}_z^{qd}(z)]}{dz} : e(z)e(z). \quad (16) $$

with $\mathcal{B}_z^{qd}$ a 3x3x3 matrix representing $\mathcal{B}^{qd}$ with the first element set to the z-axis due to the projection onto the spatial derivative in Eq. 12. The associated quadrupolar effective susceptibility $\chi^{\text{eff},q}(z)$ is given similarly by:

$$ \chi^{\text{eff},q}(z) = -\frac{1}{4} \frac{d[\rho(z)\mathcal{B}_z^{qd}(z)]}{dz} : \mathcal{T}(z)\mathcal{T}(z). \quad (17) $$

We will compare the different contributions in part IV B.

### D. Toward second harmonic detection: the polarization sheet

Once the mesoscopic polarization oscillating at the second harmonic frequency is obtained for a given perturbation, one has to compute the harmonic electric field that is emitted and then transmitted towards the detector.

The radiative equations at the second harmonic frequency are solved at the mesoscopic level between the bulk phase, the gas phase and a 'polarization sheet', which acts as the source of the second harmonic signal. We thus use the framework of the common "three-layer model"$^{28,29}$. Briefly, the mesoscopic polarization of the sheet is the source of a harmonic field in the interfacial layer. Then, the harmonic electric field in the vapor phase is computed through Fresnel factor coefficients at the harmonic frequency to describe transmission across the interface. The intensities for a P or S polarization of the second-harmonic light are then computed as a function of the polarization angle of the incident light.

The main parameters for this three layer model are the optical index of the vapor phase ($N_1 = 1$), the liquid water ($N_2 = 1.34$) the polarization sheet $n_m$ ($n_m = 1.19$) at the second-harmonic frequency (corresponding wavelength is 405 nm in this study). We will discuss the approximation made in this step as well as possible ways of improvement in Section IV D.

### III. Material and methods

In this section, we describe the numerical and experimental procedures used regarding the second harmonic generation at the liquid gas interface.

### A. Numerical procedures

Following our previous works$^{22,25,30}$, a multi-scale procedure has been used to compute the second harmonic molecular properties of the water molecule at a liquid-gas interface. To do so, we rely on a classical Molecular Dynamics (MD) calculation, followed by a Quantum Mechanics Molecular Model (QM/MM) procedure using our homemade open-access software FROG$^{31}$.

#### 1. MD details

The trajectory used is the same as in our first work on the liquid-gas interface, see Ref 22 and its Supplementary information for more details. Briefly, the LAMMPS$^{32}$ software, V.11.08.2017, has been used to perform the MD simulation along with the rigid TIP4P/2005 water force field$^{33}$. Using equilibrated bulk configurations, 9000 rigid TIP4P/2005 water molecules are placed in a simulation box ($5 \times 5 \times 40 \text{ nm}^3$) to form a water film, about 10 nm thick. The thickness of this system is large enough so that the middle part of the fluid film effectively simulates bulk water. The target temperature is 300 K using a Nose-Hoover thermostat with the period of the temperature fluctuations set to $\tau = 0.4$ ps. A first 3 ns of NVT equilibration is performed before the 20 ns NVT production run with a time step of 2 fs. Note that the center of mass of the whole system is fixed at every timestep. Both electrostatic and Lennard-Jones intermolecular interactions are computed using the long-range PPPM formalism$^{34,35}$.

#### 2. QM/MM details

The molecular hyperpolarizabilities were computed and analyzed with FROG$^{31}$. Briefly, we select a given number of water molecules in each MD snapshot and build their specific electrostatic environment made by other water molecules. Then, its response tensor for SHG in this electrostatic environment is computed thanks to the Polarizable Environment scheme$^{36,37}$ of the quantum chemical software DALTON$^{38}$.

We use the PE long-range embedding already described in our previous work$^{25}$. The direct part is considered up to 10 Å whereas the implicit part up to 20 Å: the electrostatic environment is thus included in the QM calculation up to 20 Å. The quantum box contains only one molecule per calculation with an extra electrostatic field describing the electrostatic environment. We use DFT with the CAM-B3LYP functional and the d-aug-cc-pVTZ basis set alongside the response approach provided in DALTON to compute the hyperpolarizability tensors at the optical (for a fundamental of 800 nm) and static frequency. We have shown in previous work$^{22}$ that it leads to a good comparison with CCSD results for the dipolar hyperpolarizability, and in Supp. Mat. Sect. C for the quadrupolar terms in the gas phase$^{39}$.

We use a maximum of 160 MD frames separated by 10 ns, to ensure time-uncorrelation, to perform these QM/MM calcu-
lations. The molecules for which QM/MM level calculations are performed are selected using the floating interface layer selection from FROG: the first 4 molecular layers are treated based upon the PYTIM package$^{40}$. The molecules involved in these 4 layers spread about 8 Å in the direction z, which corresponds to a total of 97625 configurations.

The interface position is defined as the Gibbs dividing altitude: the density profile is fitted by the usual tanh function and $z = 0$ is the point where the density is equal to half of the liquid density. The space is discretized in the z-direction with a 1Å step. Then, the reported data are averaged on a rectangular slab measuring [50x50x1]Å$^3$. The number of QM/MM calculation per altitude is reported in Supp. Mat. Sect. D. One can see that more than 12 000 calculations were done for molecules in the bulk. The sampling continuously decreases along the interface to reach fewer than 1000 calculations at $z = 2.5$ Å. In that case, the sampling is too small to give reliable results and the density is less than 0.05 kg/L.

Two post-treatments are noticeable: (1) Susceptibility tensor components that are expected to be equal because of symmetry reasons are averaged to increase accuracy. (2) We have fitted the evolution of the different components of $\langle \rho \mathcal{B}^{qd} : ee \rangle$ as a function of z with smooth functions to estimate its z-derivative involved in the quadrupolar polarization (see Eq.(16)).

## B. Experimental setup

Standard polarization angle-resolved S-SHG intensity measurements were realized.

The experimental setup was already described in previous studies$^{27,41}$. Briefly, the incident light is generated by a Ti:Sapphire femtosecond laser beam, which delivers pulses at a wavelength of 810 nm, a temporal length of 100 fs, at a repetition rate of 80 MHz. The laser beam is linearly polarized and the angle of polarization $\alpha$ is controlled with a motorized half-wave plate. The incident light is gently focused on pure air-water interface by a 5 cm-lens and the generated light is collected by a 10 cm-lens. After removing the reflected fundamental light with a blue filter, the harmonic light is analyzed in polarization by a half-wave plate followed by a polarized-cube to fix the outlet polarization in P (in the incidence plane) or S (perpendicular to the incidence plane). The harmonic light intensity is collected through a spectrometer (that allows removing any artifact notably due to fluorescence) and a cooled CCD camera. The incidence angle is fixed at 70° and the numerical aperture of the incident beam induced a distribution of 2° around this averaged value.

A typical polarization plot for the neat air-water interface is represented in Fig. 3.

These plots can be analyzed using the standard form of the SHG intensity in the electric dipole approximation as a function of the input polarization angle $\alpha$ from the so-called three-layers model$^{27–29}$:

$$I_P(\alpha) \propto |(a_2\chi_{XXZ} + a_3\chi_{ZXX} + a_4\chi_{ZZZ}) \cos^2\alpha + a_5\chi_{ZXX} \sin^2\alpha|^2 \quad \text{frame with :}$$

$$I_S(\alpha) \propto |a_1\chi_{XXZ} \sin(2\alpha)|^2 \quad (18)$$

<figure><img src="image_4.png" /><figcaption>FIG. 3. Normalized SHG intensity measured on pure air/water interface as a function of the inlet polarization angle. The incident angle is fixed at 70°±2°. The black continuous curves correspond to a fit following Eq. (18) (see main text for fit parameter values).</figcaption></figure>

where $a_i$ typically contains the incident angle and Fresnel coefficients whereas $\chi_{ijk}$ is the $ijk$ component of the second order susceptibility tensor (see Supp. Mat. Sect. E for more details).

Fitting the experimental measurements with this expression yields $\chi_{XXZ}/\chi_{ZZZ} = 0.37 \pm 0.01$ and $\chi_{ZXX}/\chi_{ZZZ} = 0.09 \pm 0.01$, in close agreement with results obtained in the literature (see for example Refs. 19,42 and references therein).

# IV. Results and discussions

## A. Hyperpolarizabilities

We first report the evolution of the different hyperpolarizabilities as a function of the altitude.

### a. Dipolar hyperpolarizabilities

The evolution of dipolar first order hyperpolarizability was already reported in our previous work$^{22}$. Since then, however, we have proved the considerable impact of the electrostatic electric field generated by the neighbors on this hyperpolarizability $\beta^{dd}$ in bulk liquid water$^{25}$. We also developed a methodology to improve the solvation model for $\beta^{dd}$ in bulk water by including a long-range correction taking into account electrostatic effects$^{25}$. Briefly, we compute the hyperpolarizability modeling neighboring molecules as point charges up to a radius $R_c$, and add a long-range correction proportional to the second order hyperpolarizability of water in bulk liquid water $\gamma$. This leads to $\beta^{dd} = \beta^{dd}[R_c] + \gamma: \Delta E^{DC}[R_c]$, where $\Delta E^{DC}[R_c]$ is the error made on the electrostatic field if only the neighbors up to $R_c$ are taken into account: $\Delta E^{DC}[R_c] = E^{DC}[R_c \to +\infty] - E^{DC}[R_c]$.

Here, the same correction is proposed but in the laboratory

$$\mathcal{B}^{dd} = \mathcal{B}^{dd}[R_c] + \Gamma: \Delta E^{DC}[R_c] \quad (19)$$
where $\mathcal{B}^{dd}[R_c]$ corresponds to the hyperpolarizability computed for a finished environment size $R_c$, and $\Gamma$ the second order hyperpolarizability in the laboratory frame.

The value of $E^{DC}[R_c \to \infty]$ is obtained by an extrapolation of the value of $E^{DC}[R_c]$ at $R_c$ within the limit of our MD simulation box, and fitting with a $1/R_c$ function (as expected from ordered dipole on a 2D interface):

$$E^{DC}[R_c] = E^{DC}[R_c \to \infty] - \frac{B}{R_c} \qquad (20)$$

<figure><img src="image_5.png" /><figcaption>FIG. 4. Electrostatic field as a function of $1/R_c$ at different altitudes $z$. The continuous lines are fits based on Eq. 20. Inset: evolution of the extrapolated value $E^{DC}[R_c \to \infty]$ (value when $1/R_c = 0$ in the main figure) with the altitude.</figcaption></figure>

The fit were made at different altitudes (see Fig. 4), and the evolution of the total electric field $E^{DC}[R_c \to \infty]$ is represented in the inset of Fig. 4.

We also tested for the available values of $R_c$ whether the two non-vanishing and independent second order hyperpolarizability components ($\Gamma_{ZZZZ}$ and $\Gamma_{ZXXZ}$) appearing in Eq. 19 are independent on the altitude and on $R_c$ (see Supp. Mat. Sect. F. It turns out that the second order hyperpolarizability does not evolve near the interface as it is hardly modified by the electrostatic environment, as already shown for bulk water$^{25}$. We thus consider a constant value for this tensor with $\Gamma_{ZXXZ} = 580$ a.u. and $\Gamma_{ZZZZ} = 1740$ a.u.

### b. Quadrupolar hyperpolarizabilities

The hyperpolarizabilities involved in the quadrupolar terms are also computed in the static limit and at 800nm, both in the molecular frame and the laboratory one. Figures 5 and 6 represent the evolution with the altitude of some elements of the dipolar-quadrupolar hyperpolarizability respectively in the molecular and laboratory frame.

The values computed at $z = -6\AA$ are comparable to the ones we have computed for bulk liquid water, see values in Supp. Mat. Sect. C, indicating that the bulk properties are well described at this altitude.

It appears that both $\beta^{qd}$ and its standard deviation hardly evolve through the interface. As expected, the quadrupolar-dipolar hyperpolarizabilities follow the same evolution: at this

<figure><img src="image_6.png" /><figcaption>FIG. 5. Evolution of some non-vanishing components of the dipolar hyperpolarizability in the molecular frame (see Supp.Mat., Fig. 1 for axis definition) at the quadratic order as a function of the altitude. Data are computed for a wavelength of 800 nm, with a spatial step of $1\AA$. The continuous lines represent the average value, the shadows the standard deviations of the distribution, the dots the values in the bulk, and the dashed line the values in vacuum.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>FIG. 6. Evolution of some non-vanishing components of the dipolar hyperpolarizability in the laboratory frame at the quadratic order as a function of the altitude. Data are computed for a wavelength of 800 nm, with a spatial step of $1\AA$. The continuous lines represent the average value, the shadows the standard deviations of the distribution, and the dots the values in the bulk.</figcaption></figure>

non-resonant optical frequency, the quadrupolar hyperpolarizabilities almost respect the Kleinman symmetry.

To conclude, neglecting the variation of the quadrupolar molecular hyperpolarizabilities is a valid approximation, both in the bulk and at the liquid-gas interface.
## B. Evolution of the effective susceptibility tensor with the altitude

Fig. 7 reports the evolution of the relevant components of the three effective susceptibility tensors $\chi^{\text{eff},d}$, $\chi^{\text{eff},g}$, and $\chi^{\text{eff},q}$ (see sec.II C), as a function of the altitude. We recall that these effective susceptibility tensors are the factor between the incident electric field in vacuum $E$ and their respective mesoscopic polarization $P$. They are products of the interface susceptibility tensor and the corrections to compute the local field effects at the fundamental frequency.

The dipolar components depicted in Fig.7(a) are of the order of magnitude of $10^{-2}$a.u. The XXZ and XZX components appear to be similar, in disagreement with the experimental observations. Omitting the ZZZ component, the dipolar contribution is limited to the first 5Å from the interface.

Fig. 7(b) shows that the gradient term is of the order of $10^{-1}$a.u., which is 10 times larger than the dipolar term. By definition (see eq.(14)), this contribution is governed by the gradient of the electric field across the interface. All the contributions are non-vanishing on a restricted zone around the interface, typically 6Å. This contribution has different components along XXZ and XZX.

Most importantly, the quadrupolar term plotted on Fig. 7(c) is of the order of the a.u., i.e. 50 to 100 times the magnitude of the dipolar term. This term is dictated by the density evolution across the interface, which dominates the response. For this gradient contribution, the XXZ and XZX components are similar.

In agreement with literature¹⁹, we observe that the quadrupolar contribution is dominant. Moreover, the difference between the XXZ and XZX component is due to the electric field gradient at the interface.

Summing the dipolar, gradient and quadrupolar contributions at each altitude and integrating this total $\chi^{\text{eff}}(z)$ along the interface leads to the ratios: $\chi_{XXZ}^{\text{eff}}/\chi_{ZZZ}^{\text{eff}} = 0.300$ and $\chi_{ZXX}^{\text{eff}}/\chi_{ZZZ}^{\text{eff}} = 0.105$. One can notice that the ratio between the XXZ and the ZZZ components of the effective susceptibility is under-evaluated at this step (≈ -20%) compared to the experimental one. On the contrary, the ratio between the ZXX and ZZZ components is surprisingly close to the experimental value considering the experimental uncertainty. However, these ratios cannot be directly compared to the experimental ones since they also contain the local field factor. In the following, we shall rather compare the polarization-resolved SHG intensity obtained numerically to the experimental ones.

## C. Computation of the polarization-resolved S-SHG intensity

The SHG intensity is finally computed taking into account all the elements previously described: the dipolar, gradient and quadrupolar terms are considered, with an angle of incidence of 70°. The evolution of this intensity as a function of the incident polarization angle $\alpha$ is represented as continuous blue lines when the outlet polarization is $P$ and dashed red lines when the outlet polarization in $S$ on Fig. 8 and compared

<figure><img src="image_8.png" /><figcaption>FIG. 7. Evolution of the ZZZ (green rounds), XXZ (cyan triangles) and XZX (magenta stars) components of the effective susceptibility tensor for the (a) dipolar $\chi^{\text{eff},d}$; (b) gradient $\chi^{\text{eff},g}$; and (c) quadrupolar $\chi^{\text{eff},q}$ contribution as a function of the altitude. Note that the y-scales are largely different between graphs. Color curves correspond to the different components and symbols are guides to identify curves.</figcaption></figure>

to experimental data (dots). Numerically, we study the effect of the variation of the angle of incidence from 68° to 72° to reproduce the distribution on this angle coming from the numerical aperture of the lens.

The numerical description reproduces well the experimental curve with a correct shape, and a good prediction of amplitude of the S-polarized curve. The smaller the incident angle,
<figure><img src="image_9.png" /><figcaption>FIG. 8. Comparison between experimental (left, markers) and numerical (right, lines) polarization resolved SHG intensity computed with all the contributions at an incidence of 70°. The shadows represent the effect of an uncertainty on the angle of incidence of +/-2°.</figcaption></figure>

the larger the maximum intensity in S-outlet polarization: an angle of 68° allows the experimental S-polarized curve to be reproduced well. However, the minimum of the P-curve is not well reproduced. These results are very similar to the ones by Foucaud et al.¹⁹.

To go further, we fit the numerical curves with Eq. (18) (the $a_i$ are computed with an incident angle of 70° and an optical index of the interface sheet $n_m = 1.19$) to extract the susceptibility ratio with the same procedure as for the experimental data. We obtain $\chi_{xxz}/\chi_{zzz} = 0.362 \pm 0.006$ and $\chi_{zxx}/\chi_{zzz} = 0.173 \pm 0.003$.

The predicted $\chi_{xxz}/\chi_{zzz}$ is close to the experimental one, with an error smaller than 1%, whereas the predicted $\chi_{zxx}/\chi_{zzz}$ is almost twice the experimental result. We will discuss the possible origins of this discrepancy in the following.

### D. Possible origins of the discrepancies between experimental and numerical results

We discuss here the four main limits of our approach, and possible improvements for the predicted value of $\chi_{zxx}/\chi_{zzz}$.

#### a. Local field at $\omega$ in the interfacial zone

The description of the field felt by the molecules at $\omega$ in step A is expected to take into account the collective response of the molecular assembly via the definition of a dielectric constant that varies with the altitude. The underlying assumption is that the molecules are distributed at random positions with random orientations in a volume⁴³. Equation (4) for the transfer matrix $\mathcal{T}(z)$, albeit physically sound, is eventually obtained heuristically. Besides the fundamental interest in deriving rigorously an expression for $\mathcal{T}(z)$, the model could be improved by a more careful treatment of the configurational variability of the molecular assembly, including the spatial

correlations in the molecule arrangement and the preferential orientation of the molecules that both vary with the altitude $z$. This would effectively lead to an optical anisotropy that varies with $z$ and a modification of the field felt by the molecules, which may be particularly significant at grazing angles where polarization effects are important.

#### b. Radiative intermolecular coupling at $2\omega$

Equations (6) and (7) in step B, which lead to the definition of a mesoscopic polarization in step C, neglect the fact that the dipole and quadrupole moments induced in a molecule $n$ can be affected by the field radiated by all other molecules $m \neq n$ at the second harmonic frequency. This effect may yet be expected because molecules are polarizable objects also at $2\omega$. This issue might be addressed numerically using multiple light scattering techniques (in the spirit of T-matrix methods⁴⁴), leading to the computation of “dressed” hyperpolarizabilities, or analytically (within the dipolar approximation) using a classical Green’s tensor formalism⁴⁵, although the proper consideration of the spatial arrangement and preferential orientation of the molecules would be a challenge. Nevertheless, we feel that taking the radiative coupling between molecules at $2\omega$ into account could significantly affect the mesoscopic polarization in the interfacial zone.

#### c. Second-harmonic light generation and the three-layer model

Light emission in general is known to be very sensitive to the local environment⁴⁶ and we feel that some efforts are required to address the problem of second-harmonic light generation more carefully. For instance, in part C, we present the polarization-resolved S-SHG intensity with an optical index of the interfacial layer $n_m = 1.19$. Varying this optical index from 1.15 to 1.25 ($\pm 5\%$) leads to a strong variation of the S-polarized curve without large change of the minima P-polarized curve (see Supp. Matt. Sect. G) as illustrated by the susceptibility tensor ratios $\chi_{xxz}/\chi_{zzz} = 0.361 \pm 0.03$ (10% of variation) and $\chi_{zxx}/\chi_{zzz} = 0.173 \pm 0.002$ (1% of variation). Some refinement of the Fresnel factor⁴⁷,⁴⁸ were suggested in literature, with a correction of the one involved in the normal component of the electric field. The value of the effective susceptibility tensor ratios computed in Part IV B can be used in the future to discuss this asymmetric correction.

Perhaps more importantly, the generation of the second-harmonic signal in step D, which relies on the classical three-layer model, is currently not fully consistent (albeit physically sound) with the definition of a mesoscopic polarization at $2\omega$ in step C, where the polarization explicitly depends on the altitude $z$. Rigorous light scattering theory⁴⁵ shows that the effective material parameter (dielectric constant or relative permittivity) of thin composite layers are tensorial quantities in general, implying that the local, effective refractive index of the molecular assembly (in which the emitted field propagates before reaching the far field) should be described by a $3 \times 3$ matrix with unequal diagonal terms rather than by a scalar. This is all the more true considering the preferential orientation of the molecules at certain altitudes $z$, as discussed previously. A careful treatment of the effective material parameters at $2\omega$ could affect the SHG signal. The problem at hand may be tackled starting from an analytical derivation of the effec-
tive material parameter of infinitesimal layers of the interfacial zone. The number of layers should be sufficiently large to describe the smooth profile of the interfacial zone. Thin-film stack computations would then lead to a prediction of the second harmonic signal.

### d. DC electric field at the interface

Our QM/MM calculations are performed for an ideal liquid/vapor interface of pure water. This model is simplified relative to the experimental setup using a real air/liquid interface. Notably, the real surface potential is very complicated to measure or compute¹⁶. We therefore question the amplitude of the DC field computed based on our MD simulations, and its effect on the polarization-resolved S-SHG intensity.

To illustrate the impact of this approximation on the polarization curves, we tune the DC electric field at the interface. As a first try, we artificially change the amplitude of the long range DC electric field $E^{DC}$ computed in Sec.IV A without changing its spatial expansion. The effect of a DC field varying from 0 to 10 times is illustrated on Fig. 9 for a angle of incidence $\theta = 70°$. Increasing the static electric field amplitude

<figure><img src="image_10.png" /><figcaption>FIG. 9. Comparison between experimental (left, markers) and numerical (right, lines) polarization resolved SHG intensity computed with all the contributions at an incidence of 70°. The shadows represent the variation induced when the amplitude of the DC electric field at the interface varies from 0 to 10 times the one emerging from the MD ($E^{DC}[R_c \to \infty]$) of Fig. 4).</figcaption></figure>

by this 10-factor strongly affects the minimum value of the P-resolved curve, and the ratio $\chi_{XXZ}/\chi_{ZZZ}$ is as low as 0.354, close to the experimental value of 0.368. The amplitude of the S-resolved curve is not much impacted by the static field. The other ratio, $\chi_{ZXX}/\chi_{ZZZZ} = 0.106$, remains close to the experimental one at 0.091.

These results hint that our current numerical description fails to describe accurately the electrostatic field near the interface. This might originate in the limits of the MD model, but also if a surface charge is present experimentally and not accounted for in the numerical model – for instance ions or pollutants. We will explore further this effect both numerically and experimentally.

Further numerical work must be done to compute a more realistic electrostatic field close the interface. For stronger static field, and with a deeper effect in the bulk phase longer range, the phase effects in the Debye layer may become important. Experimentally, this hypothesis will be assessed using phase-sensitive SHG experiments⁴⁸ or by studying the evolution of S-SHG intensity as a function of the surface charge (e.g. addition of surfactants⁴¹).

## V. Conclusion

We presented in this work a numerical description at the microscopic level of the polarization-resolved S-SHG intensity of vapor/pure water liquid interface. The QM/MM approach already used to describe Second Harmonic Scattering on pure water in bulk³⁰ was employed to compute the first hyperpolarizabilities at the dipolar and quadrupolar orders as a function of the distance to the interface. We also computed the evolution of the second hyperpolarizability of water in the condensed phase as a function of the altitude, confirming that it hardly depends on the electrostatic environment²⁵. The local fundamental electric field was taken into account to compute the mesoscopic polarizability as a function of the altitude, whereas the harmonic intensity was computed through the more classical 3-layer model. The numerical polarization-resolved S-SHG intensity obtained with our approach reproduces well the experimental results even if there is still a discrepancy to describe the minimum intensity of the P-polarized curve. We discussed the possible causes of this discrepancy, leaving this for future work. In particular, preliminary results suggest that a more realistic description of the static electric field close to the interface may improve the accuracy of the numerical prediction. The generation of the second-harmonic field from the mesoscopic polarization shall also have to be addressed rigorously.

Our results raise good hope that some contributions to the S-SHG intensities may be modeled analytically, without the need for QM/MM calculations, and pave the way to describe more complex systems such as liquid/liquid interfaces. In particular, we have shown that the quadrupolar first hyperpolarizabilities and second dipolar hyperpolarizability tensors hardly depend on the electrostatic environment and that they may be assumed constant for all the water molecules. Along with mean-field models for the fluid densities, surface electrostatic field, and fundamental electric field evolution, one can have a decent estimation of the induced second harmonic polarization. However, the extraction of useful information about aqueous surface from S-SHG intensities is based on the quantification of the dipolar contribution *relative* to the (dominant) quadrupolar one. Some interesting surface observables such as molecular orientation will be contained in small variations of the signal due to the dipolar contribution only. Other observables such as surface charges may influence both contributions. Coupling experiments and molecular modeling will help to go beyond their respective shortcomings.
# Supplementary Materials

The supplementary material details the following points:

A. Reminders of the response theory

B. Separation of the polarization density into different terms

C. Quadrupolar hyperpolarizability in the liquid phase

D. Number of QM/MM calculations per altitude

E. *ai* coefficients to fit polarization-resolved S-SHG experiments

F. Variation of the dipole second hyperpolarizability components in the laboratory frame

G. Change in polarization-resolved SHG intensity with the optical index of the interface

# Acknowledgements

The authors thank Pr. P-F Brevet for fruitful discussions. They gratefully acknowledge support from the NanOpTec center for experimental supports and the PSMN (Pôle Scientifique de Modélisation Numérique) of the ENS de Lyon for the computing resources. O.B thanks the financial support of the French National Agency for Research under project SOLSTICE (ANR-21-CE30-0007).

# Data

The data that support the findings of this study are available from the corresponding author upon reasonable request.

# References

<footnote n="1">¹Y. R. Shen, "Surface properties probed by second-harmonic and sum-frequency generation," Nature **337**, 519-525 (1989).</footnote>

<footnote n="2">²K. B. Eisenthal, "Liquid interfaces probed by second-harmonic and sum-frequency spectroscopy," Chemical Reviews **96**, 1343-1360 (1996).</footnote>

<footnote n="3">³F. M. Geiger, "Second harmonic generation, sum frequency generation, and $\chi^{(3)}$: Dissecting environmental interfaces with a nonlinear optical swiss army knife," Annual Review of Physical Chemistry **60**, 61-83 (2009).</footnote>

<footnote n="4">⁴P.-F. Brevet, Surface Second Harmonic Generation (PPUR presses poly-techniques, 1997).</footnote>

<footnote n="5">⁵R. J. Tran, K. L. Sly, and J. C. Conboy, "Applications of surface second harmonic generation in biological sensing," Annual Review of Analytical Chemistry **10**, 387-414 (2017).</footnote>

<footnote n="6">⁶R. Shen, "Second harmonic and sum-frequency spectroscopy," (World Scientific Publishing Company, 2023).</footnote>

<footnote n="7">⁷P. Guyot-Sionnest, W. Chen, and Y. R. Shen, "General considerations on optical second-harmonic generation from surfaces and interfaces," **33**, 15 (1986).</footnote>

<footnote n="8">⁸Y. Chen, H. I. Okur, N. Gomopoulos, C. Macias-Romero, P. S. Cremer, P. B. Petersen, G. Tocci, D. M. Wilkins, C. Liang, M. Ceriotti, and S. Roke, "Electrolytes induce long-range orientational order and free energy changes in the h-bond network of bulk water," Science Advances **2**, 1-9 (2016).</footnote>

<footnote n="9">⁹S. Pezzotti, D. R. Galimberti, and M. P. Gaigeot, "2d h-bond network as the topmost skin to the air-water interface," Journal of Physical Chemistry Letters **8**, 3133-3141 (2017).</footnote>

<footnote n="10">¹⁰G. Gonella, C. Lütgebaucks, A. G. F. de Beer, and S. Roke, "Second harmonic and sum-frequency generation from aqueous interfaces is modulated by interference," The Journal of Physical Chemistry C **120**, 9165-9173 (2016).</footnote>

<footnote n="11">¹¹Y. Litman, K. Y. Chiang, T. Seki, Y. Nagata, and M. Bonn, "Surface stratification determines the interfacial water structure of simple electrolyte solutions," Nature Chemistry **16**, 644-650 (2024).</footnote>

<footnote n="12">¹²M. D. Boamah, P. E. Ohno, F. M. Geiger, and K. B. Eisenthal, "Relative permittivity in the electrical double layer from nonlinear optics," Journal of Chemical Physics **148** (2018), 10.1063/1.5011977.</footnote>

<footnote n="13">¹³P. Guyot-Sionnest and Y. R. Shen, "Bulk contribution in surface second-harmonic generation," Physical Review B **38**, 7985 (1988).</footnote>

<footnote n="14">¹⁴S. Sun, J. Schaefer, E. H. Backus, and M. Bonn, "How surface-specific is 2nd-order non-linear spectroscopy?" Journal of Chemical Physics **151** (2019), 10.1063/1.5129108.</footnote>

<footnote n="15">¹⁵M. E. Didier, O. B. Tarun, P. Jourdain, P. Magistretti, and S. Roke, "Membrane water for probing neuronal membrane potentials and ionic fluxes at the single cell level," Nature Communications **9**, 1-7 (2018).</footnote>

<footnote n="16">¹⁶R. Hartkamp, A. L. Biance, L. Fu, J. F. Dufrêche, O. Bonhomme, and L. Joly, "Measuring surface charge: Why experimental characterization and molecular modeling should be coupled," Current Opinion in Colloid and Interface Science **37**, 101-114 (2018).</footnote>

<footnote n="17">¹⁷P. E. Ohno, H. Chang, A. P. Spencer, Y. Liu, M. D. Boamah, H.-f. Wang, and F. M. Geiger, "Beyond the gouy-chapman model with heterodyne-detected second harmonic generation," The Journal of Physical Chemistry Letters **10**, 2328-2334 (2019).</footnote>

<footnote n="18">¹⁸K. Shiratori, S. Yamaguchi, T. Tahara, and A. Morita, "Computational analysis of the quadrupole contribution in the second-harmonic generation spectroscopy for the water/vapor interface," The Journal of chemical physics **138** (2013).</footnote>

<footnote n="19">¹⁹Y. Foucaud, J. F. Dufrêche, B. Siboulet, M. Duvail, A. Jonchère, O. Diat, and R. Vuilleumier, "Why local and non-local terms are essential for second harmonic generation simulation?" Physical Chemistry Chemical Physics **24**, 12961-12973 (2022).</footnote>

<footnote n="20">²⁰C. Liang, G. Tocci, D. M. Wilkins, A. Grisafi, S. Roke, and M. Ceriotti, "Solvent fluctuations and nuclear quantum effects modulate the molecular hyperpolarizability of water," Phys. Rev. B **96**, 1-6 (2017), arXiv:1705.01649.</footnote>

<footnote n="21">²¹J. Kongsted, A. Osted, K. V. Mikkelsen, and O. Christiansen, "Second harmonic generation second hyperpolarizability of water calculated using the combined coupled cluster dielectric continuum or different molecular mechanics methods," J. Comp. Phys. **120**, 3787-3798 (2004).</footnote>

<footnote n="22">²²G. Le Breton, O. Bonhomme, P.-F. Brevet, E. Benichou, and C. Loison, "First Hyperpolarizability of Water at the Air-Vapor Interface: a QM/MM Study Questions Standard Experimental Approximations," Physical Chemistry Chemical Physics **23**, 24932-24941 (2021).</footnote>

<footnote n="23">²³T. N. Ramos and B. Champagne, "Investigation of the Second Harmonic Generation at the Water-Vacuum Interface by Using Multi-Scale Modeling Methods," ChemistryOpen **12**, e202200045 (2023).</footnote>

<footnote n="24">²⁴Y. Foucaud, B. Siboulet, M. Duvail, A. Jonchère, O. Diat, R. Vuilleumier, and J. F. Dufrêche, "Deciphering second harmonic generation signals," Chemical Science **12**, 15134-15142 (2021).</footnote>

<footnote n="25">²⁵G. Le Breton, O. Bonhomme, E. Benichou, and C. Loison, "First hyperpolarizability of water in bulk liquid phase: long-range electrostatic effects included via the second hyperpolarizability," Physical Chemistry Chemical Physics **24**, 19463-19472 (2022).</footnote>

<footnote n="26">²⁶A. Morita, Theory of Sum Frequency Generation Spectroscopy (2018).</footnote>

<footnote n="27">²⁷G. Martin-Gassin, E. Benichou, G. Bachelier, I. Russier-Antoine, C. Jonin, and P. Brevet, "Compression induced chirality in dense molecular films at the air-water interface probed by second harmonic generation," Journal of Physical Chemistry C **112**, 12958-12965 (2008).</footnote>

<footnote n="28">²⁸P. F. Brevet, "Phenomenological three-layer model for surface second-harmonic generation at the interface between two centrosymmetric media," J. Chem. Soc., Faraday Trans. **92**, 4547-4554 (1996).</footnote>

<footnote n="29">²⁹V. Mizrahi and J. E. Sipe, "Phenomenological treatment of surface second-harmonic generation," J. Opt. Soc. Am. B **5**, 660-667 (1988).</footnote>

<footnote n="30">³⁰G. Le Breton, O. Bonhomme, E. Benichou, and C. Loison, "Liquid water: when hyperpolarizability fluctuations boost and reshape the second harmonic scattering intensities," The Journal of Physical Chemistry Letters **14**, 4158-4163 (2023).</footnote>

<footnote n="31">³¹G. Le Breton, O. Bonhomme, E. Benichou, and C. Loison, "FROG: Exploiting all-atom molecular dynamics trajectories to calculate linear and</footnote>
non-linear optical responses of molecular liquids within Dalton's QM/MM polarizable embedding scheme," Journal of Chemical Physics **160** (2024), 10.1063/5.0203424.

<footnote n="32"><sup>32</sup>S. Plimpton, "Fast Parallel Algorithms for Short – Range Molecular Dynamics," J. Comput. Phys. **117**, 1–19 (1995).</footnote>

<footnote n="33"><sup>33</sup>J. L. Abascal and C. Vega, "A general purpose model for the condensed phases of water: Tip4p/2005," J. Chem. Phys. **123**, 234505 (2005).</footnote>

<footnote n="34"><sup>34</sup>R. E. Isele-Holder, W. Mitchell, and A. E. Ismail, "Development and application of a particle-particle particle-mesh Ewald method for dispersion interactions," J. Chem. Phys. **137** (2012).</footnote>

<footnote n="35"><sup>35</sup>R. E. Isele-Holder, W. Mitchell, J. R. Hammond, A. Kohlmeyer, and A. E. Ismail, "Reconsidering dispersion potentials: Reduced cutoffs in mesh-based ewald solvers can be faster than truncation," J. Chem. Theory Comput. **9**, 5412–5420 (2013).</footnote>

<footnote n="36"><sup>36</sup>J. M. Olsen, K. Aidas, and J. Kongsted, "Excited states in solution through polarizable embedding," J. Chem. Theory Comput. **6**, 3721–3734 (2010).</footnote>

<footnote n="37"><sup>37</sup>C. Steinmann, P. Reinholdt, M. S. Nørby, J. Kongsted, and J. M. H. Olsen, "Response properties of embedded molecules through the polarizable embedding model," Int. J. Quantum Chem. **119**, e25717 (2019).</footnote>

<footnote n="38"><sup>38</sup>K. Aidas, C. Angeli, K. L. Bak, V. Bakken, R. Bast, L. Boman, O. Christiansen, R. Cimiraglia, S. Coriani, P. Dahle, et al., "The dalton quantum chemistry program system," Wiley Interdisciplinary Reviews: Computational Molecular Science **4**, 269–284 (2014).</footnote>

<footnote n="39"><sup>39</sup>H. Chen, M. Liu, and T. Yan, "Molecular multipoles and (hyper) polarizabilities of water by ab initio calculations," Chemical Physics Letters **752**, 137555 (2020).</footnote>

<footnote n="40"><sup>40</sup>M. Sega, G. Hantal, B. Fábián, and P. Jedlovszky, "Pytim: A python package for the interfacial analysis of molecular simulations," **39**, 2118–2125 (2018).</footnote>

<footnote n="41"><sup>41</sup>E. Benichou, A. Bruyère, E. Forel, O. Bonhomme, and P. Brevet, "Molecular organization and phase transition at the air-water interface investigated by second-harmonic generation," in *Organic Photonic Materials and Devices XVII*, Vol. 9360 (2015) p. 93600T.</footnote>

<footnote n="42"><sup>42</sup>T. T. Pham, A. Jonchère, J.-f. F. Dufrêche, P.-f. F. Brevet, and O. Diat, "Analysis of the second harmonic generation signal from a liquid/air and liquid/liquid interface," Journal of Chemical Physics **146**, 144701 (2017).</footnote>

<footnote n="43"><sup>43</sup>K. Shiratori and A. Morita, "Molecular theory on dielectric constant at interfaces: A molecular dynamics study of the water/vapor interface," The Journal of chemical physics **134** (2011).</footnote>

<footnote n="44"><sup>44</sup>D. W. Mackowski and M. I. Mishchenko, "Calculation of the t matrix and the scattering matrix for ensembles of spheres," J. Opt. Soc. Am. A **13**, 2266–2278 (1996).</footnote>

<footnote n="45"><sup>45</sup>L. Tsang, J. A. Kong, and K.-H. Ding, *Scattering of electromagnetic waves: theories and applications*, Vol. 15 (John Wiley & Sons, 2000).</footnote>

<footnote n="46"><sup>46</sup>L. Novotny and B. Hecht, *Principles of nano-optics* (Cambridge university press, 2012).</footnote>

<footnote n="47"><sup>47</sup>X. Yu, K.-Y. Chiang, C.-C. Yu, M. Bonn, and Y. Nagata, "On the Fresnel factor correction of sum-frequency generation spectra of interfacial water," The Journal of Chemical Physics **158**, 044701 (2023).</footnote>

<footnote n="48"><sup>48</sup>L. Dalstein, K. Y. Chiang, and Y. C. Wen, "Direct Quantification of Water Surface Charge by Phase-Sensitive Second Harmonic Spectroscopy," Journal of Physical Chemistry Letters **10**, 5200–5205 (2019).</footnote>