Linear response theory of the viscosity of suspensions of spherical brownian particles

Physica 146A (1987) 417-432 North-Holland, Amsterdam

LINEAR RESPONSE THEORY OF THE VISCOSITY OF SUSPENSIONS OF SPHERICAL BROWNIAN PARTICLES B.U. FELDERHOF Institut fiir Theoretische Physik A, R.W.T.H. Aachen, Templergraben 55, 5100 Aachen, Fed. Rep. Germany

and R.B. JONES University of London, Queen Mary College, Department of Physics, Mile End Road, London E1 4NS, England

Received 15 June 1987

We study the frequency-dependentviscosityof a suspension of spherical particles on the basis of linear response theory applied to the generalized Smoluchowski equation. Hydrodynamic interactions are fully taken into account. We derive expressions for the frequency-dependent viscosity based on a cluster expansion of the linear response.

1. Introduction In this article we set ourselves the aim of deriving general expressions for the viscosity of interacting, spherical Brownian particles based on linear response theory applied to the generalized Smoluchowski equation. The fluid in which the particles are immersed is treated as a continuum background giving rise to hydrodynamic interactions between the particles. The effect of Brownian motion on the viscosity of a suspension of hard spheres was first studied by Batchelor 1). The interplay of Brownian motion and creeping flow was discussed within the framework of the generalized Smoluchowski equation by one of us2). A continuation of this study led to an expression for the stress tensor including the effects of potential and hydrodynamic interactions, as well as Brownian motion3). Here we show how linear response theory may be used to evaluate the average stress tensor of a weakly perturbed suspension. By relating the average stress tensor to the mean velocity gradient we obtain expressions for the frequency-dependent viscosity. 0378-4371/87/$03.50 (~ Elsevier Science Publishers B.V. (North-Holland Physics Publishing Division)

418

B.U. FELDERHOF AND R.B. JONES

The long range of the hydrodynamic interactions implies that the thermodynamic limit must be taken with care. The average flow velocity in a sample of finite size has a pattern which depends on the shape of the sample. Mathematically this shows up in the appearance of integrals over particle coordinates which are not absolutely convergent. However, the viscosity is a local transport coefficient which does have a proper thermodynamic limit. We employ the method of cluster expansion developed in the theory of dielectrics 4) to derive expressions for the viscosity involving absolutely convergent integrals. The absolute convergence of the integrals in the high frequency limit has been established by SchmitzS). In section 2 of this article we recall the generalized Smoluchowski equation. In section 3 we derive an expression for the force density acting on the fluid, taking account of Brownian motion. In section 4 we apply linear response theory to a suspension perturbed by an incident flow field. The discussion is analogous to our earlier study of sedimentation and diffusion6). In section 5 we formulate the cluster expansion and employ this in section 6 to derive an expression for the susceptibility kernel relating average force density and average flow velocity. In section 7 we derive the expressions for the frequencydependent viscosity. The article is concluded with a discussion.

2. Generalized Smoluchowski equation We consider a suspension of N identical spherical particles of radius a immersed in an incompressible fluid described by the linear Navier-Stokes equations for creeping flow. The particles are confined by a wall potential to a volume V. The fluid may be of infinite extent or may be contained in the volume V, in which case hydrodynamic stick or slip boundary conditions are specified at the walls. For example, in the case of shear flow the suspension may be confined by two parallel flat plates moving in opposite directions at which the fluid satisfies stick boundary conditions. The linear Navier-Stokes equations read rtoVZv-Vp =

-F(r),

V- v = O,

(2.1)

where ~/0 is the shear viscosity, v(r) is the flow velocity, p(r) is the pressure, and F(r) is the force density exerted by the particles on the fluid. The solution of (2.1) in the absence of particles is denoted by (v0, P0). Due to conditions at infinity or at the walls this solution may have a prescribed explicit timedependence. We assume that the particles at all times rotate freely and neglect rotational

VISCOSITY OF SUSPENSIONS

419

Brownian motion. Then it suffices to describe the configuration of the entire suspension by the 3N-dimensional vector X = ( R 1 , . . . , RN). The motion of the particles is influenced by the incident flow (v0, P0) and by applied forces (E/(t)}. In the absence of applied forces and interaction forces the particles are carried along by the flow. From the linearity of the flow problem it follows that the translational velocities of the purely convective motion are given by

Uc/= f c,/(r,X).vo(r,t)dr,

j= l . . . . , N ,

(2.2)

where Cj.(r) is the translational convection kernel. The applied forces impose an additional translational motion

UE/ =

~

p.jk(X)"

k

Ek(t),

(2.3)

where ~(X) is the translational mobility matrix. The same matrix gives the motion due to the interaction forces,

UKj = ~k ~/k(X)" Kk "

(2.4)

We assume that the interaction forces, as well as the force exerted by the walls, are derived from a potential according to K = -Vth. In the theory of Brownian motion the statistics of configurations is described by a probability distribution P(X, t) which is assumed to satisfy the generalized Smoluchowski equation. We may derive this from the equation of continuity OP

d---7+ V. (UP) = 0

(2.5)

by supposing that the total velocity U may be decomposed as U ( X , t) = U c + U K + U E + U B ,

(2.6)

where the Brownian velocity U B is related by V B = ~(X)"

KB

(2.7)

to the entropic forces

KBj = - k B T V / In P(X, t),

(2.8)

420

B.U. FELDERHOF AND R.B. JONES

where k B is Boltzmann's constant and T the absolute temperature. Hence one finds the generalized Smoluchowski equation

OP + ~r. [C. Ot

--

voP ] +

~7. [1~- E P ] = @ P

(2.9)

where @P = V . D(X). [VP +/3(Vq~)P],

(2.10)

with the diffusion matrix D(X) = k s r l x ( X )

(2.11)

and with/3 = 1/k s T. In the absence of an incident flow and of applied forces the Smoluchowski equation has the equilibrium solution Peq(X) = e ~*(X)/z,

(2.12)

where Z = [ dX exp(-/3q~) is the configurational integral. d

3. Force density

The hydrodynamic flow problem for freely moving particles in the presence of an incident flow (v0, P0) leads to an induced force density acting on the fluid which is linearly related to the incident flow velocity

Fc(r , X) = - f Z=(r, r', X). Vo(r' ) d r ' ,

(3.1)

where Z= is the convective friction kernelS'7). We have added the subscript oo to emphasize that the response is instantaneous. The kernel has the property that the forces and torques exerted on the fluid vanish. This implies in b r a - k e t notation

= 0,

(pjl2

= O,

(3.2)

where the coordinate representation of the bra-tensors is (~'jl = 0j(r) 1 ,

(Pjl~¢ = - %~Oj(r)(r - Rj)~ ,

(3.3)

VISCOSITY OF SUSPENSIONS

421

with Oj(r) = O(a - [r - Rjl ). The convective friction kernel may be shown to be symmetric,

2.,oo(r, r') = 2~,~.(r',

(3.4)

r),

so that also

2=l,j) = 0 ,

ILIp3 = 0 ,

(3.5)

where the kets are identical to the expressions given in (3.3). In the presence of applied forces a force density on each of the spheres of the form

Fei(r, X) = ~] IIgjk(r, X) " E k

(3.6)

k

is induced with the property that the force J ; e i = ( r j l F e j ) = E j and with vanishing torques ff'Ej = (ajlgEj) = 0. Here ~ k ( r , X) is a transfer kernel related to the convection kernel Cj(r, X) introduced in (2.2). The latter may be decomposed into contributions with the source point r located in each of the spheres, Cjk(r) = Cj(r)Ok ,

~(r)

= ~

Cjk(r ) .

(3.7)

k

The convection kernels and transfer kernels are related by the symmetry 7)

Cj~(r, X) =

W~j(r, X).

(3.8)

The forces (Kk} similarly give rise to the force density Fto(r ,

X) = ~ ~k(r, X). K k

.

(3.9)

k

Finally there is a Brownian contribution from the entropic forces (2.8)

FBj(r, X, t) = - ~

It~jk(r, X)" k B T V k In P(X, t).

(3.10)

k

The average of this quantity over the probability distribution P(X, t) gives the mean contribution to the force density acting on the fluid due to Brownian motion. An integration by parts shows that alternatively the sum of the contributions (3.9) and (3.10) may be obtained as an average of the quantity

422

B.U. FELDERHOF AND R.B. JONES

F;(r, X) = ~ [kBTV k + Kk]" Oky(r, X ) ,

(3.11)

k

where we have used the symmetry (3.8). Clearly the equilibrium average of this force density vanishes. Summing over ] we obtain

F'(r, X) = ~] F~(r, X) = ~] [kBTV] + Ki]. C](r, X ) , y

y

(3.12)

where we have used (3.7).

4. Linear response

In this section we evaluate the mean force density exerted by the spheres on the fluid in the absence of applied forces {Ej} to linear order in the incident flow velocity Vo(r). The mean force density is given by the average of the sum

F(r, X) = Fc(r, X) + F ' ( r , X ) ,

(4.1)

where Fc(r, X) is given by (3.1) and F'(r, X) by (3.12). The system is assumed to be in equilibrium in the infinite past, so that the distribution function P(X, t) is given by the equilibrium distribution Peq(X) for t--+--~. In linear response theory we write the distribution function as e ( x , t) = eeq(X) + ~P(X, t)

(4.2)

and calculate 3P(X, t) to linear order in v 0. To linear order ~P satisfies the linearized Smoluchowski equation O ~3P -- @ 5P = -V" [C" l)0Peq ] • at

(4.3)

The solution with the initial condition 8P = 0 for t = - ~ is given by

~e = - f e~¢'-r)v • [C-

Vo(t')PJ dt'.

(4.4)

It follows from (3.12) that V. [C. Vo(t)Peq ] = flPeq(X) f F ' ( r , X)- Vo(r , t) d r .

(4.5)

VISCOSITY OF SUSPENSIONS

423

Introducing the adjoint operator @*, denoted by ~f, (4.6)

= @* = ( V + / 3 K ) . D . V

and using the property (4.7)

~Peq(X)A(X) = P~q(X)LfA(X) we can cast ~P in the form

~P(X,t)=-t~Peq(X)f dr f dt' eZe(t-r)F'(r,X).eo(r,t').

(4.8)

For the average of an observable A(X) we find to linear order in o 0

A(X)'=(A(X))-~ f /

(A(X)

e:e(t-")F'(r ', X))" Vo(r', t') dr' dt' , (4.9)

where the angle brackets denote an average over the equilibrium distribution Peq(X). In particular this yields for the average force density

F(r,

x ) t = - f (Z.~(r, r', X ) ) . Vo(r', t) dr'

-fl f f (F'(r,X)e~e(t-C)F'(r',X)).Vo(r',t')dr'

dt',

(4.10)

where we have used that the equilibrium average of F ' vanishes. Clearly the first term in (4.10) gives an instantaneous response, whereas the second term describes the retarded response due to the changing distribution function. The expressions take a somewhat simpler form after Fourier transformation in time. The velocity field is transformed as :¢

Vo(r, t) =

f Vo,o(r) --c¢

e -i'°` dto,

Voo,(r) =

if

~

oo

v0(r, t) e i'°t dt.

(4.11)

--~

We denote the corresponding Fourier transform of a non-equilibrium average by

424

B.U. FELDERHOF AND R.B. JONES

A(X) , =

J

A(X) ~, e -i°'' d w ,

A(X) ~

1 2"rr

--~

J

A ( X ) ' e TM dt

(4.12)

--oe

Fourier transformation of the average force density in (4.10) yields F(r)

= -

f

r, r', w ) . vo~o(r' ) d r '

(4.13)

with the kernel :~(r, r', to) = {Z(r, r', w)) ,

(4.14)

where Z(r, r', w) = 2~(r, r') + Z'(r, r', to)

(4.15)

with the frequency-dependent contribution Z'(r, r', co) = - ~ S F ' ( r , X)[ioJ + ~ ] - l F ' ( r ' , X).

(4.16)

Evidently at high frequency only the hydrodynamic interaction term Z~ survives.

5. Cluster expansion The linear response relation (4.13) between force density and incident flow velocity is highly nonlocal due to the long range of the hydrodynamic propagator. In this and the next section we consider instead the relation between force density and average flow velocity. This is described by a transport coefficient which depends only on the local properties of the suspension. We derive the local relation with the aid of a cluster expansion, first developed in the theory of dielectrics4). We may take advantage of the fact that the particles are identical and consider instead of the total force density F ( r ) ~ the contribution F~(r) o, from a selected particle, labeled 1. The average force densities are simply related by F ( r ) ~ = N F I ( r ) ~ . Correspondingly we introduce the operator Zl(X) : 0(1)Z(X) : 0(1)2(1, N ) , where 2¢"is the set of labels 2 . . . . .

(5.1)

N. We decompose this operator with the aid

VISCOSITY OF SUSPENSIONS

425

of a rooted cluster expansion -0(1)Z(1, og) = ~

M(1, ~ ) ,

(5.2)

d4CN

where the sum is over all subsets of 2¢'. The average force density is then given by the expansion N

F ~ = , = ~ (s -1)!

"'"

, s)M(1; 2 . . . . . s)- Vow ,

dRa'''dRsn(l'''"

(5.3) where the kernel M ( 1 ; 2 , . . . , s) depends on frequency and consists of two contributions M ( 1 ; 2 , . . . s) = M=(1;2 . . . . ,s) + M ' ( 1 ; 2 , . . . , s ) ,

(5.4)

the first term corresponding to the rooted cluster expansion of -0(1)2= and the second corresponding to the rooted cluster expansion of -0(1)~'(w). The integrands in (5.3) are weighted with the partial distribution functions n ( 1 , . . . , s) corresponding to the equilibrium distribution Peq(X). Thus we have obtained a cluster expansion of the frequency-dependent response of the force density to the incident flow field. In the next section we shall evaluate also the average flow velocity and express the average force density in terms of it.

6. Susceptibility kernel

In each configuration the flow velocity of the fluid is related to the induced force density by

vow(r,X) = Vo~(r) + f G(r, r'). F~(r', X) dr'

,

(6.1)

where the Green function follows from the linear Navier-Stokes equations (2.1) for the chosen geometry and boundary conditions. For an unbounded fluid the Green function is G(r, r') = T(r- r') with the Oseen tensor 1 1+~ T(r) = - 8¢r~o r

In linear response theory the force density in (6.1) may be expressed as

(6.2)

426

B.U. FELDERHOF AND R.B. JONES £

F,o (r, X) = - J Z(r, r ,' X, o~) -Vo~ L "r'") a~r' ,

(6.3)

where the kernel is given by (4.15). The average flow velocity is given by

v(r) =

H(r, r', X, w)). Voo(r' ) dr',

(6.4)

with the kernel

(6.5)

H(X, w) = I- GZ(X, oo). The kernel may be cluster expanded as

H(N) = ~

(6.6)

L(~),

where 2( is a set of sphere labels and the sum is over all subsets of ?¢', including the empty set. Inserting the cluster expansion (6.6) in (6.4) we find the average flow velocity

v(r)

=

~

".

(6.7)

d R l ' . . d R , n ( 1 , . . . , s ) L ( 1 , . . . , s)-voo~ .

s = 0

Solving for v0~, and substituting in (5.3) we finally find a linear relation between the average force density ff o, and the average flow velocity ~°~ of the form

F(r)

=

$(r, r', o~). v(r')

dr',

(6.8)

where S(w) is a frequency-dependent susceptibility kernel. This may be cast in the form

S(O))

= s=l

(S

- l)!

(6.9)

'

where Ss(w ) involves integrals over s particle positions and is given by S, = ~ S(B),

(6.10)

(B)

where the sum is over all ordered partitions of the labels 1, 2 . . . . . condition that the label 1 be in the first subset and where

s with the

VISCOSITY OF SUSPENSIONS

S(B)

= f "" f d R , " "dR, b(B)C(B).

427 (6.11)

Here b(B) is a block distribution function which is expressed in terms of the partial distribution functions n and C(B) is a chain operator which is expressed in terms of the cluster operators L and M. The detailed expressions are given in ref. 4.

7. Viscosity

At this stage it is possible to take the thermodynamic limit N---~o% volume V--~ ~, with N/V = n constant. In this limit the equilibrium distribution functions are isotropic as well as homogeneous. The susceptibility kernel N(r, r', to) becomes translationally invariant and depends on r - r' only. It is of fairly short range and is determined by the local properties of the suspension. Because of the translational invariance it is convenient to use a Fourier transformation to express the linear response. The Fourier transforms of force density and flow velocity

= e Fq dq

,

o(r) =

. . dq e iq. .Vq

(7.1)

are related by - - t o

Fq

to

=

S( q, tO) "Wqq

(7.2)

Here the transport coefficient s(q, elements of the susceptibility kernel

to)

follows from the plane wave matrix

(qlS(to)lq') = 8~r3s(q, t o ) 8 ( q - q ' ) ,

(7.3)

with the definition

(qlSIq')= f f e-iq'rS(r,r')eiq"r'drdr'.

(7.4)

Using (6.11) and the translational invariance of the distribution functions we may write a cluster expansion directly for the transport coefficient

S(q, to)=s=l ~

( s l l ) '.

~)

"'" dR2""dR'b(B)(qlC(B)lq)'

(7.5)

428

B.U. FELDERHOF AND R.B. JONES

where in the integrand R~ may be taken to be at the origin and the thermodynamic limit is understood. Because of isotropy the transport coefficient must have the form s ( q, w) = se( q, o~)~ + st(q, w)(1 - ~ ) .

(7.6)

On account of the condition of incompressibility q. Vq ~ = 0 the first term has no effect in (7.2) and this equation may be rewritten as rq

= - s t ( q , w)q X (q X Uqq ).

(7.7)

Comparing with (2.1) we see that the frequency-dependent viscosity is given by rt(oJ) = rt0 + art(w) with art(w) = - l i m st(q, w) / q 2 . q--~0

(7.8)

An alternative expression for the viscosity may be obtained by considering a situation where the average flow is linear in the region of interest and of the co form v(r) = g~ "r, where the tensor g,o is symmetric and traceless. The average force density may in general be written as a multipole expansion F(r) o~=

2 (_ 1),_1V,_~ ~F(,)(r )

(7.9)

rt=l

The monopole force density may be put equal to F (1) = F ~ ) + F (~)'

(7.10)

as in (4.1). The convective part F ~.~)vanishes by definition and the remainder is given by r(~)'(r) = Z [k, TVj + Kj]6(r - R j ) . J

(7.11)

For field points r away from the walls this may be written in the form V(~)'(r) = - k B T V n ( r , X) + V. o'~ (r, X ) ,

(7.a2)

where o'+(r, X) is the potential contribution to the stress tensor. The averaged equations (2.1) may be cast in the form ~r-~'° = 0,

~7.~'° = 0,

(7.13)

VISCOSITY OF SUSPENSIONS

429

where the total stress tensor is given by the sum 3) o" = ~hyd + ~ r + °'6 + °rc + o r ,

(7.14)

with the hydrodynamic stress tensor

~hyd,~ = ~70(O~V~ + O~V~) -- p 6 ~ ,

(7.15)

the kinetic stress tensor ~rr = - k B Tnl, and the last two terms given by Orc= ~ ( - 1 ) " - ' v " - Z i F ~ ) ( r ) , (7.16)

.=2

=X n=2

In the region of interest of a large system where the average flow velocity is linear the higher order multipole densities will be uniform on average. Because of the spatial derivatives in (7.16) only the terms with n = 2 will contribute. In order to find the viscosity we must therefore evaluate the average of the potential contribution to the stress tensor o-, and the average of the force dipole density F ~2). In analogy to (4.13) and (4.16) we find for the first contribution tr,(r)

= -

~Z~(r, r', to). Vo~(r' ) d r ' ,

(7.17)

where the kernel is given by

~¢,(r, r', to)= - f l (or0(r, X)[ito + ~ ] - i F ' ( r ' , X ) ) .

(7.18)

The average may be evaluated by cluster expansion and leads to

~-~'~=s=l ~ (sll) ' f""

f dRl"'dR'n(1

,...

, s)M6(1;2,...

,s).

Vo~ ,

(7.19)

as in (5.3). Eliminating the incident flow field we finally obtain

~r+(r) =

f

$~,(r, r', w ) . v ( r ' )

dr',

(7.20)

where the kernel is given by a cluster expansion as in (6.9) and (6.10) with

430

B.U. FELDERHOF AND R.B. JONES

cluster integrals S,~(B)

: f ... f d R l ' . . d R , b(B)C+(B).

(7.21)

In order to find the potential contribution to the viscosity it suffices to consider the tensor ~'4,.,~o.a(B) . . . .

dRa'"dRs

b(B)(IIC,~,~,~,(B)

[xa)

I,,,=o,

(7.22)

where the matrix element is defined as (7.23) and where , )t/x indicates that the symmetric and traceless part of the tensor with indices A,/z is to be taken. In the same manner we must consider the tensor r'-n

f

f

h a/3~,6 ~2' tB'~ ~, ! . . . . . J

,a~

,

J d R 2 " " d R , b(B)(x~,lC.~,(B)lx~)

~3,6

I<=o

(7.24)

The average potential stress tensor is related to the average flow velocity by -

-

IO.~ ¢p.a/3()r

to

,aft

i __oj,T

: s~2= - - ( s1 1)[ ZB ~ o'~t3v~(B) Orv,(r)

6

,

(7.25)

and the average force dipole density is related to the average flow velocity by F(2)" , ~o ~;~ tr~

=

s=l (s

1)T.

'~ (2) tB3 ~e~a~ ,

--'~ O~va(r)

(7.26)

The tensors in (7.22) and (7.24) are both proportional to the fourth rank tensor

A ¢r~ = ½(6~6~a + 6,~6~,) -- ½6~6~ .

(7.27)

The viscosity is given by T/(w) = 7/o + A~/(w) with A n ( w ) ~-- AT~c~(OJ) -k- A T ~ ( Z ) ( w )

with the two contributions found from

(7.28)

VISCOSITY OF SUSPENSIONS

- E An*(t°)A~t3vn = 2 s=2 (s -1)! B 1 ® 1 = 2 (s - 1 )

431

(7.29) (2)

The first contribution vanishes at high frequency as may be seen from (7.18). The second contribution does not vanish at high frequency due to our assumption of nonretarded hydrodynamic interactions. For hard spheres with stick boundary conditions the term with s = 1 gives just the Einstein contribution to the viscosity8). To conclude this section we note that with neglect of hydrodynamic interactions the expressions simplify considerably. In this approximation one puts the hydrodynamic propagator G(r, r') equal to zero. As a consequence the induced flow v(r, X) is equal to the incident flow Vo(r ). Furthermore the convection kernel Cj(r, X) is approximated by 1 6 ( r - Rj), so that the force density F'(r, X) defined in (3.12) is approximated by

F'o(r ,

X) -- ~

[ k B T V k + Kk]t3(r --

Rk)

k

= - k B T V n + V. a%,

(7.30)

where the second equality is valid away from the walls. Substituting this in (7.18) and performing an integration by parts we see that (7.17) becomes

tre~(r) =

H~(r, r', to) : V'v0o,(r' ) dr' ,

(7.31)

where the kernel is given by the stress-stress correlation function

H ~ (r, r', to) = - fl ( try(r, X)[iw + ~P]-lo-6 (r', X ) ) .

(7.32)

In the convective contribution (7.26) only the single particle term s = 1 survives. This is the Einstein contribution to the viscosity. If hydrodynamic interactions are taken into account, then in the expressions (7.29) for the viscosity the hydrodynamic propagator G(r, r') may be replaced by the Oseen tensor (6.2). The reason is that the viscosity is a local transport quantity and in the thermodynamic limit the boundary effects, which in principle may be present in G(r, r'), disappear.

432

B.U. FELDERHOF AND R.B. JONES

8. Discussion

Within the f r a m e w o r k of linear response theory applied to the generalized Smoluchowski equation we have derived expressions for the frequency-dependent viscosity of a suspension of spherical particles. Full account has been taken of hydrodynamic interactions. The explicit evaluation of the viscosity remains a formidable task. As a special application we intend to perform the calculation for a suspension of spherical particles to second order in the volume fraction. This will allow us to c o m p a r e with Batchelor's work for hard spheres 1). In this article we have neglected rotational Brownian motion. For suspensions of anisotropic particles it will be essential to take this into account. The present theory may be extended to include the rotational degrees of freedom. Expressions for the viscosity analogous to (7.29) may be derived.

References

1) 2) 3) 4) 5) 6) 7) 8)

G.K. Batchelor, J. Fluid Mech 83 (1977) 97. B.U. Felderhof, Physica l18h (1983) 69. B.U. Felderhof, Physica 147A (1987) 203. B.U. Felderhof, G.W. Ford and E.G.D. Cohen, J. Stat. Phys. 28 (1982) 135. R. Schmitz, Dissertation R.W.T.H. Aachen 1981. B.U. Felderhof and R.B. Jones, Physica II9A (1983) 591. B.U. Felderhof and R.B. Jones, Physica 146A (1987) 404, this volume. A. Einstein, Ann. der Physik 19 (1906) 289, 34 (1911) 591.