Heterogeneous dynamics of ionic liquids: A four-point time correlation function approach

It is no surprise that the use of ionic liquids (ILs)—salts that exist in a molten state below 100 °C—has been tested and shows promise in a broad range of applications.^1–5 Unusual and often desirable properties, such as chemical and thermal stability, high effective polarity, nonvolatility, wide electrochemical window, etc., make them ideal for a variety of specialized tasks. Furthermore, by judiciously engineering functional groups or exchanging ions entirely, the possibility arises for generating customizable ILs.

Parallel to a material’s applicability comes interest in the underlying molecular behavior. Many experimental^6–23 and computational^23–40 studies have reported the dynamic behavior of different ILs, revealing an array of interesting features, indicating their dynamics are heterogeneous. The examples include spatially correlated dynamics, sub-diffusive translational motions, non-diffusive reorientational dynamics, non-exponential solvent relaxation, etc. These dynamic behaviors, except for spatially correlated dynamics, have been analyzed using theoretical approaches based on two-point time correlation functions, such as van Hove functions.

In this paper, we extend previous computational efforts to gain further insight into heterogeneity of IL dynamics by using a time-dependent four-point density correlation function.⁴¹ This approach enables, for example, the investigation of the cooperative IL dynamics and the related dynamic length scale for spatial correlation, not accessible via two-point time correlation functions. While four-point time correlation functions have been used to characterize the dynamic behavior of structurally correlated species in glass-like systems,^41–50 their application to ILs has been extremely limited.⁵¹ 

Here, we study the time evolution of correlated motions in choline acetate ([Cho][Ac]) and 1-butyl-3-methylimidazolium acetate ([BMI][Ac]) (Fig. 1) as a protype of hydroxylic and non-hydroxylic ILs, respectively. For comparison, acetonitrile is considered as a standard nonglassy solvent.

Details of the simulation systems are compiled in Table I. Three trajectories of each IL system listed there and one trajectory of acetonitrile were simulated with the GROMACS simulation package.^52–55 As in Ref. 38, we used CHARMM force field parameters for ILs, generated via the ParamChem v.0.9.7 interface⁵⁸ that automatically assigns atom types, force field parameters, and charges on molecules based on the 2b8 release of the CHARMM general force field.^59–61 (See supplementary material for the complete list of force field parameters for BMI employed in the present study.) We note that the density of [BMI][Ac] determined with these parameters shows a very good agreement with the available experimental data (Table I). For acetonitrile, the potential model developed in Ref. 62 was employed. Electrostatic interactions were computed using the particle mesh Ewald (PME) method with a 1.4 nm cutoff. Nosé-Hoover temperature coupling and Parrinello-Rahman pressure coupling were used to maintain constant temperature and pressure of the system, respectively. For [Cho][Ac], initial configurations for each simulation were generated by annealing the system in the isothermal-isobaric (NPT) ensemble from 700 K and then equilibrating at production temperatures for 6 ns for 400 K and 350 K simulations and 2 ns for 600 K in the canonical (NVT) ensemble. The initial structures of [BMI][Ac] were prepared by annealing the system from 800 K in the NPT ensemble, followed by equilibration at production temperatures for 10 ns. The acetonitrile simulation was preceded by 0.2 ns of equilibration in the NVT ensemble at production temperatures. All production runs were performed in the NVT ensemble.

We start by considering the static structure of the system in the form of the radial distribution function (RDF), g_AB(r), i.e.,

where i and j represent individual particles and A and B denote species. The MD results for RDFs are displayed at 400 K for [Cho][Ac] and [BMI][Ac] and at 300 K for acetonitrile in Fig. 2. Higher temperature RDFs resemble those displayed in Fig. 2 but may show minor variations in peak positions and dimensions. For reference, the position (r′) and the height (g(r′)) of the first peak in each RDF are displayed in Table II.

The ILs show a prominent oscillatory behavior beyond 1 nm. This is generally attributed to an alternating-charge shell structure, i.e., cations tend to be surrounded by anions which, in turn, tend to be surrounded by cations.⁶³ This oscillatory behavior can be found in acetonitrile to a slightly reduced degree—the amplitude of oscillations is smaller and they extend to shorter distances.

To obtain insight into the anisotropic character of relative distributions of ions, we considered three-dimensional spatial distribution functions (SDFs). The results for cation distributions around the acetate anion are presented as isodensity surfaces in Fig. 3 as an example. The isodensity values are chosen to highlight the first solvation shells.

The SDF results indicate that cations—more precisely, their polar groups—are distributed preferentially around the C–O bonds of acetate. Though not shown here, the remaining solvation shell in the direction of the Ac methyl group is filled mainly by nonpolar groups of ions, viz., the butyl chain of BMI cations and the methyl group of other Ac anions. It is worthwhile to note that the absence of long aliphatic chains in choline reduces its amphiphilic nature considerably, compared to, e.g., BMI. As a result, the nonpolar hemisphere of the first solvation shell of acetate in [Cho][Ac] is filled predominantly by the methyl group of other Ac ions.³⁸ The double peaks in Fig. 2(b) and the corresponding double shells in Fig. 3(b) are due to the interactions of C^BMI of BMI cations with Ac anions in different directions, as shown in Fig. 4.

A time-dependent order parameter to characterize the system “glassiness” for species A can be defined as

As discussed in Ref. 44, local vibrational dynamics, if present, can relax Q^p(t) rapidly. Since our main focus is on correlations associated with systemic structural deformations as opposed to weakly correlated localized vibrations, we remove the latter by employing a coarse-graining function,^42–44 

in place of the delta function in Eq. (2). This gives the species an effective volume of $V0=43πa3$ and leads to

For IL systems, we employed a = 0.4 σ. For acetonitrile, a values for its carbon and nitrogen atoms, C^AN and N^AN, were 0.5 σ and 0.6 σ, respectively. Numerical values of a used for our analysis here can be found in Table III. These a values were selected to optimize the susceptibility associated with Q(t) (see below).⁴⁴ It is important to note that Q(t) does not distinguish between particles of like species. This means that there is no contribution to the decay of Q(t) when, for example, two particles of species A exchange positions. Because of this lack of distinction between particles of like species and the finite values of a, we find that $Q(∞)/N=N0$⁠, where N₀ ≡ V₀n and n is the number density of the species; n values are listed in Table I, and N₀ values are listed in Table III.

The MD results for $Q(t)/N$ are displayed in Fig. 5. The relaxation of $Q(t)/N$ was found to be well fitted by the sum of multiple exponential decays, i.e.,

Values for a, N₀, $τiQ$⁠, and $AiQ$ can be found in Table III. [Cho][Ac] and [BMI][Ac] were fitted with three and four exponential functions, respectively. By contrast, a single exponential function well describes $Q(t)/N$ of acetonitrile, indicating that underlying dynamics are Gaussian. As is expected, $Q(t)/N$ decays faster at higher temperatures; all values of $τiQ′$s are smaller and the contribution of the fastest decay component, viz., $A1Q$⁠, to the overall relaxation increases. At 350 K and 400 K, $Q(t)/N$ for [Cho][Ac]—the more viscous of the two ILs—decays more slowly than $Q(t)/N$ for [BMI][Ac]; at all simulated temperatures, both ILs decay much more slowly than acetonitrile.

To understand the time evolution of the dynamic correlation between species A and B, we analyzed the four-point overlap correlation function $g4ol(r,t)$ between the two species

For simplicity, we consider a radially averaged correlation function, $g4ol(r,t)$⁠, where $r=r→$⁠. We note that $g4ol(r,t)$ is a direct extension of the full van Hove function; it describes the time evolution of correlation of density fluctuations occurring simultaneously at two different positions, located a distance of r apart. As such, it gauges the four-particle conditional probability, viz., given a particle of species A at position 1 and particle of species B at position 2 at the initial time t = 0, it measures (up to a constant) the probability that the same or another particle of A and the same or another particle of B will be found at positions very close to the original locations 1 and 2, respectively, at a later time. Like the order parameter described above, $g4ol(r,t)$ does not distinguish between particles of the same species. This means that at time t = 0, $g4ol(r,0)$ becomes the RDF, i.e., $g4ol(r,0)=g(r)$⁠. We also note that $g4ol(∞,t)=QA(t)/NAQB(t)/NB$ because particles are of finite size and become spatially uncorrelated at large r. In other words, $g4ol(r′,t)$ does not decay to zero due to the finite value of a.

Analogous to the fitting for $Q(t)/N$ in Eq. (5), we were able to fit $g4ol(r,t)$ for ILs accurately with three exponential functions as

We found that $g4ol(r,t)$ for acetonitrile is well described by a single exponential function. As mentioned above, b(r) in Eq. (7) becomes $QA(∞)/NAQB(∞)/NB(=NA,0NB,0)$ at large r. Values for $τiol$⁠, $Aiol$⁠, and b for r = r′, i.e., positions of primary RDF peaks, are presented in Table IV (see Table II for values of r′ and g(r′)). The corresponding relaxation behaviors of $g4ol(r′,t)$ for the systems considered here are compared in Fig. 6. At a given temperature, [BMI][Ac] shows a faster decay of $g4ol(r,t)$ than [Cho][Ac] just like the $Q(t)/N$ case. The $g4ol(r,t)$ relaxation in acetonitrile at 300 K is faster than that in [BMI][Ac] and [Cho][Ac] at 350 K by ∼4 orders of magnitude.

For further insight, we have examined the r-dependence of $τiol(r)$ and their normalized contributions $Aiol(r)/∑i=13Aiol(r)$⁠. The results for [Cho][Ac] at 400 K are shown in Fig. 7. The different $τiol(r)$ values correspond to different relaxation time-domains of $g4ol(r,t)$⁠. The fastest decay time, $τ1ol(r)$⁠, likely due to initial motions of ions within their cages, exhibits essentially no dependence on r [Fig. 7(a)]. We think that such motions are fairly localized and thus rapid density fluctuations arising from them are not affected by the motions of other ions. The slower decay times, especially $τ3ol(r)$⁠, on the other hand, show considerable r-dependence for r ≲ 1.3 nm. This means that slow density fluctuations of ions are likely influenced by and therefore correlated with the density fluctuations of other ions located within a distance of ∼1.3 nm. The results in Fig. 7(b) show that all three decay components make a significant contribution to the overall relaxation of $g4ol(r,t)$⁠. There are several noteworthy features. First, the normalized amplitudes of the two slow components, $A2,3ol(r)/∑i=13Aiol(r)$⁠, are strongly correlated while they show near perfect anticorrleation with the normalized amplitude of the fastest component. Second, the r-behavior of the amplitudes of the slow components is very similar to the corresponding RDF (cf. Fig. 2). We note that [BMI][Ac] shows behaviors similar to [Cho][Ac] and thus is not shown here.

The susceptibility χ₄(t), which describes the fluctuations of the order parameter of the system, is defined as

with β = 1/k_BT. χ₄(t) measures cooperative motions of ions at two different locations, averaged over the system. Figure 8 displays χ₄(t) for C^AC(C)⋯C^AC(C) at 400 K determined with different values of a. As mentioned above where a was first introduced, the a value with the greatest susceptibility is used for calculations throughout this work.⁴⁴ 

Figure 9 shows the MD results for the susceptibility for each system at the selected a value (Table III). As is expected from the temporal behaviors of both the order parameters and the four-point density time correlation functions, we find χ₄(t) peaks at longer times for lower temperatures. Furthermore, χ₄(t) for more viscous [Cho][Ac] is slower than that for [BMI][Cl]. We ascribe the molecular origin of this difference to the formation of hydrogen bond in the former.³⁸ Interestingly, the susceptibility calculated between like species is greater than the susceptibility calculated between unlike species.

In many prior studies,^{25,29,30,38,50,51,64,65} the deviation of the IL transport behavior from Gaussian dynamics is measured by

where Δr(t) is the displacement of ions at time t from their initial positions. Comparison of the results in Fig. 10 with those in Fig. 9 reveals that temporal behaviors of α₂(t) and χ₄(t), including their peak positions, are very similar for all temperatures and systems we considered. This indicates that the non-Gaussian behavior of single particle dynamics is directly related to time-dependent correlated motions of two particles, viz., cooperative dynamics.

Finally, we proceed to the dynamic length scale associated with cooperative motions in ILs. We employ the approach widely used in the analysis of glass systems,^44–48,50 viz., a small wave vector expansion of the four-point structure factor

The inclusion up to the q² terms in the expansion à la Ornstein-Zernike (OZ) theory yields^44,45

with α = 2. The corresponding r-space behavior is given by

While a more general expression with α ≠ 2 is often used in the analysis of glass systems,^46–48,50 we employ here the original OZ form,^44,45 i.e., α = 2 in Eq. (11), mainly due to the limited number of data points available for the numerical analysis. This limitation arises from the small system size we used to speed up the simulations because long trajectories are needed for proper sampling of highly viscous and thus slowly relaxing ILs. We computed the radially averaged structure factor, $S4ol(q,t)$ (⁠$q=|q→|$⁠), by including in the calculations only $q→$ vectors that satisfy the periodic boundary conditions of the simulation system and fitted the resulting $S4ol(q,t)$ with Eq. (11) (α = 2) at the lowest five q values using the Newton method with $S4ol(0,t)$ and ξ₄(t) treated as fitting parameters.

The results for ξ₄(t) for like ions (cation⋯cation and anion⋯anion) are shown in Fig. 11. As is expected, increasing temperature leads to shorter correlation lengths, i.e., lower peak heights of ξ₄(t). This means that the spatial range of cooperative dynamics becomes shorter as the temperature rises. Another important feature is that the spatial range of cooperative dynamics varies with time. This can be understood as follows: A dynamic event occurring at, say, position 1 cannot be influenced by another dynamic event at a different location, position 2, if the two events occur at the same time because the effect of a dynamic event cannot propagate instantaneously. However, a dynamic event occurring at a later time at position 1 can be influenced by a dynamic event at an earlier time at position 2. We expect that with time, the propagation and strength of the effect of a dynamic event grows and diminishes, respectively. As a result, the corresponding spatial range of the influence of a dynamic event will exhibit a non-monotonic behavior in t, reaching a maximum at some point, followed by a decrease. χ₄(t) measures this time-dependence via Eq. (12). Comparison with Fig. 9 shows that the peak positions of ξ₄(t) are close to those of χ₄(t).

The correlation length scale can be compared with the sizes of ions shown in RDFs, Fig. 2. For [Cho][Ac], the maximum ξ₄(t) at 600 K is smaller than the ionic radii. As the temperature decreases, the correlation length approaches and goes slightly above the ionic radii for both [Cho][Ac] and [BMI][Ac]. The length scale of acetonitrile is much less than its molecular size, indicating no correlated motions. It also shows that the correlation length depends on the dynamics across the system rather than the size of the species.

It is worthwhile to note that the temperatures of our simulations are well above the reported glass transition temperatures of both systems. We expect that if simulations were performed at lower temperatures, the length scale would further increase, which implies that both systems would show a stronger glassy behavior at room temperature.

Spatially correlated dynamics in ILs are brought about, in part, by distinct spatial clustering of “mobile” and “immobile” ions. Studies of glass systems with similar behavior implement a four-point time-correlation function to simultaneously gauge the length-scale of correlated regions and the time associated with the correlation.

In this paper, we studied the time-dependent four-point density correlation function (⁠$g4ol(r,t)$⁠), its related susceptibility (χ₄(t)), and the dynamic length scale of correlation (ξ₄(t)) for the ILs [Cho][Ac] and [BMI][Ac]. We found that $g4ol(r,t)$ exhibits r-dependent relaxation, revealing the influence of spatial correlation on collective motions of ions. ξ₄(t) determined by fitting the four-point structure factor $S4ol(q,t)$ in the low q region shows that the spatial range of collective motions initially increases with time but decreases after reaching a maximum. Compared with the traditional solvent, acetonitrile, correlation extends to greater distances in the studied ILs and is retained for significantly longer times.

Our results indicate that four-point time correlation functions can shed new light on IL dynamics, not accessible via two-point correlation functions, and thus provide a useful theoretical tool to study IL dynamics, especially their cooperative aspect. Therefore it would be worthwhile in the future to study structure-dynamics relationships in ILs, such as the influence of structural heterogeneity on collective dynamics, using the four-point time correlation functions and their variants.

See supplementary material for the force field parameters of BMI employed in the present study.

This work was supported in part by the National Science Foundation through NSF Grant No. CHE-1223988. H.J.K. thanks Professor Biman Bagchi and Professor Shinji Saito for suggesting four-point time correlation function analysis of ionic liquid dynamics.