The transition between the gas-, supercritical-, and liquid-phase behavior is a fascinating topic, which still lacks molecular-level understanding. Recent ultrafast two-dimensional infrared spectroscopy experiments suggested that the vibrational spectroscopy of N2O embedded in xenon and SF6 as solvents provides an avenue to characterize the transitions between different phases as the concentration (or density) of the solvent increases. The present work demonstrates that classical molecular dynamics (MD) simulations together with accurate interaction potentials allows us to (semi-)quantitatively describe the transition in rotational vibrational infrared spectra from the P-/R-branch line shape for the stretch vibrations of N2O at low solvent densities to the Q-branch-like line shapes at high densities. The results are interpreted within the classical theory of rigid-body rotation in more/less constraining environments at high/low solvent densities or based on phenomenological models for the orientational relaxation of rotational motion. It is concluded that classical MD simulations provide a powerful approach to characterize and interpret the ultrafast motion of solutes in low to high density solvents at a molecular level.

## INTRODUCTION

Solvent–solute interactions and the coupled dynamics between solute molecules embedded in a solvent are central for understanding processes ranging from rotational and vibrational energy relaxation to chemical reactivity in solution.^{1,2} Vibrational spectroscopy is a particularly suitable technique for following the structural dynamics of “hot” solutes in their electronic ground state interacting with a solvent environment. Experimentally, 1D- and 2D-infrared (IR) spectroscopies provide measures of the strength and time scales of solvent–solute interactions that couple to the resonant rovibrational excitation.^{3} These effects can also be probed directly by molecular dynamics (MD) simulations.^{4} The comparison of computations and experiments is a sensitive test for the quality of the intermolecular interaction potentials and provides a molecular-level understanding of energy transfer mechanisms. Furthermore, the entire structural dynamics at molecular-level detail is contained in the time-dependent motion of all solution species involved from the MD trajectories and is available for further analysis.

The understanding of solute equilibration in solvents including high-density gases and supercritical fluids (SCFs) is also important from a more practical perspective. For example, controlling the outcome of combustion reactions that are often performed in the supercritical regime requires knowledge of energy transfer rates and mechanisms in high temperature and pressure solutions.^{5,6} For reactions, it has been recognized that characterizing and eventually tuning the physicochemical properties of the solvent can be as important as determining the best catalyst.^{6} Hence, gaining molecular-level insight into the properties of solute–solvent dynamics is central not only for understanding fundamental solution interactions but also for industrial processes where solvents such as ionic, supercritical, or eutectic liquids are employed. On the other hand, the possibility to specifically manipulate the dynamical properties in a solvent system has been already successfully exploited in a wide range of applications.^{7,8}

Linear and nonlinear optical and, in particular, infrared (IR) spectroscopy is a powerful means to characterize the structural and dynamical properties of condensed-phase systems. As an example, recent spectroscopic experiments in the terahertz and infrared region combined with the MD simulations provided an atomistically resolved picture of the structural properties of a eutectic mixture.^{9} More broadly speaking, the structural and dynamical properties of solutions can and have been studied at a molecular level using two-dimensional IR (2DIR) spectroscopy through vibrational energy transfer.^{10–18} The energy transfer rate for exchange of vibrational energy is expected to follow a sixth-power, distance dependent law,^{17} similar to Nuclear Overhauser Enhancement Spectroscopy (NOESY) in nuclear magnetic resonance (NMR) spectroscopy,^{19,20} or Förster energy transfer between electronic chromophores.^{21} As an example, the presence of cross-peaks in a 2D spectrum has been directly related to the formation of aggregated structures.^{11} However, interpretation and more in-depth understanding of the solvent structure requires additional information that can be obtained, for example, from MD simulations.

Recent 2DIR experiments on gas and supercritical phase solutions have shown how rotational energy returns to equilibrium following rovibrational excitation.^{22–24} The dynamics of N_{2}O as the spectroscopic probe surrounded by xenon atoms and SF_{6} molecules has provided valuable information about the density-dependent change in the solvent structure and energy transfer dynamics as the fluid approaches and passes through the near-critical point density region. The 1D and 2D spectra of the N_{2}O asymmetric stretch (*ν*_{3}) in Xe and SF_{6} exhibit a significant dependence on solvent density.^{23} At low density, corresponding to SF_{6} and Xe in the gas phase, the FTIR (1D) band shape of the asymmetric stretch vibration is that of gas-phase N_{2}O with clearly resolved P- and R-band structures, whereas at high solvent density of the liquid Xe or SF_{6} solvent, only a Q-branch-like absorption feature peaked at the pure vibrational transition frequency was observed. Because the solvent density *ρ* can be changed in a continuous fashion between gas and liquid densities, along a near-critical isotherm (*T* ∼ 1.01*T*_{c}), it is also possible to probe the intermediate, supercritical regime of the solvent. 2DIR spectral shapes as a function of waiting time and solvent density exhibit perfectly anti-correlated features that report on the *J*–scrambling or rotational energy relaxation rates.^{22–24}

Given the molecular-level detail provided by quantitative MD simulations, the present work focuses on changes in the IR spectroscopy of N_{2}O embedded in SF_{6} and xenon and their interpretation at a structural level. For this, an accurate kernel-based representation of high-level electronic structure calculations for the spectroscopic reporter (N_{2}O) is combined with state-of-the-art treatment of electrostatic interactions employing the minimally distributed charge model (MDCM),^{25–27} which captures multipolar interactions. Extensive MD simulations for densities corresponding to those used in the recent 1D and 2D experiments were carried out.

The current work is structured as follows: first, the methods are presented. Next, the quality of the interaction potentials is discussed. Then, the density-dependent FTIR spectra are compared with the experiments, the vibrational frequency fluctuation correlation functions are characterized, and the organization of the solvent is analyzed for xenon and SF_{6} as the solvents. Finally, the findings are discussed in a broader context.

## METHODS

### Potential energy surfaces

The intramolecular potential energy surface (PES) of N_{2}O in its electronic ground state (^{1}A′) is provided by a machine-learned representation using the reproducing kernel Hilbert space (RKHS) method.^{28,29} Reference energies at the coupled-cluster level of theory CCSD(T)-F12/aug-cc-pVQZ including triple excitation perturbation were determined on a grid of Jacobi coordinates (*R*, *r*, *θ*) with *r* being the N–N separation, *R* being the distance between the center of mass of the diatom and the oxygen atom, and *θ* being the angle between the two distance vectors. The grid in *r* contained 15 points between 1.6 a_{0} and 3.1 a_{0}, 15 points between 1.8 a_{0} and 4.75 a_{0} for *R*, and a ten-point Gauss–Legendre quadrature between 0° and 90° for *θ*. All calculations were carried out using the MOLPRO package.^{30}

^{28}Reciprocal power decay kernels

*R*and

*r*). Here,

*x*

_{>}and

*x*

_{<}are the larger and smaller values of

*x*and

*x*′, respectively. For large separations, such kernels approach zero according to $\u221d1xn$ (here,

*n*= 6), which gives the correct long-range behavior for neutral atom–diatom type interactions. For the angular coordinate, a Taylor spline kernel $k[2](z,z\u2032)=1+z<z>+2z<2z>\u221223z<3$ is used, where $z=1\u2212cos\u2061\theta 2$ and

*z*

_{>}and

*z*

_{<}are again the larger and smaller values of

*z*and

*z*′, respectively, defined similarly to

*x*

_{>}and

*x*

_{>}.

Intramolecular and intermolecular force field parameters for SF_{6} are those from the work of Dellis and Samios.^{31} Intermolecular interactions are based on Lennard-Jones potentials only, and the parameters are optimized such that the MD simulations of pure SF_{6} reproduce the experimentally observed *PVT* state points for liquid and gas SF_{6}, as well as the states of liquid–vapor coexistence below and supercritical fluid above the critical temperature *T*_{c}, respectively. For xenon, the parameterization for a Lennard-Jones potential from Aziz and Slaman was used that reproduces dilute gas macroscopic properties such as virial coefficient, viscosity, and thermal conductivity over a wide temperature range but not specifically supercritical fluid properties.^{32} As discussed further below, the critical concentration for xenon as determined from the simulations (see Fig. S1) was found to be 5.19 M compared with the experimentally reported value^{33,34} of 8.45 M and variations for the critical temperature *T*_{c} for xenon compared with the experiments can also be anticipated. In comparison, the critical concentration of $\u223c5$ M for SF_{6} from the present work (see below) matches the experimental value of 5.06 M^{33,34} and is close to 5.12 M (obtained from the reported critical density of 0.74 g/ml) from the parameterization study.^{31}

Electrostatic interactions are computed based on a minimally distributed charge model that correctly describes higher-order multipole moments of a molecule.^{25–27} For parameterization, a reference electrostatic potential (ESP) of N_{2}O in the linear equilibrium conformation is computed at the CCSD/aug-cc-pVTZ level using the Gaussian program package.^{35} The optimized MDCM fit reproduces the ESP with a root mean squared error (RMSE) of 0.31 kcal/mol. For SF_{6} in its octahedral equilibrium conformation, the ESP is computed at the MP2/aug-cc-pVTZ level of theory using the Gaussian program, and the RMSE of the fitted ESP from MDCM is 0.11 kcal/mol. Recently,^{27} non-iterative polarization has also been included in MDCM, and this is also used here for N_{2}O, SF_{6}, and Xe. The polarizability of linear N_{2}O computed at the CCSD/aug-cc-pVTZ level is 2.85 Å^{3} (with each atom contributing $\u223c0.95$ Å^{3} per atom), compared with 2.998 Å^{3} from the experiment.^{36} For Xe at the CCSD/aug-cc-pVTZ, the computed value of 2.96 Å^{3} compares with 4.005 Å^{3} from the experiment.^{36} For SF_{6}, the experimentally measured polarizability of 4.49 Å^{3} was used and evenly distributed over the fluorine atoms (0.74 Å^{3} per F atom).^{37}

The atomic van der Waals (vdW) parameters (*ϵ*_{i}, *R*_{min,i}) of N_{2}O were individually optimized by least-squares fitting using the trust region reflective algorithm^{38} to best describe the nonbonded interactions with Xe and SF_{6} from comparing with energies from CCSD/aug-cc-pVTZ calculations for a number of N_{2}O–Xe and N_{2}O–SF_{6} heterodimer structures. Lorenz–Berthelot combination rules of atomic vdW parameters between atom types *i* (solute) and *j* (solvent) ($\u03f5ij=\u03f5i\u03f5j$ at an atom–atom distance of *R*_{min,ij} = *R*_{min,i}/2 + *R*_{min,j}/2) were used throughout. For the reference electronic structure calculations, the grid was defined by center of mass distances between N_{2}O and the solvent with range *r* = [2.0, 10.0] Å, and angles *α* = [0, 180]° with step size of 30° between the N_{2}O bond axis and the N_{2}O-solvent center of mass direction. Solute–solvent dimer structures with interaction energies lower than 5 kcal/mol above the dimer minimum structure were considered, which led to 102 conformations for N_{2}O–Xe, whereas for N_{2}O–SF_{6}, there were 203 structures. Reference interaction energies for fitting the vdW parameters were determined as follows: the total energy of the N_{2}O–solvent (Xe or SF_{6}) pair with the largest separation was considered to be the zero of energy, and energies for all other dimer structures were referred to this reference. Interaction energies for the equilibrium pair structure were found with −0.84 kcal/mol for N_{2}O–Xe and −1.46 kcal/mol for N_{2}O–SF_{6}. Correcting for the basis set superposition error (BSSE) through counterpoise correction (CPC)^{39} reduces the interaction energies by up to $\u223c50%$ to −0.36 and −0.85 kcal/mol, respectively. The alternative “chemical Hamiltonian approach”^{40} was found to yield similar results^{41} as the CPC, which, however, is not recommended for correlated wave function methods such as CCSD.^{42,43} In addition, the corrections due to BSSE are of a similar magnitude as the error in fitting the van der Waals parameters. Therefore, it was decided to not correct the interaction energies for BSSE. Then, the vdW parameters for each atom *i* of N_{2}O were optimized to match the interaction energy predicted by the Chemistry at “Harvard Molecular Mechanics” (CHARMM) program with the reference interaction energies from the electronic structure calculations for the respective dimer conformations. The optimized vdW parameters are given in Table S2.

### Molecular dynamics simulations

Molecular dynamics simulations were performed with the CHARMM program package.^{44} Each system (N_{2}O in Xe and N_{2}O in SF_{6} at given temperature and solvent concentration) was initially heated and equilibrated for 100 ps each, followed by the 10 ns production simulations in the *NVT* ensemble using a time step of 1 fs for the leapfrog integration scheme. The N_{2}O/Xe systems were simulated at a temperature of 291.2 K, and for N_{2}O/SF_{6}, the temperature was 321.9 K, both of which are slightly above the experimental critical temperatures for condensation of xenon and SF_{6}, respectively [*T*_{c}(Xe) = 289.74 K, *T*_{c}(SF_{6}) = 318.76 K].^{23,33,34} A Langevin thermostat (coupling 0.1 ps^{−1}) was used to maintain the temperature constant but was applied only to the solvent (Xe and SF_{6}) atoms. Positions and velocities of snapshots of the simulations were stored every 1 fs for analysis. As the intermolecular vibrational energy transfer is slow,^{24} the structure of N_{2}O was optimized, and new velocities from a Boltzmann distribution at the simulation temperature were assigned to N_{2}O after the heating step. This ensures that the kinetic energies along the asymmetric, symmetric, and bending modes match the thermal energy with respect to the target simulation temperature.

The different simulation systems were prepared according to the conditions used in the experiments.^{22–24} Table S1 summarizes the concentration *c*(N_{2}O) of N_{2}O, molar volumes *V*_{m}, and critical density ratio *ρ** = *ρ*/*ρ*_{c}. The experimentally determined critical densities are *ρ*_{c} = 1.11 g/ml for xenon and *ρ*_{c} = 0.74 g/ml for SF_{6} from which critical concentrations of 8.45 and 5.06 M for xenon and SF_{6} are obtained, respectively.^{33,34} In all setups, the simulation box contains one N_{2}O molecule and 600 Xe atoms or 343 SF_{6} molecules, which corresponds to similar simulation box volumes for similar relative density ratios of the two solvents. In the original parameterization study, a simulation box containing 343 SF_{6} molecules was used to fit temperature–pressure properties.^{31}

In the MD simulations for N_{2}O in SF_{6}, electrostatic and polarization interactions were only computed between the N_{2}O solute and SF_{6} solvent. Electrostatic and polarization contributions to the SF_{6} solvent–solvent interactions were neglected. Such a procedure ensures that the pure (liquid, gas) properties of the solvent are unaltered. All force field parameters are listed in Table S2 in the supplementary material.

### Analysis

*I*(

*ω*) of the IR spectra for N

_{2}O in different solvent densities is obtained via the Fourier transform of the dipole–dipole correlation function from the dipole sequence of the single N

_{2}O molecule,

*Q*(

*ω*) = tanh(

*βℏω*/2) was applied to the results of the Fourier transform.

^{45}This procedure yields line shapes but not absolute intensities. For direct comparison, individual spectra are, thus, multiplied with a suitable scaling factor to bring intensities of all spectra on comparable scales.

The response of the solute on the solvent structure and dynamics was evaluated by determining the frequencies of the quenched normal modes (QNMs) of N_{2}O for frames every 5 fs of the simulation. For the QNM, a steepest descent geometry optimization for N_{2}O, within a frozen solvent conformation for either 100 steps or until a gradient root mean square of $\u226410\u22124$ kcal/mol/Å was reached, was carried out. The 9 × 9 mass-weighted Hessian matrix of the N_{2}O solute in the solvent with regard to the solute atom displacements in Cartesian coordinates was determined and diagonalized to obtain time series of all nine normal modes including the effect of solvent.^{46,47} Instantaneous normal mode (INM) analyses were also applied on the solute by diagonalizing its mass-weighted Hessian matrix without prior geometry optimization on N_{2}O to obtain insights on the impact of the solvent structure on the translational and rotational modes of the N_{2}O solute.^{46}

The normalized vibrational frequency–frequency correlation function (FFCF) $\nu as(t)\u22c5\nu as(0)$/$\nu as(0)\u22c5\nu as(0)$ was computed for the time series of the asymmetric stretch frequency *ν*_{as} gathered from the QNM analyses. The amplitude *A*, lifetimes *τ*_{i}, and offset Δ of a bi-exponential function $c(t)=Ae\u2212t/\tau 1+(1\u2212A)e\u2212t/\tau 2+\Delta $ were optimized to fit the normalized FFCF for *t* ∈ [0.1, 2.0] ps.

*g*(

*r*) for solvent–solvent and solvent–solute pairs were determined from the average number of molecules ⟨d

*n*

_{r}⟩ of type B within the shell in the range of

*r*− Δ

*r*/2 <

*r*≤

*r*+ Δ

*r*/2 around molecules of type A,

*ρ*

_{local}⟩ is the local density of compound B around compound A within a range that is half the simulation box edge length, and the shell width was Δ

*r*= 0.1 Å. The average coordination number ⟨

*N*(

*r*′)⟩ of compounds B within the range

*r*′ around compound A can be obtained from

*g*(

*r*) according to

## RESULTS

### Validation of the interaction potentials

First, the quality of the intramolecular PES for N_{2}O is discussed, followed by a description of the van der Waals parameters for the N_{2}O solute fit to the *ab initio* reference calculations.

The RKHS model provides a full-dimensional, intramolecular PES for N_{2}O, which was originally developed for investigating the N + NO collision reaction dynamics.^{29} The Pearson coefficient *R*^{2} of the RKHS representation and the full set of reference values is 0.999 83, and the root mean squared error (RMSE) between the RKHS and reference energies up to 20 kcal/mol above the equilibrium structure (78 reference energies) is 0.13 kcal/mol. To establish the spectroscopic quality of the PES, the evaluation by the discrete variable representation (DVR) method using the DVR3D^{48} package yields a fundamental asymmetric stretch frequency of N_{2}O of 2229 cm^{−1} compared with 2224 cm^{−1} from the experiments in the gas phase.^{49–51} The bending and symmetric stretch frequencies obtained from the DVR3D calculations are 598 and 1291 cm^{−1}, respectively, and the overtone of the bending mode lies at 1184 cm^{−1}. Experimentally, the bending and symmetric stretches are found at 589 and 1285 cm^{−1}, respectively, while the overtone of the bending frequency is at 1168 cm^{−1}.^{49–51}

Using MDCMs for the electrostatics of N_{2}O in the two solvent environments (Xe and SF_{6}) requires a readjustment of the solute vdW parameters. To keep the solvent–solvent interaction unchanged, only the vdW parameters of N_{2}O were optimized with respect to interaction energies from reference electronic structure calculations.

Figure 1 reports the correlation between the reference and model interaction energies with optimized vdW parameters for N_{2}O for the N_{2}O–Xe [panel (A)] and N_{2}O–SF_{6} [panel (B)] complexes. The reference interaction energies Δ*E*_{CCSD} = *E*_{pair} − *E*_{solute} − *E*_{solvent} are the differences between the respective total *ab initio* energy of the pair and the energies of the solute *E*_{solute} and solvent *E*_{solvent} fragment in a fixed minimum energy conformation. The force field energies Δ*E*_{FF} = *E*_{elec} + *E*_{LJ} include electrostatic and polarization contributions and the Lennard-Jones potential between the solute and solvent at the same conformation as in the reference dataset. The energies sample attractive (negative energy) and repulsive (positive energy) ranges of the N_{2}O–solvent energies, and the comparison is restricted to energies within 5 kcal/mol of the separated molecules, which is the zero of energy. For N_{2}O–Xe, the RMSE is 0.25 kcal/mol, and for N_{2}O–SF_{6}, it is 0.29 kcal/mol. The most stable structure is stabilized by −1.41 kcal/mol (−1.37 kcal/mol from the fitted energy function) for N_{2}O–SF_{6} and −0.80 kcal/mol (−0.75 kcal/mol) for N_{2}O–Xe.

The critical concentration at which transition to a supercritical fluid occurs is another relevant property for the present work. Earlier work showed that this transition can be correlated with a pronounced increase in the local solvent reorganization lifetime *τ*_{ρ}.^{52,53} The quantity *τ*_{ρ} is the integral over the local-density autocorrelation function $\tau \rho =\u222b0\u221eC\rho (t)dt$ and characterizes the time required for the local environment around a reference particle, e.g., the solute N_{2}O in the present case, to change substantially. Hence, *τ*_{ρ} can also be considered a local-density reorganization time.^{53}

Figures S1 and S2 report the local solvent reorganization lifetimes *τ*_{ρ} for pure xenon and SF_{6}, respectively, from MD simulations. The solvent environments are defined by cutoff radii, which correspond approximately to the first and second minima of the solvent–solvent RDF; see Fig. S4. For SF_{6}, the force field was parameterized^{31} to reproduce the experimentally measured critical density with corresponding concentration *c*_{crit}(SF_{6}) = 5.06 M at the critical temperature.^{33,34} The present simulations yield a peak for *τ*_{ρ} at *c*(SF_{6}) = 5.02M, in close agreement with the experimental critical concentration at the critical temperature and the concentration with the local peak in *τ*_{2} in Fig. 4(D). For xenon, however, the parameterization of the PES^{32} did not include phase transition properties to the supercritical regime. Figure S1 shows that the maximum in *τ*_{ρ} occurs at *c*(Xe) = 5.19 M. This compares with a critical concentration of 8.45 M at *T*_{c} from the experiment.^{33,34}

### Infrared spectroscopy

The computed IR spectra for the N_{2}O asymmetric (*ν*_{as}) stretch in xenon and SF_{6} solution as a function of solvent density for near-critical isotherms allow for comparison between the experimental and simulation results. The present work focuses mainly on the change of the IR line shape at different solvent concentrations especially the P-, Q-, and R-branch structures; see Figs. 2 and 3. The P- and R-branches are the IR spectral features at lower and higher wavenumber from the pure vibrational transition frequency for the excitation of mode *ν*. The *ν*_{as}, *ν*_{s}, and *ν*_{δ} band shapes arise due to conservation of angular momentum during a vibrational excitation upon photon absorption.^{49} The Q-branch is the absorption band feature at the vibrational transition frequency. In addition to the asymmetric stretch *ν*_{as} ($\u223c2220$ cm^{−1}, symmetry A_{1}/Σ^{+}), the symmetric stretch *ν*_{s} ($\u223c1305$ cm^{−1}, A_{1}/Σ^{+}) and the bending vibration *ν*_{δ} ($\u223c615$ cm^{−1}, E_{1}/Π) fundamentals from the MD simulations are reported and analyzed as a function of solvent density.

From a quantum mechanical perspective, the rotational structure of an IR-active vibrational mode for a linear molecule (such as N_{2}O) in the gas phase leads to P- and R-branches. Selection rules dictate the change of the rotational and vibrational quantum states, *j* and *ν*, satisfying Δ*ν* = ±1 and Δ*j* = ±1 for vibrational modes of A_{1}/Σ^{+} symmetry and Δ*j* = 0, ±1 for vibrational modes of E_{1}/Π symmetry, respectively. For N_{2}O, the asymmetric *ν*_{as} and symmetric stretch *ν*_{s} are vibrational modes with a parallel vibrational transition dipole moment to the bond axis and of A_{1}/Σ^{+} representation and the bending mode *ν*_{δ} with a perpendicular vibrational transition dipole moment of E_{1}/Π representation. These rotational selection rules are only strictly valid in the absence of perturbations (exact free rotor) and break down with increasing deviation from a free rotor model, which, for example, can be due to embedding the rotor into a solvent of different density, which is a means to tune the strength of the perturbation. As a consequence of the perturbation, a Q-branch-like spectral feature emerges for parallel bands and becomes dominant with a band maximum at the wavenumber of the corresponding pure vibrational transition energy at sufficiently high solvent density.

For N_{2}O in xenon, the IR band structure of *ν*_{as} in Fig. 2(C) ranges from resolved P- and R-branches in gaseous xenon up to a single near-Lorentzian-shaped band in liquid xenon in these classical MD simulations. From the simulations, at xenon concentrations higher than 5.19 M or relative density of *ρ** = *ρ*/*ρ*_{c} = 0.62, the xenon solvent is a supercritical fluid. In the simulated absorption spectra at this density and above, the *ν*_{as} P- and R-band structure is not resolved, and a centrally peaked band at the pure vibrational frequency dominates the band shape. The black dashed lines in Fig. 2(C) are the experimentally determined spectra at the given solvent concentration.^{22–24}

The agreement between experimental and computed N_{2}O in Xe absorption band shapes is good and captures the dramatic gas phase-to-condensed phase line shape change as a function of solvent density. However, the frequency position of the *ν*_{as} band needs to be blue-shifted by a small frequency shift *ω* to achieve the best overlap with the experimental IR line shape. The density-dependent shifts are indicated in the panel and range from +6 to +22 cm^{−1} with no discernible trend. The shift originates from different effects, including insufficient sampling of the amount and distribution of internal energy within the N_{2}O solutes vibrational degrees of freedom, remaining small inaccuracies in the intermolecular interactions, neglecting many-body contributions and slightly underestimating the anharmonicity in the PES along the relevant coordinates.

Figure 2(B) shows the corresponding band shape for *ν*_{s} of N_{2}O in Xe around 1305 cm^{−1} in comparison to experimentally measured 1291 cm^{−1} for N_{2}O in the gas phase. Similar to *ν*_{as}, it also displays resolved P- and R-branches at low solvent concentrations and changes into a Q-branch-like dominated structure in supercritical and liquid xenon. The IR band of the N_{2}O bending mode *ν*_{δ} appears at around 615 cm^{−1} shown in Fig. 2(A), compared with 589 cm^{−1} from the experiment in the gas phase.^{49} At low solvent concentrations, the band structure of *ν*_{δ} from the MD simulations is a sharp Q-branch with weaker P- and R-branch side bands consistent with the quantum mechanical selection rules. The intensity of these bands is no longer evident relative to the central bending feature at higher solvent concentration. The first overtone of the *ν*_{δ} vibration is also detected between 1220 and 1230 cm^{−1} in Fig. 2(B) but with low intensity. The computed bending overtone (2*ν*_{b}) exhibits the same P- and R-branches at low density and the change to a single peaked band at high solvent concentration as observed for the *ν*_{s} and *ν*_{as} band structures, see Figs. S5 and S6, in agreement with the experimental results.^{54}

The IR band structure of *ν*_{as} of N_{2}O in SF_{6} is shown in Fig. 3(C). At the lowest solvent concentrations 0.20 and 0.51 M (*ρ** = {0.04, 0.10}), the simulated distinct P- and R-branch structure is evident but gradually disappears for [SF_{6}] > 1.52 M (*ρ** > 0.30) for corresponding experimental state points. As the solvent concentration of SF_{6} increases, the band shape changes into a single-peaked band in agreement with the experimental IR spectra for [SF_{6}] = 7.65 and 9.47 M (*ρ** = {1.51, 1.87}). The computed IR spectra for 0.81 ≤ [SF_{6}] ≤ 5.93 M (0.16 ≤ *ρ** ≤ 1.17) do not exhibit the double-peak structure or are narrower than the experimental absorption spectra. The frequency position of the *ν*_{as} band is also adjusted by a small frequency blue shift *ω* within the range of +8 to +27 cm^{−1} to maximize the overlap with the experimental IR line shape.

The line shapes for the *ν*_{δ} and *ν*_{s} vibrational modes are shown in Figs. 3(A) and 3(B) as a function of SF_{6} concentrations, respectively. The *ν*_{s} mode changes from P-/R-branches at low solvent concentrations into a single featured structure in supercritical and liquid SF_{6}. The calculated IR band of *ν*_{δ} around 615 cm^{−1} exhibits a sharp Q-branch with resolvable weak P- and R-branch satellites only at low SF_{6} concentration. The first overtone of the perpendicularly polarized *ν*_{δ} mode is again detected between 1220 and 1230 cm^{−1} with low intensity and a band shape as observed for the parallel polarized fundamental *ν*_{s} and *ν*_{as} band structures.

### Vibrational frequency–frequency correlation function

The vibrational frequency fluctuation correlation function (FFCF), which probes the coupling of the solute vibrational modes to the solvent environment, can be determined from 2DIR experiments. Of particular interest are the time scales of the lifetimes *τ*_{i} and amplitudes *A*_{i} with which the FFCF decays, which are shown in Fig. S7. In the condensed phase, the FFCF is also related to the change of the center-line slope at different 2DIR waiting times.^{3}

Figures 4(A) and 4(B) report the FFCFs for the asymmetric stretch *ν*_{as} mode of N_{2}O in xenon and SF_{6} from QNM, respectively.^{1,55} The fit of the FFCFs using a bi-exponential decay yields two times: a rapid inertial component (*τ*_{1} ∼ 0.05 ps) and a longer spectral diffusion time scale (*τ*_{2} ∼ 0.2 … 0.5 ps). It is not possible to directly compare the computed vibrational FFCFs from INM with the experimentally measured 2DIR spectra as they are overwhelmingly dominated by the contribution of rotational energy relaxation dynamics.^{22–24} Nevertheless, these time scales can be compared with measured and computed time scales [*τ*_{1}, *τ*_{2}, *τ*_{3}] for CN^{−} ([0.2, 2.9, n.a.] and [0.04, 0.87, 9.2] ps) and for $N3\u2212$ ([n.a., n.a., ∼1.2] and [0.04, 0.23, 1.2] ps) in H_{2}O.^{55–58} Here, “n.a.” refers to the time scales that were not determined or available from the experiments but clearly appeared in the analysis of the simulations. A similar behavior of the vibrational FFCF as the one found in the present work—a first rapid relaxation time on the 0.05–0.1 ps time scale, followed by a slower time scale of 0.5–1.0 ps—was reported in an earlier MD simulation study for liquid HCl.^{59}

It is important to note that experimentally decay times on the several 10 fs time scale are difficult to determine with confidence from vibrational FFCFs. Therefore, only those longer than that are considered in the following. It is noted that for both ions (CN^{−} and $N3\u2212$) in water, mentioned above, the longer decay times are 1 ps or longer compared with $\u223c0.5$ ps for N_{2}O in SF_{6}. This is consistent with the fact that ion–water interactions are considerably stronger than N_{2}O–SF_{6} interactions. Furthermore, the N_{2}O–Xe interaction is weakest among all those considered here, which leads to the rather short *τ*_{2} value for this system, even relative to SF_{6}. Given that for the present system the intermolecular interactions are weak and the decay times are rapid, it is anticipated that only *τ*_{2} for N_{2}O in SF_{6} would be amenable to the 2DIR experiments.

Figures 4(B) and 4(D) show the dependence of the short (*τ*_{1}) and long (*τ*_{2}) time scales as a function of the xenon and SF_{6} solvents concentrations, respectively. In both solvents, the longer lifetime slows down with increasing concentration, and a peak in the lifetime at 5.19 M in Xe and 5.02 M in SF_{6} is observed. The pronounced increase in *τ*_{2} is likely to be related to approaching the critical point as the local peak in *τ*_{2} at a SF_{6} solvent concentration of 5.02 M matches the critical concentration for pure SF_{6}. Similarly, for N_{2}O in Xe, the peak at 5.19 M is consistent with a lengthening of the *τ*_{2} timescale in this solution at the critical concentration for pure xenon which was obtained from the local solvent reorganization lifetimes *τ*_{ρ}; see Fig. S1. Thus, both solvents exhibit critical slowing (lengthening of *τ*_{2}) at the critical concentration of the respective pure solvent.

It is also noted that the slope for *τ*_{2}(*c*) for N_{2}O in SF_{6} is considerably steeper than for Xe as the solvent. This is most likely also due to the increased solvent–solute interaction strength between N_{2}O and SF_{6} compared with Xe. The results for *τ*_{2} indicate that interesting dynamical effects, such as critical slowing, can be expected to develop around the critical point of the solvent. It is worthwhile to note that at the critical concentration *c*(Xe) = 5.19 M for the transition to a SCF, the spectroscopy of the asymmetric stretch changes from P-/R-branches to Q-branch-like (Fig. 2), and the decay time *τ*_{2} from the FFCF of the asymmetric stretch [see Fig. 4(B)] features a pronounced maximum. The fact that for xenon the critical density from the simulations is underestimated (5.19 vs 8.45 M from experiments), but the concentration-dependence of the spectroscopy agrees with the experiment suggests that orientational relaxation and inhomogeneous effects play the dominant role for the density dependence of the dipole correlation function.

A direct comparison of the vibrational FFCFs from INM with the prior experimental FFCFs is not meaningful at present because of the overwhelming contribution to the reported FFCF due to the rotational energy relaxation dynamics. The 2DIR spectral features corresponding to the vibrational FFCF overlap the rotational energy relaxation but may be accessible in the future experimental studies that employ higher spectral resolution detection than in Refs. 22–24. Furthermore, unlike typical 2DIR spectra of condensed phase systems, the available experimental FFCFs report only on the quasi-free rotor members of the ensemble contrary to the computed vibrational FFCFs.

### Radial distribution functions and size of the solvent shells

A structural characterization of the solvent environment is afforded by considering the radial distribution functions *g*(*r*) between the center of mass of N_{2}O and the central sulfur atom and the Xe atom, which are shown in Figs. S3A and S3C, respectively. They show a wide first peak corresponding to the first solvation shell around the N_{2}O solute within $\u223c4.3$ Å [panel (A)] for the xenon and $\u223c4.6$ Å [panel (C)] for the SF_{6} solvent, followed by considerably weaker features around 8 Å and beyond depending on the density.

A more concise comparison of how the solvent environment depends on solvent density can be obtained by considering the average coordination number ⟨*N*(*r*′)⟩ of the Xe and SF_{6} solvents within the first solvation shell [Figs. 5(A) and 5(B)] that are computed from the respective *g*(*r*), see Fig. S3. At the respective lowest sampled solvent concentration, ⟨*N*(*r*′)⟩ yields an average of 0.31 Xe atoms [panel (A) at 0.34 M] and 0.29 SF_{6} molecules [panel (B) at 0.20 M] within the first solvation shell. The coordination number rises monotonically with increasing solvent concentration but shows a decrease in the slope within the SCF regime as the solvation shell fills. It is interesting to note that both solvents exhibit two slopes for *N*(*r*) depending on *c*: a first, steeper one for the gas-phase-like solvent into the SCF regime, and a second, flatter one, leading into the liquid-like regime. The slopes for the gas-phase-like regime is steeper (1.40 M^{−1}) for SF_{6} than for Xe (0.94 M^{−1}), which may be related to the solvent–solute interaction strength.

It is also of potential interest to consider the present results in light of the independent binary collision (IBC) model^{60–62} although the model has also been criticized to be potentially oversimplified when applied to liquids.^{63,64} In the simplest implementation of the IBC model—where the contact distance is associated with the position of the first maximum of the solvent–solute radial distribution function, which only insignificantly changes with density—the rate for a given process, e.g., vibrational relaxation, depends linearly on (a) the collision rate and (b) a density independent probability that a given collision is effective; see Figs. S3A and S3C. As a consequence, rates adequately described by an IBC model depend linearly on density, which itself is proportional to the occupancy of the first solvation shell. This is what is found in the present work for low densities; see Fig. 5. However, toward higher densities, another linear dependence develops with lower slope for both solvents. The breakdown of the IBC model occurs at a slightly higher solvent concentration for Xe than for SF_{6} (*c*(Xe) ∼ 3.5 M vs *c*(SF_{6}) ∼ 2.9 M, see Fig. 5) compared with values of $\u223c3.4$ M and 4.0 M from 2DIR experimentally determined rotational energy relaxation, respectively.^{24} Just as found for comparison of the simulated and experimental absorption spectra, further refinement of the N_{2}O–SF_{6} interactions will provide even better quantitative agreement with the experiment, for example, within a PES morphing approach.^{65}

## DISCUSSION

Classical MD simulations have been previously used to determine and analyze rotational vibrational spectra of simple solutes, including CO or HCl in a solvent (argon)^{66,67} or for water.^{68} These calculations demonstrated that classical MD simulations are capable of realistically capturing P- and R-branch structures in moderate to high density gases and even in liquid water.^{68} Along similar lines, the present MD simulations of N_{2}O in xenon or SF_{6} as the solvent reproduce well the experimentally observed splitting in P-, Q-, and R-branches and the overall line shapes in the IR spectra of N_{2}O depending on the solvent concentration (see Figs. 2 and 3). The agreement in the computed line shapes of the FTIR *ν*_{as} spectrum of N_{2}O in xenon does match with the experimental spectra as a function of solvent concentration. The separation between the maxima of the P- and R-branches is ∼[27, 25, 29, 22] cm^{−1} at [0.34, 0.84, 1.26, 3.10] M from the simulations, compared with ∼[25, 25, 22] cm^{−1} at [0.84, 1.26, 3.10] M from the experiments. The IR signals for the *ν*_{as} vibration of N_{2}O in SF_{6} show a larger contribution of the Q-branch-like feature at lower *ρ** than that observed in the experimental spectra. For simulations in both solvents, the computed spectra are shifted to the red compared with the results from experiments for *ν*_{as}.

The computed IR spectra from the classical MD simulations of N_{2}O/Xe in Fig. 2 agree with the relative intensity and band width of the *ν*_{as} vibration signal at the corresponding solvent concentration. The IR spectra of N_{2}O in SF_{6} in Fig. 3 deviate more from the corresponding experimental line shapes by overestimating solvent–solute interactions, leading to a somewhat stronger contribution of the Q-branch-like absorption character at lower solvent concentration. The experimental *ν*_{as} signal at the lowest measured solvent density of 0.81 M overlaps better with the computed signal at lower solvent densities of 0.51 and 0.20 M. However, the line shape in highly dense supercritical and liquid SF_{6} agrees well with the experiment where the Q-band-like feature dominates the spectral line shape and resolvable P- and R-branch contributions are not evident. At the same absolute density, the asymmetric stretch of N_{2}O is indicative of a more liquid-like character of SF_{6} as the solvent compared with xenon; see Figs. 2 and 3(C).

Splitting of the IR spectrum into P- and R-branches correlates with free rotor properties of the N_{2}O solute molecule. Based on this, the results in Fig. 3 indicate an overestimation of the solute–solvent interaction in SF_{6} because a merged line shape for *ν*_{as} arises between 0.81 and 1.52 M from the computations, whereas experimentally, this occurs for concentrations between 3.39 and 4.36 M. Consequently, the rotation of N_{2}O in the simulations starts to become hindered at lower solvent densities compared to the experiment. This observation correlates with the long-range and stronger electrostatic interaction between N_{2}O and polar SF_{6} rather than xenon as the solvent.

As a comparison, the electrostatic contribution to the interaction energy of a N_{2}O/SF_{6} pair in its equilibrium conformation with −1.42 kcal/mol is considerably stronger then the polarization contribution in an N_{2}O/Xe pair with −0.84 kcal/mol. As the electrostatic MDCM and vdW parameter fit are only applied to single molecules or molecule pairs, respectively, the energy function might insufficiently capture many-body interactions between the N_{2}O solute and multiple SF_{6} solvent molecules.

An overtone of the *ν*_{δ} vibration (615 cm^{−1}) of N_{2}O is found at about twice the frequency between 1220 and 1230 cm^{−1} in both solvents. The splitting into P- and R-branches at low solvent densities originates from excitation of simultaneously two bending mode quanta (2*ν*_{δ}) and the change of the rotational state by Δ*J* = ±1.^{69,70} Consequently, as the splitting of the 2*ν*_{δ} IR band depends on the free rotor properties such as for the *ν*_{as} and *ν*_{s} IR band structure, the band shape changes toward a single Q-branch peak at higher solvent concentrations. In other words, the rotational “selection rules” followed by the bending overtone are the same as for the asymmetric stretch vibration, both of which have the same symmetry. This quantum mechanical result is seen in the present classical MD simulations.

The perturbation of the rotational modes by the solvent is visualized by the INM in Figs. S8A and S8B for N_{2}O in SF_{6} and xenon at different densities, respectively. The INM histogram for N_{2}O in SF_{6} shows two peaks around the zero frequency line for low positive and imaginary frequencies (positive and negative second derivative of the potential projected along normal mode displacement, respectively). The emerging shoulder(s) in the INM histogram toward higher solvent densities (Fig. S8) are due to low-frequency modes involving solute–solvent interactions. This leads to larger deviations from a free rotor movement of the N_{2}O, akin to motion in a density-dependent rotationally constraining solvent potential.^{23} Visual inspection of the normal modes identifies these as librational modes or frustrated rotations. The INM histogram for N_{2}O in xenon in Fig. S8B shows broader peaks and lower shoulders for the gaseous solvent in comparison with the INM distribution of N_{2}O in SF_{6} in Fig. S8A. The higher mode density around low frequencies for the INMs of N_{2}O in xenon relative to SF_{6} also indicates lower perturbation of the free rotor character in gaseous and supercritical xenon solvent than SF_{6}.

The change of the IR spectra for linear molecules in gaseous, supercritical up to liquid solvents can be compared with the experimental and theoretical work on hydrogen cyanide in an electric field.^{71,72} With increasing field strength, the molecules evolve from a free rotor to those trapped in pendular states, which is one way to control the population of rotational states in molecular beam experiments. The electric field lifts the degeneracy of the quantum number *m*_{j} (the projection of the angular momentum on the electric field direction) for each rotational state *j* and leads to a manifold of different transitions that correspond to the IR spectra. Pendular states are linear combinations of field-free (*j*, *m*_{j}) states covering a range of *j*–values but sharing the same *m*_{j} quantum number, i.e., *j*–mixing occurs.^{71} With increasing electric field strength, the field-free selection rule Δ*j* = ±1 breaks down. For excitation along dipole transition vectors parallel to the electric field, the P- and R-branches change into a single Lorentzian-shaped Q-branch at large electric field strength.^{71}

The effects found for polar molecules in external electrical fields can be compared with the motion of N_{2}O in a supercritical or liquid solvent. Here, the solvent molecules form a cavity around the solute where the direction of the angular momentum is not energetically degenerate, and N_{2}O behaves like in a pendular state. In analogy to a molecule in an electric field, the P- and R-branches collapse into a single Q-branch-like band structure for vibrational modes with transition dipole moments perpendicular to the angular momentum of the molecule upon breakdown of the usual selection rules.

*ϕ*=

*ϕ*

_{vib}·

*ϕ*

_{rot}a product ansatz of a vibrational

*ϕ*

_{vib}and a rotational component

*ϕ*

_{rot}.

^{68}The composition of the IR band by P/R-branch and Q-branch depends on the angle

*θ*between the vibrational transition vector and the rotational axis. Assigning the vibration and rotation with an angular frequency

*ω*

_{vib}and

*ω*

_{rot}, the time-dependence of the correlation function

*ϕ*(

*t*) for an idealized classical rotor (without damping) has been expressed as

^{68}The Fourier transform of

*ϕ*(

*t*) in Eq. (6) leads to one Q-branch signal at

*ω*

_{vib}for parallel alignment of the rotational axis and vibrational transition dipole vector (cos

*θ*= 1) and two P- and R-branch signals at

*ω*

_{vib}−

*ω*

_{rot}and

*ω*

_{vib}+

*ω*

_{rot}for perpendicularly aligned rotational axis and vibrational transition dipole vector (cos

*θ*= 0). For other cases or an unstable rotation axis, signals for all P-, Q-, and R-branches arise in the IR spectra.

Figure 6 shows computed IR spectra from the Fourier transform of the dipole–dipole correlation function from the MD simulations of a single N_{2}O molecule with different initial conditions. For all cases (A)–(D), intramolecular energy is distributed along all vibrational normal modes with respect to the thermal energy at 321.9 K. Without any rotational energy assigned, the IR spectra in Fig. 6(A) show no split into P- and R-branches. With rotational energy assigned around the rotational axis *parallel* to the transition dipole vector, the IR spectra in Fig. 6(B) show separation into P- and R-branches by 28 cm^{−1} for the *ν*_{as} stretch that correlates for a rotational constant of *B*_{j=16} = 0.43 cm^{−1}. The IR signal of the bending mode primarily consists of the Q-branch with two small P- and R-branch satellites separated by ± 14 cm^{−1} each from the Q-branch. For rotational energy assigned around a rotational axis *perpendicular* to both stretch and bending transition dipole moments, all IR signals in Fig. 6(C) split in P- and R-branches as well (*B*_{j=16} = 0.42 cm^{−1}), but no Q-branch is observed. The rotational energy corresponding to the rotational quantum state *j* = 16 is the state with the highest probability according to the Maxwell–Boltzmann distribution function *P*_{j}(*T*) for N_{2}O at 321.9 K.

At lower solvent density, the appearance of P-, Q-, and R-branches in the model IR spectra in Fig. 6(B) can be explained by the mechanical instability of a rotation around the transition dipole vector of the bending mode *e*_{y}. For bent N_{2}O, the moment of inertia are ordered by *I*_{x} > *I*_{y} ≫ *I*_{z}. According to the tennis racket theorem,^{73} the rotation around the axis *e*_{y} with the middle moment of inertia is unstable. Thus, the angle between rotation axis or angular momentum vector and the transition dipole vector of the bending mode varies in time and gives rise to P-, Q-, and R-branch signals.

A superposition of model IR spectra at varying rotational energy states and randomly sampled rotation axis but perpendicular to the N_{2}O bond axis yields a line shape close to the computed spectra of N_{2}O in the gaseous solvent at lower concentrations. The contribution of the single IR spectra at rotational state *j* to the model spectra in Fig. 6(D) iss scaled by the probability value according to *P*_{j}(*T* = 321.9 K). An average rotational constant *B* = 0.42 cm^{−1} is computed from the spectra with rotational states *j* = [2, 26], which is in perfect agreement with the experimentally determined rotational constant of *B*_{exp} = 0.419 cm^{−1}.^{74} Hence, overall two regimes can be distinguished: for low densities, the rotational motion of the solute is oscillatory and guided by the mechanical instability around the middle rotational axis (tennis racket theorem), whereas at higher densities, the constraining effect of the solvent due to tighter packing leads to a more diffusive rotational motion of the solute.

For additional clarification for how the transition between P-/R-features and Q-branch-like line shapes occurs, two different models are briefly considered. The first one operates directly on the dipole–dipole correlation function of a single N_{2}O molecule, whereas the second one uses a Langevin-type simulation of one N_{2}O molecule in the gas phase. The influence of a solvent can be heuristically included by multiplication of the dipole–dipole correlation function *ϕ*(*t*), see Eq. (6), with an exponential damping, i.e., considering *ϕ*(*t*) × exp(−*t*/*τ*). The damping with characteristic time *τ* models all sources of the dipole correlation function decay, primarily the solvent density dependence on orientational relaxation due to collisions between N_{2}O and the solvent in the present case. With increasing density of the solvent, for example, *τ* decreases, which leads to loss of knowledge about its rotational motion and broadening of the computed IR signals; see Fig. 7(A). Consequently, the P-/R-features collapse into a single Q-branch-like peak for relaxation decays 0.1 ≤ *τ* ≤ 0.2 ps, which are considerably shorter than the rotational periods.

Alternatively, interaction with a solvent can be modeled from the MD simulations with a stochastic thermostat, such as in Langevin dynamics, which applies random forces on the atoms to simulate orientational interactions with the heat bath. For this, 10 ns Langevin simulations with varying coupling strengths *γ* were run for a single N_{2}O molecule in the gas phase using CHARMM, and the dipole–dipole correlation function was determined. Figure 7(B) shows the computed IR spectra. At low coupling strength *γ*_{i} = 10^{−4} ps^{−1}, the IR spectra still yield a split of the asymmetric stretch into P- and R-branches (yellow trace). With increasing friction *γ*_{i} = 0.01 ps^{−1}, still multiple peaks are visible in the asymmetric stretch absorption line shape, but the spectral Q-branch-like region is getting filled in (pink trace). At yet higher *γ*_{i} ≥ 0.1 ps^{−1} (green, red, and blue traces), the P-/R-features wash out entirely.

Finally, it is of some interest to discuss the results on the spectroscopy (Figs. 2 and 3) and the solvent structure (Fig. 5) from a broader perspective. The interaction between N_{2}O and Xe is weaker by about 50% than that between N_{2}O and SF_{6}. The type of interaction also differs: pure van der Waals (in the present treatment) for Xe vs van der Waals and electrostatics for SF_{6}. It is interesting to note that the transition between the P-/R- and the Q-band-like spectra in Xe occurs at *ρ** ∼ 0.6 ($\u223c5.2$ M) compared with *ρ** ∼ 0.2 ($\u223c0.9$ M) for SF_{6}; see Figs. 2 and 3. This is consistent with the increased interaction strength for N_{2}O–SF_{6} as fewer solvent molecules are required to sufficiently perturb the system to illicit the change in the spectroscopy. Similarly, the slope of *N*(*r*) vs concentration can also serve as a qualitative measure for the interaction strength. To stabilize a solvent shell of a certain size, a smaller number (lower concentration) of solvent molecules is required if the solvent–solute interaction is sufficiently strong. These qualitative relationships are also reflected in the longer *τ*_{2} values of the FFCFs for N_{2}O in SF_{6} compared with Xe as the solvent; see Fig. 4. For large polyatomic molecules in highly compressible fluids, “attractive” solutes have been found to recruit increased numbers of solvent molecules, which leads to local density augmentation.^{75} While the amount of augmentation was found to correlate well with the free energy of solvation of the solute from simulations, computations and experiments were reported to follow different correlations.

## CONCLUSION

In conclusion, the present work establishes that the rotational structure of the asymmetric stretch band is quantitatively described for N_{2}O in xenon and is qualitatively captured in SF_{6} compared with the experiment. Specifically, the change from P-/R-branches in the gas phase environment changes to a Q-branch-like structure, eventually assuming a Lorentzian-like line shape once the solvent passes through the supercritical fluid to the liquid state. Transition between the gas- and SCF-like solvent is reflected in a steep increase in the longer correlation time *τ*_{2} in the FFCF for both solvents. As the density of the solvent increases, additional solvent shells appear in *g*(*r*). The present work suggests that atomistic simulations together with machine learned and accurate electrostatic interactions yield a quantitative understanding of the spectroscopy and structural dynamics^{4} from very low density to the liquid environment with the possibility to characterize transition to the SCF.

## SUPPLEMENTARY MATERIAL

See the supplementary material for additional data and figures referred to in the article that further support the findings of this work.

## ACKNOWLEDGMENTS

This work was supported by the Swiss National Science Foundation, Grant Nos. 200021-117810 and 200020-188724, the NCCR MUST, and the University of Basel (to M.M.) and by European Union’s Horizon 2020 research and innovation program under the Marie Skłodowska-Curie Grant Agreement No. 801459—FP-RESOMUS, which is gratefully acknowledged. The support of the National Science Foundation, Grant No. CHE-2102427 (L.D.Z.), and the Boston University Photonics Center is gratefully acknowledged. We thank James L. Skinner for valuable correspondence on Ref. 68.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Kai Töpfer**: Conceptualization (equal); Data curation (equal); Formal analysis (lead); Investigation (lead); Methodology (equal); Software (lead); Validation (equal); Visualization (lead); Writing – original draft (lead); Writing – review & editing (equal). **Debasish Koner**: Investigation (supporting); Methodology (supporting); Software (supporting); Validation (supporting); Writing – review & editing (equal). **Shyamsunder Erramilli**: Conceptualization (equal); Methodology (equal); Writing – review & editing (equal). **Lawrence D. Ziegler**: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Supervision (equal); Writing – review & editing (equal). **Markus Meuwly**: Conceptualization (lead); Formal analysis (lead); Investigation (equal); Methodology (equal); Project administration (lead); Resources (lead); Supervision (lead); Writing – original draft (lead); Writing – review & editing (equal).

## DATA AVAILABILITY

Raw data were generated at the University of Basel large scale facility. The derived data supporting the findings of this study are available from the corresponding author upon reasonable request. The data and source code that allow for reproducing the findings of this study are openly available at https://github.com/MMunibas/N2O_SF6_Xe. The original release of the N2O model potential is available at https://github.com/MMunibas/N2O-PESs.

## REFERENCES

*Concepts and Methods of 2D Infrared Spectroscopy*

^{+}, Na

^{+}, K

^{+}, Cs

^{+}) thiocyanate (SCN

^{−}) aqueous solutions

*J*scrambling and perfectly anticorrelated cross peaks

*J*-scrambling, vibrational relaxation, and the onset of liquid character

_{6}solutions

_{2}O and dynamics for the N + NO ↔ O + N

_{2}and N

_{2}+ O → 2N + O reactions

*Handbook of the Thermodynamics of Organic Compounds; Section on Vapor-Liquid Critical Constants of Fluids*

*CRC Handbook of Chemistry and Physics; CRC Handbook of Chemistry and Physics*

*Infrared and Raman Spectra of Polyatomic Molecules*

_{2}O) below 1.2

*μ*

_{2}O

^{−}in water

*ab initio*potentials: A systematic study of Ne–HF

_{2}O trapped in argon matrices: Comparison with O

_{3}and CO

_{2}

_{2}O)

_{2}and Ar–N

_{2}O complexes in the 2

*ν*

_{2}overtone region of N

_{2}O

_{2}O