THE JOURNAL OF CHEMICAL PHYSICS 141, 214109 (2014)
Markovian master equation for a classical particle coupled with arbitrary strength to a harmonic bath Maxim F. Gelin Department of Chemistry, Technische Universität München, D-85747 Garching, Germany
(Received 20 March 2014; accepted 13 November 2014; published online 3 December 2014) We consider a classical point particle bilinearly coupled to a harmonic bath. Assuming that the evolution of the particle is monitored on a timescale which is longer than the characteristic bath correlation time, we derive the Markovian master equation for the probability density of the particle. The relaxation operator of this master equation is evaluated analytically, without invoking the perturbation theory and the approximation of weak system-bath coupling. When the bath correlation time tends to zero, the Fokker-Planck equation is recovered. For a finite bath correlation time, the relaxation operator contains contributions of all orders in the system-bath coupling. © 2014 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4902438] I. INTRODUCTION
II. NAKAJIMA-ZWANZIG MASTER EQUATION
The system-bath approach is one of the key paradigms in statistical mechanics. The idea is to derive the (irreversible) master equation for the probability density function (density matrix in quantum case) ρ(t) of the system starting from the combined system+bath (reversible) Hamiltonian dynamics. If the characteristic time of the bath-induced correlations is short on the timescale of the system dynamics, one can hope to derive a Markovian master equation1–8 of the form
Let us consider a classical particle bilinearly coupled to a classical harmonic bath. The total Hamiltonian H can be partitioned into the system (S) Hamiltonian, the bath (B) Hamiltonian, and their coupling,
∂t ρ = (−iLS + R)ρ. Here, the Liouville operator LS describes the time evolution of the system alone, and R is responsible for the bath-induced relaxation. Even if we restrict our consideration to classical mechanics, it is not an easy task to obtain the explicit form of the relaxation operator R. In general, one has to resort to additional approximations and/or assumptions, those like weak systembath coupling or delta-correlated Brownian forces (FokkerPlanck equations), low-density approximation (Boltzmann equations), mean-field approximation (Vlasov equations).8 The aim of the present work is to show that the relaxation operator R can be analytically evaluated without additional approximations in the particular, but important case of a classical particle bilinearly coupled to a harmonic bath. The paper is organized as follows. To setup the stage, we introduce the Nakajima-Zwanzig (NZ) master equation in Sec. II and discuss its Markovian limit in Sec. III. The explicit form of the relaxation operator R in the absence of an external force is derived in Sec. IV. General properties of the obtained master equation are discussed in Sec. V. Section VI outlines a generalization of the master equation. The predictions of the generalized Markovian master equation are numerically tested against the exact results in Sec. VII, and its possible applications are discussed in Sec. VIII. Section IX briefly summarizes our main findings. 0021-9606/2014/141(21)/214109/8/$30.00
H = HS + HB + HSB .
(1)
In our case, HS =
p2 + U (x), 2m
(2)
where x, p, and m are the position, momentum, and mass of the particle, and U(x) is an external potential. The system-bath coupling and the bath are described by 2 1 pi2 ci 2 + mi ωi xi − x , (3) HSB + HB = 2 i mi mi ωi2 where xi , pi , mi , and ωi are the positions, momenta, masses, and frequencies of the bath oscillators, and ci are the systembath coupling coefficients. The influence of the bath on the dynamics of the particle is determined by the bath correlation function2, 3 g(t) =
1 ci2 cos{ωi t}. m i mi ωi2
(4)
The Liouville equation for the total system+bath probability density function ρ tot (t) reads ∂t ρtot (t) = −iLρtot (t), where the total Liouvillian is defined as p i iL ≡ iLS + F∂p + ∂x − F i ∂p . i mi i i
(5)
(6)
Here,
141, 214109-1
iLS ≡
p ∂ + F (x)∂p m x
(7)
© 2014 AIP Publishing LLC
214109-2
Maxim F. Gelin
J. Chem. Phys. 141, 214109 (2014)
is the Liouvillian corresponding to the Hamiltonian HS , F(x) ≡ −dU(x)/dx is the external force, and ci ci Fi , Fi = xi − x , (8) F= mi ωi2 i
which guarantees that the particle attains the Boltzmann equilibrium distribution H (17) ρS = ZS−1 exp − S , kB T
where F is the total force exerted on the system by the bath oscillators. The initial condition to Eq. (5) is considered in the form
ZS being the partition function. We also define the reduced probability density function ∂t ρ(p, t) ≡ dxρ(x, p, t), (18)
ρtot (0) = ρ0 ρeq .
(9)
Here, ρ 0 is a certain distribution in the system phase space and H + HB −1 (10) exp − SB ρeq = Zeq kB T is the equilibrium bath distribution at the presence of the system, Zeq is the partition function, kB is the Boltzmann constant, and T is the bath temperature. Note that ρ eq of Eq. (10) is x-dependent and accounts for the initial system-bath correlations. We prefer to work with this distribution due to the following two reasons. First, it is more physical than the (commonly used) product system-bath distribution ρ 0 ρ eq (x = 0). A general discussion of the importance of system-bath correlations can be found in Refs. 9–11. Second, the correlated initial distribution (9) guarantees that the stochastic force f(t) generated by the Hamiltonian (1) in the generalized Langevin equation (Sec. IV C) possesses zero mean, f(t) = 0. See Chap. 3 of the monograph in Ref. 3 and Ref. 12 for the discussion of these issues. We introduce the operator9, 10, 13 P ≡ ρeq dpdx (11) (p and x denote collectively the momenta and positions of the bath oscillators) which projects the total system+bath probability density function into the correlated initial system-bath distribution (10). The application of the NZ formalism2, 3, 9 to Eq. (5) yields then a formally exact non-Markovian master equation for the probability density function of the particle, t dt R(t − t )ρ(x, p, t ). ∂t ρ(x, p, t) = −iLS ρ(x, p, t) + 0
which describes relaxation of the particle in momentum space. If the external potential is zero (F(x) = 0), the integration of Eq. (12) over x yields an exact NZ master equation ∂t ρ(p, t) =
0
where
t
dt Rp (t − t )ρ(p, t ),
(19)
Rp (t) ≡
dxR(t)|F (x)=0 .
(20)
III. NZ MASTER EQUATION IN THE MARKOVIAN LIMIT
Let us now assume that the characteristic decay time τ B of the bath correlation function (4) is short on the timescale of the system dynamics, τ S . If we are interested in the evolution of the particle for t > τ B , we can extend the time integration in Eq. (12) to infinity and expand the probability density in the Taylor series, ρ(p, t ) = ρ(p, t) + ∂t ρ(p, t)(t − t) + · · · ,
(21)
τ B /τ S being an effective small parameter of the expansion.14 In the leading order in τ B /τ S , we replace ρ(p, t ) by ρ(p, t) and arrive at the Markovian master equation ∂t ρ(x, p, t) = (−iLS + R)ρ(x, p, t).
(22)
Similarly, the Markovian analogue of Eq. (19) reads ∂t ρ(p, t) = Rp ρ(p, t). Here,
∞
R≡ 0
dt R(t ), Rp ≡
∞ 0
dt Rp (t )
(23)
(24)
(12) The relaxation operator in Eq. (12) can be written as a generalized Fokker-Planck operator,9 p . (13) R(t) = ∂p G(t) ∂p + mkB T
are the Markovian relaxation operators. In general, the evaluation of the operator R requires additional approximations (see Refs. 9, 15, 16, and Sec. VI). Rp , on the other hand, can be calculated analytically (see below).
Here,
IV. EVALUATION OF THE REDUCED RELAXATION OPERATOR
G(t) =
dpdxFe−i(1−P)Lt Fρeq
(14)
is the force-force correlation operator. According to Eq. (13), the relaxation operator obeys normalization dpdxR(t) = 0, (15) which insures the conservation of probability ( dpdxρ(x, p, t) ≡ 1) and detailed balance R(t)ρS = 0,
(16)
The explicit evaluation of Rp by Eqs. (13), (20) and (24) looks difficult. Fortunately, there exists an alternative way for achieving this goal. The derivation involves three steps. A. Step No. 1
Let A(p) and B(p) be arbitrary functions of the momentum of the particle. The correlation function CAB (t) = A(p) exp(Rp t)B(p)
(25)
214109-3
Maxim F. Gelin
J. Chem. Phys. 141, 214109 (2014)
(. . . denotes equilibrium averaging) can be calculated through the reduced probability density function (18) as (26) CAB (t) = dpA(p)ρ(p, t),
(F(x) = 0) it reads d p(t) = − dt
with ρ0 (p) =
1 2π mkB T
exp −
p2 2mkB T
(27) .
The integral relaxation time of CAB (t) is defined as ∞ dtCAB (t). τAB ≡
t
dt g(t − t )p(t ) + f (t),
0
(28)
Here, f(t) is a Gaussian colored noise with zero mean (f(t) = 0), and the memory function g(t) is defined by Eq. (4). The generalized Langevin equation (34), in turn, is equivalent to the non-Markovian Fokker-Planck equation19, 20 p ρ(p, t). (35) ∂t ρ(p, t) = a(t)∂p ∂p + mkB T Here,
(29)
a(t) = −
0
Obviously, the CAB (t) evaluated through the NZ master equation (20) and the Markovian master equation (23) are different. However, their integral relaxation times τ AB are identical. This is immediately proven by the direct time integration of Eqs. (20) and (23) and the use of Eq. (24). B. Step No. 2
The relaxation operator R(t) obeys the detailed balance (16), which simplifies to Rp ρ0 (p) = 0. We assume that Rp is Hermitian, and its eigenvalues are all negative, except the eigenfunction ρ 0 (p) possessing zero eigenvalue.17 For compactness, we adopt the bra-ket notation Rp |α = −να |α , α = 0, 1, 2, . . .
(30)
(ν 0 = 0 ensures the detailed balance, ν α > 0 for α = 1, 2, . . . ). The Hermiticity of Rp is difficult to prove in full mathematical rigor (cf. Ref. 18), but the requirement is reasonable from physical point of view and will be demonstrated a posteriori in Sec. IV C. To reveal the physical meaning of ν α , we specify the correlation function α(0) |α(t) ≡ α| exp(Rp t) |α = exp(−να t).
(31)
The integral relaxation time of this correlation function is defined as ∞ dt α(0) |α(t) = να−1 , α = 1, 2, . . . . (32) Tα ≡
∞ α=1
|α
1 α| Tα
(36)
The eigenfunctions are proportional to the Hermite polynomials (see, e.g., Ref. 14) 1 p p p|n = ρ0 (p)Hn , k|p = Hk n! 2mkB T 2mkB T (38) and are orthonormal, with the scalar product defined as k |n ≡ dp k |p p |n = δkn . (39) We then specify the fundamental correlation function t n(0) |n(t) ≡ Cn (t) = exp{−nb(t)}, b(t) = dt a(t ) 0
(40) (p(t) is a Gaussian process, hence Cn (t) = (C1 (t))n ). The integral relaxation time of this correlation function reads ∞ dtCn (t) (41) τn = 0
(the term with n = 0, which corresponds to the equilibrium distribution (28), yields τ 0 = ∞). Therefore, according to Eq. (33), Rp ≡ −
Hence, we can write
dC1 (t)/dt C1 (t)
is the time-dependent friction and C1 (t) ≡ p(0)p(t) is the momentum correlation function. As in Sec. IV B, we adopt the bra-ket notation and introduce the eigenfunctions of the Fokker-Planck operator, p |n = −n |n . (37) ∂p ∂p + mkB T
0
Rp ≡ −
(34)
f (t)f (t ) = g(t − t ).
where ρ(p, t) obeys Eq. (23) with the initial condition ρ(p, 0) = ρ0 (p)B(p)
∞
|n
n=1
1 n| . τn
(42)
(33)
(T0−1 = 0 due to the detailed balance).
V. GENERAL PROPERTIES OF THE MASTER EQUATION
C. Step No. 3
The relaxation operator (42) can be represented in several equivalent forms. Let us introduce the rate ν and the relaxation cross-sections σ n according to the expression
The representation (33) of the operator Rp can be easily constructed. The Hamiltonian (1)–(3) is equivalent to the generalized Langevin equation.2, 3 If the external potential is zero
1 ≡ ν(1 − σn ), τn
(43)
214109-4
Maxim F. Gelin
J. Chem. Phys. 141, 214109 (2014)
where σ 0 ≡ 1 (the decomposition (43) is not unique, of course). Using the definition of the bra and ket eigenfunctions (38) and the orthonormality relation (39), we can explicitly write Rp as an integral operator and recast the master equation (23) in the standard integro-differential form (44) ∂t ρ(p, t) = −νρ(p, t) + ν dp T (p|p )ρ(p , t). The relaxation kernel is defined as ∞ T (p|p ) ≡ σn p|nn|p .
(45)
n=0
Interestingly, the master equation (44) describes a Markovian, but non-Gaussian process, in contrast to the non-Markovian Fokker-Planck equation (35) and the generalized Langevin equation (34), which describe a Gaussian, but non-Markovian process. Alternatively but equivalently, the Markovian master equation (44) can be written as a Kramers-Moyal expansion21 ∂t ρ(p, t) =
∞ n=1
Mn (p) = τn−1
∂pn Mn (p)ρ(p, t),
B. ∂ x -dependence of R
Physically, ∂ x accounts for the change of the position of the particle during the bath correlation time τ B . Neglecting the ∂ x -dependence is tantamount to assuming that the particle cannot significantly move during τ B (see the discussion in Ref. 24). The assumption is reasonable, and its validity is similar to that of the Markovian approximation, τ B /τ S 1. If the ∂ x -dependence of R cannot be neglected, the Markovian approximation itself is questionable. C. Approximate master equation
(46) dp (p − p )n p |n n|p .
Equation (46) is obtained from the master equation (44) by the expansion of ρ(p , t) in Taylor series and integrations by parts.21 In the limit τ B → 0, only two first terms in the Kramers-Moyal expansion, M1 (p) = ξ (mkB T)−1 p and M2 (p) = ξ , remain nonzero (g(t) → δ(t), τ n → 1/(ξ n)). Then Eq. (46) reduces to the Markovian Fokker-Planck equation p ∂t ρ(p, t) = ξ ∂p ∂p + ρ(p, t). (47) mkB T For finite τ B , all Mn (p) = 0. Thus, the Markovian master equation (44) can be viewed as a tool for summing up the entire Kramers-Moyal series (46).
Above we have presented certain physical arguments supporting the idea that the relaxation operator R should not significantly depend on Fx and ∂ x . More solid arguments can be obtained by expanding the total Liouvillian L in Eq. (14) in powers of LS (cf. Ref. 9). It is straightforward to show that τ B /τ S is an effective small parameter of this expansion, since the force-force operator G(t) of Eq. (14) is proportional to the bath correlation function g(t) of Eq. (56) in the leading (zero) order in τ B /τ S . Then both the F(x)- and the ∂ x -dependence of the relaxation operator can be dropped and we can identify R with Rp R(F (x ); ∂x , p, ∂p ) ≈ R(0; 0, p, ∂p ) ≡ dxR(0; ∂x , p, ∂p ) ≡ Rp .
(49)
Adopting the approximation of Eq. (49), we can write the Markovian master equation (22) as follows:
VI. MASTER EQUATION FOR ρ(x, p, t)
The reduced probability density ρ(p, t) has a rather limited range of applications. One normally deals with the master equation (22) for the full probability density ρ(x, p, t) of the particle in an external potential U(x). According to Eqs. (13) and (14), the relaxation operator R in the master equation (22) is, in general, a functional of the external force F(x), and also depends on the operators ∂ x , p, ∂ p R ≡ R(F (x ); ∂x , p, ∂p ).
on the timescale τ B , the dependence of R on F(x) can be neglected. This should be a reasonable approximation in many practically relevant cases, since the short-time dynamics of the particle is predominantly determined by the (strong) interaction with its nearest neighbors.22, 23 We emphasize that the F(x)-independence of R does not mean that the influence of the external force is neglected, since F(x) enters the Markovian master equation through the system Liouvillian LS .
(48)
Due to the isotropy of space, R cannot depend on x explicitly, but only through the external force F(x) (this follows immediately from Eqs. (8) and (14)). Let us consider now possible approximations to Eq. (48). A. F(x)-dependence of R
According to Eq. (14), F(x) enters R through the projected total Liouvillian (1 − P)L. If the impact of the external force F(x) is smaller than that of the bath-induced forces F
∂t ρ(x, p, t) = (−iLS + Rp )ρ(x, p, t).
(50)
As is argued above, the validity of Eq. (49) is closely related to the validity of the Markovian approximation itself. It is thus not surprising that virtually all Markovian master equations available in the literature are of the form of Eq. (50), with the relaxation operators acting Rp acting in the momentum space only.1–6, 8, 25 If necessary, the corrections due to F(x ) and ∂ x can be evaluated perturbatively (cf. Refs. 9 and 15). VII. TESTING THE APPROXIMATE MASTER EQUATION
The Hamiltonian (1)-(3) generates the nonlinear generalized Langevin equations d x(t) = p(t)/m, dt t d p(t) = F (x(t))/m − dt g(t − t )p(t ) + f (t). dt 0
(51) (52)
214109-5
Maxim F. Gelin
J. Chem. Phys. 141, 214109 (2014)
If the memory function g(t) can be represented as a linear combination of N exponentials, Eqs. (51) and (52) can exactly be mapped into the Markovian Fokker-Planck equation involving x, p, and N additional auxiliary variables.26–28 Hence, the simplest (N = 1) exact Fokker-Planck equation depends on three variables. Below we study, numerically, whether the master equation (50) delivers a good approximation to the exact three-variable Fokker-Planck equation.
where the bra and ket eigenfunctions are proportional to the Hermite polynomials √ 1 ω |n = √ exp{−ω2 /2}Hn (ω/ 2), 2π n! (60) √ n |ω = Hn (ω/ 2). The relaxation times are specified by Eqs. (40) and (41). For the exponential memory function (56), they can be evaluated analytically
A. Formulation of the problem
We consider one-dimensional rotator (specifying by the angle φ and the angular velocity ω) in a N-well potential
τn =
U (φ) = {1 − cos(N φ)}. (53) N The rotator is subject to stochastic torques, so that its motion is described by the rotational analogue of the generalized Langevin equations (51) and (52)
δ≡
d φ(t) = ω(t), (54) dt t d d ω(t) = − U (φ(t)) − dt g(t − t )ω(t ) + f (t). dt dφ 0 (55) In this section, we adopt the dimensionless variables: Time, angular velocity, and energy are expressed in units of τ r , τr−1 , and kB T, respectively. Here, τr = I /(kB T ) is the period of free rotation, and I is the moment of inertia of the rotator. Our choice of the system specified by Eqs. (54) and (55) is stipulated by two reasons. Fundamentally, Eqs. (54) and (55) define a paradigmatic model of rotational relaxation in liquids, which mimics rotation of a molecule in the cage formed by its nearest neighbors (see, e.g., Refs. 15 and 29). Technically, Eqs. (54) and (55) are essentially nonlinear, and a proper description of the ensuing dynamics is a significant challenge for any approximate method. We restrict ourselves to the consideration of the exponential bath correlation function ξ exp{−t/τB }. (56) g(t) = τB Here, ξ and τ B are the strength and the correlation time of the Gaussian torque f(t). B. Rotational master equation
Equations (53)-(56) generate the rotational counterpart of the Markovian master equation (50)
∂t ρ(φ, ω, t) = − iLr + RME ρ(φ, ω, t) (57) r (the subscript “r” stands for “rotator”). The bath-free evolution of the rotator is described by the Liouvillian iLr = ω∂φ −
dU (φ) ∂ . dφ ω
(58)
According to Eq. (42), the relaxation operator is defined as ≡− RME r
∞ n=1
|n
1 n| , τn
(59)
n 2τB n (δ + 1)m (δ − 1)n−m Cm , (2δ)n m=0 n + (n − 2m)δ
1 − 4ξ τB , Cmn ≡
n! . m!(n − m)!
(61) (62)
In the limit τ B → 0, we recover the Markovian Fokker-Planck result τ n = 1/(ξ n). Explicitly, the corresponding FokkerPlanck equation is written as
(63) ∂t ρ(φ, ω, t) = − iLr + RFP r ρ(φ, ω, t), where
2 RFP r = ξ ∂ω + ∂ω ω .
(64)
For finite τ B , the integral relaxation times (61) are the ratios of polynomials in ξ τ B ; the degree of these polynomials increases with n. The perturbative evaluation of τ n with ξ τ B as a small parameter would require tedious calculations. C. Exact three-variable Fokker-Planck equation
The generalized Langevin equation (55) with the exponential memory function (56) can equivalently be represented as26–28 d d (65) ω(t) = − U (φ(t)) + ξ/τB q(t), dt dφ d q(t) = − ξ/τB ω − τB−1 q + fq (t), (66) dt where q(t) is an auxiliary variable and fq (t) is a deltacorrelated Gaussian process, fq (t)fq (t ) = 2τB−1 δ(t − t ).
(67)
The stochastic equations (55), (65) and (66) are equivalent to the Fokker-Planck equation
ρ(φ, ω, q, t) (68) ∂t ρ(φ, ω, q, t) = − iLr + REX r (the superscript “EX” stands for “exact”). Here, Lr is defined by Eq. (58) and
ξ/τB (ω∂q − q∂ω ) + τB−1 ∂q2 + ∂q q . (69) REX r =
D. Computational details
To solve Eq. (57), we expand the probability density into the series, ω |n e−ikφ ρnk (t) (70) ρ(φ, ω, t) = n,k
Maxim F. Gelin
J. Chem. Phys. 141, 214109 (2014)
(ω|n defined per Eq. (60) are proportional to the Hermit polynomials). Similarly, the solution of the exact threevariable Fokker-Planck equation (68) can be written as ω |n q |m e−ikφ ρnmk (t). (71) ρ(φ, ω, t) = n,m,k
1.1
The expansions (70) and (71) generate the matrix representation of the evolution operators
− iLr + RME and − iLr + REX (72) r r nk,n k nmk,n m k
n m k
facilitating and substantially speeding up the numerics (see, e.g., Refs. 30 and 31 for more details). E. Energy relaxation
There exist several interesting aspects of rotational and orientational relaxation which can be studied by invoking Eqs. (57) and (68). Here, we concentrate on the rotational energy relaxation. First, this problem has recently received considerable attention in the context of the rotational relaxation under non-equilibrium conditions and breakdown of the linear response description.18, 32–35 Second, relaxation of the rotational energy reveals importance of the bath memory effects, which cannot be grasped by the standard Markovian Fokker-Planck equation (63). Fig. 1 shows the equilibrium energy correlation function C(t) = E(0)E(t) , E = ω2 /2 + U (φ)
(73)
(. . . means equilibrium averaging). The system-bath coupling is moderate (ξ = 1). Panel (a) corresponds to a moderate bath memory (τ B = 1), while panel (b) corresponds to a long bath memory (τ B = 10). C(t) calculated via the exact three-mode Fokker-Plank equation (68), the Markovian master equation (58), and the Markovian Fokker-Planck equation (63) correspond to green, red, and black lines, respectively. We prefer to consider unnormalized correlation functions (vide infra), hence C(0) = 1. If the bath memory is moderate, all three descriptions grasp the essentials of the energy relaxation (panel (a)). The Markovian master equation is superior over the Markovian Fokker-Planck equation at short times t < 2, but no qualitative differences between the three approaches occurs. Panel (b) offers a qualitatively different picture. The Markovian Fokker-Planck equation is significantly off, predicting an order of magnitude shorter relaxation time (cf. black and green
1 0.9 0.8 0.7 0.6 0
1
3
2 t (b)
1.2 1.1 C(t)
and convert the master equations (57) and (68) into the linear matrix first-order differential equations for the “vectors” ρ nk (t) and ρ nmk (t), respectively. These matrix equations are solved by using the fourth-order Runge-Kutta integrator. In our calculations, we consider a two-well potential of Eq. (53) with N = 2 and = 1, so that the torque −dU/dφ = − sin (2φ). In this case, all matrix elements of Eq. (72) are evaluated analytically in closed form and in elementary functions. Furthermore, the structure of the matrices of Eq. (72) permits analytical evaluation of all matrix-vectors multiplication of the kind
− iLr + REX ρ (t), r nmk,n m k n m k
(a)
1.2
C(t)
214109-6
1 0.9 0.8 0.7 0.6 0
10
20
30
t
FIG. 1. Equilibrium energy correlation functions C(t) evaluated for a moderate system-bath coupling ξ = 1. Panel (a) corresponds to a moderate bath memory τ B = 1, while panel (b) corresponds to a long bath memory τ B = 10. Green lines: exact results; black lines: Markovian Fokker-Planck equation; red lines: Markovian master equation.
curves). The Markovian master equation performs much better. It cannot reproduce the oscillatory behavior of the exact C(t) (see the discussion below), but delivers a smoothened version of C(t) and correctly grasps its relaxation timescale (cf. red and green curves). Fig. 2 monitors equilibration of the rotational energy after the non-equilibrium preparation. It reproduces the mean energy evolution S(t) = E(t)ne .
(74)
Here, . . . ne means that the rotator at t = 0 possesses equilibrium Boltzmann distribution over its angular velocity ω, while the orientation angle is fixed at φ = π , which corresponds to the minimum of one of the potential wells. As in Fig. 1, panels (a) and (b) correspond to a moderate (τ B = 1) and long (τ B = 10) bath memory. Green, red, and black lines correspond to the exact three-mode Fokker-Plank equation (68), the Markovian master equation (58), and the Markovian Fokker-Planck equation (63), respectively. All three curves in Fig. 2 start at S(0) = 1/2, since ω2 /2ne = 1/2, but Une = 0 because U(π ) = 0. As in Fig. 1, all three approaches deliver similar S(t) for a moderate bath memory (panel (a)), but the Markovian Fokker-Planck equation is significantly off for the long-memory bath (panel (b)). In this latter case, the exact S(t) exhibits quite irregular oscillations, which are caused by the bath memory and inertial effects. Fig. 3 offers a more detailed view on the energy relaxation in a long-memory bath. It provides the partitioning of
214109-7
Maxim F. Gelin
J. Chem. Phys. 141, 214109 (2014)
0.55
1
(a)
(a) 0.9
0.5 0.45
kin
(t)
S(t)
0.8
S
0.7
0.4
0.6
0.35 0.5 0
2
6
4
0.3 0
t
10
30
20 t
1
0.5
(b)
(b)
0.9
0.4 0.3
pot
(t)
S(t)
0.8
S
0.7
0.2
0.6
0.1 0.5 0
10
20
30
t
0 0
20
10
30
t FIG. 2. Non-equilibrium energy relaxation functions S(t) evaluated for a moderate system-bath coupling ξ = 1. Panel (a): τ B = 1. Panel (b): τ B = 10. Green lines: exact results; black lines: Markovian Fokker-Planck equation; red lines: Markovian master equation.
S(t) depicted in Fig. 2(b) into the kinetic energy part Skin (t) = ω(t)2 /2ne (panel (a)) and the potential energy part Spot (t) = U(φ(t))ne (panel (b)). Both panels confirm that the Markovian master equation reproduces the timescale of the energy relaxation (note that Spot (t) is somewhat faster than for Skin (t)), while the Markovian Fokker-Planck equation does not. Parenthetically, panel (a) exhibits an interesting effect. Despite the rotator has equilibrium initial distribution over its angular velocities (hence Skin (0) = 1/2), it loses its kinetic energy at t < 1 and regains it at t > 1, reaching again the equilibrium value 1/2 at t → ∞. This is a manifestation of the energy exchange with the bath in the course of non-equilibrium relaxation. This behavior cannot be caused by the transformation of the rotator kinetic energy into the potential energy. The latter increases (on the average), starting from Spot (0) = 0 (Fig. 3(b)). VIII. APPLICATION OF THE APPROXIMATE MASTER EQUATION
It is appropriate to discuss here perspectives of possible applications of the Markovian master equation (50). The master equation fulfills all necessary formal requirements, in particularly normalization and detailed balance. In addition, it is exact in several particular cases. In the limit τ B → 0, it yields the Markovian Fokker-Planck equation. The same Fokker-Planck equation is recovered in the case of weak system-bath coupling for any g(t). The master equation (50) can certainly be applied beyond these particular cases. Evaluating its performance, we
FIG. 3. Non-equilibrium kinetic energy relaxation function Skin (t) (panel (a)) and potential energy relaxation function Spot (t) (panel (b)) calculated for a moderate system-bath coupling ξ = 1 and the bath correlation time τ B = 10. Green lines: exact results; black lines: Markovian Fokker-Planck equation; red lines: Markovian master equation.
have to realize that Eq. (50) possesses a set of relaxation rates (30), which produce exponential decay of the momentum correlation functions (31). Hence, the Markovian master equation (50), as any Markovian master equation, can yield non-monotonous (possibly oscillatory) time correlation functions in the case of weak system-bath coupling only, owing to the contribution of the system Liouvillian LS (7). If the system-bath coupling is moderate or strong, Markovian master equations usually yield exponentially decaying correlation functions. On the other hand, bath memory effects can produce oscillatory correlation functions even in the case of strong system-bath coupling. Obviously, the Markovian master equation (50) cannot reproduce oscillatory correlation functions in the regime where τ B > τ S (see, e.g., Figs. 1(b) and 2(b)). However, it can grasp the timescale of the oscillatory relaxation, delivering the correlation functions which are smoothened over oscillations. The point is that the relaxation rates of the Markovian master equation depend significantly on the bath memory. This substantially extends the domain of validity of the Markovian master equation in comparison with the Markovian Fokker-Planck equation, whose relaxation rates are solely determined by the friction.36, 37 IX. CONCLUSIONS
We have considered a classical particle bilinearly coupled to a harmonic bath. Assuming that the evolution of the
214109-8
Maxim F. Gelin
particle is monitored on the timescale which is longer than the characteristic bath correlation time, we have derived a Markovian master equation for the probability density of the particle. The relaxation operator of this master equation has been evaluated analytically without invoking the weak system-bath coupling approximation. We have found a simple method which avoids the Kramers-Moyal-like expansions of the relaxation operator (which require tedious calculations9, 15 ) and is valid for strong system-bath coupling. When the bath correlation time tends to zero or the weak system-bath coupling is considered, the Markovian Fokker-Planck equation is recovered. For a finite bath correlation time, the relaxation operator contains contributions of all orders in the system-bath coupling. The master equation (50) requires the correlation function C1 (t) of the momenta of the particle as an input. It thus permits the evaluation of various transport coefficients38 and reaction rates39, 40 starting from some first principles. One can proceed as follows: (i) retrieve C1 (t) from, e.g., molecular/Brownian dynamics simulations; (ii) compute the corresponding integral relaxation times τ n via Eq. (41); (iii) use the τ n for the construction of the master equation (50). The approach developed in the present work can straightforwardly be generalized to the construction of quantum Markovian master equations in the Wigner representation. ACKNOWLEDGMENTS
This work has been supported by the Deutsche Forschungsgemeinschaft (DFG) through a research grant and through the DFG Cluster of Excellence “Munich Centre for Advanced Photonics” (www.munich-photonics.de). I am grateful to Wolfgang Domcke and Danny Kosov for stimulating discussions. 1 I.
Prigogine, Non-equilibrium Statistical Mechanics (Interscience Publishers, New York, 1962). 2 R. Zwanzig, Nonequilibrium Statistical Mechanics (Oxford University Press, Oxford, 2001). 3 U. Weiss, Quantum Dissipative Systems (World Scientific, Singapore, 1999).
J. Chem. Phys. 141, 214109 (2014) 4 Y.
Tanimura, J. Phys. Soc. Jpn. 75, 082001 (2006). Rivas, A. D. K. Plato, S. F. Huelga, and M. B. Plenio, New J. Phys. 12, 113032 (2010). 6 Y. J. Yan and R. X. Xu, Annu. Rev. Phys. Chem. 56, 187 (2005). 7 M. F. Gelin, D. Egorova, and W. Domcke, Phys. Rev. E 84, 041139 (2011). 8 H. Spohn, Rev. Mod. Phys. 53, 569 (1980). 9 G. T. Evans, Mol. Phys. 36, 49 (1978). 10 V. Romero-Rochin and I. Oppenheim, J. Stat. Phys. 53, 307 (1988). 11 A. G. Dijkstra and Y. Tanimura, Philos. Trans. R. Soc. A 370, 3658 (2012). 12 S. Banerjee and A. Dhar, Phys. Rev. E 73, 067104 (2006). 13 P of Eq. (11) involves the system coordinate x and differs from its “purely bath” counterpart PB = P|x=0 normally used in the literature. 14 P. Resibois and M. de Leener, Classical Kinetic Theory of Fluids (John Wiley and Sons, New York, 1977). 15 W. Coffey, M. Ewans, and P. Grigolini, Molecular Diffusion and Spectra (John Wiley and Sons, New York, 1984). 16 A. P. Blokhin and M. F. Gelin, Physica A 251, 469 (1998). 17 If at least one of ν < 0, then ρ(p, t) would diverge at long times. α 18 M. F. Gelin and D. S. Kosov, J. Chem. Phys. 125, 224502 (2006). 19 S. A. Adelman, J. Chem. Phys. 64, 124 (1976). 20 P. Hänggi, Z. Phys. B 31, 407 (1978). 21 H. Risken, The Fokker-Planck Equation: Methods of Solution and Applications (Springer Verlag, Berlin, 1996). 22 G. Goodyear and R. M. Stratt, J. Chem. Phys. 105, 10050 (1996). 23 G. Goodyear and R. M. Stratt, J. Chem. Phys. 107, 3098 (1997). 24 M. F. Gelin and D. S. Kosov, J. Chem. Phys. 130, 134502 (2009). 25 M. R. Hoare, Adv. Chem. Phys. 20, 135 (1971). 26 F. Marchesoni and P. Grigolini, J. Chem. Phys. 78, 6287 (1983). 27 R. Zwanzig, J. Chem. Phys. 86, 5801 (1987). 28 P. Siegle, I. Goychuk, P. Talkner, and P. Hänggi, Phys. Rev. E 81, 011136 (2010). 29 E. Praestgaard and N. G. van Kampen, Mol. Phys. 43, 33 (1981). 30 W. Domcke and G. Stock, Adv. Chem. Phys. 100, 1 (1997). 31 J. Krˇ cmáˇr, M. F. Gelin, and W. Domcke, Chem. Phys. 422, 53 (2013). 32 A. S. Moskun, A. E. Jailaubekov, S. E. Bradforth, G. Tao, and R. M. Stratt, Science 311, 1907 (2006). 33 G. Tao and R. M. Stratt, J. Chem. Phys. 125, 114501 (2006). 34 G. Tao and R. M. Stratt, J. Phys. Chem. B 112, 369 (2008). 35 M. F. Gelin and D. S. Kosov, J. Chem. Phys. 126, 144511 (2007). 36 The so-called Fokker-Planck equation with memory describes nonMarkovian effects.37 Unfortunately, the equation is lacking microscopic justification. 37 A. P. Blokhin and M. F. Gelin, Physica A 229, 501 (1996). 38 M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids (Clarendon Press, Oxford, 1992). 39 P. Hänggi, P. Talkner, and M. Borkovec, Rev. Mod. Phys. 62, 251 (1990). 40 A. I. Burshtein, Adv. Chem. Phys. 129, 105 (2004). 5 A.
The Journal of Chemical Physics is copyrighted by the American Institute of Physics (AIP). Redistribution of journal material is subject to the AIP online journal license and/or AIP copyright. For more information, see http://ojps.aip.org/jcpo/jcpcr/jsp
Markovian master equation for a classical particle coupled with arbitrary strength to a harmonic bath Maxim F. Gelin Department of Chemistry, Technische Universität München, D-85747 Garching, Germany
(Received 20 March 2014; accepted 13 November 2014; published online 3 December 2014) We consider a classical point particle bilinearly coupled to a harmonic bath. Assuming that the evolution of the particle is monitored on a timescale which is longer than the characteristic bath correlation time, we derive the Markovian master equation for the probability density of the particle. The relaxation operator of this master equation is evaluated analytically, without invoking the perturbation theory and the approximation of weak system-bath coupling. When the bath correlation time tends to zero, the Fokker-Planck equation is recovered. For a finite bath correlation time, the relaxation operator contains contributions of all orders in the system-bath coupling. © 2014 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4902438] I. INTRODUCTION
II. NAKAJIMA-ZWANZIG MASTER EQUATION
The system-bath approach is one of the key paradigms in statistical mechanics. The idea is to derive the (irreversible) master equation for the probability density function (density matrix in quantum case) ρ(t) of the system starting from the combined system+bath (reversible) Hamiltonian dynamics. If the characteristic time of the bath-induced correlations is short on the timescale of the system dynamics, one can hope to derive a Markovian master equation1–8 of the form
Let us consider a classical particle bilinearly coupled to a classical harmonic bath. The total Hamiltonian H can be partitioned into the system (S) Hamiltonian, the bath (B) Hamiltonian, and their coupling,
∂t ρ = (−iLS + R)ρ. Here, the Liouville operator LS describes the time evolution of the system alone, and R is responsible for the bath-induced relaxation. Even if we restrict our consideration to classical mechanics, it is not an easy task to obtain the explicit form of the relaxation operator R. In general, one has to resort to additional approximations and/or assumptions, those like weak systembath coupling or delta-correlated Brownian forces (FokkerPlanck equations), low-density approximation (Boltzmann equations), mean-field approximation (Vlasov equations).8 The aim of the present work is to show that the relaxation operator R can be analytically evaluated without additional approximations in the particular, but important case of a classical particle bilinearly coupled to a harmonic bath. The paper is organized as follows. To setup the stage, we introduce the Nakajima-Zwanzig (NZ) master equation in Sec. II and discuss its Markovian limit in Sec. III. The explicit form of the relaxation operator R in the absence of an external force is derived in Sec. IV. General properties of the obtained master equation are discussed in Sec. V. Section VI outlines a generalization of the master equation. The predictions of the generalized Markovian master equation are numerically tested against the exact results in Sec. VII, and its possible applications are discussed in Sec. VIII. Section IX briefly summarizes our main findings. 0021-9606/2014/141(21)/214109/8/$30.00
H = HS + HB + HSB .
(1)
In our case, HS =
p2 + U (x), 2m
(2)
where x, p, and m are the position, momentum, and mass of the particle, and U(x) is an external potential. The system-bath coupling and the bath are described by 2 1 pi2 ci 2 + mi ωi xi − x , (3) HSB + HB = 2 i mi mi ωi2 where xi , pi , mi , and ωi are the positions, momenta, masses, and frequencies of the bath oscillators, and ci are the systembath coupling coefficients. The influence of the bath on the dynamics of the particle is determined by the bath correlation function2, 3 g(t) =
1 ci2 cos{ωi t}. m i mi ωi2
(4)
The Liouville equation for the total system+bath probability density function ρ tot (t) reads ∂t ρtot (t) = −iLρtot (t), where the total Liouvillian is defined as p i iL ≡ iLS + F∂p + ∂x − F i ∂p . i mi i i
(5)
(6)
Here,
141, 214109-1
iLS ≡
p ∂ + F (x)∂p m x
(7)
© 2014 AIP Publishing LLC
214109-2
Maxim F. Gelin
J. Chem. Phys. 141, 214109 (2014)
is the Liouvillian corresponding to the Hamiltonian HS , F(x) ≡ −dU(x)/dx is the external force, and ci ci Fi , Fi = xi − x , (8) F= mi ωi2 i
which guarantees that the particle attains the Boltzmann equilibrium distribution H (17) ρS = ZS−1 exp − S , kB T
where F is the total force exerted on the system by the bath oscillators. The initial condition to Eq. (5) is considered in the form
ZS being the partition function. We also define the reduced probability density function ∂t ρ(p, t) ≡ dxρ(x, p, t), (18)
ρtot (0) = ρ0 ρeq .
(9)
Here, ρ 0 is a certain distribution in the system phase space and H + HB −1 (10) exp − SB ρeq = Zeq kB T is the equilibrium bath distribution at the presence of the system, Zeq is the partition function, kB is the Boltzmann constant, and T is the bath temperature. Note that ρ eq of Eq. (10) is x-dependent and accounts for the initial system-bath correlations. We prefer to work with this distribution due to the following two reasons. First, it is more physical than the (commonly used) product system-bath distribution ρ 0 ρ eq (x = 0). A general discussion of the importance of system-bath correlations can be found in Refs. 9–11. Second, the correlated initial distribution (9) guarantees that the stochastic force f(t) generated by the Hamiltonian (1) in the generalized Langevin equation (Sec. IV C) possesses zero mean, f(t) = 0. See Chap. 3 of the monograph in Ref. 3 and Ref. 12 for the discussion of these issues. We introduce the operator9, 10, 13 P ≡ ρeq dpdx (11) (p and x denote collectively the momenta and positions of the bath oscillators) which projects the total system+bath probability density function into the correlated initial system-bath distribution (10). The application of the NZ formalism2, 3, 9 to Eq. (5) yields then a formally exact non-Markovian master equation for the probability density function of the particle, t dt R(t − t )ρ(x, p, t ). ∂t ρ(x, p, t) = −iLS ρ(x, p, t) + 0
which describes relaxation of the particle in momentum space. If the external potential is zero (F(x) = 0), the integration of Eq. (12) over x yields an exact NZ master equation ∂t ρ(p, t) =
0
where
t
dt Rp (t − t )ρ(p, t ),
(19)
Rp (t) ≡
dxR(t)|F (x)=0 .
(20)
III. NZ MASTER EQUATION IN THE MARKOVIAN LIMIT
Let us now assume that the characteristic decay time τ B of the bath correlation function (4) is short on the timescale of the system dynamics, τ S . If we are interested in the evolution of the particle for t > τ B , we can extend the time integration in Eq. (12) to infinity and expand the probability density in the Taylor series, ρ(p, t ) = ρ(p, t) + ∂t ρ(p, t)(t − t) + · · · ,
(21)
τ B /τ S being an effective small parameter of the expansion.14 In the leading order in τ B /τ S , we replace ρ(p, t ) by ρ(p, t) and arrive at the Markovian master equation ∂t ρ(x, p, t) = (−iLS + R)ρ(x, p, t).
(22)
Similarly, the Markovian analogue of Eq. (19) reads ∂t ρ(p, t) = Rp ρ(p, t). Here,
∞
R≡ 0
dt R(t ), Rp ≡
∞ 0
dt Rp (t )
(23)
(24)
(12) The relaxation operator in Eq. (12) can be written as a generalized Fokker-Planck operator,9 p . (13) R(t) = ∂p G(t) ∂p + mkB T
are the Markovian relaxation operators. In general, the evaluation of the operator R requires additional approximations (see Refs. 9, 15, 16, and Sec. VI). Rp , on the other hand, can be calculated analytically (see below).
Here,
IV. EVALUATION OF THE REDUCED RELAXATION OPERATOR
G(t) =
dpdxFe−i(1−P)Lt Fρeq
(14)
is the force-force correlation operator. According to Eq. (13), the relaxation operator obeys normalization dpdxR(t) = 0, (15) which insures the conservation of probability ( dpdxρ(x, p, t) ≡ 1) and detailed balance R(t)ρS = 0,
(16)
The explicit evaluation of Rp by Eqs. (13), (20) and (24) looks difficult. Fortunately, there exists an alternative way for achieving this goal. The derivation involves three steps. A. Step No. 1
Let A(p) and B(p) be arbitrary functions of the momentum of the particle. The correlation function CAB (t) = A(p) exp(Rp t)B(p)
(25)
214109-3
Maxim F. Gelin
J. Chem. Phys. 141, 214109 (2014)
(. . . denotes equilibrium averaging) can be calculated through the reduced probability density function (18) as (26) CAB (t) = dpA(p)ρ(p, t),
(F(x) = 0) it reads d p(t) = − dt
with ρ0 (p) =
1 2π mkB T
exp −
p2 2mkB T
(27) .
The integral relaxation time of CAB (t) is defined as ∞ dtCAB (t). τAB ≡
t
dt g(t − t )p(t ) + f (t),
0
(28)
Here, f(t) is a Gaussian colored noise with zero mean (f(t) = 0), and the memory function g(t) is defined by Eq. (4). The generalized Langevin equation (34), in turn, is equivalent to the non-Markovian Fokker-Planck equation19, 20 p ρ(p, t). (35) ∂t ρ(p, t) = a(t)∂p ∂p + mkB T Here,
(29)
a(t) = −
0
Obviously, the CAB (t) evaluated through the NZ master equation (20) and the Markovian master equation (23) are different. However, their integral relaxation times τ AB are identical. This is immediately proven by the direct time integration of Eqs. (20) and (23) and the use of Eq. (24). B. Step No. 2
The relaxation operator R(t) obeys the detailed balance (16), which simplifies to Rp ρ0 (p) = 0. We assume that Rp is Hermitian, and its eigenvalues are all negative, except the eigenfunction ρ 0 (p) possessing zero eigenvalue.17 For compactness, we adopt the bra-ket notation Rp |α = −να |α , α = 0, 1, 2, . . .
(30)
(ν 0 = 0 ensures the detailed balance, ν α > 0 for α = 1, 2, . . . ). The Hermiticity of Rp is difficult to prove in full mathematical rigor (cf. Ref. 18), but the requirement is reasonable from physical point of view and will be demonstrated a posteriori in Sec. IV C. To reveal the physical meaning of ν α , we specify the correlation function α(0) |α(t) ≡ α| exp(Rp t) |α = exp(−να t).
(31)
The integral relaxation time of this correlation function is defined as ∞ dt α(0) |α(t) = να−1 , α = 1, 2, . . . . (32) Tα ≡
∞ α=1
|α
1 α| Tα
(36)
The eigenfunctions are proportional to the Hermite polynomials (see, e.g., Ref. 14) 1 p p p|n = ρ0 (p)Hn , k|p = Hk n! 2mkB T 2mkB T (38) and are orthonormal, with the scalar product defined as k |n ≡ dp k |p p |n = δkn . (39) We then specify the fundamental correlation function t n(0) |n(t) ≡ Cn (t) = exp{−nb(t)}, b(t) = dt a(t ) 0
(40) (p(t) is a Gaussian process, hence Cn (t) = (C1 (t))n ). The integral relaxation time of this correlation function reads ∞ dtCn (t) (41) τn = 0
(the term with n = 0, which corresponds to the equilibrium distribution (28), yields τ 0 = ∞). Therefore, according to Eq. (33), Rp ≡ −
Hence, we can write
dC1 (t)/dt C1 (t)
is the time-dependent friction and C1 (t) ≡ p(0)p(t) is the momentum correlation function. As in Sec. IV B, we adopt the bra-ket notation and introduce the eigenfunctions of the Fokker-Planck operator, p |n = −n |n . (37) ∂p ∂p + mkB T
0
Rp ≡ −
(34)
f (t)f (t ) = g(t − t ).
where ρ(p, t) obeys Eq. (23) with the initial condition ρ(p, 0) = ρ0 (p)B(p)
∞
|n
n=1
1 n| . τn
(42)
(33)
(T0−1 = 0 due to the detailed balance).
V. GENERAL PROPERTIES OF THE MASTER EQUATION
C. Step No. 3
The relaxation operator (42) can be represented in several equivalent forms. Let us introduce the rate ν and the relaxation cross-sections σ n according to the expression
The representation (33) of the operator Rp can be easily constructed. The Hamiltonian (1)–(3) is equivalent to the generalized Langevin equation.2, 3 If the external potential is zero
1 ≡ ν(1 − σn ), τn
(43)
214109-4
Maxim F. Gelin
J. Chem. Phys. 141, 214109 (2014)
where σ 0 ≡ 1 (the decomposition (43) is not unique, of course). Using the definition of the bra and ket eigenfunctions (38) and the orthonormality relation (39), we can explicitly write Rp as an integral operator and recast the master equation (23) in the standard integro-differential form (44) ∂t ρ(p, t) = −νρ(p, t) + ν dp T (p|p )ρ(p , t). The relaxation kernel is defined as ∞ T (p|p ) ≡ σn p|nn|p .
(45)
n=0
Interestingly, the master equation (44) describes a Markovian, but non-Gaussian process, in contrast to the non-Markovian Fokker-Planck equation (35) and the generalized Langevin equation (34), which describe a Gaussian, but non-Markovian process. Alternatively but equivalently, the Markovian master equation (44) can be written as a Kramers-Moyal expansion21 ∂t ρ(p, t) =
∞ n=1
Mn (p) = τn−1
∂pn Mn (p)ρ(p, t),
B. ∂ x -dependence of R
Physically, ∂ x accounts for the change of the position of the particle during the bath correlation time τ B . Neglecting the ∂ x -dependence is tantamount to assuming that the particle cannot significantly move during τ B (see the discussion in Ref. 24). The assumption is reasonable, and its validity is similar to that of the Markovian approximation, τ B /τ S 1. If the ∂ x -dependence of R cannot be neglected, the Markovian approximation itself is questionable. C. Approximate master equation
(46) dp (p − p )n p |n n|p .
Equation (46) is obtained from the master equation (44) by the expansion of ρ(p , t) in Taylor series and integrations by parts.21 In the limit τ B → 0, only two first terms in the Kramers-Moyal expansion, M1 (p) = ξ (mkB T)−1 p and M2 (p) = ξ , remain nonzero (g(t) → δ(t), τ n → 1/(ξ n)). Then Eq. (46) reduces to the Markovian Fokker-Planck equation p ∂t ρ(p, t) = ξ ∂p ∂p + ρ(p, t). (47) mkB T For finite τ B , all Mn (p) = 0. Thus, the Markovian master equation (44) can be viewed as a tool for summing up the entire Kramers-Moyal series (46).
Above we have presented certain physical arguments supporting the idea that the relaxation operator R should not significantly depend on Fx and ∂ x . More solid arguments can be obtained by expanding the total Liouvillian L in Eq. (14) in powers of LS (cf. Ref. 9). It is straightforward to show that τ B /τ S is an effective small parameter of this expansion, since the force-force operator G(t) of Eq. (14) is proportional to the bath correlation function g(t) of Eq. (56) in the leading (zero) order in τ B /τ S . Then both the F(x)- and the ∂ x -dependence of the relaxation operator can be dropped and we can identify R with Rp R(F (x ); ∂x , p, ∂p ) ≈ R(0; 0, p, ∂p ) ≡ dxR(0; ∂x , p, ∂p ) ≡ Rp .
(49)
Adopting the approximation of Eq. (49), we can write the Markovian master equation (22) as follows:
VI. MASTER EQUATION FOR ρ(x, p, t)
The reduced probability density ρ(p, t) has a rather limited range of applications. One normally deals with the master equation (22) for the full probability density ρ(x, p, t) of the particle in an external potential U(x). According to Eqs. (13) and (14), the relaxation operator R in the master equation (22) is, in general, a functional of the external force F(x), and also depends on the operators ∂ x , p, ∂ p R ≡ R(F (x ); ∂x , p, ∂p ).
on the timescale τ B , the dependence of R on F(x) can be neglected. This should be a reasonable approximation in many practically relevant cases, since the short-time dynamics of the particle is predominantly determined by the (strong) interaction with its nearest neighbors.22, 23 We emphasize that the F(x)-independence of R does not mean that the influence of the external force is neglected, since F(x) enters the Markovian master equation through the system Liouvillian LS .
(48)
Due to the isotropy of space, R cannot depend on x explicitly, but only through the external force F(x) (this follows immediately from Eqs. (8) and (14)). Let us consider now possible approximations to Eq. (48). A. F(x)-dependence of R
According to Eq. (14), F(x) enters R through the projected total Liouvillian (1 − P)L. If the impact of the external force F(x) is smaller than that of the bath-induced forces F
∂t ρ(x, p, t) = (−iLS + Rp )ρ(x, p, t).
(50)
As is argued above, the validity of Eq. (49) is closely related to the validity of the Markovian approximation itself. It is thus not surprising that virtually all Markovian master equations available in the literature are of the form of Eq. (50), with the relaxation operators acting Rp acting in the momentum space only.1–6, 8, 25 If necessary, the corrections due to F(x ) and ∂ x can be evaluated perturbatively (cf. Refs. 9 and 15). VII. TESTING THE APPROXIMATE MASTER EQUATION
The Hamiltonian (1)-(3) generates the nonlinear generalized Langevin equations d x(t) = p(t)/m, dt t d p(t) = F (x(t))/m − dt g(t − t )p(t ) + f (t). dt 0
(51) (52)
214109-5
Maxim F. Gelin
J. Chem. Phys. 141, 214109 (2014)
If the memory function g(t) can be represented as a linear combination of N exponentials, Eqs. (51) and (52) can exactly be mapped into the Markovian Fokker-Planck equation involving x, p, and N additional auxiliary variables.26–28 Hence, the simplest (N = 1) exact Fokker-Planck equation depends on three variables. Below we study, numerically, whether the master equation (50) delivers a good approximation to the exact three-variable Fokker-Planck equation.
where the bra and ket eigenfunctions are proportional to the Hermite polynomials √ 1 ω |n = √ exp{−ω2 /2}Hn (ω/ 2), 2π n! (60) √ n |ω = Hn (ω/ 2). The relaxation times are specified by Eqs. (40) and (41). For the exponential memory function (56), they can be evaluated analytically
A. Formulation of the problem
We consider one-dimensional rotator (specifying by the angle φ and the angular velocity ω) in a N-well potential
τn =
U (φ) = {1 − cos(N φ)}. (53) N The rotator is subject to stochastic torques, so that its motion is described by the rotational analogue of the generalized Langevin equations (51) and (52)
δ≡
d φ(t) = ω(t), (54) dt t d d ω(t) = − U (φ(t)) − dt g(t − t )ω(t ) + f (t). dt dφ 0 (55) In this section, we adopt the dimensionless variables: Time, angular velocity, and energy are expressed in units of τ r , τr−1 , and kB T, respectively. Here, τr = I /(kB T ) is the period of free rotation, and I is the moment of inertia of the rotator. Our choice of the system specified by Eqs. (54) and (55) is stipulated by two reasons. Fundamentally, Eqs. (54) and (55) define a paradigmatic model of rotational relaxation in liquids, which mimics rotation of a molecule in the cage formed by its nearest neighbors (see, e.g., Refs. 15 and 29). Technically, Eqs. (54) and (55) are essentially nonlinear, and a proper description of the ensuing dynamics is a significant challenge for any approximate method. We restrict ourselves to the consideration of the exponential bath correlation function ξ exp{−t/τB }. (56) g(t) = τB Here, ξ and τ B are the strength and the correlation time of the Gaussian torque f(t). B. Rotational master equation
Equations (53)-(56) generate the rotational counterpart of the Markovian master equation (50)
∂t ρ(φ, ω, t) = − iLr + RME ρ(φ, ω, t) (57) r (the subscript “r” stands for “rotator”). The bath-free evolution of the rotator is described by the Liouvillian iLr = ω∂φ −
dU (φ) ∂ . dφ ω
(58)
According to Eq. (42), the relaxation operator is defined as ≡− RME r
∞ n=1
|n
1 n| , τn
(59)
n 2τB n (δ + 1)m (δ − 1)n−m Cm , (2δ)n m=0 n + (n − 2m)δ
1 − 4ξ τB , Cmn ≡
n! . m!(n − m)!
(61) (62)
In the limit τ B → 0, we recover the Markovian Fokker-Planck result τ n = 1/(ξ n). Explicitly, the corresponding FokkerPlanck equation is written as
(63) ∂t ρ(φ, ω, t) = − iLr + RFP r ρ(φ, ω, t), where
2 RFP r = ξ ∂ω + ∂ω ω .
(64)
For finite τ B , the integral relaxation times (61) are the ratios of polynomials in ξ τ B ; the degree of these polynomials increases with n. The perturbative evaluation of τ n with ξ τ B as a small parameter would require tedious calculations. C. Exact three-variable Fokker-Planck equation
The generalized Langevin equation (55) with the exponential memory function (56) can equivalently be represented as26–28 d d (65) ω(t) = − U (φ(t)) + ξ/τB q(t), dt dφ d q(t) = − ξ/τB ω − τB−1 q + fq (t), (66) dt where q(t) is an auxiliary variable and fq (t) is a deltacorrelated Gaussian process, fq (t)fq (t ) = 2τB−1 δ(t − t ).
(67)
The stochastic equations (55), (65) and (66) are equivalent to the Fokker-Planck equation
ρ(φ, ω, q, t) (68) ∂t ρ(φ, ω, q, t) = − iLr + REX r (the superscript “EX” stands for “exact”). Here, Lr is defined by Eq. (58) and
ξ/τB (ω∂q − q∂ω ) + τB−1 ∂q2 + ∂q q . (69) REX r =
D. Computational details
To solve Eq. (57), we expand the probability density into the series, ω |n e−ikφ ρnk (t) (70) ρ(φ, ω, t) = n,k
Maxim F. Gelin
J. Chem. Phys. 141, 214109 (2014)
(ω|n defined per Eq. (60) are proportional to the Hermit polynomials). Similarly, the solution of the exact threevariable Fokker-Planck equation (68) can be written as ω |n q |m e−ikφ ρnmk (t). (71) ρ(φ, ω, t) = n,m,k
1.1
The expansions (70) and (71) generate the matrix representation of the evolution operators
− iLr + RME and − iLr + REX (72) r r nk,n k nmk,n m k
n m k
facilitating and substantially speeding up the numerics (see, e.g., Refs. 30 and 31 for more details). E. Energy relaxation
There exist several interesting aspects of rotational and orientational relaxation which can be studied by invoking Eqs. (57) and (68). Here, we concentrate on the rotational energy relaxation. First, this problem has recently received considerable attention in the context of the rotational relaxation under non-equilibrium conditions and breakdown of the linear response description.18, 32–35 Second, relaxation of the rotational energy reveals importance of the bath memory effects, which cannot be grasped by the standard Markovian Fokker-Planck equation (63). Fig. 1 shows the equilibrium energy correlation function C(t) = E(0)E(t) , E = ω2 /2 + U (φ)
(73)
(. . . means equilibrium averaging). The system-bath coupling is moderate (ξ = 1). Panel (a) corresponds to a moderate bath memory (τ B = 1), while panel (b) corresponds to a long bath memory (τ B = 10). C(t) calculated via the exact three-mode Fokker-Plank equation (68), the Markovian master equation (58), and the Markovian Fokker-Planck equation (63) correspond to green, red, and black lines, respectively. We prefer to consider unnormalized correlation functions (vide infra), hence C(0) = 1. If the bath memory is moderate, all three descriptions grasp the essentials of the energy relaxation (panel (a)). The Markovian master equation is superior over the Markovian Fokker-Planck equation at short times t < 2, but no qualitative differences between the three approaches occurs. Panel (b) offers a qualitatively different picture. The Markovian Fokker-Planck equation is significantly off, predicting an order of magnitude shorter relaxation time (cf. black and green
1 0.9 0.8 0.7 0.6 0
1
3
2 t (b)
1.2 1.1 C(t)
and convert the master equations (57) and (68) into the linear matrix first-order differential equations for the “vectors” ρ nk (t) and ρ nmk (t), respectively. These matrix equations are solved by using the fourth-order Runge-Kutta integrator. In our calculations, we consider a two-well potential of Eq. (53) with N = 2 and = 1, so that the torque −dU/dφ = − sin (2φ). In this case, all matrix elements of Eq. (72) are evaluated analytically in closed form and in elementary functions. Furthermore, the structure of the matrices of Eq. (72) permits analytical evaluation of all matrix-vectors multiplication of the kind
− iLr + REX ρ (t), r nmk,n m k n m k
(a)
1.2
C(t)
214109-6
1 0.9 0.8 0.7 0.6 0
10
20
30
t
FIG. 1. Equilibrium energy correlation functions C(t) evaluated for a moderate system-bath coupling ξ = 1. Panel (a) corresponds to a moderate bath memory τ B = 1, while panel (b) corresponds to a long bath memory τ B = 10. Green lines: exact results; black lines: Markovian Fokker-Planck equation; red lines: Markovian master equation.
curves). The Markovian master equation performs much better. It cannot reproduce the oscillatory behavior of the exact C(t) (see the discussion below), but delivers a smoothened version of C(t) and correctly grasps its relaxation timescale (cf. red and green curves). Fig. 2 monitors equilibration of the rotational energy after the non-equilibrium preparation. It reproduces the mean energy evolution S(t) = E(t)ne .
(74)
Here, . . . ne means that the rotator at t = 0 possesses equilibrium Boltzmann distribution over its angular velocity ω, while the orientation angle is fixed at φ = π , which corresponds to the minimum of one of the potential wells. As in Fig. 1, panels (a) and (b) correspond to a moderate (τ B = 1) and long (τ B = 10) bath memory. Green, red, and black lines correspond to the exact three-mode Fokker-Plank equation (68), the Markovian master equation (58), and the Markovian Fokker-Planck equation (63), respectively. All three curves in Fig. 2 start at S(0) = 1/2, since ω2 /2ne = 1/2, but Une = 0 because U(π ) = 0. As in Fig. 1, all three approaches deliver similar S(t) for a moderate bath memory (panel (a)), but the Markovian Fokker-Planck equation is significantly off for the long-memory bath (panel (b)). In this latter case, the exact S(t) exhibits quite irregular oscillations, which are caused by the bath memory and inertial effects. Fig. 3 offers a more detailed view on the energy relaxation in a long-memory bath. It provides the partitioning of
214109-7
Maxim F. Gelin
J. Chem. Phys. 141, 214109 (2014)
0.55
1
(a)
(a) 0.9
0.5 0.45
kin
(t)
S(t)
0.8
S
0.7
0.4
0.6
0.35 0.5 0
2
6
4
0.3 0
t
10
30
20 t
1
0.5
(b)
(b)
0.9
0.4 0.3
pot
(t)
S(t)
0.8
S
0.7
0.2
0.6
0.1 0.5 0
10
20
30
t
0 0
20
10
30
t FIG. 2. Non-equilibrium energy relaxation functions S(t) evaluated for a moderate system-bath coupling ξ = 1. Panel (a): τ B = 1. Panel (b): τ B = 10. Green lines: exact results; black lines: Markovian Fokker-Planck equation; red lines: Markovian master equation.
S(t) depicted in Fig. 2(b) into the kinetic energy part Skin (t) = ω(t)2 /2ne (panel (a)) and the potential energy part Spot (t) = U(φ(t))ne (panel (b)). Both panels confirm that the Markovian master equation reproduces the timescale of the energy relaxation (note that Spot (t) is somewhat faster than for Skin (t)), while the Markovian Fokker-Planck equation does not. Parenthetically, panel (a) exhibits an interesting effect. Despite the rotator has equilibrium initial distribution over its angular velocities (hence Skin (0) = 1/2), it loses its kinetic energy at t < 1 and regains it at t > 1, reaching again the equilibrium value 1/2 at t → ∞. This is a manifestation of the energy exchange with the bath in the course of non-equilibrium relaxation. This behavior cannot be caused by the transformation of the rotator kinetic energy into the potential energy. The latter increases (on the average), starting from Spot (0) = 0 (Fig. 3(b)). VIII. APPLICATION OF THE APPROXIMATE MASTER EQUATION
It is appropriate to discuss here perspectives of possible applications of the Markovian master equation (50). The master equation fulfills all necessary formal requirements, in particularly normalization and detailed balance. In addition, it is exact in several particular cases. In the limit τ B → 0, it yields the Markovian Fokker-Planck equation. The same Fokker-Planck equation is recovered in the case of weak system-bath coupling for any g(t). The master equation (50) can certainly be applied beyond these particular cases. Evaluating its performance, we
FIG. 3. Non-equilibrium kinetic energy relaxation function Skin (t) (panel (a)) and potential energy relaxation function Spot (t) (panel (b)) calculated for a moderate system-bath coupling ξ = 1 and the bath correlation time τ B = 10. Green lines: exact results; black lines: Markovian Fokker-Planck equation; red lines: Markovian master equation.
have to realize that Eq. (50) possesses a set of relaxation rates (30), which produce exponential decay of the momentum correlation functions (31). Hence, the Markovian master equation (50), as any Markovian master equation, can yield non-monotonous (possibly oscillatory) time correlation functions in the case of weak system-bath coupling only, owing to the contribution of the system Liouvillian LS (7). If the system-bath coupling is moderate or strong, Markovian master equations usually yield exponentially decaying correlation functions. On the other hand, bath memory effects can produce oscillatory correlation functions even in the case of strong system-bath coupling. Obviously, the Markovian master equation (50) cannot reproduce oscillatory correlation functions in the regime where τ B > τ S (see, e.g., Figs. 1(b) and 2(b)). However, it can grasp the timescale of the oscillatory relaxation, delivering the correlation functions which are smoothened over oscillations. The point is that the relaxation rates of the Markovian master equation depend significantly on the bath memory. This substantially extends the domain of validity of the Markovian master equation in comparison with the Markovian Fokker-Planck equation, whose relaxation rates are solely determined by the friction.36, 37 IX. CONCLUSIONS
We have considered a classical particle bilinearly coupled to a harmonic bath. Assuming that the evolution of the
214109-8
Maxim F. Gelin
particle is monitored on the timescale which is longer than the characteristic bath correlation time, we have derived a Markovian master equation for the probability density of the particle. The relaxation operator of this master equation has been evaluated analytically without invoking the weak system-bath coupling approximation. We have found a simple method which avoids the Kramers-Moyal-like expansions of the relaxation operator (which require tedious calculations9, 15 ) and is valid for strong system-bath coupling. When the bath correlation time tends to zero or the weak system-bath coupling is considered, the Markovian Fokker-Planck equation is recovered. For a finite bath correlation time, the relaxation operator contains contributions of all orders in the system-bath coupling. The master equation (50) requires the correlation function C1 (t) of the momenta of the particle as an input. It thus permits the evaluation of various transport coefficients38 and reaction rates39, 40 starting from some first principles. One can proceed as follows: (i) retrieve C1 (t) from, e.g., molecular/Brownian dynamics simulations; (ii) compute the corresponding integral relaxation times τ n via Eq. (41); (iii) use the τ n for the construction of the master equation (50). The approach developed in the present work can straightforwardly be generalized to the construction of quantum Markovian master equations in the Wigner representation. ACKNOWLEDGMENTS
This work has been supported by the Deutsche Forschungsgemeinschaft (DFG) through a research grant and through the DFG Cluster of Excellence “Munich Centre for Advanced Photonics” (www.munich-photonics.de). I am grateful to Wolfgang Domcke and Danny Kosov for stimulating discussions. 1 I.
Prigogine, Non-equilibrium Statistical Mechanics (Interscience Publishers, New York, 1962). 2 R. Zwanzig, Nonequilibrium Statistical Mechanics (Oxford University Press, Oxford, 2001). 3 U. Weiss, Quantum Dissipative Systems (World Scientific, Singapore, 1999).
J. Chem. Phys. 141, 214109 (2014) 4 Y.
Tanimura, J. Phys. Soc. Jpn. 75, 082001 (2006). Rivas, A. D. K. Plato, S. F. Huelga, and M. B. Plenio, New J. Phys. 12, 113032 (2010). 6 Y. J. Yan and R. X. Xu, Annu. Rev. Phys. Chem. 56, 187 (2005). 7 M. F. Gelin, D. Egorova, and W. Domcke, Phys. Rev. E 84, 041139 (2011). 8 H. Spohn, Rev. Mod. Phys. 53, 569 (1980). 9 G. T. Evans, Mol. Phys. 36, 49 (1978). 10 V. Romero-Rochin and I. Oppenheim, J. Stat. Phys. 53, 307 (1988). 11 A. G. Dijkstra and Y. Tanimura, Philos. Trans. R. Soc. A 370, 3658 (2012). 12 S. Banerjee and A. Dhar, Phys. Rev. E 73, 067104 (2006). 13 P of Eq. (11) involves the system coordinate x and differs from its “purely bath” counterpart PB = P|x=0 normally used in the literature. 14 P. Resibois and M. de Leener, Classical Kinetic Theory of Fluids (John Wiley and Sons, New York, 1977). 15 W. Coffey, M. Ewans, and P. Grigolini, Molecular Diffusion and Spectra (John Wiley and Sons, New York, 1984). 16 A. P. Blokhin and M. F. Gelin, Physica A 251, 469 (1998). 17 If at least one of ν < 0, then ρ(p, t) would diverge at long times. α 18 M. F. Gelin and D. S. Kosov, J. Chem. Phys. 125, 224502 (2006). 19 S. A. Adelman, J. Chem. Phys. 64, 124 (1976). 20 P. Hänggi, Z. Phys. B 31, 407 (1978). 21 H. Risken, The Fokker-Planck Equation: Methods of Solution and Applications (Springer Verlag, Berlin, 1996). 22 G. Goodyear and R. M. Stratt, J. Chem. Phys. 105, 10050 (1996). 23 G. Goodyear and R. M. Stratt, J. Chem. Phys. 107, 3098 (1997). 24 M. F. Gelin and D. S. Kosov, J. Chem. Phys. 130, 134502 (2009). 25 M. R. Hoare, Adv. Chem. Phys. 20, 135 (1971). 26 F. Marchesoni and P. Grigolini, J. Chem. Phys. 78, 6287 (1983). 27 R. Zwanzig, J. Chem. Phys. 86, 5801 (1987). 28 P. Siegle, I. Goychuk, P. Talkner, and P. Hänggi, Phys. Rev. E 81, 011136 (2010). 29 E. Praestgaard and N. G. van Kampen, Mol. Phys. 43, 33 (1981). 30 W. Domcke and G. Stock, Adv. Chem. Phys. 100, 1 (1997). 31 J. Krˇ cmáˇr, M. F. Gelin, and W. Domcke, Chem. Phys. 422, 53 (2013). 32 A. S. Moskun, A. E. Jailaubekov, S. E. Bradforth, G. Tao, and R. M. Stratt, Science 311, 1907 (2006). 33 G. Tao and R. M. Stratt, J. Chem. Phys. 125, 114501 (2006). 34 G. Tao and R. M. Stratt, J. Phys. Chem. B 112, 369 (2008). 35 M. F. Gelin and D. S. Kosov, J. Chem. Phys. 126, 144511 (2007). 36 The so-called Fokker-Planck equation with memory describes nonMarkovian effects.37 Unfortunately, the equation is lacking microscopic justification. 37 A. P. Blokhin and M. F. Gelin, Physica A 229, 501 (1996). 38 M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids (Clarendon Press, Oxford, 1992). 39 P. Hänggi, P. Talkner, and M. Borkovec, Rev. Mod. Phys. 62, 251 (1990). 40 A. I. Burshtein, Adv. Chem. Phys. 129, 105 (2004). 5 A.
The Journal of Chemical Physics is copyrighted by the American Institute of Physics (AIP). Redistribution of journal material is subject to the AIP online journal license and/or AIP copyright. For more information, see http://ojps.aip.org/jcpo/jcpcr/jsp
