J. Phys. Chem. B 2001, 105, 6675-6683
6675
Stochastic Dynamics in Near-Critical Supercritical Fluids† Alexander N. Drozdov‡ and Susan C. Tucker* Department of Chemistry, UniVersity of California, DaVis, California 95616 ReceiVed: January 26, 2001
An extension of the generalized Langevin equation (GLE) is introduced in order to model diffusion processes in near-critical fluids. Unlike the standard GLE, this new approach allows a flexible description of the influence of slowly fluctuating local density distributions that occur in a fluid as its liquid-vapor critical point is approached. We find that, given an initial nonequilibrium ensemble, diffusion coefficients differing substantially from the expected bulk value can be observed at intermediate times, that is, at times longer than the velocity relaxation time but shorter than the local environment interconversion time.
I. Introduction The development of efficient technological processes such as supercritical separation and extraction1 has recently directed a great deal of theoretical and experimental attention to the transport properties of supercritical fluids (SCFs).2 There has also been a growing interest in using SCF solvents to promote reaction processes. Numerous investigations have demonstrated that chemistry in SCFs frequently yields unexpected reaction products, mechanisms, and rates.3-6 From the technological point of view, the most appealing feature of SCF solvents is that their density, F, can be varied continuously over extremely broad ranges simply by changing the temperature T or pressure p. Because the transport and solvating properties of a fluid depend strongly on its density, an ability to tune the solvent density translates into an ability to manipulate the efficiency of separation processes and reaction outcomes.4-6 However, one obstacle that hinders realization of this potentially efficient thermodynamic control is our lack of understanding (and resultant inability to predict) the behavior of SCF solvents in the vicinity of the critical point where the desired control over fluid properties is maximized. Product formation in chemical reactors is mainly controlled by elementary rate mechanisms such as diffusion and thermal activation, and it is therefore of considerable importance to be able to predict these processes in fluids under supercritical conditions. Herein, we shall focus on the simplest of the diffusion mechanisms, self-diffusion. Although molecular dynamics (MD) simulations7,8 have historically provided the most detailed information about diffusion processes in fluids, 9-11 the applicability of this method to study diffusion processes in nearcritical SCF solvents seems rather doubtful, due to the presence of critical phenomena. In particular, the long-length scale density inhomogeneities, which generate the divergence of the compressibility at the critical point, fluctuate only very slowly, such that the times required to adequately sample the equilibrium system may become quite long. As a result, the fact that the typical simulation times are of the order 10-10 s imposes a severe obstacle to applying the MD method in the critical region.12-14 Similarly, typical finite system sizes also pose a problem, as the divergence of the correlation length under near-critical †
Part of the special issue “Bruce Berne Festschrift”. Permanent address: Institute for High Temperatures, 13/19 Izhorskaya Street, 127412 Moscow, Russia. ‡
conditions clearly cannot be captured. This means that the simulation of dynamic processes near the critical point, while feasible on the fastest, parallel supercomputers, is not viable for routine investigations, and we therefore explore less computationally intensive approaches. Of the other theoretical approaches to modeling solute dynamics in liquids, the method that has been most successfully and extensively employed in the past is the theory of Brownian motion.15,16 The basic idea underlying this theory is the observation that the velocity of a (Brownian) solute particle moving through a fluid medium varies by a large number of small jumps that originate from the random collisions between the particle and the molecules of the medium. As a consequence, the Brownian particle experiences two forces: one is a velocity proportional friction force and the other is a random force, f(t). In many practical situations the autocorrelation time, τf, of the random force, f(t), i.e.,
τf )
∫0∞ dt
〈f(t)f(0)〉 〈f2(0)〉
(1.1)
is not vanishingly small compared to the time scale, τv, on which the velocity, x˘ (t), of the diffusing particle becomes uncorrelated,
τv )
∫0∞ dt Cv(t)
(1.2)
where the velocity autocorrelation function Cv(t) is defined by
Cv(t) )
〈x˘ (t)x˘ (0)〉 〈x˘ 2(0)〉
(1.3)
The fact that τf ≈ τv gives rise to correlations between successive changes of the velocity, or, in other words, to a retarded frictional force. The corresponding dynamics are then governed by the generalized Langevin equation (GLE) which, in the simplest one-dimensional case, has the form17
mx˘ ) -m
∫0t dt′ γ(t - t′)x˘ (t′) + f(t)
(1.4)
Here m is the mass of the particle, x the particle coordinate, and the dot denotes the time derivative. The frictional memory kernel γ(t) entering eq 1.4 represents the response of the solvent from the past to the present. This kernel is directly related to
10.1021/jp010354s CCC: $20.00 © 2001 American Chemical Society Published on Web 06/14/2001
6676 J. Phys. Chem. B, Vol. 105, No. 28, 2001
Drozdov and Tucker
the correlation function of the random force f(t) by the second fluctuation dissipation theorem reading
〈f(t)f(t′)〉 ) kBTmγ(t - t′)
(1.5)
where kBT is the Boltzmann constant multiplied by absolute temperature. Last, the static friction coefficient, η, which characterizes the overall magnitude of the solvent damping, is defined by the zero-frequency Laplace transform of γ(t)
η)
∫0∞ dt γ(t)
(1.6)
Both the friction kernel γ(t) and the static friction coefficient η are typically determined either from experiment18,19 or by analysis of full MD simulations.9-11 In theoretical calculations an exponential form for γ(t) is often assumed20
γ(t) ) (η/τf)e-t/τf
(1.7)
with τf being the correlation time of the random force, eq 1.1. Although the Brownian picture has been quite successful in describing diffusion in liquids,16 we find that near-critical SCFs pose yet an additional complication, which is not describable by even the GLE. This complication arises because the slowly fluctuating, long-length scale inhomogeneities generate a very broad distribution, P(Fl), of local solvent environments, ranging from very low to very high local densities, which interconvert only very slowly,13,21 on a time scale τF. We assume that the local environment around an atom can be well characterized by the local density, Fl, surrounding that atom. We therefore define an instantaneous local density around a tagged particle, Fl(t), as the number of solvent particles (excluding the tagged particle) falling within a pre-specified local radius rl (typically taken to be the range of 1 or 2 solvent shells) of the specified particle at time t, divided by the volume of this local region (see ref 21). To understand why this situation of a broad, but slowly interconverting distribution of local densities is not describable by the usual GLE, first recognize that the net solvent force acting on the “solute” (i.e., on the diffusing particle, be it a solute or simply a tagged solvent particle) at any time t, f(t), will be determined by the locally surrounding solvent molecules.22 In particular, the magnitude of f(t) will, on average, reflect the number of such local solvent molecules present at time t, and hence will be sensitive to the instantaneous local solvent density, Fl(t). Additionally, near the critical point, the time scale on which these local densities interconvert, τF, becomes very large, such that τF . τf, the time scale on which the solvent forces f(t) remain correlated. Hence, we see immediately that, in this limit, the force correlations, and thus γ(t), will have decayed to zero long before the local environments interconvert. Physically, then, for a given time period, the friction γ(t), which decays to zero on the time scale τf, will reflect only the current local density, Fl(t), such that γ(t) will be the friction expected in a homogeneous fluid of bulk density F ) Fl(t), i.e., γFl(t). Note that since γFl(t) decays to zero before Fl(t) changes, it can in no way reflect the solvent’s local density dynamics occurring on the time scale τF. Moreover, as is well known, the magnitude of γFl(t) will be strongly dependent on the density, Fl, such that over times τF, the function γ(t) appearing in the GLE ought to vary in magnitude. But, because the GLE assumes that γ(t) is a function independent of time origin, it does not allow for the presence of such an additional slow solvent time scale, τF, which exceeds the longest correlation time τf. Specifically, the GLE formalism would require us to use 〈γFl(t)〉, where the average is over the
distribution of local solvent densities, P (Fl), and would thus provide a mean field approximation to the true dynamics. That is, the GLE would require one to assume that the solute dynamics are controlled by an average solvent friction, even though the time scale over which the true ensemble of frictions would be sampled is very long, thus making the usual GLE a poor method for approximating the true dynamics. In section II of the present paper, we introduce a multipletime-scale generalized Langevin equation (MTSGLE), which enables the inclusion of the solvent dynamics occurring on the time scale τF, while, at the same time, retaining the solventforce-relaxation time scale τf. This formalism is seen to be a variant of the irreversible GLE (iGLE) recently introduced by Hernandez and co-workers23 for the description of nonstationary dissipative systems. In section III, we discuss the connection between the usual GLE and the diffusion coefficient, and, additionally, we note how these relations may be altered in an inhomogeneous fluid that must be described by an MTSGLE. In section IV, we present our results, the diffusion dynamics for both initial equilibrium and nonequilibrium ensembles of diffusing particles in this inhomogeneous medium, as predicted by the MTSGLE. Conclusions follow in section V. II. Multiple-Time-Scale Generalized Langevin Equation (MTSGLE) A. Formalism. Our approach begins by recognizing that at any time, t, the static friction coefficient, η, which will control the motion of the diffusing particle, will be strongly dependent on the local density, Fl(t), which surrounds the diffusing particle at that time.24 In contrast, the time dependence of the force correlations, τf, and the functional form of the friction kernel’s decay with time are relatively insensitive to variations in the density.24 We therefore follow Hernandez and co-workers23 and split the friction kernel γ(t,t′) into a time-dependent static friction contribution xη(t)η(t′), which will depend on the local densities Fl(t) and Fl(t′), and a stationary friction kernel, γ0(t - t′), which will be independent of Fl(t). This split then gives
γ(t, t′) ) xη(t)η(t′)γ0(t - t′)
(2.1)
where, by construction, γ0(t) is normalized to unity,
∫0∞ dt γ0(t) ) 1
(2.2)
and the static friction prefactor may vary on times τF even though γ0(t) decays on times τf , τF. With this division, one may write the stochastic equation of motion in a general form as
mx¨ ) -
∂U(x) -m ∂x
∫0t dt′ xη(t)η(t′)γ0(t - t′)x˘ (t′) + xkBTmη(t)f0(t)
(2.3)
where U(x) is a potential of mean force along a chosen intramolecular (solute) coordinate x, and, for the diffusion problem, this term is absent. In the present application, the normalized stationary friction kernel γ0(t) and the normalized random force f0(t) should be taken to be the same as those experienced by a diffusing particle in a homogeneous fluid at the given temperature (as noted above,24 these quantities are rather insensitive to density). Once γ0(t) has been determined, f0(t) may be generated according to the second fluctuation dissipation theorem,
Stochastic Dynamics in Near-Critical SCFs
J. Phys. Chem. B, Vol. 105, No. 28, 2001 6677
〈f0(t)f0(t′)〉 ) γ0(t - t′)
(2.4)
Notice that the slow dynamic properties of the near-critical solvent, occurring on the time scale τF, will enter eq 2.3 through the static friction η(t). Before proceeding with the formulation used to generate η(t), we note that using eq 2.1 for γ(t,t′) and introducing a new function f(t) ) xkBTη(t)f0(t) casts eq 2.3 into the standard irreversible form23
mx¨ ) -
∂U(x) -m ∂x
∫0t dt′ γ(t, t′)x˘ (t′) + f(t)
(2.5)
where, by construction,23
〈f(t)f(t′)〉 ) kBTmγ(t,t′)
(2.6)
which is a natural extension of the fluctuation-dissipation relation to a nonstationary or stochastic friction kernel. The latter property ensures that the equilibrium distribution will be of the standard Maxwell-Boltzmann form. What distinguishes the MTSGLE, which represents a stationary equilibrium system involving a slow reorganizational time scale, from the iGLE, which represents nonstationary processes, is the nature of the time-dependence of the static friction coefficient η(t). In the iGLE, η(t) is a nonstationary variable which thus depends explicitly on the absolute value of the time, whereas in the MTSGLE, η(t) is a stationary stochastic variable having fluctuations that remain correlated over a long time scale, here τ F. In the present problem, the fluctuations in the static friction reflect fluctuations in the local density around the diffusing particle. Hence, we choose to use the local density, Fl(t), as the underlying stochastic variable, and model η(t) as η[Fl(t)], where the dependence of η on Fl will be extracted from MD simulations10 (see section II.B). To model the correlated fluctuations of the stochastic variable Fl(t), we assume the standard form of an Ornstein-Uhlenbeck process for which the relevant Langevin equation reads
F˘l(t) ) -τ-1 F [Fl(t) - 〈Fl〉] + ξ(t)
(2.7)
where ξ(t) is a Gaussian white noise normalized to
〈ξ(t)〉 ) 0, 〈ξ(t)ξ(0)〉 ) (2φ2/τF)δ(t)
(2.9)
and has the standard Gaussian form
[Fl - Fm(t)]
Peq(Fl) ) (2πφ )
}
2
2ϑ2(t)
(2.10)
where
Fm(t) ) 〈Fl〉 + (F′l - 〈Fl〉) exp(-t/τF)
(2.11)
ϑ2(t) ) φ2[1 - exp(-2t/τF)]
(2.12)
and
Finally, the equilibrium probability, Peq(Fl) ) P(Fl, t f ∞|F′l ),
[
]
(Fl - 〈Fl〉)2
exp -
2φ2
(2.13)
thus defining the parameters φ and 〈Fl〉 as the width and mean of the equilibrium distribution of local densities, respectively. Note that, because we have taken Peq(Fl) to be symmetric around 〈Fl〉, this prescription would become problematic at low bulk densities, or more precisely, when the mean local density, 〈Fl〉, is small relative to the width of the distribution, φ, for in that case a significant number of negative local densities might be generated by eq 2.7. In the case considered herein, only a small fraction of the generated densities were negative. Hence we simply set all densities Fl(t) < 0 to 0. B. Parameters. The parameters necessary for simulating a near-critical SCF with the MTSGLE are (i) the normalized stationary friction kernel γ0(t), (ii) the parameters 〈Fl〉, φ, and τF, and (iii) the relation η[Fl(t)]. First, for simplicity, we let
γ0(t) ) 2δ(t)
(2.14)
Second, the parameters 〈Fl〉, φ, and τF can be determined from MD simulation. Specifically, the mean value 〈Fl〉 and the width φ can be derived from MD data for the equilibrium local density distribution, P(Fl), where Fl(t) is defined as in ref 21, by making use of the standard relations
〈Fl〉 )
∫0∞ dFlP(Fl)Fl
(2.15)
and φ ) x〈Fl2〉-〈Fl〉2. Since, in the critical region, the average local density around a tagged particle does not necessarily coincide with the bulk density, F, it will be convenient to introduce the shift ∆F ) 〈Fl〉 - F, where it is recognized that the Gaussian approximation Peq(Fl), eq 2.13, is centered at 〈Fl〉 ) F + ∆F. Similarly, the lifetime of local density correlations, τF (also called the local density reorganization time), corresponds to the time required for the local environment around a tagged particle to change substantially. It can be determined in terms of the local density autocorrelation function,13
(2.8)
P(Fl,0|F′l) ) δ(Fl - F′l)
{
2 -1/2
CF(t) )
and where the parameters 〈Fl〉, φ, and τF will reflect the SCF of interest (see section II.B). The resulting transition probability, P(Fl,t|F′l ), is required to satisfy the initial condition
P(Fl,t|F′l) ) [2πϑ2(t)]-1/2 exp -
reads
〈δFl(0)δFl(t)〉 〈[δFl(0)]2〉
δFl(t) ) Fl(t) - 〈Fl〉
(2.16)
by picking out the lowest frequency component of the relaxation via the zero-frequency Laplace transform:13
τF )
∫0∞ dt CF(t)
(2.17)
An extensive study of the behavior of atom-centered local environments in a neat two-dimensional Lennard-Jones fluid sufficiently close to its critical point has been recently made by Tucker and co-workers13,21 (see also papers by Yoshii and Okazaki12 and Kajimoto et al.25). The calculation was performed at a wide range of state points including the critical region. A way to relate the corresponding properties of the threedimensional Lennard-Jones fluid to the local density data from the two-dimensional fluid, which is based upon equating state points having the same correlation length, has been proposed in a previous paper.26 We therefore use this method, along with the two-dimensional data of Tucker and co-workers,13,21 to extract the local density parameters, ∆F, φ, and τF, for the three-
6678 J. Phys. Chem. B, Vol. 105, No. 28, 2001
Drozdov and Tucker TABLE 1: Parameters Used in the Calculation
Figure 1. Bulk density dependence of the shift, ∆F (solid circles), and the width, φ (open circles) of the local density distributions on a near-critical isotherm, T ) 1.236. Solid and dashed lines provide a guide for the eye.
Figure 2. Same as in Figure 1 but for the local density relaxation time.
dimensional Lennard-Jones fluid. Note that in doing so, we exclude the limit of low densities, F e 0.1 σ-3, from consideration, since the Brownian picture is not fully adequate in such a case. The correspondence method also uses the critical parameters of the three-dimensional Lennard-Jones fluid, which can be found in a paper by Panagiotopoulos;27 these are Tc ) 1.176 ( 0.008 /kB and Fc ) 0.33 ( 0.01 σ-3. Additionally, to simplify notations we set ) kB ) σ ) m )1 from here on. Our analysis26 of the MD data available for the local density dynamics13,21 shows that the latter can be reasonably approximated by an Ornstein-Uhlenbeck process, as assumed above. The shifts, ∆F, and widths, φ, extracted from P(Fl) are shown in Figure 1. As evidenced by this figure, a tagged solvent molecule experiences a mean local density which is enhanced relative to the bulk value, whereas the distribution becomes considerably broadened in the critical region. In Figure 2, the local density relaxation times computed from MD data for CF(t) are shown as a function of the bulk density along the T ) 1.24 near-critical isotherm. This quantity is seen to attain its maximum at the critical density when the critical fluctuations are also maximized. Finally, although the present prescription provides the necessary parameters as a function of density, in this paper we consider only a single near-critical state point for which ∆F, φ, and τF are given in Table 1. The remaining ingredient in our model of the fluctuating static friction coefficient, η(t), is the dependence of η on the local density Fl(t). To establish η[Fl(t)], we first recognize that,
T
F
〈Fl〉
φ
τF
ηeq
Deq
1.24
0.33
0.40
0.06
25
3.4
0.38
because the local density determines the friction, η[Fl] ) η(F), that is, the instantaneous static friction experienced by a particle in an inhomogeneous fluid in a region of local density Fl will be the same as that experienced by a particle in a homogeneous fluid of bulk density F ) Fl (at the same temperature). Hence, knowledge of η(F) at T for the homogeneous fluid (where F is the bulk density) provides an excellent approximation to η[Fl] at this temperature. In a homogeneous fluid, however, the static friction is largely nonfluctuating, and it may therefore be related to the diffusion coefficient, D(F,T), by the standard relation (m ) 1 and kB ) 1)
η(F,T) ) T/D(F,T)
(2.18)
such that knowledge of the diffusion coefficient may be used to generate η(F,T) and thus η[Fl] at T. Conveniently, Rowley and Painter10 have derived a formula for D(F,T) for the threedimensional Lennard-Jones fluid by fitting an extensive set of MD simulated diffusion data generated at points in the F,T plane, which cover nearly the entire phase diagram. Since these authors avoided the near-critical region, their fitted formula, even when extrapolated through this region, can be assumed to apply to the homogeneous fluid case. Hence, in our MTSGLE model, we evaluate η(t) by using eqs 2.7 and 2.8 to generate Fl(t) and eq 2.18 to generate η[Fl] at T, where D(F,T) is given by the expression of Rowley and Painter.10 Finally, we note that the procedure we have followedsfirst obtaining model parameters from MD simulation and then using the MTSGLE to examine the diffusion process in the compressible regimesis, computationally, much more efficient than is a procedure in which the diffusion coefficient is extracted directly from an MD simulation in the compressible regime. This enhanced efficiency arises because direct simulation of the diffusion coefficient requires simulation times which far exceed the slowest relaxation time in the system, which, near the critical point, is extremely long (τF). In contrast, the diffusion coefficients used as input to the MTSGLE theory are all evaluated well away from the critical region where relaxation times are fast and simulation times are correspondingly short. The parameters controlling the local density distribution and its time evolution are more problematic. In principle they can be evaluated via simulations that are shorter than required for accurate evaluation of the diffusion coefficient, but they are still limited by the need to equilibrate the slowly relaxing nearcritical system. Consequently, it is most efficient to estimate these parameters by other means. Herein we extracted these parameters from direct simulations of near-critical twodimensional systems, which are much less costly than are simulations of the corresponding three-dimensional systems. However, it is reasonable to expect that in the future we will learn both how to extract these local environmental parameters from spectroscopic observables5,28,29 and how to extrapolate them from state points of smaller correlation length to those of larger correlation length. III. Diffusion A. Standard Relations. As a preliminary, we briefly outline the main relationships following from the standard GLE, eq 1.4. Because this equation is linear in x˘ , it can be easily solved
Stochastic Dynamics in Near-Critical SCFs
J. Phys. Chem. B, Vol. 105, No. 28, 2001 6679
exactly to yield
Cv(t) ) 〈exp[-
x˘ (t) ) Cv(t)x˘ (0) +
∫0t dt′ Cv(t - t′)f(t′)
(3.1)
where the velocity autocorrelation function Cv(t) is defined by eq 1.3. The Laplace transform of this function,
C ˆ v(z) )
∫0t dt e-tzCv(t)
(3.2)
is determined by the Laplace transform of the time-dependent friction, γˆ (z), through the relation
C ˆ v(z) ) [z + γˆ (z)]-1
(3.3)
such that, by eq 1.6, C ˆ v(0) ) η-1. The quantity of interest that is much easier to measure in experiments18,19 than the velocity correlation function is the diffusion coefficient D. It is related to the velocity correlation lifetime τv, eq 1.2, by the Green-Kubo7 equation
D ) 〈x˘ 2(0)〉τv
(3.4)
Using eqs 3.3 and 1.6 and assuming that the initial velocity has an equilibrium distribution, i.e., 〈x˘ 2(0)〉 ) T, leads immediately to the relation
D ) Tτv ) T/η
(3.5)
which is the well-known Einstein result for the diffusion coefficient. Yet another important quantity, which is often used to monitor the Brownian motion and which is much easier to measure than the velocity autocorrelation function, is the mean-square displacement ∆(t). If we assume that the Brownian particle starts at time t ) 0 at x ) x(0) with the velocity x˘ (0), its mean-square displacement is given by15
∆(t) ) 〈[x(t) - x(0)]2〉 )
∫0t dt′ ∫0t dt′′ 〈x˘ (t′)x˘ (t′′)〉
(3.6)
Unlike the velocity autocorrelation function, which decays quickly and thus characterizes the short time behavior of the diffusing particle, the mean-square displacement describes its long-time limit. Using eqs 3.1 and 3.5, one can easily show that for ηt . 1, the leading term in the mean-square displacement is ∆(t) ) 2Dt, or, in other words, that the time derivative of ∆(t) approaches a limiting value
lim∆˙ (t) ) 2D tf∞
(3.7)
where D is given by eq 3.5. It should be noted that eq 3.7 holds true for any initial conditions x(0) and x˘ (0). B. MTSGLE Relations. For the more complicated dynamics associated with the MTSGLE, which is required to model inhomogeneous, near-critical SCFs, the standard relations, eqs 3.4, 3.5, and 3.7, are no longer valid. To calculate the velocity autocorrelation function Cv(t), eq 1.3, from the MTSGLE, we have to find a particular solution of eq 2.3, which at t ) 0 is equal to x˘ (0). In the case of free Brownian motion [U(x) ) 0 and γ0(t) defined by eq 2.14], the resulting expression for Cv(t) reads
) N-1
∫0t dt′ η(t′)]〉
N
∑ exp{-∫0 dt′ η[Fl,n(t′),T]} n)1 t
(3.8)
where n is the realization index. Since the static friction η is a nonlinear function of a stochastic variable, Fl,10 numerical schemes would be required to evaluate Cv(t), thus making straightforward evaluation via eqs 2.7 and 3.8 most effective. However, in some limiting situations this evaluation can be done analytically. Indeed, when the local density relaxation time, τF, is small relative to typical relaxation times of Cv(t), τF , τv ∼ η-1 eq , where
ηeq )
∫0∞ dFl η(Fl, T) Peq(Fl)
(3.9)
a negligible error will be introduced if one replaces in eq 3.8 η(t) by ηeq. This gives the expression
Cv(t)|τFf0 ) exp(-ηeqt)
(3.10)
which holds true in the entire time domain not only for equilibrium but also for nonequilibrium initial distribution of local densities, P0(Fl). In the opposite limit, τF . τv ∼ η-1 eq , one may also safely neglect the time dependence of η(t) in eq 3.8 by setting η(t) ) η(0). This immediately yields
Cv(t)|τFf∞ ) exp(-η0t)
(3.11)
where η0 is the static friction averaged over the initial distribution of local densities, P0(Fl),
η0 )
∫0∞ dFl η(Fl, T) P0(Fl)
(3.12)
It is thus seen that in both limits of small and large local density relaxation times the velocity autocorrelation function is characterized by single-exponential decay, eqs 3.10 and 3.11. In either of these two limits, then, τv, eq 1.2, represents an exponential decay constant. Since we are interested in nearcritical SCFs, though, we are most interested in the latter limit in which τF . τv. In this case, Cv decays via a single exponential because the velocities of the diffusing particles become uncorrelated well before the local density of the particles’ surroundings changes. Yet, on times of order τF, this environment will change. As a result, η will change and Cv measured from a later initial time may decay at a different rate than Cv measured from an earlier initial time. It is thus important to remember that, dependent upon the initial ensemble, Cv is not necessarily a stationary process, i.e., Cv may depend on the absolute time as Cv(t, t + t′)
Cv(t, t + t′) )
〈x˘ (t + t′)x˘ (t)〉 〈x˘ 2(t)〉
(3.13)
Consequently, the velocity correlation time τv may also depend on the absolute time, τv(t). That is, because η(t) and, in fact, also 〈η(t)〉t, where 〈 〉t represents an ensemble average at time t,
〈η(t)〉t )
∫0∞ dFl η(Fl) Pt(Fl)
(3.14)
6680 J. Phys. Chem. B, Vol. 105, No. 28, 2001
Drozdov and Tucker
with
Pt(Fl) )
∫0∞ dF′l P(Fl,t|F′l)P0(F′l)
(3.15)
realizations (see eq 3.8). In deriving eq 3.19, we assumed that the initial velocity x˘ (0) has an equilibrium distribution, i.e., 〈x˘ 2(0)〉 ) T. IV. Results
may change over long times, the definition of τv must incorporate this possible time dependence as
τv(t) )
∫0∞ dt′ Cv(t, t + t′)
(3.16)
[Note that this nonstationary behavior can be observed only if one begins with a nonequilibrium ensemble P0(Fl), see below]. And, clearly, if τF and τv are no longer well separated, Cv(t) will not decay exponentially, although τv may still be uniquely defined via eq 3.16. It is important to recognize that, for the MTSGLE, the simple relations between D, τv, and η, eq 3.5, no longer need rigorously apply in all cases. Indeed, since η is a fluctuating quantity, it clearly cannot be inversely proportional to a diffusion constant, although possibly such a relation might hold for 〈η〉. Similarly, the possible time dependence in both τv and 〈η〉 brings into question the relationship of these quantities to a traditional diffusion coefficient. This is especially true given that τv may change over times τF . τv during which one would expect to observe true diffusion motion for which D is constant. We therefore use the mean square displacement, eq 3.6, which is well-defined even for the MTSGLE, and its long-time limit [defined by t > τv and a linear time dependence in ∆(t)], eq 3.7, to define the diffusion coefficient in the inhomogeneous, near-critical SCF. Thus, in the usual fashion, we examine
1 D(t) ) ∆˙ (t) 2
(3.17)
where the diffusion constant itself is defined by the t > τv limit of D(t). By analogy with ηeq, eq 3.18, one can also define an equilibrium diffusion coefficient
Deq )
∫F∞
l min
dFl D(Fl,T)Peq(Fl)
(3.18)
Ideally, the lower limit of integration in eq 3.18 would be set to zero, Fl min ) 0. However, this integral typically diverges at low local densities, Fl f 0, because the zero local density distribution Peq(0), although very small, tends to be finite at vanishing Fl, because locally isolated particles exist, whereas the diffusion coefficient D(F,T) diverges in the low-density limit as F-1. To overcome this problem, we have used the requirement that at low bulk densities, where the fluid tends to become homogeneous, the mean diffusion coefficient computed from eq 3.18 should coincide with the homogeneous bulk fluid result, i.e, 〈D(Fl)〉Ff0,T ) D(F f 0,T), to fix Fl min. Finally, note that for the MTSGLE describing free Brownian motion, eq 3.8 with U(x) ) 0 and γ0(t) given by eq 2.14, an explicit expression for the mean square displacement reads
∫0t dt′ exp[-∫0t′ dt′′η(t′′)]}2〉 + t t t′′ 2T〈∫0 dt′ η(t′) {∫t′ dt′′ exp[-∫t′ dt′′′ η(t′′′)]}2〉
∆(t) ) T〈{
(3.19)
where the average is again to be performed over local density
We now use the MTSGLE to numerically examine the behavior of the self-diffusion coefficient in a near-critical SCF. The computation was carried out for a near-critical state point which is specified in Table 1, along with some other parameters that characterize the fluctuations of the local density Fl(t) at this state point. For each initial distribution of local densities, P0(Fl), an ensemble of N ) 10 000 realizations, {Fl,n(t)} (n ) 1, ..., N), was generated on a grid, ti ) i∆t (∆t ) 0.05), by making use of the transition probability P(Fl,∆t|F′l ), eq 2.10. Last, the time integrals entering eqs 3.8 and 3.19 were evaluated by means of the rectangular rule. In addition to the parameters controlling the thermodynamic state and η(t) shown in Table 1, we also have the flexibility to choose the initial ensemble of diffusing particles. That is, each realization of the MTSGLE requires an initial choice for the local density surrounding the diffusing particle at t ) 0, Fl(t ) 0). Thus, the initial distribution of diffusing particles within the fluid is characterized by the ensemble of initial local densities, P0(Fl), we use to generate successive realizations of the MTSGLE. The most obvious initial distribution of local densities to consider is the equilibrium distribution, P0(Fl) ) Peq(Fl), eq 2.13. A simple thought experiment, however, indicates that other, nonequilibrium initial ensembles might also be of interest in inhomogeneous, near-critical SCFs. To begin, let the molecules comprising the neat SCF under consideration have a spectroscopically accessible excited state, and let the energy of this excited state be sensitive to the local solvent density. Then, in the near-critical SCF, wherein the distribution of local densities experienced by the solvent molecules becomes quite broad, the spectral peak associated with the density-dependent excited state will be substantially broadened. Indeed, such broadened spectra have been observed experimentally in near-critical N2.28 It would then be possible, with a narrow band light source, to pick out just a subset of this broad distribution, i.e., to excite only those molecules that are initially in very low or very high local density regions. One could then follow the distribution of this nonequilibrium ensemble of excited solvent molecules.30 For simplicity, we shall assume that the spectroscopic excitation picks out only a single initial local density, F′l, such that
P0(Fl) ) δ(Fl - F′l)
(4.1)
We here consider four initial local density distributions: the equilibrium distribution P0(Fl) ) Peq(Fl), and three nonequilibrium distributions, eq 4.1 with F′l equal to a low local density, F′l ) 0.2, a high local density, F′l ) 0.6, and the mean local density, F′l ) 〈Fl〉 ) 0.4. As it turns out, both the equilibrium distribution and the mean local density distribution, F′l ) 〈Fl〉, predict nearly identical diffusion behavior, so we report only the results for F′l ) 〈Fl〉. In Figure 3 we present the initial time velocity correlation function Cv(t ) 0, t′), eq 3.13, for the three initial local density ensembles, eq 4.1 with F′l ) 0.2, 0.4, and 0.6. We see immediately that Cv decays on very short times, although, as should be expected, the time scale for this decay is strongly dependent on the local density. That is, the velocities become uncorrelated much more slowly, τv(t ) 0) ) 0.67, for atoms initially in low local density regions, F′l ) 0.2, than for atoms initially residing high local density regions, F′l ) 0.6, for which τv(t ) 0) ) 0.15. (Recall that all times are given in Lennard-
Stochastic Dynamics in Near-Critical SCFs
J. Phys. Chem. B, Vol. 105, No. 28, 2001 6681
Figure 3. Initial time velocity correlation function Cv(t ) 0,t′), eq 3.13, versus time. The initial local density distribution is given by eq 4.1 with F′l ) 0.2 (solid line), 0.4 (dashed line), and 0.6 (dotted-dashed line). The rest of the parameters are given in Table 1
Figure 5. Same as in Figure 3 but for the mean square displacement, ∆(t), eq 3.6 (top panel), and its time derivative, D(t), eq 3.17 (bottom panel).
Figure 4. Same as in Figure 3 but for the mean static friction, 〈η(t)〉t, eq 3.14 (top panel), and the velocity relaxation time, τv(t), eq 3.16 (bottom panel).
Jones units; 1 τLJ corresponds to 2 ps for Ar.) This behavior is easily understood from Figure 4, which shows the expected result that the mean static friction 〈η(t)〉t, eq 3.14, is much smaller for atoms in low local density than in high local density environments. Also shown in Figure 4 is τv(t), eq 3.16. As evidenced by the figure, both 〈η(t)〉t and τv(t) remain nearly constant on the time scale on which Cv(t ) 0, t′) decays, i.e., for times t , 3. This near constancy results because the values of 〈η(t)〉t and τv(t) are controlled by the distribution of local densities at time t. Since Fl changes on times of order τF, which is here τF ) 25 . 3, little evolution of the initial density ensembles, P0(Fl), will have occurred on the short time scale over which Cv decays. In Figure 5 we show the mean square displacements of the particles in the different ensembles. We see immediately that,
after short time transients on times of order τv(t ) 0), linear behavior is observed in ∆(t) (Figure 5, top panel), indicating the onset of diffusive motion, as expected. The onset of this diffusive behavior is also reflected in the nearly constant values of the diffusion coefficient (Figure 5, bottom panel) for times greater than ∼2. The most striking result from Figure 5, however, is the fact that the particles initially in low local density regions of the fluid (F′l ) 0.2) diffuse with a diffusion coefficient that is nearly four times as large as that describing the diffusion of particles that were initially in high local density regions of the fluid. Clearly, even though diffusion is a long time scale process, i.e., it occurs on times longer than τv, it is controlled by the local solvent environments experienced by the diffusing particle. Again, keep in mind that on the times considered, the distributions of local densities Pt(Fl) have not changed much from their initial distributions P0(Fl). Finally, note that resultant diffusion “constants” for each ensemble are reasonably well approximated by Tτv(0), given by the dotted lines, as might have been expected on the basis of eq 3.5. For the average local density ensemble, F′l ) 〈Fl〉 ) 0.4, we also find that the diffusion coefficient is nearly constant and equal to Deq, the equilibrium diffusion coefficient defined by eq 3.18. The same result is achieved if one starts with an initial equilibrium distribution of local densities, P0(Fl) ) Peq(Fl). We now turn our attention to the diffusion behavior on very long times, t > τF. Over such long times, the diffusing particles will sample many different local density environments, and, more importantly, the initial nonequilibrium local density ensembles will evolve into the equilibrium distribution, i.e., P0(Fl) ) δ(Fl - F′l) f Peq(Fl). As a consequence, the mean t.τF static friction experienced by the ensemble of diffusing particles, 〈η(t)〉t, eq 3.14 (see Figure 6, top panel), will evolve from its initial values, controlled by the nonequilibrium ensembles, to its mean value in the equilibrium ensemble, ηeq, eq 3.9. Notice
6682 J. Phys. Chem. B, Vol. 105, No. 28, 2001
Drozdov and Tucker
Figure 7. Local density relaxation time dependence of D(t), eq 3.17. Solid line: τF ) 25; dashed line: τF ) 4; dotted line: τF ) 1; dasheddotted line: τF ) 0.04. The initial local density distribution is given by eq 4.1 with F′l ) 0.2 (top panel) and 0.6 (bottom panel). The rest of the parameters are the same as in Table 1.
Figure 6. Same as in Figures 4 and 5 but for a larger time interval.
that, for the initial ensemble F′l ) 〈Fl〉 ) 0.4, 〈η(t)〉t)0 ≈ ηeq, such that the mean frictional forces on the diffusing particles remain unchanged as P0 f Peq, thus explaining why F′l ) 〈Fl〉 ) 0.4 and P0(Fl) ) Peq(Fl) generate the same diffusive behavior. On such long time scales, t . τF, we see that the mean square displacement (Figure 6 middle panel) again behaves in a linear fashion, and, also, that deviations from the linear behavior seen on the extremely short time scales become evident at more intermediate times. This behavior is more readily ascertained by examining the time derivative of ∆(t), D(t), eq 3.17, which is shown on a semilog plot in Figure 6 (bottom panel). There we see that (neglecting short time transients on times t j 1), while the diffusion constant for the F′l ) 〈Fl〉 ) 0.4 distribution is indeed nearly constant over all times, the diffusion coefficients for the initially low and high local density ensembles change as the local density distributions evolve in time. Indeed, a short plateau region of apparent constancy of D is observed (see also Figure 5, bottom panel) on times τv < t < τF for which Pt(Fl) ≈ P0(Fl), but this is followed by evolution of D(t) to the equilibrium ensemble value Deq, eq 3.18, as Pt(Fl) evolves to Peq(Fl). Hence, at long times the diffusion constant becomes independent of the initial ensemble considered, as it must for an equilibrium system as is considered here. Note, however,
that although the diffusion coefficient becomes independent of the initial ensemble at long times, the mean square displacement (Figure 6, middle panel) still reflects the initial ensemble differences, with particles initially in low density regions having traveled further, on average, than particles initially residing in high-density regions (see also ref 31). Last, in Figure 7 we artificially explore the effect of changing τF while keeping all other parameters of the MTSGLE constant, a variation which elucidates the importance of τF but which does not correspond to any realistic SCF state points.13,21 It is immediately evident that, as we have suggested, the rate of onset of the equilibrium diffusion behavior, Deq, is controlled by τF, the rate at which the local solvent densities surrounding the diffusing particle evolve. In the long time limit, τF . τv, intermediate behavior is indicated by a region of relatively invariant DFl′ )0.2 or 0.6 * Deq (see also Figure 5, bottom panel), whereas when the solvent environments interconvert rapidly, τF ) 0.4 , τv, DFl′ ) Deq immediately upon cessation of the short time transients. More complex intermediate behavior is, of course, observed when τF and τv are comparable. V. Conclusions We have proposed a microscopic model for diffusion in nearcritical supercritical solvents. The model, which we call the multiple time scale generalized Langevin equation (MTSGLE), is based on the generalized Langevin equation, but is further extended to include slow changes in the solvent response that are caused by the underlying density inhomogeneities. The theory uses as an input both equilibrium [∆F(F, T) and φ(F, T)] and dynamical [τF(F, T)] information about the local solvent environments that surround the diffusing particle over time. This information and the other ingredients entering the model [D(F,
Stochastic Dynamics in Near-Critical SCFs T)] were determined from MD simulations. The random force and the stochastic friction kernel of the resulting MTSGLE satisfy a generalization of the fluctuation-dissipation theorem, thus ensuring that it is consistent with the usual GLE in the limit of the homogeneous solvent. The present theory should provide a useful alternative to MD simulation for determining the velocity autocorrelation function and the mean-square displacement in such slowly evolving, complex fluids, since the direct simulational determination of these quantities in the vicinity of the liquid-vapor critical point is quite difficult.32 Application of the MTSGLE to diffusion in an inhomogeneous, near-critical SCF uncovered a rich and unusual diffusive behavior. In particular, it was observed that the time-scale separation between the local density relaxation time, τF, and the velocity correlation time, τv, i.e., τF . τv, present in nearcritical SCFs generates a diffusion behavior which may evolve over long times, τF. If, for example, one spectroscopically generates nonequilibrium initial local density ensembles comprised of particles residing in either low or high local density regions of the fluid, these particles will initially diffuse according to a diffusion coefficient which is larger or smaller than the equilibrium value, respectively. Over times t > τF, however, the initial local density ensemble will evolve into the equilibrium ensemble, causing the ensemble’s diffusion coefficient to approach the equilibrium value. At very long times t . τF, normal diffusion with D ) Deq will be observed, despite the large inhomogeneities present in the fluid. Finally, it should be remembered that even though the long time diffusion behavior is quite normal, the diffusion coefficient Deq ≈ D(〈Fl〉) * D(F) under near-critical conditions.26 Although we restricted our application to the case of free Brownian motion, the theory also offers a new paradigm for describing activated rate processes in near-critical solvents. It combines, in a natural way, the many-body structure of the density inhomogeneities with the intramolecular solute dynamics that is the cornerstone of the celebrated Kramers theory of chemical reactions.15,16 For example, the distribution of local solvent densities around a solute particle is an inherently manybody quantity that cannot be obtained from the pair distribution function. On the other hand, according to our theory, reactive events can be described by considering stochastic motions along a chosen solute coordinate. During these events, the solute particle will typically sample many local densities with different diffusion rates that change slowly over time as compared to the lifetime of random force correlations. The most appealing feature of the present theory is that it allows one to explicitly incorporate both time scales, a feature liable to make it applicable to many other complex stochastic processes, such as the diffusion of particles in a crystal containing inhomogeneous domains. This is in contrast to the standard GLE approach, which assumes that the medium is stationary and homogeneous, with its response taking place on the time scale defined only by the random force correlations. Acknowledgment. We thank Rigoberto Hernandez for helpful and insightful discussions. This work has been supported by the National Science Foundation under Grant No. CHE9727361, and S.C.T. acknowledges a Camille Dreyfuss TeacherScholar award. References and Notes (1) McHugh, M. A.; Krukonis, V. J. Supercritical Fluids Extraction: Principles and Practice, 2nd ed.; Butterworth: Boston, 1994. (2) Sengers, J. V.; Luettmer-Strathmann, J. In Transport Properties of Fluids: Their Correlation, Prediction and Estimation; Millat, J., Dymond, J. H., Nieto de Castro, C. A., Eds.; IUPAC, Cambridge University Press: New York, 1996; pp 113-137.
J. Phys. Chem. B, Vol. 105, No. 28, 2001 6683 (3) Subramaniam, B.; McHugh, M. A. Ind. Eng. Chem. Res. 1986, 25, 1. (4) Johnston, K. P.; Haynes, C. AIChE J. 1987, 33, 2017. (5) Wu, B. C.; Klein, M. T.; Sandler, S. I. Ind. Eng. Chem. Res. 1991, 30, 822. (6) Bennecke, J. F.; Chateauneuf, J. E. Chem. ReV. 1999, 99, 433. (7) Savage, P. E.; Gopalan, S.; Mizan, T. I.; Martino, C. J.; Brock, E. E. AIChE J. 1995, 41, 1723. Brennecke, J. F. Chem. Ind. 1996, 831. (8) Kajimoto, O. Chem. ReV. 1999, 99, 355. (9) Tucker, S. C. Chem. ReV. 1999, 99, 391. Tucker, S. C.; Maddox, M. W. J. Phys. Chem. B 1998, 102, 2437. (10) Hoover, W. G. Molecular Dynamics; Springer: Berlin, 1986. Simulation of Liquids and Solids; Ciccotti, G., Frenkel, D., McDonald, I. R., Eds.; Noth-Holland: Amsterdam, 1987. Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Clarendon Press: Oxford, 1989. Haile, J. M. Molecular Dynamics Simulation: Elementary Methods; Wiley: New York, 1992. Rowley, R. L. Statistical Mechanics for Thermophysical Property Calculations; Prentice Hall: Englewood Cliffs, NJ, 1994. (11) For a recent review on MD, see: Tuckerman, M. E.; Martyna, G. J. J. Phys. Chem. B 2000, 104, 159. (12) Heyes, D. M. Phys. ReV. B 1988, 37, 5677. Heyes, D. M. Chem. Phys. Lett. 1988, 153, 319. Stoker, J. M.; Rowley, R. L. J. Chem. Phys. 1989, 91, 3670. Castillo, R.; Villaverde, A.; Orozco, J. Mol. Phys. 1991, 74, 1315. Iwai, Y.; Higashi, H.; Uchida, H.; Arai, Y. Fluid Phase Equilibria 1997, 127, 251. (13) Rowley, R. L.; Painter, M. M. Int. J. Thermophys. 1997, 18, 1109. (14) Mishra, B.; Berne, B. J. J. Chem. Phys. 1995, 103, 1160. Kurnikova, M. G.; Waldeck, D. H.; Coalson, R. D. J. Chem. Phys. 1996, 105, 628. Biswas, R.; Bagchi, B. J. Chem. Phys. 1996, 105, 7543. (15) Yoshii, N.; Okazaki, S. J. Chem. Phys. 1997, 107, 2020. (16) Maddox, M. W.; Goodyear, G.; Tucker, S. C. J. Phys. Chem. B 2000, 104, 6266. (17) Martinez, H. L.; Ravi, R.; Tucker, S. C. J. Chem. Phys. 1996, 104, 1067. Goodyear, G.; Maddox, M. W.; Tucker, S. C. J. Chem. Phys. 2000, 112, 10327. (18) Risken, H. The Fokker-Planck Equation, Methods of Solution and Applications, 2nd ed.; Springer: New York, 1989. (19) Ha¨nggi, P.; Talkner, P.; Borkovec, M. ReV. Mod. Phys. 1990, 62, 251. Mel’nikov, V. I. Phys. Rep. 1991, 209, 1. ActiVated Barrier Crossing; Ha¨nggi, P., Fleming, G., Eds.; World Scientific: Singapore, 1992. New Trends in Kramers’ Reaction Rate Theory; Talkner, P., Ha¨nggi, P., Eds.; Kluwer Academic: Dordrecht, 1995. (20) Zwanzig, R. J. Stat. Phys. 1973, 9, 215. (21) Trappeniers, N. J.; Oosting, P. H. Phys. Lett. 1966, 23, 445. Martynets, V. G.; Matizen, EÄ . V. SoV. Phys. JETP 1970, 31, 228; 1975, 40, 507. Tsekhanskaya, Yu. V. Russ. J. Phys. Chem. 1971, 45, 744. (22) Liong, K. K.; Wells, P. A.; Foster, N. R. J. Supercrit. Fluids 1991, 4, 91. Eaton, A. P.; Akgerman, A. Ind. Eng. Chem. Res. 1997, 36, 923. He, C.-H. AIChE J. 1997, 43, 2944. Liu, H.; Silva, C. M.; Macedo, E. A. Ind. Eng. Chem. Res. 1997, 36, 246. (23) Straub, J. E.; Borkovec, M.; Berne, B. J. J. Chem. Phys. 1985, 83, 3172; 1986, 84, 1788. (24) Tucker, S. C.; Tuckerman, M. E.; Berne, B. J.; Pollak, E. J. Chem. Phys. 1991, 95, 5809. (25) Reese, S. K.; Tucker, S. C.; Schenter, G. K. J. Chem. Phys. 1995, 102, 104. (26) Maddox, M. W.; Goodyear, G.; Tucker, S. C. J. Phys. Chem. 2000, 104, 6248. (27) We consider only the case where the intermolecular interaction potential is short ranged, such as for the Lennard-Jones potential. (28) Hernandez, R.; Somer, F. L. J. Phys. Chem. B 1999, 103, 1064. Hernandez, R. J. Chem. Phys. 1999, 111, 7701. (29) The validity of this assumption is based on MD simulations which show that the autocorrelation time of the random force, τf, is rather insensitive to both bulk and local densities (in homogeneous and inhomogeneous fluids, respectively), whereas, on the other hand, the magnitude, η, is very sensitive to density changes. (30) Kajimoto, O.; Sekiguchi, A. K.; Nayuki, T.; Kobayashi, T. Ber. Bunsen-Ges. Phys. Chem. 1997, 101, 600. (31) Drozdov, A. N.; Tucker, S. C. J. Chem. Phys. 2001, 114, 4912. (32) Panagiotopoulos, A. Z. Int. J. Thermophys. 1994, 15, 1057. (33) Clouter, M. J.; Kiefte, H.; Deacon, C. G. Phys. ReV. A 1986, 33, 2749. (34) Song, W.; Biswas, R.; Maroncelli, M. J. Phys. Chem. A 2000, 104, 6924. (35) Technically, one would have to assume a fast relaxation down to a state having the same intermolecular interactions as does the ground state for our subsequent analysis, which assumes that the ground state Peq(Fl) controls the MTSGLE parameters, to apply to this thought experiment. (36) Cao, J. Phys. ReV. E 2001, 63, No. 4. (37) Kitao, O.; Tanabe, K.; Ono, S.; Kumakura, S.; Nakanishi, K. Fluid Phase Equilib. 1998, 144, 279.
