An OHD-RIKES and simulation study comparing a benzylmethylimidazolium ionic liquid with an equimolar mixture of dimethylimidazolium and benzene

The principal difference between 1-benzyl-3-methyl-imidazolium triflimide [BzC1im][NTf2] and an equimolar mixture of benzene and dimethylimidazolium triflimide [C1C1im][NTf2] is that in the former the benzene moieties are tied to the imidazolium ring, while in the latter they move independently. We use femtosecond optical heterodyne-detected Raman-induced Kerr effect spectroscopy (OHD-RIKES) and molecular simulations to explore some properties of these two systems. The Kerr spectra show small differences in the spectral densities; the simulations also show very similar environments for both the imidazolium rings and the phenyl or benzene parts of the molecules. The low frequency vibrational densities of states are also similar in the model systems. In order to perform the simulations we developed a model for the [BzC1im] + cation and found that the barriers to rotation of the two parts of the molecule are low.


Introduction
OHD-RIKES (optical heterodyne-detected Raman induced Kerr effect spectroscopy) is a nonlinear optical time-domain technique that gives information about the low frequency fluctuations in the polarizability of a liquid and is complementary to dielectric spectroscopy. [1][2][3] Because of its ease of use and high quality of the data, it is one the most common methods for studying the low frequency intermolecular modes of liquids [1][2][3] and in particular ionic liquids.  By use of a Fouriertransform-deconvolution procedure, 33,34 the OHD-RIKES timedomain data can be converted to a spectral density or Kerr spectrum, which is directly related to the depolarized Rayleigh/ Raman spectrum of the liquid. 35 The low-frequency part of a Kerr spectrum is dominated by intermolecular motions, but is unable to discriminate between different types of motion such as librational (hindered rotation) and translational. Further one does not obtain any direct information about the local environments of different molecular species. Molecular simulation can complement these spectroscopic results by giving a molecular level picture of the local environments of molecules and the local molecular dynamics. Recently Quitevis and co-workers have compared the Kerr spectra of the neat ionic liquid 1,3-dimethylimidazolium triflimide, [C 1 C 1 im][NTf 2 ], with its mixtures with benzene and with neat benzene. 36 Simulations were used to show that the benzene response was at a lower frequency than that of the ionic liquid and that in a 1 : 1 molar mixture of this ionic liquid with benzene the benzene response moves to a higher frequency while the [C 1 C 1 im] response moves to a lower frequency. This is not unexpected as it implies that the ions move more freely on dilution with benzene while the benzene motion is impeded on dilution with an ionic liquid. Translational motion tends to be at a lower frequency than the librational motion. These results were based on an analysis of different contributions to the vibrational density of states.
This work on ionic liquid-benzene mixtures led Quitevis and coworkers to ask what would happen if the benzene and imidazolium were tied together as 1-benzyl-3-methyl-imidazolium [BzC 1 im] + as shown in Fig. 1. Shirota et al. 37 performed a study of the lowfrequency Kerr spectra and physical properties of ionic liquids functionalized with phenyl moieties, including the [BzC 1 im][NTf 2 ]. In this paper we describe the local molecular environments found in simulation and Kerr spectra of [BzC 1 im][NTf 2 ]. The current study goes beyond the work of Shirota et al. 37 in that we compare the results with those from a 1 : 1 mixture of C 6 H 6 and [C 1 C 1 im][NTf 2 ]. In order to do this, a force field for benzylimidazolium ions was needed; Section 3.1 describes how this was developed. The resulting local environments of the phenyl groups and the imidazolium rings are compared with those for the mixture in Section 5.1 while in Section 5.2 we compare the vibrational densities of states and the Kerr spectra.

Sample preparation
The syntheses of [C 1  The mixture was prepared by weighing a given amount of [C 1 C 1 im][NTf 2 ] in an argon-filled glovebox and adding the required weight of benzene to make up an equimolar mixture. A syringe fitted with a 0.1 mm filter was used to deliver aliquots of the pure liquids and the mixtures, while still in the argonfilled glovebox, into 2 mm path-length quartz cuvettes fitted with vacuum stopcock valves to insure airtight sealing of the samples during OHD-RIKES measurements.

OHD-RIKES apparatus and procedures
The apparatus has been described previously. 40,41 Briefly, a labbuilt, self-modelocked Ti:sapphire laser, pumped by a Coherent Verdi V6 laser, generates a train of 38 AE 3 fs pulses at a repetition rate of 82 MHz. A (01/451/À451) polarization configuration is used for the pump, probe and analyzer polarizers, with a l/4-plate inserted in the probe beam before the sample. Stray birefringences, such as from the windows of the sample cell, are nulled by adjusting the l/4-plate. A local oscillator (LO) field E LO is introduced for the signal field E S by slightly misaligning (AEd1) the probe polarizer. 42 The two signals obtained in this way are given by where I + is the signal with the in-phase (+d) local oscillator and I À is the signal with the out-of-phase (Àd) local oscillator. In these studies the heterodyne contribution <|E LO E S *|, which is linear in the material response, is obtained from the difference of I + and I À .
The temperature of the sample is kept constant at 295 K by placing the sample cell in a lab-built, thermostated copper cell holder, with temperature control from a thermoelectric heatercooler system. A data set is obtained by averaging at least 15 scans, with each scan being performed in 10 fs steps for time delays between À1 and 4 ps and in 50 fs steps for time delays greater than 4 ps to save the data collection time.

Analysis of the OHD-RIKES data
The OHD-RIKES signals in the 0.4-10 ps time range are fitted by the empirical decay function where the t 1 and t 2 terms are sub-picosecond components, the t 3 term is a picosecond component, and B is a constant that accounts for components in the decay function relaxing on a time scale much longer than the time-range of the measurements. OHD-OKE time-domain measurements 6,20,25 over a much wider range of timescales (from sub-ps to sub-ns) than that of our measurements clearly indicate ionic liquids are slowly relaxing systems. Maroncelli and coworkers 43 have shown through molecular dynamics simulations that the slow dynamics are attributed to the tendency of the charge-ordered structure of an ionic liquid to resist reorganization. Moreover, the reorientational dynamics are consistent with that of a complex fluid as reflected by an intermediate power law at short times (few ps to hundreds of ps) followed by a Schweidler power law and a temperature dependent exponential. In the frequency domain, these regimes correspond to low-frequency modes associated with aand b-relaxations. Clearly in the time window of our measurements, the fit of the empirical decay function to the data in Fig. 6 to an OHD-RIKES signal (e.g., Section 4.1 below) captures the intermolecular dynamics and the early part of the intermediate power region of the reorientational dynamics. By convoluting the part of the decay given by the t 3 term and the constant B with the pulse-intensity autocorrelation, a reorientational response is generated, which is subtracted from the OHD-RIKES signal to yield a ''reduced'' response consisting only of the electronic and the sub-picosecond nuclear responses. The Fourier-transform-deconvolution procedure is applied to the reduced response and results in the RSD (Reduced Spectral Density) corresponding to the part of the OKE spectrum associated with intramolecular and subpicosecond intermolecular modes of the liquid. The use of a window function, as described previously, 44 reduces the noise in the low-frequency band in the 0-200 cm À1 region of the RSD without affecting its line shape. In principle, the relative intensities of RSDs can be determined if the OHD-RIKES signals of each of the samples are measured under the same experimental conditions. In practice this is difficult to achieve because of having to make adjustments of the laser and tweak the alignment of the optics in the apparatus between runs. In our previous study 36 a referencing technique was used. However in the current study, because the OHD-RIKES signals for all the systems were obtained in a single run in such a way that constant experimental conditions were maintained for all the data sets, the referencing technique was not used.

Force field development for simulations
We use a standard type of force field with terms describing bond stretching, angle bending and dihedral motion within each molecular species together with partial charges and Lennard-Jones terms on each atomic site to describe intermolecular interactions. In earlier work 45 on the adsorption of ethylene and ethane in various ionic liquids the experimental work was complemented by some simulations by Padua and co-workers on ionic liquids including [BzC 1 im][NTf 2 ]. They used standard force fields and did not reconsider the charge or dihedral parameters. As we are here particularly concerned with the low frequency motions of the liquid, we examined the changes in energy with the dihedral rotations about the methylene group linking the phenyl and imidazolium rings as well as the partial charges.
The partial charges of some of the configurations generated in the determination of the dihedral potential were determined by iterated stockholder analysis (ISA). 46,47 In this method a calculated electronic charge density, (here obtained using density functional theory with PBE0 functional and asymptotic correction, and the Sadlej basis 48 ), is partitioned between atoms using a spherical weight function for each atom, and the resulting density for each atom is spherically averaged to obtain a revised weight function. This procedure is iterated to obtain an optimum partition into atomic densities that are as nearly spherical as possible. Lillestolen and Wheatley 46 showed that the method is guaranteed to converge, but their original grid-based method converged very slowly, and the basisfunction-based method of Misquitta et al. 47 converges very much more quickly. The method yields atomic multipole moments as well as charges, but the charges alone give a good account of the electrostatic potential. Fig. 2 shows the difference between the electrostatic potential from the full wave function and that generated at twice the van der Waals radii. It can be seen that the fit is good even though only site charges and not higher multipoles are included. The r.m.s. error over the surface is 0.03 V, about 1% of the potential, which varies between about 2 V and 3 V over the surface. The values of the partial charges vary by a small amount, of the order of 0.01e, with molecular conformation. For a classical simple point charge model one needs constant charges; we used an average of the charges obtained at the two geometries with energy minima (0 and 1351) truncated to two decimal places and adjusted to give a total charge of +1e.
The bond stretching and angle bending terms for the methyl imidazolium moiety were taken from the CLAPS force field 49,50 which is based on the OPLS force field. 51,52 As the lowest vibrational frequency of benzene is 467 cm À1 , the internal motions of the phenyl group are unlikely to affect the vibrational spectrum below 300 cm À1 , and the phenyl ring was treated as a rigid body. The Lennard-Jones terms were initially taken from these force fields, but finally were scaled by 0.95 to obtain the correct pressure of the liquid at 300 K and the experimental density.
In order to obtain the variation of molecular energy with the dihedral angle, f 1 , which describes the rotation of the imidazolium ring with respect to the bridging methylene group (see Fig. 3), conformations with fixed values of this angle were energy-minimised with respect to all other degrees of freedom. This was repeated for values of f 1 from 0 to 1801 in steps of 151 in a relaxed scan. The quantum chemistry calculations were carried out at the MP2 level using a cc-pVTZ basis set as implemented in the Gaussian 09 software package 53 Further details are given in the ESI. † The dihedral energy of this rotation is shown in Fig. 4 as red circles connected by a full red line. There is a global minimum at f 1 = 01. Rotation of the imidazolium group leads to a global maximum of 6.8 kJ mol À1 near 801, a secondary minimum of 2.8 kJ mol À1 near f 1 E 1401 and finally a secondary local maximum of about 3.8 kJ mol À1 at f 1 = 1801. As kT E 2.5 kJ mol À1 at 300 K this barrier would be easily surmounted at room temperature for an isolated ion. This is what we found in the simulations (see below).
The second dihedral angle of interest, f 2 , describes the rotation of the phenyl ring relative to the bridging methylene group (see Fig. 3). It was investigated in a similar way although the energy minimisations were carried out with the constraint that f 1 = 901. This was done to avoid steric clashes between the phenyl and imidazolium rings. The results of this second set of calculations are shown in Fig. 5 as red circles connected by a full line. The dihedral energy for rotation of the phenyl group has a minimum at f 2 = AE901 and a maximum at f 2 = AE01. The barrier is about 10 kJ mol À1 or 4kT at 300 K. The fitted function (open red circles) has no explicit dihedral terms connecting the phenyl group and the methylene group; the observed variation of energy with f 2 is solely due to other terms in the force field including electrostatic and dispersion interactions between the phenyl group and the rest of the molecule.
The dihedral terms in the force field for f 1 between the unique imidazolium carbon site and first phenyl carbon site (see Fig. 3) was expressed as the sum of three cosine terms and the coefficients A 1 , A 2 and A 3 were found by fitting the total force field energy to the quantum chemical results. The resulting energies are shown in Fig. 4 by the red circles connected by a dashed line.
As the dihedral term between the phenyl group and the imidazolium nitrogen site was found to be described reasonably well with only a scaling of the 1 : 4 electrostatic and dispersion interactions by 0.5, no further dihedral term was included. The complete force field is given in the ESI. † The distribution of conformers in the bulk depends on local intermolecular interactions as well as on the intramolecular potential. In simulations of bulk liquid we find that the observed free energies of both dihedral angles have considerably lower barriers than the single molecule energy results. These free energies A(f) (determined from the observed probability distributions p(f) by eqn (4)) are shown by full blue lines in Fig. 4 and 5.

Simulations and analysis
Bulk simulations of liquid [BzC 1 im][NTf 2 ] were then carried out using these potentials for the benzylimidazolium cation and the potential due to Ludwig and co-workers 54,55 for the triflimide anion [NTf 2 ] or (N(SO 2 CF 3 ) 2 ) À . This latter potential has a united CF 3 group. These simulations were carried out using a modified version of DL_POLY. 56 128 ion pairs in a periodically repeated box with sides of 4.014 nm were used and the electrostatics was treated by Ewald summation with a real space cutoff of 1.2 nm and convergence parameter = 2.3 nm À1 . This corresponds to a precision of 10 À5 in the energy. A time step of 2 fs was used and the temperature was 300 K. Simulations of the 1 : 1 mixture under similar conditions were described in our previous paper 36 and the results from that work used here. The vibrational densities of states are found by Fourier transforming the velocity autocorrelation functions of individual atomic sites. Here we report the density of states of the imidazolium ring and the phenyl group separately. These were calculated by taking the sum of the velocity autocorrelation functions for the 5 heavy atoms in the ring for the imidazolium density of states and the sum for the 6 carbon atoms in the     2 ] there is also a deep oscillation at t E 0.33 ps, with width Dt E 0.17 ps. If we assume this feature to be the first period of an underdamped oscillation of an intramolecular mode, then we estimated its frequency to be E100 cm À1 .
The solid curves in Fig. 6 are fits of eqn (2) to the data in the 0.4-10 ps time range, with fit parameters listed in Table 1. The dashed lines correspond to the slow t 3 term or t 2 + t 3 terms (where the t values are greater than 1 ps) plus the constant B. In addition to the fit parameters in eqn (2), the values of the average relaxation time for the fast part of the decay, defined by ht s i = (A 1 t 1 + A 2 t 2 )/(A 1 + A 2 ), are also listed in Table 1. From Table 1 it can be seen that ht s i of the benzene-[C 1 C 1 im][NTf 2 ] mixture has a much larger value (ht s i = 0.61 ps) than that of [BzC 1 im][NTf 2 ] (ht s i = 0.15 ps) which indicates that the OHD-RIKES signal decays more slowly in the mixture.
The value of t 3 = 6.4 AE 0.2 ps for the mixture is associated with the intermediate power-law region of the reorientational dynamics, which cannot be attributed solely to either C 6 H 6 or the IL. However the value of t 3 of mixture is much larger than the value for [BzC 1 im][NTf 2 ] (t 3 = 3.7 AE 0.2 ps).   Table 1.   2 ] mixture have been specifically discussed in our recent published paper. 36 To see more clearly how the intermolecular spectrum of [BzC 1 im][NTf 2 ] differs from that of C 6 H 6 -[C 1 C 1 im][NTf 2 ], the RSDs of these systems in Fig. 7 are overlapped with each other and compared in Fig. 8. As can be seen in the Fig. 8a

Multicomponent line-shape analysis
The featureless broad low frequency RSD spectrum contains contributions that arise from both intramolecular and intermolecular vibrational motions. The inter-and intramolecular dynamics overlap in this spectral region and cannot be distinguished unambiguously, so the RSDs in the 0-250 cm À1 region were fitted to an empirical multicomponent line shape model in order to gain further insight into these spectra. As described previously, in this model the lowest component is given by the Bucaro-Litovitz (BL) function  and higher frequency components by the antisymmerized Gaussian (AG) function This analysis allows us to separate intermolecular and intramolecular contributions in the RSD in the 0-250 cm À1 region. This multicomponent line shape analysis gives a quantitative description of the intensity of RSD in this frequency region, but there is no a priori reason that the component bands correspond to actual modes of the liquids. Fig. 9 shows multicomponent line shape fits of RSD in the 0-250 cm À1 region. The spectral parameters for the main band (first spectral moment M 1 and full-width-at-half maximum Do) and component bands (peak frequency, o pk ; peak width, Do; and fractional area f area ) that underlie the main band are shown in Table 2.
The spectral parameters of C 6 H 6 , [C 1 C 1 im][NTf 2 ] and the C 6 H 6 -[C 1 C 1 im][NTf 2 ] mixture have been discussed in detail previously. 36 The  Table 3 gives the frequencies, and relative intensities of the component bands, whereas the ESI † gives fit parameters corresponding to eqn (5) and (6). In both cases there is a high probability of finding either an imidazolium ring or another phenyl group in the volumes above and below the phenyl group or benzene molecule. The cutoffs are given as ratios to the average molecular number densities of each component. These are equal to 1.98 nm À3 and 1.81 nm À3 for [BzC 1 im][NTf 2 ] and the 1 : 1 mixture respectively. Fig. 11 shows the corresponding distributions around an imidazolium group. Here there is more of a difference between  Table 2.  the united ion and the mixture. In the mixture the distribution is symmetrical with evidence of stacking of benzene and [C 1 C 1 im] ions in the z or polar direction (perpendicular to the imidazolium ring) with the anions distributed around the equator in three broad bands. This distribution is clearly perturbed by the benzyl substituent. The local distributions around both benzene and imidazolium ions have been seen in earlier simulations 36,58 and in quantum mechanical ion pair studies. 59,60

Simulation results for densities of states
The lower panel of Fig. 12 shows the computed densities of states for dimethyl imidazolium in a 1 : 1 mixture and for the imidazolium ring in [BzC 1 im][NTf 2 ]. Below 100 cm À1 there is little difference between the densities of states, while at 150 cm À1 and above one sees low frequency internal modes in [BzC 1 im] + . The strong similarity of these two spectra can be correlated with the similarity in the local environment shown in Fig. 10. The upper  . The peak of the benzene density of states in the mixture is at a slightly higher frequency than the peak in the phenyl spectrum. This implies that the motion of the phenyl group is less impeded in [BzC 1 im][NTf 2 ] than in the mixture. This is an unexpected result, which suggests that the packing around a benzene molecule in the mixture is tighter than that around the unwieldy [BzC 1 im] + cation.

Discussion
The results shown in Fig. 7  are remarkably similar. The detailed comparison of the RSD plots in Fig. 8 and 9 show that the main difference is a decrease in the low frequency side of the main response compared to the high frequency side which shifts the maximum in [BzC 1 im][NTf 2 ] to higher frequency. There is also an increased response around 160-250 cm À1 which is due to internal modes involving the relative motion of the phenyl and imidazolium groups which are absent in the mixture. The Kerr signal in the time domain is due to the time correlation function of the anisotropic polarisability of the liquid as a whole. This is made up of contributions from individual molecules, interactions between molecules (such as the dipole-induced dipole terms), and cross terms between these two parts. At low frequencies these terms are affected by motions of molecules relative to each other. Hindered molecular rotation (libration) modulates the contributions from individual molecules, while hindered translation affects the interaction between molecules. All these motions are captured in the vibrational densities of states obtained from velocity correlation functions of individual atoms. These are easily calculated from a simulation. It is possible to simulate the Kerr response directly 18,31,61,62 but we have not done so. As the Kerr response is a collective property of the liquid or solution it is noisier than the single particle vibrational densities of states and so more expensive to calculate. Further it would require a polarisable model which we do not use. Rather we calculate the vibrational densities of states in the liquid. The advantage of calculating the single particle vibrational densities of states is that different contributions can easily be distinguished; the disadvantage is that cross correlations are missing and that the strength of different contributions cannot be estimated. Thus the simulation results in Fig. 12 cannot be directly compared to the observed results in Fig. 8. For example it is possible that differences in the relative weighting of the phenyl and imidazolium rings in the mixture and in [BzC 1 im] + causes the change in shape of the RSDs in Fig. 8. Shirota and co-workers 37 speculate that coupling of the phenyl librations in [BzC 1 im] + to the intraionic modes complicates the interpretation of the Kerr spectra. We can confirm the identification of the signal between 170 and 250 cm À1 in [BzC 1 im][NTf 2 ] as being due to floppy internal modes which show in the phenyl, the CH 2 and most strongly in the imidazolium response. It is surprising that the phenyl density of states in [BzC 1 im] + peaks at a slightly lower frequency than that for C 6 H 6 in the mixture, showing that the motion is less impeded when the phenyl group is tied to the imidazolium ring.
One can understand the similarities between the C 6 H 6 -[C 1 C 1 im][NTf 2 ] mixture and [BzC 1 im][NTf 2 ] to some extent by looking at the strong similarities in the spatial distribution functions in three dimensions shown in Fig. 10 and 11. The distributions around the benzene in the mixture and the phenyl group in [BzC 1 im][NTf 2 ] both show stacking of imidazolium and benzene rings above and below a benzene moiety. Although the presence of the phenyl group in the latter compound does distort the stacking to some extent and changes the distribution of anions in the equatorial planes of the phenyl moieties, the changes are small. Thus it is not surprising that the local intermolecular motions are similar. The distributions around the imidazolium rings in the mixture and in [BzC 1 im][NTf 2 ] are again similar, but the changes are perhaps more pronounced. Earlier we commented on the low values of the barriers to twisting the phenyl groups and imidazolium rings in this ion. It seems that although the phenyl and imidazolium parts of the ion are tied through a methyl group, they behave as if they are almost free.
It is difficult to interpret the difference between the Kerr spectra of the mixture and [BzC 1 im] + in terms of these simulation results. In both cases the differences are small but the experimental and simulation results show opposite trends. In the multicomponent analyses of the experimental line shapes shown in Fig. 9 the difference between the mixture and [BzC 1 im] + spectra arises from the relative strengths of the components rather than changes in their frequency. It is not possible to determine the relative strengths of different contributions from the simulations, so the changes in densities of states could be due to changes in the relative weighting. It is also possible that cross terms are important and also that the use of non-polarisable potentials is not adequate for this comparison.

Conclusions
We have performed OHD-RIKES experiments and simulations of [BzC 1 im][NTf 2 ] and compared the results with an equimolar mixture of [C 1 C 1 im][NTf 2 ] and benzene. The difference between these liquids is that in the first the benzene is tied to the imidazolium ring and in the second the benzene molecules are free to move relative to the cation. We found that the relative rotation of the two moieties of the [BzC 1 im] + ion is nearly free with only small barriers. The local environments of both imidazolium rings and phenyl groups are similar as are the computed densities of states for intermolecular motion. The experimental Kerr spectra are more similar to each other than to either neat [C 1 C 1 im][NTf 2 ] or neat C 6 H 6 , but that of [BzC 1 im][NTf 2 ] has a stronger high frequency component, possibly due to internal modes involving relative motion of the imidazolium ring and the phenyl group.