Unified approach to the calculation of inertial corrections in diffusion

Physica 125A (1984) 1-24 North-Holland, Amsterdam

UNIFIED APPROACH TO THE CALCULATION OF INERTIAL CORRECTIONS IN DIFFUSION Ulrich R. STEIGERt Department of Chemistry, University of California, Davis, California 95616, USA

Received 20 December 1982 Revised 22 August 1983

Inertial effects in diffusion due to the nonlinear equations of motion are investigated. The starting point is a Fokker-Planck equation for a general mechanical system with Lagrange function L =~gij(q l. . . . . qn)qiqj_ U(qt . . . . . qn) under the influence of random forces. The reduction to a description in position space is achieved by using the time ordered cumulant expansion and the boson operator representation. The first term of the expansion gives the Smoluchowski type diffusion equation. The next term leads to inertial corrections. This theory can be applied, for example, to the diffusion of N asymmetric molecules immersed in a fluid undergoing coupled translational and rotational diffusion and interacting via intermolecular and hydrodynamic forces. Comparison with the literature is presented.

I. Introduction Brownian motion can be formulated on different levels. A detailed description, which was first given by Klein 1) and Kramers2), uses a F o k k e r - P l a n c k equation to describe the time evolution o f the distribution function in the full phase space. According to Smoluchowski3), it is sufficient to consider the reduced distribution function on the position space only. The relationship between these two descriptions was unclear for a long time. During the last few years systematic procedures have been developed in order to derive the Smoluchowski equation, and generalizations o f it, from the F o k k e r - P l a n c k equation. Translational Brownian motion was investigated, for instance, by Wilemski4), TitulaerS), Steiger and Fox6). Wilemski used a hierarchy of equations in coordinate space for successively higher m o m e n t s of the velocity. Titulaer used an expansion procedure o f the C h a p m a n - E n s k o g type with the inverse friction coefficient playing the role of the expansion parameter. Steiger and Fox used the time ordered cumulant expansion combined with the boson operator representation. The advantage o f the last method lies in the fact that it permits writing down explicit expressions for all terms o f the expansion. These expressions become increasingly complicated the t Present address: FASELEC AG, Postfach, CH-8045 Ziirich, Switzerland. 0378-4371/84/$03.00 © Elsevier Science Publishers B.V. (North-Holland Physics Publishing Division)

2

U.R. STEIGER

further the expansion is pushed. One can show, by reordering some of the differential operators which appear, that many terms nearly cancel if the potential is sufficiently smooth. The problem of convergence of the time ordered cumulant expansion is related to similar considerations concerning the Mayer-series for the grand canonical partition function for the real gasT). Fox has shown that the time ordered cumulant expansion satisfies a time-cluster property which is analogous to the cluster property of the UrseU functions. The cluster property shows in particular that the higher order cumulants are small if the correlation time of the velocity is small. All higher order terms, except the second order one (the Smoluchowski term), vanish in the limit as ~--*0 where z denotes the correlation time of the velocity. The Smoluchowski equation is exact in the limit as z --*0. Corrections have been calculated which arise due to the fact that z > 0 in a physical situation. The different authors used various procedures; the most important ones are mentioned above. The results are in agreement, although the first method and the Chapman-Enskog procedure cannot be applied to short times t < z. The time ordered cumulant expansion is used in this work because valid results can be obtained for short and for large times and because all terms of this expansion are well defined. The cluster-property guarantees that the higher order cumulants are small if the correlation time is small and if the potential is sufficiently smooth. The same methods can be applied to rotational Brownian motion. The first theory of rotational Brownian motion was developed by Debye s) who investigated the response behavior of a dielectric material to an external electric field with circular frequency 09. Debye's theory did not agree with experimental results because he did not include inertial effects. These become observable at high frequency 09 with to >~z -~. During the last few years investigations by Ford et al.9), Hubbard t°) and McCorlnell ~L12)succeeded in including inertial effects for rotating Brownian particles of any shape in a systematic way. In a general theory translational and rotational Brownian motion cannot be treated independently, since translational and rotational motion of a large molecule moving in a dense fluid are coupled~3). Steiger and Fox ~4) have derived a generalized Smoluchowski equation for this case which applies to short and to large times. Deutsch et al. ~s) have investigated Brownian motion of several particles including hydrodynamic interactions. They calculated corrections of the diffusion constants of particles interacting via hydrodynamic forces. Higher order corrections were given by Titulaer~6). Hydrodynamic interactions are important for applications of Brownian motion in physical chemistry, for instance, for the study of the polymer dynamics and for the theory of diffusion controlled reactions. The purpose of this work is to present a formalism which includes all the corrections mentioned above. Such a unification is suggested by the fact that these

INERTIAL CORRECTIONS IN DIFFUSION

3

corrections can be viewed as inertial effects, as explained in the next section. The merit of an unified theory is twofold. Firstly, it leads to a better, more complete understanding of the interplay of Hamilton's dynamics and Brownian motion. Secondly, the unified theory permits studying new effects arising due to the simultaneous influence of external forces, the coupling between translations and rotations, inertial effects of the Brownian particles, and hydrodynamic interactions. The theory developed here is even more general. It can be viewed as a generalization of rotational Brownian motion which applies to any Lie group G with finite dimensional Lie algebra. Rotational Brownian motion is obtained in the case G = SO(3). If we choose for G the Euclidean group E we obtain translational and rotational Brownian motion. Hydrodynamic interactions are incorporated via a position dependent friction tensor.

2. Qualitative description of the four inertial effects (examples) The Smoluchowski equation can be viewed as the conservation law of a current which consists of a systematic part proportional to the external force and a dissipative part proportional to the gradient of the distribution function in coordinate space P(q), j = - D ( g r a d P + P grad U/kT).

(1)

D is the diffusion constant, T the temperature, k the Boltzmann constant and U = U(q) is the external potential. The equation of motion is -- P + divj = 0. Ot

(2)

2.1. Short time behavior It is impossible to write down a general diffusion equation in position space which would apply to short times t ~< ¢ because the time evolution depends also on the initial velocity distribution. The diffusion coefficient becomes time dependent in a corrected Smoluchowski equation. One can show that

D(t) = D®(1 - e -'#)

(3)

if the average velocity vanishes initially6). 2.2. Anisotropic diffusion coefficients Diffusion is not isotropic anymore in the presence of an external potential. It has been shown 6) that the diffusion constant D may be replaced by a matrix D U.

4

U.R. S T E I G E R

We have D+ ~ D 60 +

m OqiOqJ]

.

(4)

This anisotropy can be removed to a certain extent by means of the coordinate transformation T2

Q - q + VoZ + a ~ ,

(5)

with a being the acceleration of a particle with mass m due to the potential U; a = - m ~d U/dq. We have assumed that the average initial velocity vanished and set consequently v0 = 0. Eq. (5) is the orbit of a classical particle in a constant force field. This motion leads to a distortion of the Euclidean coordinate frame as described by the anisotropic diffusion matrix (4). If the diffusion equation is expressed in the new coordinates Q one gets back equations of the original form, eqs. (1) and (2), with isotropic diffusion coefficients. (This correspondence holds in first order of I1Q - q II/llqLI)The correction described by eq. (4) can be neglected if Q ~ q on a macroscopic length scale. 2.3. Rotational Brownian motion We consider rotational Brownian motion of a sphere with radius R and moment of inertia I in order to illustrate a type of inertial corrections which does not occur in the case of translational diffusion. We introduce the dimensionless quantity V, kTl

? = ¢~,

(6)

where ~R is the rotational friction constant ~l~ = 8•R3r/13), t/is the viscosity of the fluid. One finds experimentally that ~ is less than a few percent, in general. Ford et al. showed that there are non-vanishing corrections of Einstein's relation kT

DR = ¢-~

(7)

even in the absence of an external potentialH). The first correction is kT

+o,,, = ~ (l + ½y + . . . ) .

(8)

General expressions for particles of arbitrary shape have been calculated in ref. 11. Because rotations and translations are in general coupled these inertial corrections lead also to corrections of the translational diffusion coefficients. To our knowledge, nobody has investigated these effects.

INERTIAL CORRECTIONS IN DIFFUSION

2.4. H y d r o d y n a m i c

5

interactions

We consider a system of N interacting particles for which the hydrodynamic interactions are included through the presence of a position dependent interparticle friction tensor ~ = ~(qt . . . . . qN). Modified Smoluchowski equations describing these effects by a position dependent diffusion tensor have been given by several authorslS'16"17). Titulaer calculated higher order corrections for two spheres interacting via a hydrodynamic force and an intermolecular force. Hydrodynamic interactions are long ranged forces due to inertial effects of the fluid. The higher the density the larger the spatial correlations of the velocity field of the fluid. 3. Basic equations and definition of the coordinate frame

The time evolution of the system considered here consists of a systematic part and a random part. The first one is described by classical mechanics. We have n degrees of freedom with coordinates q = (ql, q2. . . . . q") and the Lagrange function (9)

L = ½guv~v j - U ,

where g,7 = gij(q) and U = U ( q ) . The Euler-Lagrange equations are equivalent to (10)

f? + F ~j~vJvk = - g U ( e y U ) ,

where the base vectors are defined by 0

(1 1)

ei = Z / a q j ,

with arbitrary coefficients X/= z / ( q ) 19.20,21).The velocity components are defined by q t = v jZj,,

(12)

where 0 ~ stands for the time derivative of the coordinate qt. The commutation coefficients c~ and the connection coefficients F~k are given by (1 3)

[ei, ej] = Cijkek , Fi

_

1

il,,

jk = i g tgtk, j + gtj,k - - gjk,l + Clkj + Ctjk - -

Cjkl)

,

(14)

where ctkj - c~k"g.j

(15)

g~gk:=

(16)

and 6/.

6

U.R. STEIGER

For vanishing potential U, eq. (10) is reduced to the geodesic equation in a noncoordinate base2°). The metric gij is constant in the following. Hence, we have as a consequence the symmetry relation F ok = -- F j,k ,

(17)

where F/jk = g~lF tjk. We are mainly interested in the translational and rotational Brownian motion of rigid particles, for which the infinitesimal generators of the Eucledian group are needed. We denote with 7"1, T2, T3 the infinitesimal translations along the three body fixed axis x, y, z, respectively. The infinitesimal rotations about the same axes are M1, M2, M3. These operators are given by Tk = Rjk(CP, O, 6 ) Oq j ,

(18)

where ql, q2, q3 denote the position of the center o f mass. ~b = q4, 0 --= q5, 6 -= q6 are the Eulerian angles and R(tp, 0, 6 ) is the rotation matrix which gives a vector in the body fixed frame when R(tp, 0, 6 ) is applied to the same vector in the laboratory fixed frame (compare ref. 21, page 242),

~ -- ctg 0 sin 6 ~a, M1 = cos 6 ~ + sin 6 sin1 0 0~ d 1 8 8 M2 = sin 6 ~-~ -- cos 6 sin~-0~--7 + ctg 0 cos 6 ~7z_,_, vql tpc

(19)

M3 =~-~. These operators fulfill the cummutation relations 2°'21'22)

[T,, r A = 0 , [Mi, Tj] = --E,jkT,,

(20)

[M,, Mj] = - c~kMk,

which can be verified by direct calculation. E~jkis the completely antisymmetric Levi-Citiva tensor with E~23= 1. The base vectors el to e6 are defined by el = Tl, e2 = T2, e3 = 7"3, e4 = M1, e5 = M2, e6 = M3.

(21)

The structure constants c~ can be written in terms of six 6 x 6 matrices defined by (Ck)ij = Ckij. We have

I N E R T I A L C O R R E C T I O N S IN D I F F U S I O N

C1=

0

0

0

'

(72=

0

0

-1

° ,°1' ° OlO

-1

/

o

0

0

--1

0

1

0

0

0

0

0

0

1

0

0

0

0

0

-1 0

0

0

0

0

'

C 3 =

oo

0

o

'

0

0

0

0

0

-1

1

0

1

0

-10

0 0

0 0

0

C 5=

--

7

(22)

0

0_10 c~=

1

0

0

00 0

0

0

0

--1

1

0

0

0

0

"

The Lagrange function of the rigid body is obtained by setting gu = Mu,

(23)

where

M is the generalized inertia matrix which contains the mass m and the tensor of inertia with respect to the center of mass L Eq. (10) becomes now M i k l ) k "1- ClkJMjiv iv k = __ elU.

The vector v contains the translational velocity 4 velocity in the body fixed coordinate frame, sin 0 sin ~k + 6I cos ~, o92 = q~ sin 0 cos ~k - 6I sin v = (¢,o). The equation above contains Euler's

(24) = (ql, q2, q3) and the angular o=(o91, o92,o93) with o91= ~b and o93 = ff + q~ cos 0 is). equation of the rigid body

IdJ + ta x I ~ = O.

4. The F o k k e r - P l a n c k

equation, boson operator representation

4.1. G e n e r a l i z e d L i o u v i l l e e q u a t i o n

The volume element in phase space is conserved by a Hamiltonian flow according to Liouville's theorem19). Using the coordinates qi and the canonical

8

U.R. STEIGER

conjugate momenta p; = dT/t3?l ~= ~ i l j this volume element can be written d# = d q l . . , dqn dpl • • • dpn = det(g/j) d q l . . , dqn dt~l.., dq".

(25)

Using the transformation (12) this equation can be written

d# = h(q) dq~ . . . dq" dvl . . . dv n.

(26)

The weight function h(q) is h (q) = det(x~t),

(27)

X,, = Z,%,.

(28)

h (q) defines the invariant measure of the group G, the Haar measure2°). In the case considered before, the motion of a rigid body, the transformation (12) with eqs. (27) and (28) gives (up to a constant factor)

h(q) = sin 0,

(29)

which is the Haar measure of the Euclidean group. We consider now distributions on f2, where t2 is the phase space consisting of the points (q~ . . . . . qn, v ~ , . . . , vn). P(t, S) denotes the probability of finding a particle in the subset S c t2 at time t. The distribution function f ( t , q~ . . . . . q", v ~. . . . . v") is defined by the following equation which holds for all Sct2: (30)

P(t, S) = f d # f ( t ) . t/ S

Because the probability is conserved as well as the volume element d#, the total time derivative along an orbit defined by the equation of motion (10) vanishes d

~ f (t, ql(t ) . . . . . qn(t ), v'(t) . . . . . v"(t)) = 0.

(31)

Thus, the generalized Liouville equation can be written

- ~ f + vieif - {F~kVJVk + (ejU)g ji} ~ i f --O.

(32)

4.2. The Fokker-Planck equation Klein and Karmers modified the Liouville equation in order to describe stochastic influences too. This is achieved by adding an extra term to eq. (32). The new term is ~'2)

K f = ~ ¢iJ(q)(kT t3 + gjtv')f. ~v ~

\

dvJ

(33)

INERTIAL CORRECTIONS IN DIFFUSION

9

The friction tensor ~ # depends on the position if hydrodynamic interactions are included. The friction tensor is the ensemble average of the random forces acting on the Brownian particle23,24), cO

~/J : ~T | dt (if/(0, 0)ffy(t,q ) ) ,

(34)

~d 0

(35)

~tm = ~qgilgym .

ff~(t, q) denotes the ith component of the random force at time t and position q. The Fokker-Planck equation is (36)

fit f = (L + K ) f , where the operator L is defined by eq. Kramers-Liouville equation, in the following.

(32).

Eq. (36) is called the

4.3. Dimensionless variables, boson operators One knows from equilibrium statistical mechanics that the velocity distribution relaxes to the Maxwell distribution

gM(v I ....

,

detl/2(giJ) exp( -- gijvivJ/2k T) v") = (2nkT),/2

(37)

The standard deviation of the velocity distribution of a point particle in thermal equilibrium is a = x/rk-T-/m. All expectations which can be calculated are a function of tr only because the velocity distribution is Gaussian with vanishing mean. The dimensionless variable v* = via is normally distributed indeed. Generally, the covariance of two velocity components with distribution (37) is H) cov(v i, v j) = kTg # •

(38)

We introduce dimensionless variables v *~. . . . . v*" by imposing the condition cov(v*i, v'J) = 6 U.

(39)

Thus, the components v */ are independent, identical, normally, distributed random variables. The matrix gU = (M-l)u is symmetric; hence one can define a square root of the matrix M -l, denoted by M -~/2, simply by diagonalizing this matrix and calculating the square roots of the diagonal elements. The new base vectors are

e *i = (kT) ]/2M- 1/2ijej

(40)

10

U.R. STEIGER

and the dimensionless velocity components are (41)

v ,i = (kT)-l/2Ml/2ifl)J.

Using these definitions the Lagrange function (9) can be written (42)

L" = ½6o.v*iv *j -- U /k T ,

where L ' = L / k T . The new metric is just the Kroenecker symbol 6 U.The equations (13), (14), (15), (40) and (42) determine the new commutation coefficients c*~ and the new connection coefficients F*~k, which can be obtained from %k and F~k, respectively, by a linear transformation. It is not necessary anymore to distinguish between upper and lower indices because the new metric is the identity matrix, C * / j k ~--- C,jk, *

F*~k=F*0k,

etc.

(43)

Thus, the Kramers-Liouville equation can be written a

~ f + v *ie * J - {r *'jkv*Jr ,k + (e *,U)/kT} c

~--~*'-fov

a

The matrix ~*° is defined by (45)

~ *ij __ ( b T ~ - 1 ~ f l / 2 ~ l r n l l A l / 2 --

\,vx

I

.rail

~

~.A mj

•

4.4. Boson operators Brownian motion can be viewed in three different ways: in terms of the Langevin equation, the corresponding Fokker-Planck equation and in terms of boson operators. Boson operators were introduced in this context independently by Titulaer 5) and Steiger et al.6). They are defined by a

a

a+ . Or*' . . a.= . Ov*

v*,

b +=e*,

b-,--e*-(e*U)/kT.

(46)

The symbol a + is a shorthand notation for a "vector" operator which has n components a + = (a T. . . . . a,t) = (a +'. . . . at").

(47)

The same remark applies also to the other operators defined by eqs. (46). These operators satisfy the following commutation relations:

[a,, a}] = 6,j,

[6~, hi] = c*kb~, [b~,bj] = --c*kbk,

[b,, b]] = c*~b]:+ Uj;. (48)

All other commutators vanish. The tensor U o is defined by UU= ( e ' e * U ) / k T and

INERTIALCORRECTIONSIN DIFFUSION

11

is not symmetric, in general, because the operators e* and e* do not have to commute. The n th order partial derivatives of the potential U can be written in terms of

* • . . e*U)/kT. Uili2. . . ik ==-(e il* e i2

(49)

The following bilinear mapping plays an important role in the calculation of the diffusion operator:

( f , g )(ql . . . q,) =_ f d61. . . d r , g ~ ( v ~ , . . . v") xf*(q l. . . . q", v 1. . . . vn)g(q I. . . . q", v 1. . . . v")

oxp/V. + x f * ( q l . . , q", D*1

. . . .

v * " ) g ( q l . . , q", v * l . . , v*").

(50)

The complex-valued functions f and g are defined on the whole phase space. The inverse of the Maxwell distribution gM, eq. (37), and the inverse of the standard normal distribution in the variables v'i, respectively, serve as the weight function. Eq. (50) can also be viewed as scalar product on the velocity space. However, it is useful to define a scalar product on the full phase space too,

( ( f , g ) ) = f dktg~l~(q, v)f*(q, v)g(q, v) U(q) =[ f d"qh(q)exp(~)l

fd"qh(q)exp(-~){f,g)(q).

(51)

gMa(q, 6) is the Maxwell-Boltzmann distribution. The operators a, a t, b, b t have several extraordinary properties. Firstly, they fulfill the commutation relations (48). Secondly, it is easy to verify that agM = agMa = bgMa = 0.

(52)

Thirdly, the operator a t is the formal adjoint of a and b t is the formal adjoint of b. We have

(af, g ) = ( f , a t g ) ,

(53)

((bf, g)) = ((f, btg)),

(54)

where f and g are functions with the property ( ( f , f ) ) < oo a n d ((g, g ) ) < oo. These equations can always be proved by integrating by parts. Eq. (54) would not be correct if one would not use the Haar measure h(q) in the definition of the scalar product ( ( f , g ) ) , eq. (51). In order to give a general proof of eq. (54) we

12

U.R. STEIGER

have to recall that the Haar measure is invariant under the group transformations. Thus, for such a transformation A, we have ~ d"qh e V / k r ( A f ) g = ~ & q h ( A - l eV/kr) f ( A - ~ g ) . We set A = e~* and calculate the derivative with respect to 2 at both sides of this equation. Eq. (54) is now obtained by setting 2 = 0. The properties of the operators a, a t, b, b + are summed up by the equations (48) and (52) to (54). These equations show that the operators atk and ak are creation and annihilation operators, respectively2S). This is an important result of the remarks above. Unfortunately, this does not apply to b and b+ because these operators do not fulfill the canonical commutation relations. The Kramers-Liouville equation, eq. (32) or (44), shows a simple structure if it is expressed in terms of the operators a, a t, b and b t. The Kramers operator can be written (55)

K = - ~*Ua~aj.

The Liouville operators can be split in two parts. We have L = L ~+) - L ~-) , (56)

L (+) = b t a + Fiykajaka * ti , L (-) = ba + + r*ka, a~a~.

Note that the ordering of the operators of the qubic terms (e.g. F~ka,a]aD is unimportant because the connection coefficients are antisymmetric in the first two indices, eq. (17), and because of the canonical commutations relations, eq. (48) and because trace Ck = 0 26). Thus, • +a ~ = Y ~ a ] a i a ~ = F qkaiaj •

t

*

t

* +a~a~, F ijkaj *

(57)

t

F uka iaja k = F ukaja ia k = F ukajaka i •

The operator L (-) is the formal adjoint operator of L (+). Hence, the Liouville operator, L = L ~+) - L ~-), is skew symmetric with respect to the scalar product (51). However, Kramers operator, eq. (55) is symmetric. The discovery of this structure of the Liouville operator, as described by eqs. (56), constitutes an important result although the proof is simple indeed. These equations are obtained in a straightforward way by applying the definitions of the operators a, a*, b and b t, eqs. (46), to eq. (44). Note that v* = - ( a + at). 4.5. H e i s e n b e r g p i c t u r e in velocity space Two different representations are used in quantum mechanics: the Schr6dinger representation and the Heisenberg representation. Similarly, the relaxation process of the velocities due to random forces can be investigated from these two

INERTIAL CORRECTIONSIN DIFFUSION

13

point of views. The Kramers equation ( O / d t - K ) f = 0 corresponds to the Schr6dinger picture. Here the states (the function f ) are time dependent. In the Heisenberg picture the operators are time dependent while the states are time independent. There is an operator O(t) corresponding to each observable. They satisfy the equation of motion d

Tt o(t) =

[o(t),/q

(58)

which replaces the Kramers equation used in the Schr6dinger picture. It is sufficient to consider operators O which can be written as a polynomial in the creation and annihilation operators a~ and ai, respectively. In this case, the time dependent operator O(t) is obtained simply by replacing the operators a ~ by the corresponding time dependent operators (where a ~ stands for either a or at). This follows from the fact that the sum and the product of two solutions of the equation of motion (58) are again a solution of this equation. One obtains

a*(t) = e¢'ta t,

a(t) = e-~*ta.

(59)

Using the commutation relations (48) one can easily verify that these equations are solutions of eq. (58) with the initial conditions a ~'(0)= a #. (Note that ~* is symmetric, compare eq. (45), thus e-~'rta = e-¢*ta.) In the next section equilibrium averages of operators acting on the velocity space are calculated in order to achieve the reduction to a description of Brownian motion in position space only. The (equilibrium) expectation value of an operator O(t) is defined by

( O ( t ) ) =
(60)

with the Maxwell distribution gM. O(t) is a polynomial in the operators a ~(t). In order to evaluate expressions of the form of eq. (60) all annihilation operators are moved to the right by using the commutation relations (48). One uses the fact that agr~ = 0. This procedure is analogous to Wick's normal ordering used to evaluate expectation values of field operators in quantum field theory25). By this means, one can show that all expectation values (O (t)) can be expressed in terms of the correlation function 2~(t).

~(t) - (a(t)a*(O)) =

e -¢''

.

(61)

The physical meaning of this expression is quite clear. £(t) is the time dependent velocity correlation matrix (expressed in the base e*).

U.R.STEIGER

14

4.6. Interaction picture NextweconsidertheLiouvilleoperatorintheinteractionpicture.Operatorsacting on the velocity space are time dependent. We define, as in ref. 6,

;

L”(t)= [L(t), zq.

(63)

The Liouville operator in the interaction picture t(t) is obtained from the expressions (56) by replacing the operators a # and b # by a “(8) and b *(t), respectively. The operators b “(1) are given by b?(t) = e-“bf

erK= b: &- t~~~~~(t~~(t).

WI

The minus sign holds for bf and the plus sign for b,, c& denotes the partial derivative with respect to the coordinate I, (8, zz eT and b #(tJ do not have to commute at unequal times although we have always [a”(t), b”(t)] =o. The Liouvifle operator

(66) in the interaction

picture can now be written

Z(t) = D+‘(t) - Z’-‘(t) f

(67)

B+)(t) = bf(l) * a(2) -f- ~~~Uj(~)ff~(~)~~~~),

631

i?‘(t)

= b(t) - at(t) i- ~~~~i~~~~~~z}~~~~~ e

05%

It also follows from eqs. (59) and (66) to (69) that ((L”‘+‘(0f, 8)) = ((St L”‘-‘( - tk>> -

(70)

Thus the operator Z’+)(t) can be obtained by calculating the formal adjoint of L”(-)(E) and inverting the direction of time.

5. General expression for the inertial corrections 5.1. Smohchawski

type d@usion equations

In order to calculate the diffusion operator G will use the same method as in ref. 6. C describes the time evolution of the reduced distribution function in

INERTIAL

CORRECTIONS

IN DIFFUSION

15

position space P ( t , q), Ot P = G P ,

(70)

G = ~ G ~") .

(71)

n=l

The time-ordered cumulant expansion o f G can be derived from the theorem on time ordered exponentials which is an extension of the theorem of exponentials widely used in renormalization group calculations. The more general version is needed because the conditional averages considered here are non-commuting operators. The cumulants G ~") are explicitly defined in ref. 6. For the following calculation the first two nonvanishing terms are needed. Because the odd cumulants vanish in our case the second and fourth cumulants are the important contributions. They are - using the notation introduced by eq. (60) l

6~2)(t) = { d s ( f f , ( t ) f f , ( s ) ) ,

(72)

i/ 0 tl a(4)(tl)

=

t2

t3

fdt2fdt3fdt4{(E(tl)E(t2)E(ts)E(t4))-(E(tOE(tz))(E(ts)E(t4)) 0 0 0 -

(ff,(tl)ff,(ts))(L(t2)ff,(t4))

- (ff,(tl)ff,(t4))(ff,(t2)ff,(t3))).

(73)

In order to calculate the equilibrium averages in the velocity space (/[(tl) • •. L(tn)) only the properties o f the boson operators described by eqs. (48), (52) and (53) are used. One obtains for instance (£(tOff,(t2)) = - (/~+~(tl)E~-)(t2)) =

-(bt(h)

"a(tOb(t2)'at(t2))

= - b*~ (a,(t,)a)(t2))bj = - b*~,o(tl - tz)bj.

(74)

%(t) is the velocity auto-correlation matrix defined by eq. (61). Note that 27 may depend also on the position if the hydrodynamic interactions do not vanish. The operators b '~ and 27 do, therefore, not commute in general. However, the first non vanishing cumulant is always G~2)(t) = - - b * . ¢*-1(1 - e-'¢')b.

(75)

16

U.R. S T E I G E R

Thus, using the equations (40), (45) and (46) the (first) diffusion equation can be written

Ot P(t, q) = e. D(1 - exp( - At))(e + eU/kT)P(t, q),

(76)

where

(O--1)i j

=

~o./k T

and

AO._ ~ij = ~Ug0 =

¢ilglj.

Eq. (76) is the generalization of the Smoluchowski equation which describes diffusion in an arbitrary mechanical system (described by the Euler-Lagrange equations (10)) under the influence of (strong) random forces. Hydrodynamic interactions are also included and the short time behavior is described correctly. This generalized Smoluchowski equation does not describe inertial effects. It describes diffusion in the limit as z tends to zero. It is assumed that the correlation time of the velocity vanishes or is very small. But this means indeed that inertial effects are neglected. All higher order cumulants vanish in the limit z--,0. 5.2. Inertial corrections In order to investigate the inertial effects one has to calculate the fourth cumulant. In the case of translational diffusion it was shown that the fourth cumulant leads only to a modification of the diffusion coefficients (compare eq. (4)) but the structure of the diffusion equation remains the same6). On the other hand, if we consider rotational Brownian motion, the fourth cumulant leads to a third order partial differential operator. Translational diffusion is a special case because the cummutation coefficients vanish. The calculation of the fourth cumulant can be made more transparent by the following graphical representation. We have

c =- <1234)c =1 2 3,_4+1 2 3 4 + 1 2 3 4 - 1

1

+

234-1

3 2 4-1

4 2 3

4

4

+

4

I

j3 +

4

(77)

1 3

I N E R T I A L C O R R E C T I O N S IN D I F F U S I O N

17

These diagrams are similar to the Feynman diagrams used in quantum field theories for the calculation of vacuum expectation values25). They follow, basically, from the properties of the creation and annihilation operators a* and a. The lines stand for the matrix Z(ti - t~) which is the analogue of the propagator of a (virtual) particle in quantum field theories. The vertex connected by one line only represents the operators b and b*, respectively (b* for an outgoing line on the right side and b for an outgoing line on the left side. We have for instance 1 2 = b~*Zq(tl - t2)bj). Therefore, the orientation of the graphs is important. The diagrams are written in a way such that the time variables increase by going from the left to the right. The vertices with three external lines correspond to the connection coefficients F'k, e.g. 2

1 ~

4 - b~Zij(tl

-- 12)btk~,kt(t2

-- ta)Zjm(t2

-- t3)['*nmt Z, nr(t3 -- t4)br.

(78)

In this equation the correlation matrix Z ( f i - t3) has been written as the product compare eq. (61). Eq. (78) must be written in this particular form because the matrix Z ( t 2 - / 3 ) may not commute with the operators b '~ Similarly, we have

~ ( / 1 - - /2)~(12 - - /3),

1~4

= b~Zu( fi - tE)b~Z,~(t2- t3)Sjk(t2- ta)bk~,m~(t3 -- ta)b,.

(79)

The first three terms of eq. (77) represent the Gaussian contributions of the fourth moment (1 2 3 4 ) = 1 2 3 4 + 1 2 3 4 + 1 2 3 4. The next three terms follow directly from the definition of the cumulant averages. The remaining terms are the non-Gaussian contributions. They are a consequence of the non-linearity of the equation of motion. They vanish if the group G is commutative such as, for instance, the translations. For translational diffusion we proved that the Liouville operator in the interaction picture fulfills the Gaussian property6). But this property no longer holds if the operators e,. do not commute. The rules outlined above for the calculation of the a v e r a g e s ( / ~ ( l l ) . • • L(ln)) can be viewed as a generalization of theorem 3 in ref. 6. Using this procedure one obtains for the fourth cumulant average

(1 2 3 4)~ = bt,S,~j(A2){b~,Y,jt(A3)Zk.(Aa)bl + b~,S,kI(A3)Zj~(A3)bt • * -- Zjk(A3)bkb~Zt,,(A3) + Fjk,Zkm(A3)Ztr(A3)F,,mr

+ FjtkZkm(A3)Z,t.(A3)F.~r + + rj*k-rk.(A3)-rt~(A3)b~ +

r~t,~k.(A3)St~.(A3)b~.

b~,~,k,(A3)Sjm(A3)r~m

+ b*kZk~A3)S,.(A3)r*~m,}S.,(a4)b. - b*~Z~(a2)Zjk(a3)rk~a,)b,b~Z.(a3)b.,

(8O)

18

U.R. S T E I G E R

where Ai = ti_l - ti. Eq. (80) gives the contribution of the fourth cumulant in the most general case. In order to obtain the inertial corrections one has to calculate the integrals over the variables t~. It is possible to write down a general expression for these integrals (all functions of time of the integrand are exponentials) but it is complicated. There is no simple result in general. Eq. (80) reduces to four terms in the case of translational diffusion with position dependent friction coefficient in one dimension. The result of this calculation is given by eqs. (89). If the matrix ~* does not depend on the position it is possible to choose a coordinate frame such that ~ * becomes diagonal (~ * = ~ .r). If the hydrodynamic interactions are considered too, then there is in general no such global coordinate frame. However, it is possible to choose a local coordinate frame such that ¢ * is diagonal at a certain point q. Such a choice would not help because one needs not only ~* but also its derivative el*. One obtains a simple expression for the fourth cumulant if the friction coefficient is constant. We have in this case ~ * = diag(21 . . . . 2,).

(81)

Note that the inertia matrix M and the friction tensor ~ given by eq. (34) are in general not diagonal in this coordinate frame. Eq. (80) can now be written (1 2 3 4)¢ = b~[b~, b,]bj exp( - 2,A2 - 2~A3- Aid 3 -- ;vA4) + b,[b~bj, b,] exp( - 2,A2 - 2,A3 - 2jA3 - 2,A4) --I- 2b t~b f

* * exp( -- 2i A 2 -- 2kA 3 --21A 3 -- 2jA 4) i~kt)F j(kt)

+ 2bjbkbfi(kj) t t • exp( - ).jA2 - 2kA3 -- 2jA3 -- 2~A4)

+ 2b~bjbkF~k)exp( -- ).,A2 - 2jA3 - 2kA3 -- 2jA,),

(82)

where Fiqk) * = !2(if'Ok • q- Fikj). • We consider now the Laplace transform of the cumulants defined by oo

d(")(s) -= [ d t exp( - st)G(")(t).

(83)

Q/ 0

The diffusion equation (71) becomes (with • denoting the convolution) s P ( s ) = [t~Z)(s) + 1~4)(s)+'" "]* P ( s ) + P ( O ) .

(84)

With eqs. (75) and (83), we find that 1~(2)(s) = - b *

1

s(s + ¢*)

b,

(85)

INERTIAL CORRECTIONS IN DIFFUSION

19

where 1/(s + ~*)=_ (s + ~*)-~, the Laplace transform of the correlation matrix X(t). Using the expressions (73), (82) and the convolution theorem for the Laplace transformation") one obtains

8(~(s) =

cgkb,b~bj + bi*Ujibj * * s(s + ,~3(s + ,~, + ~j)(s + ,~j) .q

cukb * ,*b~bj + cokbibjb * * * k + b*,~,bj s(s + ~ y ( s + ,~, + ~j)

bibj(ciktcj q (b~bkbj + bJb~b,)(c*k + c~j) * * *it "[- CilkCjkl) * * s(s + x,)(s + xj)(s + zk + ~,) s(s + ~j)(s + z,)(s + ~ + ~j)

(86)

For most applications one uses instead of the diffusion equation (70) the time-independent version 0

P ( t ) = [G(2)(oo) + G(4)(oo) + . - -]P(t).

(87)

As it turns out the cumulants G(")(oo), where we have taken the limit t--*oo, can be obtained directly from its Laplace transform (without performing the inverse transformation). We have oo

G(.)(~ ) = [-d G(.)(t ) d t = lim [sG(")(s)]. 3tu

(88)

0

(Note that G<")(0)= 0.) Using this relationship between G~")(s) and G~")(oo) one obtains the inertial corrections G~4)(oo) directly from the expression (86). The limit (88) is trivial. 5.3. Discussion and applications The Smoluchowski type diffusion equation (76) can be used if the higher order corrections are small. This implies that the potential is sufficiently smooth and that the damping by the random forces is relatively strong. For particular cases one can make more precise statements of course. For the damped harmonic oscillator one obtains the condition that E ,~ 1 where E = (mo~/c) ~ with o~ the natural frequency. For a rotating sphere undergoing Brownian motion the inertial corrections are small if ? <~ 1 (compare eqs. (6)-(8)). In order to discuss the inertial effects we consider the following typical and instructive situations. (i) Translational diffusion, one dim., with position dependent friction coefficients. (ii) Translational diffusion, three dim., ~*--const. (iii) Rotational diffusion, U -- 0, ~* -- const.

20

U.R. STEIGER

(iv) Coupled translational and rotational diffusion, two dim., U = 0 , * = const. These four cases contain simplifications which make the calculation easier but they still show important properties. Most results for (i)-(iii) have been published earlier. We use those results to check the equations (80) and (86) once more and to show that the present treatment is very general. The results obtained in case (iv) are new.

Example (i). The fourth cumulant becomes in the limit as t ~ oo (compare remarks following eq. (80))

G~'

~ fmU" = - dq \ c2

me'l-8 1) ~ Lk T -~q + U " j, j =

kT~8 U') c \ ~q + -k-T "

(89)

The prime denotes the partial derivative with respect to q. Eq. (89) is identical with Risken's result (eq. 4.50 in ref. 27) obtained by different methods. This is a third order partial differential operator. The position dependence of the friction coefficient c (c corresponds to ¢ in eq. (33)) leads to an additional correction of the diffusion constant proportional to U'.

Example (ii). Translational diffusion in three dimensions was discussed in detail in ref. 6. Those results are obtained from eq. (86) by setting C*k=O, b~ = ~ 8/dqi and b~ = - x / k T / m 8/8q~- ~ ( k T m ) d U / ~ q i for i = 1, 2 and 3. Example (iii). Rotational diffusion has many applications for the theory of dielectric materials. A great deal o f work has been done in this area. The commutation constants of the group of rotations are c~ = -EUk. We have in the absence of external forces b~ = Ci*jk =

-

-

bi = (kT)'/2I~ 1/2Mj,

rbT'd/2~ - - ~.t~ l i

1-1121-112I+1/2 % f l ~ l ai ~t [3j l ~k

(90) .

(91)

Using these expressions and evaluating (86) in the limit (88) it becomes clear that the fourth cumulant can be written 2) 3) G~) = DI~M,Mm + Dt~,,MjM,,,M,,.

(92)

The coefficients Dt~) and D / ~ can be chosen such that D/~!..,, = n(k) "-'~('l).... ('*)' for z ~ Sk

(93)

Sk is the group of permutations o f k objects. An antisymmetric part would reduce to a lower order partial differential operator due to the commutation relation MIMz - M,,MI = --Et,,~Mk. In ref. 11 the coefficients DIE) and Dt~, were calculated

INERTIAL CORRECTIONS IN DIFFUSION

21

for a special case. It was assumed that the friction tensor and the tensor of inertia commutes which permits diagonalizing both matrices at the same time. Eq. (86) agrees with McConnell's result (ref. 11, eq. 11.5.22-23). Using eq. (86), we may consider a situation where the tensor of inertia I is not diagonal (we still choose a coordinate frame such that ~* is diagonal). More general expressions for Dt~!..tk are now obtained. We have, for instance, Dt~, _

(kT)2 3

r 1-1/2I-I/2f-1/21-1/2 r-'atflz(n)'~ti "iz(I) l f l j "jr(m)

E3 ~ES

2~(2,.+ 2j)

= --1/3 ~ E~a~(,)D~(oEa~(m),

(94)

~tES3

where

D = kT~ -t and E = kT~ -1I~ -1. Example (iv). Eq. (92) describes correction of the rotational diffusion constants for nonspherical molecules. The translational diffusion constants are modified too because in this case translations and rotations are not independent. The following example illustrates this situation. We consider a rod like molecule moving on a two-dimensional surface. We obtain by specializing eqs. (18) to (22) el = cos 0

e2

=sin0 0

- sin 0 dq2, + c o s 0 dq2,

(95)

e3 = ~ - ~ .

e3 corresponds to e 6 in eq. (21). The nonvanishing commutation coefficients are C132= C321 ~ --C231 = --C312 =

1.

(96)

Using the equations (86), (96) and (40) one finds for the fourth cumulant

d2 + ~ cos

~3 + c cos 20 Oq f~q20

02 )

Oqlc3q2 COS20

O] '

(97)

22

U.R. STEIGER

with =

~22--3.,)y~

I

+

3.1"~'2'~'3(~1+ 3.2) + 3.12'32.2 '~

(98)

b - "°~"2"~3 - ~ \3.~ + 3.3 + -3-2-+ 3.3 3.2 - 3.1 c = 2 3.~(3.~ + 3.3)(22 + 3.3)"

The relaxation times 271 are determined by 21 = IT~m, 22 = IT/m and 23 = IR/I. i T are the translational friction coefficients along the /-axis in the body fixed coordinate frame. CR is the rotational friction constant. Since the rotational relaxation process is faster, in general, than the spatial evolution one can often assume a uniform angular distribution. The variables describing the orientation can be eliminated from the diffusion equation. The translational diffusion equation is then -- P =/SAP.

(99)

c~t

/3

is

the

average

translational diffusion constant. We have first approximation (in two dimensions). But eq. (98) shows that there are higher order corrections due to the coupling with the rotational motion. From eq. (97) and (98)

JO= (kT/2)[(i T + ~T)/~T~2T] in

d~kT kT'X m2(kT kT~2(

1

1

)

Consider a slightly deformed sphere. We write 13) i f = 6 r c r / ( R +ER), i T = 6n~/(R - ER), i R = 81rt/R 3 + ~(E) and I = 2mR2/5 + ~(E) assuming a uniform mass distribution. One obtains

D= kT ( 1 4 0 0 2 6rct/---R

)

- - ~ - E 7 + 60(~3) ,

(101)

with the dimensionless constant ~ = IkT/~ ~- mkT/(8Orrbl2R4). 7 never exceeds a few percent, depending on the viscosity o f the fluid or gas. The inertial corrections are larger for small molecules and they are proportional to the square of the temperature T.

6. Conclusion

Many authors have calculated inertial corrections for special situations mentioned in the introduction. A general theory can be viewed as progress because

INERTIAL CORRECTIONS IN DIFFUSION

23

most of those calculations are rather complicated and sometimes not very transparent. We derived a Smoluchowski type diffusion equation for an arbitrary finite dimensional mechanical system under the influence of strong random forces. It has been shown that the nonlinear terms in the equation of motion do not contribute to this type of diffusion equation. It does not contain any inertial corrections (except the position dependence of the diffusion tensor due to the hydrodynamic interactions, which are inertial effects of the fluid). Inertial corrections have been calculated based on the time ordered cumulant expansion. We investigated in particular the coupled translational and rotational diffusion. Applying this theory to rotational Brownian motion one obtains a certain extension of theory of McConnellH). The inertial corrections are given in terms of the commutation coefficients, the potential and its first three derivatives and the friction tensor. They are time-dependent and can be applied to small and large times. Our results agree with results obtained by various other methods such as the Chapman-Enskog procedure or the continued fraction method. This theory however represents an important generalization, and unification.

Acknowledgements The author wishes to thank Ronald F. Fox and Joel E. Keizer for helpful discussions and for reading this manuscript. This work was supported by NSF Grant CHE-8009816.

References 1) 2) 3) 4) 5) 6) 7) 8) 9) 10) 11)

O. Klein, Ark. Mat. Astron. Fys. 16 (1922) 5. H.A. Kramers, Physica 7 (1940) 284. M.V. Smoluchowski, Ann. d. Phys. 48 (1915) 1103. G. Wilemski, J. Stat. Phys. 14 (1976) 153. U.M. Titulaer, Physica 91A (1978) 321. U.R. Steiger and R.F. Fox, J. Math. Phys. 83 (1982) 1678. R.F. Fox. J. Math. Phys. 17 (1976) 1148. P. Debye, Ber. dt. phys. Ges. 15 (1913) 777. G.W. Ford, J.T. Lewis and J. McConnell, Phys. Rev. A 19 (1976) 907. P.S. Hubbard, Phys. Rev. A6 (1972) 2421. J. McConnell, Rotational Brownian Motion and Dielectric Theory (Academic Press, London, 1980). 12) J. McConnell, Physica 111A (1982) 85. 13) J. Happel and H. Brenner, Low Reynolds Number Hydrodynamics (Prentice Hall, Englewood Cliffs, N.J., 1965). 14) U.R. Steiger and R.F. Fox, J. Math. Phys. 2,3 (1982) 296.

24

U.R. STEIGER

15) J.M. Deutsch and I. Oppenheim, J. Chem. Phys. 54 (1971) 3547. B.U. Felderhof, Physica 89A (1977) 373. P. Mazur, Physica ll0A (1982) 128. 16) U.M. Titulaer, Physica 100A (1980) 251. 17) T.J. Murphy and J.L. Aguirre, J. Chem. Phys. 57 (1972) 2098. 18) H. Goldstein, Classical Mechanics (Addison-Wesley, Reading, MA, 1950). 19) V.I. Arnold, Mathematical Methods of Classical Mechanics, Springer Graduate Texts in Math. No. 60 (Springer, New York, NY, t978). 20) R. Abraham and J.E. Marsden, Foundations of Mechanics (Benjamin, Reading, MA, 1978). 21) C.W. Misner, K.S. Thorne and J.A. Wheeler, Gravitation (Freeman, San Francisco, CA, 1973). 22) M.E. Rose, Elementary Theory of Angular Momentum (Wiley, New York, NY, 1957). 23) S. Chandrasekhar, Rev. Mod. Phys. 15 (1943) 1. 24) D. Forster, Hydrodynamic Fluctuations, Broken Symmetry and Correlation Functions (Benjamin, Reading, MA., 1975). 25) S.M. Bilenky, Introduction to Feynman Diagrams (Pergamon Press, Oxford, N.Y., 1974). 26) The equation trace Ck = 0 can be verified directly for the special case given by eq. (37). However, it must hold in general because, otherwise, the Liouville operator (56) would not be skew symmetric neither w.r.t, the normal scalar product nor (because LgMB= 0) w.r.t, the scalar product (51). This would be a contradiction to Liouville's theorem. 27) H. Risken, H.D. Vollmer and M. Mrrsch, Z. Phys. B 40 (1981) 343.