A simulation study of CS₂ solutions in two related ionic liquids with dications and monocations

Dicationic ionic liquids based on pairs of imidazolium ions linked by alkyl chains are liquid over a wider temperature range and are thermally more stable than the corresponding monocationic liquids, which make them of interest for electrolyte applications. An understanding of the physical chemistry and microscopic properties of these liquids is needed. Dielectric and RIKES spectroscopies are suitable experimental methods for studying liquids, while molecular simulation provides complementary information.

CS₂ is often used as a probe for optical heterodyne-detected Raman-induced Kerr effect spectroscopy (OHD-RIKES)^1,2 as it has a strong optical Kerr signal, which, in solution, is mainly due to the effect of its librational motion on the molecule’s large anisotropic polarisability giving rise to fluctuations in the total anisotropic polarisability of the liquid. The maximum frequency and the first moment of the spectra in the dicationic ionic liquids were observed to have higher values than those in the corresponding monocationic liquid, suggesting that the motion of the CS₂ probe is more restricted in the former.³ This is not unexpected as the viscosity is higher in dicationic liquids than that in the corresponding monocationic liquids, but atomistic simulation can give further insight into the motion of molecules in the liquid and the local environment.

Recently, Gurung et al.³ have been studying the Kerr spectra of solutions of CS₂ in pairs of imidazolium ionic liquids [C_nC₁im][NTf₂] and $[C2n(C1im)2][NTf2]2$⁠, with n = 3, 4, 5. The dications are formed from two monocations joined at the methyl end groups of the alkyl tails (see Fig. 1). The anion is bis(trifluoromethylsulfonyl)imide which is also known as [NTf₂]^– or [TFSA]^–. In this paper, we present the results of atomistic simulations of CS₂ solutions in [C₄C₁im][NTf₂] and $[C8(C1im)2][NTf2]2$⁠. These show that the maximum in the low frequency vibrational density of states of CS₂ does indeed shift to higher frequency in the dicationic solution at approximately 1 atm but is sensitive to pressure (or equivalently anion number density). At the same anion number density as the monocation solution, the maximum shifts to lower frequency in the dication solution. We also compare the local structure and dynamics in dication and monocation solutions at different state points.

Several simulation studies of neat dicationic imidazolium ionic liquids have been published^4–10 and some of tricationic imidazolium ionic liquid.¹¹ Shirota and co-workers^7–9 have studied the Kerr spectra of these liquids.

The intermolecular forces comprised site-site Coulomb and Lennard-Jones terms. In the interest of speed, methyl groups and methylene groups were modeled by united atoms situated on the C sites, but the imidazolium ring protons, which are important for hydrogen bonding, were included explicitly. Bond-bending and dihedral and bond stretching parameters for the alkyl chains of the cations were taken from the CL&P¹² and OPLS (Optimised Potentials for Liquid Simulations)¹³ force fields, while site charges were taken from Zhong et al.¹⁴ The anion was modeled with united CF₃ groups with parameters taken from Ködderman et al.^15,16 CS₂ was treated as a rigid body with Lennard-Jones’ centres on each atomic site and charges on these sites and 4 supplementary sites to describe the electrostatics.¹⁷ The model used was non-polarisable, as preliminary work showed that including the realistic polarisabilities from Lynden-Bell and Stone’s work¹⁷ made a negligible difference to the local structure and energetics and came with a significant computational cost. All the parameters for the force field are given in the supplementary material.

Molecular dynamics NVT simulations were performed with a modified version of the program DL_POLY.¹⁸ The cubic simulation cell contained 108 anions, 12 CS₂ molecules, and either 54 dications or 108 monocations, i.e., a mole percent of CS₂ of 18.2% and 10%. Note that the ratio of CS₂ to anions remains the same. The [C₄C₁im][NTf₂] solutions were simulated in a cell with edge length = 3.9 nm. This gives number densities of 0.202 CS₂ molecules/nm³ and 1.82 ion pairs/nm³. Simulations of dication solutions used a number of different cell sizes. Runs were performed with time steps of 1 fs. Several independent runs were carried out in the NVT ensemble at 300 K using a Berendsen thermostat. The amplitude of fluctuations in pressure in such a constant volume simulation is about 0.5 kbar. The electrostatic energy was calculated using an Ewald sum with precision 10⁻⁵. This corresponds to the Ewald parameter α = 0.2247 Å⁻¹ and a maximum k vector equal to (8 8 8). The stability of the Ewald sum was checked by monitoring the Coulomb energy and Coulomb virial, which remained constant and equal and opposite within the required precision. The real space cutoff for both the Ewald and the Lennard-Jones terms was chosen to be 1.25 nm. Configurations were saved to trajectory files at intervals for further analysis after initial equilibration.

In this paper, we give results for three systems at 300 K. The number densities, temperatures, and average pressures are given in Table I.

The motion of molecules in liquids can be studied directly in a molecular dynamics simulation, through measurements of molecular diffusion and determination of the density of states of intermolecular vibrations. Experimentally intermolecular motion can be probed using dielectric spectroscopy, low frequency Raman spectroscopy, or OHD-RIKES (Kerr effect spectroscopy). The Kerr spectrum is a collective property of the solution, related to the fluctuations in the anisotropic polarisability tensor for the liquid. As the anisotropic polarisability of individual CS₂ molecules is large, the contribution of CS₂ is an important part of the RIKES spectrum. Xue et al.¹ introduced an additivity model to determine this contribution. The RIKES spectrum K(ω) is determined by the fluctuations in the collective anisotropic polarisability of the liquid,

where Π_ij is the i, j component of the polarisability of the sample, $F$ is the Fourier transform at frequency ω, and a is a constant, a = 1/(ℏ²k_BT). Although it is possible to obtain the RIKES spectra from simulations,^19–22 it is expensive as a large matrix must be diagonalised at each step and, being a collective property rather than a molecular property, the simulation must be run for $N$ times as long to obtain the same signal to noise ratio as for a single molecule property. Here N is the number of molecules in the molecular dynamics cell. On the other hand, it is straightforward to calculate vibrational densities of states for individual molecules for specific normal modes. The Kerr spectrum of CS₂ consists of molecular terms and interaction induced terms. The former are due to fluctuations in the sum of molecular polarisabilities. As the principal axes of the molecular polarisability of each CS₂ molecule are along and perpendicular to its molecular axis, the fluctuations in the sum of polarisabilities are due to reorientations (librational motion) of individual molecules. The collision-induced contribution, on the other hand, is affected by the motion of CS₂ molecules relative to their surrounding molecules and is affected by both librational and translational motions of the molecules. In a dilute solution, collective contributions are likely to be small; thus the dominant contribution to the CS₂ RIKES spectra is due to librations. The librational and translational vibrational density of states for CS₂ is given by

where $vlibi$ is the difference of the perpendicular velocities of the two sulfur atoms in the ith molecule and $vCi$ is the velocity of the C atom in the ith CS₂ molecule. A and B are constants. Note that as the CS₂ molecules are rigid the motion of the centre of mass and of the central C atom is identical. In solution, the dominant contribution to the CS₂ Kerr spectrum is from the librational motion, which can be correlated with the CS₂ librational density of states in the low frequency librational spectrum. We present results for the latter quantity.

Molecular diffusion constants were calculated from the slopes of mean square displacements of the C sites in CS₂, the ring centres in the imidazolium cations, and the central N sites in the anions as a function of time. Radial distribution functions, energies, and pressures were calculated on the fly in the DL_POLY program.

Figure 2 shows the vibrational densities of states for CS₂ librations and translations in dication and monocation solutions at approximately 1 atm. The maxima are at a higher frequency in the dication solution, suggesting that the motion of the CS₂ is more restricted. Using the model of Xue et al.,¹ Gurung et al.³ found that the observed first moment of the CS₂ contribution in the Kerr spectra increases from 34.8 cm⁻¹ in a 10 mol. % [C₄C₁im][NTf₂] solution to 39.8 cm⁻¹ in the equivalent 18.2 mol. % $[C8(C1im)2][NTf2]2$ solution, while the peak position moved from 20.6 cm⁻¹ to 30.3 cm⁻¹. In the simulations, the first moment of the CS₂ librational density of states increased from 44.3 cm⁻¹ to 47.3 cm⁻¹ and the maxima moved from 29.1 cm⁻¹ to 33.5 cm⁻¹ in the corresponding solutions. The first moment of a spectrum K(ω) is defined by

where the integrals are over the frequency range of the intermolecular motion and the ω_max is chosen to be low enough to exclude intramolecular vibrations. Although the simulation values are of the vibrational density of states rather than the Kerr spectra, the trends are expected to be similar in solutions.

Figure 3 shows the solvent contributions. The density of states of the anion is very similar in the two liquids, but there is a significant change in the density of states of the imidazolium ring atoms. In the monocationic liquid, this is broad with a rather flat top, while in the dicationic liquid there is a broad maximum at higher frequency. We note that as the polarizability anisotropy is mainly associated with the imidazolium rings rather than the alkyl chains, this is the correct property to compare with the Kerr spectrum. In the Kerr spectrum, Gurung et al.³ found that the part associated with the cation is broad for both liquids with a weak maximum at higher frequency for the dicationic liquid.

The maxima in these vibrational densities of states are sensitive to pressure. Figure 4 shows how the first moments [defined in Eq. (4)] of the CS₂ translational and librational densities of states of the dication solutions vary as the pressure is varied (above) or as the number density is varied (below). The first moments in the [C₄C₁im][NTf₂] solutions are equivalent to those in a $[C8(C1im)2][NTf2]2$ solution with a number density of 1.7 rings/nm³ or p = −0.3 kbar. This suggests that the confinement of the CS₂ molecules would be similar in these solutions and that the observed shift in the Kerr spectra is due to increased confinement in the dication solution at the same state point.

Table II shows the average Lennard-Jones and Coulomb contributions to the energies of interaction of solute and solvent ions with their surroundings (embedding energies) for the 3 systems. The first two rows refer to the dication and monocation systems at approximately the ambient pressure and 300 K, which correspond to the experimental conditions. It seems that the magnitude of both contributions to the interaction energy of the CS₂ with the solvent increases by about 8%-9% on going from the monocation to the dication liquid, which explains the increase in the frequencies of the maxima of the librational and translational densities of states. There is a smaller increase in the anion interaction energies. The bottom row of the table gives more insight into this as it refers to a system (dicat2) with the same number density of rings as the monocation system. The [C₄C₁im]⁺ system can be thought of as the consequence of cutting all the dications in system dicat2 in the centres of their alkyl links, giving more crowding and a higher pressure, but greater entropy. The interaction energy of CS₂ and its surroundings follows the trend seen in the first moments of the librational spectra, being greatest (most negative) for the dication solution at the same pressure as the monocation solution and least (least negative) for the dication solution at the same number density as the monocation solution. Note that the Lennard-Jones contribution to the solvation energy of CS₂ is much greater than the Coulomb (electrostatic) energy, while the interaction of individual ions with their surroundings is dominated by the electrostatic energy.

In order to understand the effects of pressure and density on a molecular scale, we can look at the local structure in the liquids. Figure 5(a) shows the radial distribution functions between the central C site of CS₂ and the unique carbon site on the ring (C₂) in three examples, namely, in the dicationic solution (dicat) at ambient pressure (dark blue line), at low density (dicat2–cyan line) and in the monocationic solution (red dashed line). Part (b) shows the corresponding changes in the g(r) between CS₂ and the central anion site (N). These show clearly the increasing confinement in the series of solutions: ambient dicationic, monocationic, and high pressure dicationic, confirming that it is this increase which causes the shifts in the vibrational density of states and the observed Kerr spectrum.

There is also a significant difference in the radial distribution functions between CS₂ and the sites in the alkyl chains. This is shown in Fig. 6. The CS₂ molecules are more likely to be close to the terminal methyl group of the butyl group than to any of the sites in the linking alkyl chain in the dication.

Another probe of molecular motion in these solutions is provided by molecular diffusion constants. While the spectroscopic and energetic results are mainly due to nearest neighbour interactions, the diffusion constant is measured over a longer time and hence a larger length scale. We measured the diffusion constants of the solutes, cations and anions, in the solutions, and the results are shown in Fig. 7. The diffusion constant of the CS₂ molecule (which interacts less strongly with the surroundings, see Table II) is an order of magnitude higher than those for the ions. Unlike the shifts in the vibrational density of states shown in Fig. 2 and the changes in the local radial distribution functions shown in Fig. 5, the diffusion constant for CS₂ is slightly higher in the monocation solution than those in both dication solutions. One may think of the small non-polar CS₂ as sneaking through the more slowly moving anions and cations and being less impeded when the imidazolium rings are not tied together. The diffusion constants for the anions are lower than those for the imidazolium rings in every case, whether the rings are tied together as in the dications or free as in the [C₄C₁im]⁺ ionic liquids. Comparing dication and monocation results at constant pressure (dicat vs [C₄C₁im]⁺), one finds that the ion diffusion constants are faster by a factor of about 3 for the ions and 2 for the CS₂ when the dication is split into a pair of [C₄C₁im]⁺ ions. However at constant imidazolium density (dicat2 vs [bmim]), the increase is less (a factor of 1.5 for the ions and only 2% for CS₂).

The aim of this work was to compare molecular motion in the solvent [C₄C₁im][NTf₂] and the corresponding dication solvent $[C8(C1im)2][NTf2]2$ using CS₂ as a probe and atomistic simulation as the technique. In particular, we wished to understand the observed differences in the Kerr spectra. The main contribution to the low-frequency Kerr spectrum is due to the librational motion of the dissolved CS₂ molecules. These simulations show that these spectra are sensitive to the number density of the imidazolium rings; if this is the same for the mono and dication solutions, the first moment of the CS₂ librational motion is lower in the dication solution, while if the pressure is the same (and the number density higher in the dication solution) then the first moment is greater as in the observed RIKES spectra. We deduce that the environment of the CS₂ molecules is more restricted in the dication solution than in the monocation solution at the same pressure, but less so at the same number density. Radial distribution functions show that the distances between CS₂ probe molecules and ionic liquid cations and anions in the first shell decrease in the dicationic solutions with increasing density, while the environment in the monocation solution is intermediate between these. The changes in the solute-solvent interaction energies follow the same trend. By contrast, the CS₂ diffusion constants are larger in the monocation solution than in either of the dication solutions (the same density or the same pressure). This difference can be attributed to the fact that the CS₂ density of states and the differences in the radial distribution functions are primarily due to interactions with the first solvation shell while the diffusion constants are measured on a longer length scale.

See supplementary material for a list of the force field parameters for the dication solutions. The same charges and Lennard-Jones parameters were used for the monocation solutions.

The authors thank Dr. Eshan Gurung for providing spectral data before publication. E.L.Q. acknowledges support from the National Science Foundation, Grant No. CHE-1153077.