J. Phys. Chem. B 2007, 111, 4705-4714
4705
On the Response of an Ionic Liquid to External Perturbations and the Calculation of Shear Viscosity† Zhonghan Hu and Claudio J. Margulis* Department of Chemistry, The UniVersity of Iowa, Iowa City, Iowa 52242 ReceiVed: October 27, 2006; In Final Form: December 11, 2006
In this article, attention is directed to three related problems: (1) the response of the ionic liquid (IL) 1-hexyl3-methylimidazolium chloride ([HMIM+][Cl-]) to different external perturbations, (2) the calculation of its shear viscosity, and (3) the investigation of the range of validity of linear response theory for these types of systems. For this purpose, we derive a set of equations linking bulk hydrodynamic predictions with microscopic simulations which are valid when linear response theory is applicable. As far as we are aware, this article reports results from the largest atomistic simulations ever performed on this liquid. Our study shows that even for systems with a box length as large as 0.03µ the viscosities computed from perturbation frequencies compatible with this box size have not yet reached the bulk hydrodynamic limit. This is in sharp contrast with the case of other solvents such as water in which the hydrodynamic limit can be achieved by using perturbations on a length scale of typical molecular dynamics simulation box sizes. In order to achieve our goals, we comprehensively investigated how the IL relaxed upon weak external perturbations at different wavenumbers. We also studied the steady-state flow created by external shear acceleration fields. The short time behavior of instantaneous velocity profiles was compared with the results of linear response theory. The short time response appears to match the prediction from linear response theory, while the long time response deviates as the external field becomes stronger. From this study, the range on which a perturbation can be considered “weak” in the linear response sense can be established. The relaxation of initial velocity profiles was also examined and correlated to the decay of the transverse-current autocorrelation function. Even though none of our calculations reached the bulk hydrodynamic limit, we are able to make predictions for the shear viscosity of the bulk system at different temperatures which qualitatively agree with experimental data.
1. Introduction The viscosity of ionic liquids (ILs) has demanded a great deal of interest in recent years, since this property determines the usability of these compounds for particular applications. Several groups have measured the coefficients of shear viscosity of different ILs.1-10 In general, the reported viscosities are 2-4 orders of magnitude larger than that of room temperature water. One of the main issues with experiments is that many ILs contain impurities such as chloride and water. These impurities cause experimental inaccuracy, since viscosity is known to significantly change in the presence of these contaminants. Computer simulations are also not without problems, since different models may produce slightly different liquid properties. A few theoretical and computational studies have already addressed in some detail the problem of viscosity in ILs.11-13 However, a complete microscopic theory of viscosity is currently not available. It is a challenging task to accurately compute the viscosity of a complex system by means of simulation methods. For a system with high viscosity, it is extremely difficult to reach the hydrodynamic limit (zero wavenumber) where the experimental data is observed. This is because, in order to reach this limit, a very large simulation box is required. Traditional simulation methods to calculate the shear viscosity of fluids fall into two categories: (1) the evaluation of the transverse-current autocorrelation function (TCAC) through †
Part of the special issue “Physical Chemistry of Ionic Liquids”. * To whom correspondence should be addressed.
equilibrium molecular dynamics (EMD) trajectories and (2) nonequilibrium molecular dynamics (NEMD) simulations that impose a periodic perturbation or Couette flow.14-21 In an interesting recent article, Hess compared most of the above methods by performing simulations of Lennard-Jones and water system.20 They concluded that the NEMD method using a periodic shear perturbation was the best choice. Recently, another novel method called momentum impulse relaxation (MIR) was developed by Maginn and co-workers.22,23 Theoretical similarities and differences between these methods based on a microscopic theory of viscosity will be the focus of a future publication in preparation.24 In general, most computational methods suffer from the fact that they depend on a wavenumber which needs to be small in order to reach the hydrodynamic limit in which the viscosity is a number, not a function. This limit is hard to reach, since the smallest possible wavenumber scales like the inverse of the simulation box size which is microscopic. For many molecular systems, this is not a problem, but this is not the case for ILs. In order to gauge the effect of size on the response of ILs upon perturbation, we have performed different simulations of the 1-hexyl-3-methylimidazolium chloride ([HMIM+][Cl-]) system varying from a very large box containing 8232 pairs of ions to a traditional simulation box of 343 pairs of ions. We have used equilibrium simulations in order to compute TCACs at different wavenumbers. Within the nonequilibrium framework, we have examined the case of externally imposed steady-state periodic drags and also the relaxation of systems
10.1021/jp067076n CCC: $37.00 © 2007 American Chemical Society Published on Web 01/19/2007
4706 J. Phys. Chem. B, Vol. 111, No. 18, 2007
Hu and Margulis
prepared with an initial momentum profile. This last approach is related to the MIR method. We have carefully examined computationally and theoretically the relationship between equilibrium and nonequilibrium methods. In this way, we present a complete description of the response of [HMIM+][Cl-] to external shear perturbations. From our equilibrium TCACs, we have obtained information about the spontaneous fluctuation of these systems at relatively low frequencies. Since we also have corresponding data from nonequilibrium simulations, we have been able to check the validity of linear response theory and the fluctuation-dissipation relationships for nonequilibrium perturbations of different amplitudes. From this, we have gained an understanding on what constitutes a “weak” perturbation in the linear response sense for this liquid. The rest of this paper is organized as follows. In section 2, we provide needed theoretical background both on hydrodynamic equations and on linear response theory. A detailed description of current available methods based on a microscopic theory will be published elsewhere.24 However, this brief section makes this work self-contained. The simulation details are provided in section 3. Section 4 presents and discusses the results of theoretical and computational analysis. In section 5, we draw conclusions from our present work. 2. Theoretical Background 2.1. Hydrodynamic Equations. In order to derive a formula for shear viscosity, it is sufficient to consider only the transverse part of the Navier-Stokes equation:25
∂ux(z,t) ∂2ux(z,t) F ) Fax(z,t) + η ∂t ∂z2
(1)
where we define the z direction as the longitudinal direction, ux(z,t) and ax(z,t) are the transverse velocity field and external acceleration, respectively, F is the mass density, and η is the coefficient of shear viscosity. This macroscopic equation is only valid at long time scales and large length scales (small wavenumbers). We are going to use this equation as the starting point to discuss three different cases: spontaneous fluctuation (SF), periodic perturbation (PP), and initial pulse perturbation (IPP). 2.1.1. Spontaneous Fluctuation. Here, we consider an equilibrium case where the external force is zero in eq 1 and define the spatial Fourier transform of the dynamical variable, u(r), as
u(k,t) )
∫-∞+∞u(r,t)eik‚r dr
(
)
(
and directly integrate both sides of eq 5 to obtain
C(k,ω)0) ) C(k,t)0)
∫0+∞exp(- ηFk2t) dt ) ηkF 2
(7)
The left-hand side of the above equation is just the area under the normalized time correlation function, C(k,t)/C(k,t)0). The coefficient of shear viscosity, η, can be calculated from this normalized area or the zero frequency value of the spectra (eq 6). This method has been used by Balucani et al.19 for a water system and Urahata and Ribeiro12 for ionic liquid systems. 2.1.2. Periodic Perturbation. In the case of PP, we assume an initial velocity and an acceleration profile as
{
ux(z,t)0) ) 0 ax(z,t) ) a0 cos(kz)
(8)
The solution to eq 1 becomes
ux(z,t) ) a0τ(1 - e-t/τ) cos(kz) ) u0 cos(kz)
(9)
where the relaxation time, τ, is defined as τ ) F/ηk2. By fitting the steady-state velocity profile at sufficiently long time (e-t/τ f 0) to a cos(kz) form and measuring the amplitude, u0, one is able to get the coefficient of shear viscosity as the following:
η)
a0F u0k2
(10)
This method called the periodic perturbation (PP) method was first developed by Gosling et al.15 Instead of imposing a cos(kz) shape acceleration, Backer and co-workers have used a step function21
{
a0 if z > 0 ax(z,t) ) -a 0 if z < 0
(11)
Correspondingly, the resulting velocity profile has to be fitted to a parabolic function instead of a cos(kz) function. Without any loss of generality, in this work, we only focus on cos(kz) shape accelerations. 2.1.3. Initial Pulse Perturbation. In the case of IPP, an initial velocity profile of the form
(3)
ux(z,t) ) w0 e-t/τ cos(kz) ) w0(t) cos(kz) (4)
)
(6)
is created by a pulse interaction and the external force is zero (a0 ) 0); the solution to eq 1 gives the decay of the velocity profile as a function of time:
evolves exponentially as25
η C(k,t) ) C(k,0) exp - k2t F
∫0+∞C(k,t)eiωt dt
ux(z,t) ) w0 cos(kz)
Therefore, the autocorrelation function
C(k,t) ) 〈u*(k,0) u(k,t)〉
C(k,ω) )
(2)
The equilibrium solution to eq 1 in k space is
η u(k,t) ) u(k,0) exp - k2t F
the Fourier transform of the time correlation function, C(t), as
(5)
where C(k,t) is only a function of the magnitude of the wavenumber vector k for isotropic systems. We further define
(12)
(13)
From this solution, we can see that the velocity profile will decay exponentially at large time scales and length scales and the relaxation time, τ ) F/ηk2, is inversely proportional to the coefficient of shear viscosity. Thus, shear viscosity can be determined by directly measuring this relaxation time. 2.2. Linear Response Theory. Although most of the simulation methods to compute viscosity were originally developed from hydrodynamic equations, it is insightful to understand them on the basis of a microscopic theory. This is because a real
Response of an IL to External Perturbations
J. Phys. Chem. B, Vol. 111, No. 18, 2007 4707
simulation usually uses a box size corresponding to finite wavenumbers (k * 0). We therefore start from linear response theory which is generally valid for weak perturbations at arbitrary wavenumbers. In general, for a system subject to an external field, Fe(t), the dynamics of the system obeys the following equation:26
R4 q ) vq + CqFe(t)
mqv3 q ) Fq + DqFe(t)
∫0tdτ 〈B(t-τ) J(t)0)〉Fe(τ)
(15)
where the dissipative flux, J(t), is defined as
∑q (-Dq‚vq + Cq‚Fq)
J(t) ≡
(16)
2.2.1. Periodic Perturbation and Spontaneous Fluctuation. For the case of PP, we consider again an external acceleration of the form of eq 8. The dynamical variable Cq is always zero here. The dynamical variable Dq and the external field can be taken as the following:
Dq(k,t) ) exmq cos(kzq)
Fe(t) ) a0h(t)
〈B(z,t)〉 )
Noticing that 〈B(z,t)〉lz/M corresponds to ux(z,t) in eq 9, replacing u0 with a0C(k,ω)0)/C(k,t)0), we have
η(k) )
{
Therefore, the dissipative flux is
(18)
The dynamical variable B we are interested in can be written as
B(z,t) ≡
∑q mqVqx(t)δ(z - zq(t))
(19)
By manipulating the integral in eq 15 (see section 6), the nonequilibrium ensemble average of the dynamical variable B(z,t) becomes
〈B(z,t)〉 )
C(k,τ)
∫0tdτ C(k,t)0)
M a cos(kz) lz 0
(20)
where M is the total mass and lz is the z direction length of the system. M/lz gives the length density of the system. The transverse-current autocorrelation function (TCAC) is defined as
C(k,t) ) 〈
∑q mqVqxe ∑p mpVpxe -ikz
ikz
〉
(21)
Note that we use momentum here instead of velocity, as was used in eq 4. Macroscopically, these two are only different by a constant. Equations 5 and 7 are still valid given the definition of eq 21, while it is only exact to use eq 21 rather than eq 4 for multiple component systems.
a0F 2
u0k
)
F C(k,t)0) k2 C(k,ω)0)
(23)
This equation is exactly consistent with eq 7. Therefore, eqs 7 and 10 are equivalent even in the case of finite wavenumbers (k * 0) as long as linear response theory applies. In order to look at the instantaneous response as a function of time only, it is straightforward to define an instantaneous velocity amplitude by eliminating the space dependence part, cos(kz), from both sides of eq 20 (see section 6)
V(t) )
∑q 2mqVqx(t) cos(kzq(t))/M
(24)
We therefore have (see section 6)
〈V(t)〉 ) a0
(17)
1 t>0 h(t) ) 0 t
