An ultrafast vibrational study of dynamical heterogeneity in the protic ionic liquid ethyl-ammonium nitrate. I. Room temperature dynamics

Protic ionic liquids (PILs) are room-temperature molten salts of a Brønsted acid and Brønsted base.^1,2 Their high thermal stability, low vapor pressure, high proton conductance, and “water-like” hydrogen bonding make them promising as proton-conducting electrolytes for next-generation fuel cells.^3,4 The exact molecular mechanism of proton transport, however, is still unclear. Proton transport, in turn, is fundamentally controlled by the hydrogen bond network—its structure and fluctuations dictate if protons move independently of the proton donor (Grotthus transport) or with the proton donor (vehicular transport). The lack of a clear mechanistic picture for hydrogen bond dynamics in PILs limits the strategic development and chemical synthesis for applications in industrial processing and proton-conducting electrolytes.

The PIL ethyl-ammonium nitrate (EAN), the first room temperature ionic liquid to be discovered,⁵ has been extensively studied ever since it was synthesized in 1914. EAN is a “good” PIL (ΔpKa = pK_a(base) − pK_a(acid) = 11.93),^6,7 meaning the acid and base dissociate completely. In addition, EAN behaves like an ideal solution of independent ions² and conducts electricity by free diffusion of the ions.⁷ 

EAN forms an extended, three-dimensional hydrogen bonding network comparable to that of H₂O(l), suggesting it may have a similar proton transfer mechanism to H₂O(l). EAN consists of three hydrogen bond donors on the cation and three hydrogen bond acceptors on the anion; the large hydrogen bond angles between the cation and anions result in a majority of bent hydrogen bonds.⁸ Far-IR spectroscopy shows that the hydrogen bond stretching and bending modes in EAN are comparable in energy to that of H₂O(l).⁹ Dielectric spectroscopy¹⁰ and combined dielectric relaxation and optical Kerr effect spectroscopies¹¹ provide a detailed picture of the temperature dependence of structural relaxation modes in the liquid. Classical molecular dynamics (MD) simulations find the average hydrogen bond lengths in EAN to be comparable to those of strong hydrogen bonds, and the timescales of molecular rotations of the ammonium cation (NH stretch) are faster in comparison to PILs with weakly interacting single hydrogen bonds.¹² Furthermore, ultrafast mid-IR spectroscopy suggests that the ammonium of EAN undergoes large angle jumps when switching hydrogen bond partners.¹³ 

EAN displays both structural and dynamical heterogeneity.^14,15 Ionic liquids (ILs), in general, have a propensity to form polar and non-polar domains with sizes on a nanometer length-scale. The polar domains are charge dense and dominated by electrostatic ordering of the cations and anions, while the non-polar domains are charge depleted and dominated by weaker van der Waals interactions between alkyl groups.^16,17 Neutron scattering experiments and modeling show that EAN has alternation of charge-enhanced and charge-depleted regions on a nanometer lengthscale.^14,15 In EAN and ethylammonium formate, ultrafast IR-pump–IR-probe experiments reveal different rotational times depending on the region under investigation; ionic sub-domains relax much more slowly than the hydrophobic sub-domains.¹⁸ In MD simulations, increasing the alkyl chain length of a primary ammonium PIL increases the dynamical heterogeneity.¹⁹ Solutes can effectively report the dynamical heterogeneity because the solute’s vibrational and rotational dynamics can differ in the soft or stiff regions.^20,21

The nitrile stretch of SCN⁻, the ν₃ mode, is sensitive to its local environment.²² In water, neutron scattering suggests that SCN⁻ accepts one hydrogen bond along with two additional “hydration bonds,” which are Coulombic in nature, within the first hydration shell.²³ Early MD simulations note few linear hydrogen bonds with water²⁴ and many arrangements of dipolar character, while more recent MD simulations of SeCN⁻²⁵ find an average of 3–4 hydrogen bonds, with one of them being an axial (linear) hydrogen bond. Ultrafast infrared spectroscopy can use the ν₃ mode of pseudo-halide anions to investigate ultrafast vibrational dynamics of water and other polar solvents,^22,26–31 concentrated ion solutions,^32–34 ILs,^35–39 supported IL membranes,⁴⁰ deep eutectic solvents,⁴¹ and colloid emulsions.⁴² Furthermore, a combination of MD simulations and ultrafast infrared spectroscopy has developed a spectroscopic map of SeCN⁻ in D₂O that describes the frequency dependence of the nitrile stretch on the hydrogen bonding environment.²⁵ 

Two-dimensional infrared spectroscopy (2D-IR) spectroscopy can report two-point correlations in vibrational frequency and orientation. The frequency fluctuation correlation function (FFCF) is revealed through a change in shape of the 2D-IR spectra. The intensities of the spectra encode the orientations of the chromophores with respect to the lab frame. Orientational relaxation can be measured through both the orientational anisotropy from pump–probe spectra, the standard technique, and also two-dimensional anisotropy from 2D-IR spectra.⁴³ When an IR chromophore’s vibrational frequency depends on its instantaneous orientation (like a Stark shift), then the frequency and orientational dynamics become coupled.^44,45 The relative angle of the molecule’s dipole with the local electric field determines the magnitude of the vibrational frequency shift. As a result, either rotation of the molecule with the respect to a fixed solvation cage or rotation of the solvent electric field relative to the molecule will produce fluctuations in the vibrational frequency. First described by Rivera et al. for vibrational sum frequency generation⁴⁶ and elaborated by Kramer et al. for 2D-IR spectroscopy,^44,45 the FFCF extracted from 2D spectra can be decomposed into reorientation-induced spectral diffusion (RISD) and structural spectral diffusion (SSD) when polarization information is available.

The dynamics of SCN⁻ in EAN at room temperature are presented in this paper in comparison to H₂O. The 2D-IR lineshapes of SCN⁻ in H₂O indicate only a single ensemble and are thus a negative control for sub-ensembles in this context. In Paper II of this series, we will compare the activation energies for the dynamics of SCN⁻ in EAN and H₂O in detail.⁴⁷ 

In this paper, we develop an assignment of the two sub-ensembles of SCN⁻ in EAN. The non-Gaussian structure of the IR absorption spectrum of SCN⁻ in EAN hints at the presence of two hydrogen bonding environments. The two sub-ensembles are indicated by the frequency-dependent rate of vibrational relaxation in the 2D-IR measurements. The SCN⁻ band shifts and narrows as a function of the population time, t₂, which strongly indicates two different vibrational relaxation rates. We also observe that the rates of decay of rotational anisotropy are frequency dependent, and the frequency fluctuations also vary somewhat across the spectra. Nevertheless, the time constants for rotational motion and frequency fluctuations are the same for the two sub-ensembles, within error. Treating the frequency and rotational dynamics as representing an average of the two ensembles, a wobbling-in-a-cone RISD model shows that rotational motion suffices to account for the observed spectral diffusion on a few picosecond timescale. The slower times in the rotational anisotropy decay (∼25 ps) are believed to represent the slower process of total orientational relaxation either through rotational diffusion or cone-jump processes. The temperature dependence and activation energies of the dynamics of SCN⁻ in EAN determined through this approach are reported in Paper II of this series.⁴⁷ 

Ammonium thiocyanate (NH₄SCN, purity >97.5%), ethylamine (67–70 wt. %), and 15.7 M HNO₃ were purchased from Sigma-Aldrich. EAN was synthesized by an acid–base titration. Before use, the PIL was dried using a vacuum pump (10 μTorr). Solutions were 100 mM of NH₄SCN added to dry EAN and H₂O, respectively. The water content in EAN was <40 mM after drying, confirmed by FTIR (supplementary material).

Samples were prepared on two CaF₂ windows with a 25 μm teflon spacer in between the two windows in a Harrick cell. A few μl of the given solution was sandwiched between the two windows with the desired spacer. The samples were prepared in a glove bag under an inert N₂ atmosphere.

The FTIR spectra were obtained using a Nicolet 6700 FTIR instrument purchased from ThermoFisher Scientific. The instrument was equipped with a N₂ purge to eliminate atmospheric bands that would interfere with IR vibrational bands under investigation. The resulting spectra were analyzed using MATLAB. The spectra were baselined to zero using the region 3950–4000 cm⁻¹.

The 2D-IR spectra were collected at the Central Laser Facility at the Rutherford Appleton Laboratory (RAL) in the United Kingdom. The SCN⁻ data were collected on the Lifetime system developed by Donaldson et al.,⁴⁸ which is a 100 kHz 2D-IR spectrometer based on Yb:KGW amplifier technology. Briefly, pump and probe pulses were generated from optical parametric amplifiers (Orpheus-1 and Orpheus-HP, Light Conversion) driven by individual 15 and 6 W Yb amplifiers (Pharos-SP and Pharos-HP, Light Conversion) operating with 300 and 180 fs pulse durations, respectively. The 2D-IR spectrometer incorporated in the system allowed for rapid acquisitions with a few μOD of noise on 5000 laser shots.⁴⁸ Purely absorptive 2D-IR spectra were acquired in the pump–probe geometry using a collinear pair of pump pulses generated by a mid-IR pulse shaper (Phasetech Spectroscopy). The focused pump and probe spot sizes were measured to be ∼60 μm full width at half maximum (FWHM). The total amount of pump and probe light incident on the sample during rapid acquisition was ∼12 mW and 2–3 mW, respectively. Sample heating was estimated to be 2–3 °C.⁴⁸ The instrument response time was determined to be ∼270 fs FWHM using optical Kerr effect spectroscopy. Shot-by-shot interferometry (time step 22 fs) with phase cycling at the repetition rate of the lasers (100 kHz) was used to collect interferograms for Fourier transformation into 2D-IR spectra.

The 2D-IR spectra were recorded simultaneously in the parallel (⟨ZZZZ⟩) and perpendicular (⟨ZZXX⟩) directions for both SCN⁻ in EAN and H₂O. This enabled polarization dependent datasets to be collected with high signal-to-noise in a matter of minutes. For this paper, each pair of simultaneously collected ⟨ZZZZ⟩ and ⟨ZZXX⟩ 2D-IR spectra were recorded with 10 s of averaging. Thus, a total of 37 waiting times per temperature meant that 23 temperatures could be collected in under 3 h. Before the sample, a wire-grid polarizer was accurately set to polarize the probe beam at 45° to the bench. The pump pulses were kept parallel to the bench. After the sample, another wire gird polarizer was used to reflect and transmit the ⟨ZZZZ⟩ and ⟨ZZXX⟩ signals, respectively, which were routed to separate 128 element MCT detectors. The ⟨ZZXX⟩ signal was corrected after collection by multiplying the amplitude of the 2D signal by a correction factor (0.89, i.e., 10% error). The correction factor was determined by adjusting the anisotropy from a negative offset attributed to a 10% transmission error in the wire-grid polarizer. The spectral resolution was 4 cm⁻¹ for ω₁ and 3 cm⁻¹ for ω₃. The population time (t₂) was varied up to 11 ps for both samples. A thermalization effect was evident for EAN at a long t₂ (>11 ps); the short vibrational lifetime of the CN stretch in H₂O (2 ps) limited the effective t₂ range (∼5.5 ps).

The timescales of molecular motions are revealed by the change in shape of the 2D-IR spectra as a function of t₂. Molecules in different local environments can have different vibrational frequencies. As a result, many Lorentzians make up a vibrational band due to the “natural” linewidth. When t₂ is short compared to the molecular motions, the molecules will have limited time to move around. The initial frequency is highly correlated with the final frequency, and the 2D peak stretches along the diagonal. As t₂ increases, the molecules experience different local environments and lose correlation to their original environment, and the 2D peak becomes round. The random motions of molecules cause a random walk in frequency space, called spectral diffusion.⁴⁹ 

Center line slope (CLS)⁵⁰ analysis is used to monitor the shape of the 2D spectra for both EAN and H₂O as a function of population time (t₂). The change in shape of the 2D spectra can be quantified by various methods—ellipticities,^51,52 CLS,^50,53,54 and nodal-slope⁵⁵ and phase-slope.⁵² All of these quantities are approximately equal to the normalized vibrational frequency fluctuation function in the limit of inhomogeneous broadening. In this work, we used CLS to characterize the frequency fluctuations; the error in the measurements was estimated by propagating the uncertainty through the peak maximum determination, the slope fit, and the correlation fitting. For the fitting, a total of 10 ω₁ slices were taken from 2040 to 2082 cm⁻¹. Each slice was fit to a double Gaussian to determine the local minimum over an ω₃ range of 2040–2082 cm⁻¹. A total of 29 CLS values were used to extract the PW-FFCF from a t₂ of 0.2–11 ps.

The linear IR spectrum of the nitrile stretching mode (ν₃) of SCN⁻ EAN provides the first hint that there may be multiple local environments of SCN⁻ in EAN (Fig. 1). In H₂O, the center frequency of the ν₃ mode of SCN⁻ is 2065 cm⁻¹ with a FWHM of ∼36 cm⁻¹. The center frequency of the ν₃ mode of SCN⁻ in EAN is 2059 cm⁻¹ with a FWHM of ∼50 cm⁻¹; the ν₃ band is also asymmetric with a shoulder at 2070 cm⁻¹. The appearance of a shoulder suggests that SCN⁻ populates two different environments in EAN, though it is by no means conclusive. The 2D-IR spectra provide a more rigorous test of this hypothesis.

The overlapping features in the FTIR spectrum are also unresolved in the 2D-IR spectra of SCN⁻ in EAN. The 2D-IR spectra of SCN⁻ in EAN are shown in comparison with SCN⁻ in H₂O [Fig. 2(a)]. In EAN, the inhomogeneous linewidth of the nitrile band is ∼40 cm⁻¹ at early t₂ times. The 2D peaks remain stretched along the diagonal at 10 ps, indicating that spectral diffusion is not complete in the timeframe of the experiment. In addition, while spectral diffusion occurs, the center frequency of the nitrile band shifts from 2065 cm⁻¹ at t₂ = 0.2 ps to ∼2075 cm⁻¹ at 10 ps and simultaneously narrows. In comparison, the inhomogeneous linewidth in H₂O is narrower than in EAN, ∼20 cm⁻¹, and the initial tilt of the spectrum is less than in EAN. By ∼1 ps, the tilt of the 2D band is nearly gone, and, beyond 5 ps, there are no more changes except the loss of signal due to vibrational relaxation. The shifting and narrowing of the band observed in EAN is not observed in H₂O.

To visualize the change in the shape of the 2D-IR spectrum of SCN⁻ in EAN, we fit the slices of the spectrum along the pump axis (ω₁) to two Gaussian peaks, one positive and one negative, to account for the negative ground state bleach (GSB) and stimulated emission (SE) and the anharmonically shifted, positive excited state absorption (ESA). These fits determine two parameters, the peak frequency, ω_center, and the peak width, σ, which are functions of the initial frequency, ω₁. The anharmonic shift is constant with frequency and waiting time (∼20 cm⁻¹).

For Markovian dynamics on a harmonic potential, one expects the centers of the peaks, ω_center, to relax to the center of the band (∼2065 cm⁻¹).^50,51,56 Vibrational modes initially at both low and high frequency will, on average, move toward the average frequency as a function of time. Additionally, one expects the widths of the peaks, σ, to begin at the homogeneous linewidth and then increase together until they reach the average linewidth.

This simple expectation, however, is not what is observed here. The center frequencies [Fig. 3(a)] initially move to toward the center (during the first 1–2 ps) and then all begin to shift to higher wavenumbers (∼2070 cm⁻¹). The widths start near the homogeneous width, ∼12 cm⁻¹, and they all increase with time. One side of the spectrum, however, broadens more than the other [Fig. 3(b)].

The vibrational relaxation rate also seems to vary across the spectrum. To isolate population relaxation from orientational relaxation, we calculate the isotropic 2D-IR spectrum, S_iso(ω₁, t₂, ω₃), from the ⟨ZZZZ⟩ and ⟨ZZXX⟩ measurements, S_ZZZZ(ω₁, t₂, ω₃) and S_ZZXX(ω₁, t₂, ω₃),

Each slice along ω₁ is fit to positive and negative Gaussian profiles representing the ESA and GSB/SE peaks in the 2D-IR spectrum. We assume single exponential loss of the ESA and SE, which refill the GSB. Fitting the integrated Gaussian representing the ESA at each ω₁ to an exponential decay, we determine the vibrational relaxation time constant, T₁ [Fig. 3(e)]. At low ω₁ (2030 cm⁻¹), the apparent T₁ time is ∼2 ps, and at high ω₁ (2090 cm⁻¹), it is ∼4 ps. The loss of signal on the low frequency side of the spectrum relative to the high frequency side of the spectrum is, therefore, due to faster vibrational relaxation. This effect causes the peak in the 2D-IR spectrum to shift as a function of t₂.

To test the hypothesis that two overlapping sub-ensembles with different vibrational relaxation dynamics cause the observed changes in the spectra, we simulate the spectra and perform a global fit of the model parameters. We calculate third-order response functions based on truncating the cumulant expansion at second order (equivalent to assuming Gaussian statistics). For each sub-ensemble, we include independent, free-fitting variables for the timescales of frequency fluctuations, orientational relaxation, and population relaxation.

To describe frequency fluctuations, we used a Kubo lineshape⁴⁹ function incorporating a fast (homogeneous) and a slow component,

where T₂ is the dephasing time, Δ is the range of frequencies sampled by the chromophore, and τ is the correlation time. To describe orientational relaxation, we assume single exponential decay of the orientational correlations, and the appropriate orientational averages and orientational relaxation were incorporated to reproduce parallel and perpendicular spectra.⁴⁹ Errors are estimated by bootstrapping.⁵⁷ The nonlinear least squares fitting is repeated on 100 synthetic datasets, in which random samples of the full dataset are chosen, with replacement, to preserve the number of data points. The observed variation in fit parameters serves as an estimate of the error. To describe population relaxation, we assume single exponential loss of the ESA and SE, which refills the GSB and is described by a time constant, T₁. In the fitting procedure, each sub-ensemble has an independent lineshape, orientational relaxation time, and vibrational lifetime. The ⟨ZZZZ⟩ and ⟨ZZXX⟩ were fit simultaneously at each respective t₂ and temperature.

The complete set of best fit parameters is given in the supplementary material (Table S1). The most notable results are that the spectra are best fit with two sub-ensembles separated by ∼15 cm⁻¹. The low-frequency sub-ensemble contributes a higher initial intensity by a factor of 1.4 but decays more rapidly (T₁ = 2.0 ps vs 5.7 ps). The lower frequency sub-ensemble has faster dynamics, both reorientational times (22 ± 1 ps vs 34 ± 1 ps) and spectral diffusion times (3.3 ± 0.5 ps vs 11.1 ± 0.5 ps).

Overall, the two-population 2D-IR simulation replicates the main features in the 2D-IR spectra of SCN⁻ in EAN data (Fig. 4). At early t₂, the GSB peak is elongated along the diagonal. As t₂ increases, the GSB visibly shifts to the higher frequency side as the lower frequency region of the GSB relaxes in a few ps. At a long t₂ (10 ps), there is a systematic error in the residuals due to a hot ground state (HGS) that is not accounted for.

Analyzing the simulated spectra in the same way as the experimental spectra (Fig. 3) shows reasonable, semi-quantitative agreement. The band centers shift toward higher frequency on the same timescale as the experiment [Fig. 3(b)]. As the lower frequency population relaxes, the longer-lived population dominates and causes the spectrum and the peak centers to converge toward higher frequency, though the magnitude of the shift is somewhat less than the experiment. The widths of the bands also increase [Fig. 3(d)], though the spread is less than in the experiment. Finally, the trend of the apparent T₁ as a function of ω₁ for the simulated data is similar to the experiment [Fig. 3(e)]. The two populations overlap, and T₁ increases smoothly across ω₁. At low ω₁ (2030–2060 cm⁻¹), the trend of the apparent T₁ is concave up, indicating the most contributing sub-ensemble in that frequency range. [Note that the T₁ times extracted from treating the simulated spectra like the experimental spectra are similar but not identical to the T₁ times in the best-fit simulations (2 and 5.7 ps).] At high ω₁ (2060–2080 cm⁻¹), the apparent T₁ goes through an inflection point revealing a new sub-ensemble. On the high frequency side, there is a moderate difference in magnitude of the apparent T₁, which we attribute to systematic error from thermalization effects (hot ground state).

The simulated spectra do not require exchange between the sub-ensembles to accurately reproduce the evolution in 2D-IR lineshapes. Two sub-ensembles with different vibrational frequencies could generate chemical exchange crosspeaks if the local environments switch within the vibrational relaxation time of the vibrational probe. In the current case, the short T₁ time of the lower frequency sub-ensemble, 2 ps, precludes the ability to detect the timescale of interconversion of these local environments.

There are some limitations to the spectral modeling. At long t₂ (>5 ps), the two-population simulated model falls short of duplicating the exact experimental trends for the center frequency and width and at high frequency for T₁. The HGS shifts the experimental center frequency artificially increasing the lineshape width. In addition, the long-lived HGS inflates T₁ at high frequencies. There may also be inaccuracies in the model lineshape due to several assumptions. The model assumes that (1) the FFCF is polarization independent, (2) the FFCF is single exponential, and (3) orientation and frequency factorize. There are indications that each of these assumptions is not rigorously correct, which will be highlighted in Secs. III D and III E. Despite these limitations, the model supports the hypothesis that two sub-ensembles with different vibrational relaxation rates are responsible for the apparent shift in the center position of the peaks in the 2D-IR spectra and the observed frequency-dependent variation of the spectral linewidth.

In Sec. III D, we show the correlation between the orientational dynamics and the initial frequency of the vibrational chromophore by examining the 2D rotational anisotropy, which, then, motivates the subsequent RISD analysis (Sec. III E 2).

Rotational anisotropy in the 2D-IR data is inconclusive regarding the two overlapping sub-ensembles in EAN. The full two-dimensional rotational anisotropy and the timescales for decay of rotational anisotropy across the 2D-IR spectrum indicate that rotational motion and spectral diffusion are coupled.

The two-dimensional anisotropy map directly shows how the anisotropy depends on the vibrational frequency [Fig. 5(a)]. As was first demonstrated for liquid water,⁴³ a two-dimensional anisotropy map can be calculated from each point in a set of 2D-IR spectra,

At early times, t₂ = 0.2 ps, the anisotropy is near 0.4 along the diagonal and decreases to below 0.3 away from the frequency diagonal. If rotations and frequency fluctuations were independent, then the value should be constant across the entire ω₁, ω₃ plane.

The 2D anisotropic signal decays at different rates for different frequencies. Timescales can be extracted by fitting each point of the 2D anisotropic signal to a biexponential function. The amplitude of the fastest motions is small (∼0.05) for all points (data not shown), but including this term increases the quality of the fits. The weighted correlation time, $τrottotal$⁠,

varies from ∼25 ps near to the ω₁–ω₃ diagonal to ∼5 ps in the off-diagonal regions [Fig. 5(b)]. This is the first important indicator of the coupling of rotations and vibrational frequency fluctuations (Sec. III E 2). This observation suggests that molecules whose vibrational frequencies remain correlated also tend to maintain rotational correlations. Along the diagonal, there are two regions of slightly slower decay of the rotational anisotropy, which are suggestive of the two sub-ensembles, but the times for each region are similar.

Rotational anisotropy for the excited state absorption decays on longer timescales, up to ∼100 ps in some regions. In the ideal case, the rotational anisotropy should decay identically for GSB/SE and ESA. Our calculations suggest that a combination of the spectral diffusion and frequency shifts due to the frequency-dependent vibrational relaxation rate leads to a slower decay of the apparent rotational anisotropy in the region of the ESA; this effect is not observed in H₂O. As a note of caution, the longest population times in the experiment are 12 ps so the rotational timescales >10 ps should be treated as estimates based on the initial slope of the anisotropy decay. The reported times are thus lower bounds, and slower processes could be present. Thermal signals do appear in the 2D-IR spectra beyond ∼10 ps, and thermal transients can be anisotropic.¹³ Nevertheless, based on the deviation of the experiment and the simulated signals (Fig. 4), we estimate the error they could cause in the reported anisotropy to be less than 15% at ∼10 ps.

In summary, the two-dimensional rotational anisotropy shows that rotational motions and vibrational frequency fluctuations of the SCN⁻ are correlated. Section III E will treat the coupled dynamics with the theory for reorientation-induced spectral diffusion.

The 2D-IR signal depends on the vibrational frequency of the SCN⁻ probe molecule, the laser polarizations, and the orientation of the molecule with respect to the laser fields. Typically, the orientation and vibrational frequency of the vibrational chromophore are assumed to be uncorrelated, in which case the third-order response function factorizes into terms that depend on orientation and vibrational frequency. In such an ideal case, the vibrational frequency fluctuation correlation function can be determined uniquely in any polarization condition. When the molecule’s vibrational frequency depends on its orientation, however, one can no longer separate the frequency fluctuations from rotational motions. Thus, the measured frequency fluctuation correlation functions depend on the polarization condition employed in the experiment. To emphasize this dependence on polarization, the extracted correlation functions will be referred to as polarization-weighted frequency fluctuation correlation functions (PW-FFCFs).

This section will describe the measured PW-FFCFs and subsequent analysis using a model of RISD and SSD.

The PW-FFCFs of SCN⁻ in H₂O and EAN depend on the laser polarization configuration, which indicates RISD (Fig. 6). Separating the RISD from the PW-FFCFs, c₂(t₂), of SCN⁻ shows the dependence of the frequency fluctuations on rotational motions. To extract the PW-FFCF, we perform a CLS analysis over both sub-ensembles (2040–2082 cm⁻¹); thus, the PW-FFCF consists of contributions from both low- and high-frequency sub-ensembles, which is reasonable, given the broad similarity of rotational times. The CLS of both the ⟨ZZZZ⟩ and ⟨ZZXX⟩ spectra begins near ∼0.8 (a large inhomogeneity) (Fig. 6). Spectral diffusion is not complete for either the ⟨ZZZZ⟩ or ⟨ZZXX⟩ within the timeframe of the experiment (∼11 ps). As the CLS of ⟨ZZXX⟩ decays, it is always less than ⟨ZZZZ⟩, which is a strong indicator of RISD.

The PW-FFCFs of SCN⁻ in H₂O also show RISD, but it is less pronounced (Fig. 6, inset). In H₂O, the initial CLS value is much less than EAN. At an early time (t₂ ∼ 0.2 ps), the initial CLS values for the ⟨ZZZZ⟩ signal (0.296 ± 0.005) are similar but systematically larger than the ⟨ZZXX⟩ signal (0.267 ± 0.008). For both polarization signals, spectral diffusion occurs within a few ps, and again, the CLS of the ⟨ZZXX⟩ signal is always lower than ⟨ZZZZ⟩.

The frequency fluctuation dynamics are clearly much slower in EAN than in H₂O (Table I). Different mathematical functions are required to extract the PW-FFCF from the CLS results for SCN⁻ in EAN and in H₂O. In EAN, a biexponential fit captures both fast and slow frequency fluctuations; in H₂O, a single exponential is sufficient to fit the decay. Since the individual time constants and amplitudes are cross-correlated due to the linear dependence of the exponential functions, we also report an integrated correlation time, $τtotalFFCF$⁠,

The integrated correlation time accounts for motions entering the motionally narrowed limit; these motions contribute less to spectral diffusion, which is observed as a decrease in the amplitude of the CLS. The $τtotalFFCF$ extracted from ⟨ZZZZ⟩ is slower than that for ⟨ZZXX⟩ in both EAN and H₂O, as expected for RISD. The timescales of frequency fluctuations are slower in EAN than H₂O by ∼10 ps. One explanation for a slower timescale is that EAN (32 cP⁵⁸) is 30 fold more viscous than H₂O. The individual τ_i vary by a factor ∼14, and the τ_total timescales differ by a factor of 31.7 ± 0.5, nearly the same ratio as the viscosities of EAN and H₂O. We explore the assignment of the EAN dynamics further in Sec. IV and in Paper II of this series.⁴⁷ 

While we fit the low and high frequency parts of the spectrum simultaneously, the CLS values and PW-FFCF fits are mildly dependent on the ω₁ range for SCN⁻ in EAN (data not shown). While the exponential fit constants are consistent with faster dynamics on the low frequency edge and slower dynamics on the high frequency side, the differences were limited to the first ∼2 ps and the 2–10 ps range decayed with very similar rates. Having few data points along ω₁ also leads to uncertainty in the slope of center lines, so the error bars in the low and high windows largely overlap with the average CLS.

We also attempted to fit the two overlapping sub-ensembles with a method proposed by Yuan and Fayer.⁵⁹ This approach also failed to resolve separate dynamics for the two sub-ensembles (supplementary material). Because the two sub-ensembles overlap in frequency (much more than in the work of Yuan and Fayer⁵⁹) and one spectral feature decays in amplitude much faster than the other, the extracted fitting parameters are strongly cross-correlated.

Neither the ⟨ZZZZ⟩ nor ⟨ZZXX⟩ PW-FFCFs alone indicate the presence of more than one sub-ensemble of SCN⁻ in EAN. To determine CLS requires fitting over a frequency range (2040–2082 cm⁻¹). We interpret the resulting fits as a weighted average of the two sub-ensembles. Nevertheless, the FFCFs decay in a nearly exponential fashion in the observed time-window, leaving no obvious signatures of differences in dynamics. We suspect that this is a result of the significant overlap of the two distributions, which makes separating the sub-ensembles in the spectrum impossible, and also of the differences in vibrational relaxation rates, which biases the extracted PW-FFCFs to the longer-lived species.

Analyzing the 2D-IR spectra directly provides an average characterization of the frequency fluctuation dynamics of SCN⁻ in EAN, but it was unable to capture robustly the differences in dynamics of the two sub-ensembles. In Sec. III E 2, further analysis based on a model of RISD demonstrates that hindered rotational motion is sufficient to explain the observed spectral diffusion of SCN⁻ in EAN.

One framework to understand the frequency fluctuations of a vibrational chromophore considers the coupling of the vibrator’s dipole to the local electric field created by the surrounding solvent, such as a Stark effect. This model naturally couples rotational motion and the vibrational frequency, which explains the observed variation of timescales and activation energies in the different polarization conditions.

The simplest qualitative picture of reorientation-induced spectral diffusion (RISD) is the Stark effect. If the vibration of a molecule is anharmonic and has a projection onto the dipole of the molecule, then the vibrational frequency of that mode will depend on the angle between the local electric field and the molecular dipole. Molecules aligned with the local electric field will have frequencies shifted in one direction, molecules aligned against the field will have frequencies shifted in the other direction, and molecules perpendicular to the field will have no net shift. The constant of proportionality for this frequency shift is the Stark tuning rate. The consequence for 2D-spectroscopy is that different polarization conditions will observe different rates of spectral diffusion. In the all parallel configuration, ⟨ZZZZ⟩, molecules that have not rotated contribute the most to the spectrum. In the crossed-polarization configuration, ⟨ZZXX⟩, molecules that have rotated contribute the most to the spectrum. Molecules that do not rotate tend to keep their orientation relative to the local solvent field, so they do not change their vibrational frequency as much as the molecules that do rotate and tend to move to a different angle relative to the local field and a new vibrational frequency. As a result, the FFCF measured in the ⟨ZZZZ⟩ polarization condition will relax to zero more slowly than the ⟨ZZXX⟩ polarization condition when there is RISD.

A theory for the effects of RISD on 2D-IR spectra is derived in the context of the Stark effect,^44,45 but the essential physical picture is the same even when the origin of the vibrational frequency shifts is not a true Stark effect. Theoretical work has shown that both electrostatics beyond the dipole approximation and chemical interactions, such as charge transfer, play important roles in determining the C≡N stretching frequency.^25,60,61 Nevertheless, if these interactions are vectorial in nature, i.e., they depend on the orientation of the molecule, then they can be incorporated into the conceptual framework of RISD. For SeCN⁻, the electric field at the position of the C projected onto the CN bond-vector does accurately predict the C≡N stretching frequency, but this is largely due to a cancellation of effects.²⁵ This phenomenological association of the local electric field with vibrational frequency still strongly supports the concept that the vibrational frequency will be sensitive to molecular orientation.

The contributions of molecular rotation and structural fluctuations of the solvent can be quantitatively separated.^44–46 As a starting ansatz, we take the PW-FFCF to be a sum of two frequency fluctuation correlation functions,

Here, δω_v is the frequency shift due to interactions that depend on molecular orientation, so they are vectorial in nature, and δω_s is the frequency shift due to orientation independent effects that are scalar in nature. The term $Cpv$ carries all of the polarization dependence and represents the relative orientation of the solvent cage. The scalar term, C_s, is polarization independent and relates to isotropic physical effects, such as density fluctuations. The polarization independent scalar term becomes necessary to include when the two PW-FFCFs (⟨ZZZZ⟩ and ⟨ZZXX⟩) converge to a non-zero value (an offset) in the system.⁶² We neglect this term, as in previously studied systems,^44,45 because the PW-FFCFs do not converge to an offset. In this approximation, the PW-FFCF can be simplified to the product of SSD, F(t), and RISD, R_p(t),

Additionally, the time dependence of the RISD, R_p(t), can be expressed solely in terms of the reorientation time for each polarization condition after the appropriate orientational averaging of the signals.^44,45 For ⟨ZZZZ⟩, the result is

while for ⟨ZZXX⟩, the result is

where Δ is a normalization factor for the total frequency fluctuation amplitude. Assuming the reorientation of the dipole’s motion to be diffusive, the lth-order Legendre polynomial orientational correlation function is⁴⁴ 

where the orientational diffusion constant, D_m, is inversely proportional to the molecular rotational timescale from an anisotropy experiment,

To emphasize the parametric dependence of R_p(t) on the reorientation time, τ^rot, we will write it as R_p(t; τ^rot).

In H₂O, a model of free rotational diffusion of SCN⁻ and a single exponential decay of SSD,

adequately describe the data (Fig. 6, inset). The total SSD timescale (⁠$τtotalSSD=aSSDτSSD$⁠) is similar to the rotational timescale (τ^rot), and therefore, both rotations and local structural changes contribute to spectral diffusion of SCN⁻ in H₂O. The best fit parameters for this model are given in Table II.

In EAN, neither a model based on a single rotational time and a single exponential SSD [Eq. (12)] (data not shown) nor biexponential SSD,

fits the ⟨ZZZZ⟩ and ⟨ZZXX⟩ data simultaneously [Fig. 7(a)]. The biexponential SSD fit falls between the two PW-FFCFs. The data suggest the presence of both faster rotational motion, which causes the observed separation of ⟨ZZZZ⟩ and ⟨ZZXX⟩, and also slower reorientation, which maintains the offset from zero.

Models that restrict or bias reorientation naturally provide these two timescales. The wobbling-in-a-cone model provides convenient analytical forms to parameterize the PW-FFCF data. Following Kramer et al.,⁴⁵ the orientational correlation function, L_l(t), can be approximated in terms of an exponential decay to an offset, S_l. The term S_l is a complicated function of the order, l, and the half-cone angle, θ. The time constant of the decay, τ_l,eff, itself depends both on the rotational diffusion constant, D_m, and the order, l,

This form of L(t) can be inserted in Eqs. (8) and (9). We found that a constant SSD, F(t) = a^SSD, was sufficient to fit the data. The constant accounts for the product of the Stark tuning rate and the variance of the electric field fluctuations due to the solvent [Ref. 45, Eq. (15)]. Our final fitting function for the PW-FFCFs was then

The best fit parameters for this model are summarized in Table II.

We were also able to model the data successfully as the sum of two rotational terms,

where $τ1rot$ and $τ2rot$ are two independent rotational timescales. This model fits the PW-FFCF data because it also provides two rotational times (decay to an offset), but the model is not consistent with the observed 2D anisotropy or frequency-dependent CLS (data not shown). We, therefore, disfavor this model and prefer the model of restricted reorientation [Eqs. (14) and (15)].

In summary, for SCN⁻ in EAN, a model of restricted reorientation fits the PW-FFCFs. Rotational motion on a 2.6 ps timescale fully accounts for the observed spectral diffusion (Table II), and no SSD time can be extracted. For SCN⁻ in H₂O, both SSD and RISD contribute to spectral diffusion with roughly similar time constants (Table II).

Our experiments suggest two sub-ensembles for SCN⁻ in EAN. In previous pump-probe measurements, heterogeneous dynamics were detected in ammonium based PILs.¹⁸ Here, population relaxation of SCN⁻ shows strong indications of dynamical heterogeneity in EAN. One sub-ensemble is characterized by a lower vibrational frequency (2059 cm⁻¹) and shorter vibrational relaxation time (∼2 ps). The other sub-ensemble has a higher vibrational frequency (2070 cm⁻¹) and longer vibrational relaxation time (∼5 ps). Based on these features, we propose a molecular description of the solvation environment around the SCN⁻ in each sub-ensemble. Based on the calculations of hydrogen bonding to SeCN⁻ in water,⁶³ the apparent vibrational frequency shifts suggest that the low frequency ensemble has more, but weaker, hydrogen bonds, while the high frequency ensemble has fewer, but stronger, hydrogen bonds. Many other experiments have used the rate of vibrational relaxation of SCN⁻ as a probe of the local solvation environment.^{22,28,36,64–66} Ren et al.³⁶ also identified sub-ensembles of SCN⁻ with different vibrational relaxation rates in an aprotic ionic liquid. In that case, the lower frequency, shorter lived sub-ensemble was identified as an effect of water in the first one or two solvation shells. In the present case, the linear spectra and 2D-IR dynamics did not change significantly over the range of water content that we could access (1 M–40 mM), which suggests, but does not completely rule out, that the effect is not due to the presence of water. In a Fermi’s golden rule picture,^64,66,67 the vibrational lifetime is proportional to the product of the square of a coupling constant and the vibrational density of states of accepting modes. Our experiments are not able to separate these two possible contributions, but at least one of the two quantities must be high enough to double the rate of vibrational relaxation.

Our proposed assignment is also consistent with the heterogeneous dynamics observed in molecular simulations of imidazolium ionic liquids and EAN. Molecular simulations^20,21 provide evidence that translational and rotational motion of a vibrational probe are slower when it is involved in a strong, specific hydrogen bond. In imidazolium ionic liquids, these structures are associated with charge dense regions of the liquid, while the faster domains are associated with relatively charge-depleted regions. EAN also consists of regions that are enhanced and depleted charge density;¹⁵ the charge-depleted regions are relatively concentrated in the alkyl chain length, and the charge dense regions are relatively concentrated in the ammonium headgroups and nitrate ions. Similarly, molecular simulations show that the cation–anion regions of EAN resulted in slower reorientation dynamics due to the strong interactions, while the alkyl chain region reoriented on a faster timescale due to a lack of strong Coulombic interaction.¹⁹ 

In EAN, the frequency fluctuations are dominated by rotational motions, while in H₂O, both SSD and RISD contribute. While dynamical heterogeneity is present with respect to T₁ in EAN, the 2D anisotropy plot reveals that there is no difference in anisotropy along ω₁. Therefore, it is reasonable to analyze the dynamics over the average of the two sub-ensembles. In EAN, RISD fully captures the frequency fluctuation dynamics because the SSD is unresolved. The RISD in EAN, however, is not biexponential, i.e., a separate rotational timescale for each sub-ensemble. Instead, a wobbling-in-a-cone RISD approach provides a single fast timescale (∼2 ps). The timescale correlates with SCN⁻ undergoing free diffusive rotations in a constrained angular space. In previous pump–probe measurements,¹³ a long-lived rotation was identified by examining the ND stretch of EAN with a large angle jump (106°) comparable to the angle extracted from the wobbling-in-a-cone approach used here (108°). The 2D anisotropy analysis suggests longer rotational timescales (20–40 ps), which could be large angle jumps to a new local environment where the SCN⁻ explores a new cone. The observed cone angles are large, which is inconsistent with a specific steric hindrance. (Considering both the S and N ends of the molecule, the SCN⁻ would sweep out a full sphere even at half-angles of 180°.) The cone could also be the effect of an energetic bias in the rotational distribution, which was also observed for SCN⁻ in an ionic liquid colloidal dispersion.^42,62

As a word of caution, the two-population model may be a simplification of additional solvation sub-ensembles, and molecular simulations are needed to test our interpretation of how EAN solvates SCN⁻. The insights into the solvation environment of EAN have stimulated a temperature-dependent 2D-IR study where we explore an interpretation of the thermodynamics of the transition state of SCN⁻ in EAN compared with H₂O in Paper II.⁴⁷ 

In this paper, we quantitatively characterize the dynamics and heterogeneity of the 3D hydrogen bonding network in the PIL EAN. The 2D-IR experimental approach used allowed us to gather the data in a timescale of hours by collecting the ⟨ZZZZ⟩ and ⟨ZZXX⟩ simultaneously at 100 kHz.⁴³ Combined with the analysis approach synthesized from the work of others,^43–45 this approach provides a powerful protocol for the further characterization of other PILs with a view to establishing 2D-IR spectroscopy as a useful tool for developing molecular models of PILs in support of tuning their properties for practical applications.

In summary, we observe dynamical heterogeneity in the solvation of the small ionic solute SCN⁻ in the 3D hydrogen bonding network of the PIL EAN. The linear FTIR spectrum reveals a high frequency shoulder on the vibrational band of the nitrile stretch. The 2D-IR spectra show complex evolution with population time due to frequency-dependent rates of vibrational relaxation. Simulations of the spectra support the hypothesis that SCN⁻ resides in at least two structural sub-ensembles with vibrational relaxation rates that are different by a factor of more than two. Vibrational energy relaxation is, itself, an important reflection of the stochastic coupling of the probe to its bath and is critical in many physical and chemical processes.⁶⁷ The observed vibrational relaxation rates result from either differences in the magnitude of the coupling of the ν₃ mode to the bath or differences in the density of states of the bath able to accept energy from the ν₃ mode. Vibrational frequency fluctuations and reorientational correlations for the two sub-ensembles, however, are the same with the experimental error. Applying a RISD analysis, we determine that the frequency fluctuations are dominated by rotational motion (i.e., are in the RISD limit), and a wobbling-in-a-cone model extracts a single fast timescale of molecular rotation. There are slower reorientations (∼30 ps timescale) responsible for the final loss of rotational anisotropy. Finally, we postulate that the two different sub-ensembles in EAN might be correlated with charge-depleted (alkyl chain) and charge dense regions (ammonium headgroups and nitrate ions). In Paper II,⁴⁷ we develop the investigation of the dynamics of hydrogen bonding of SCN⁻ in EAN and H₂O to obtain the temperature dependence and activation energies of SCN⁻ in EAN and H₂O.

Details of the spectroscopic fitting, signal-to-noise, center line slope, and ellipticity fits are provided in the supplementary material.

This work was supported by the National Science Foundation (USA) under Award Nos. CHE-1454105 and CHE-1954848, the BBSRC from the grant Alert 13 (BB/L014335/1), and the STFC Central Laser Facility. We thank Maria Grazia Izzo for championing early 2D-IR studies of the NH and ND stretches of EAN at the Central Laser Facility, Nishal Chandarana for practical assistance, and Sunayana Mitra for support and assistance.

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