- Email: [email protected]
Statistical mechanics of irreversible processes
Brout, R.
Physica XXII 263-272
Prigogine, I. 1956
STATISTICAL MECHANICS OF IRREVERSIBLE PROCESSES PART VI : THERMAL CONDUCTIVITY OF CRYSTALS b y R. B R O U T a n d I. P R I G O G I N E Facult~ des Sciences de l'Universit~ Librc de Bruxelles, Belgique
R~sum~
Nous avons montr~ ant~rieurement 1)3) que les fonctions de distribution locales tendent, par suite des interactions harmoniques, vers une forme gaussienne et cela moyenant des hypotheses assez g4n6rales sur la distribution initiale. L'~volution vers la distribution d'6quilibre d'un r4seau 16g~rement anharmonique peut d~s lors ~tre de composs6e en deux 4tapes proc6dant suivant des 6chelles de temps diff4rentes : une *chelle de temps courte li6e aux fr6quences harmoniques et une 4chelle de temps longue liSe ~t l'anharmonicit~ au cours de laquelle s'6galent les Snergies des modes normaux. Cette ,,preparation" initiale du syst~me par les seules forces harmoniques est utilis6e pour obtenir une th~orie tr~s 614mentaire de la conductibilit6 thermique dans les r4seaux. En particulier nous obtenons une expression explicite du coefficient de conductibilit6 thermique en termes de grandeurs mol4culaires. I. Introduction. In parts I I I and IV of this series it has been shown t h a t u n d e r r a t h e r general conditions (summarized in § 2 of this paper) local distribution functions in an infinite h a r m o n i c crystal t e n d a s s y m t o t i c a l l y to m u l t i v a r i a t e gaussian distributions. This result has interesting consequences which are discussed in the present paper. T h e evolution t o w a r d s equilibrium of a slightly a n h a r m o n i c crystal can be subdivided into two steps which proceed according to quite different time scales (cf. § 2). T h e first step corresponds to the a p p e a r a n c e of the m u l t i v a r i a t e gaussian distribution and is r e a c h e d in times of order 1/co (co = characteristic f r e q u e n c y of the h a r m o n i c crystal). T h e second step is the disappearance of all second order correlations (except of course autocorrelations) t h r o u g h the effect of a n h a r m o n i c forces. This p a r t corresponds to the evolution t o w a r d s e q u i p a r t i t i o n of e n e r g y b e t w e e n n o r m a l modes a n d corresponds to a " l o n g " time scale related to t h e a n h a r m o n i c i t y of t h e crystal. W e shall show t h a t if we t a k e a c c o u n t of this initial " p r e p a r a t i o n " of the s y s t e m b y h a r m o n i c forces, the molecular t h e o r y of t h e r m a l c o n d u c t i v i t y t a k e s a r e m a r k a b l y simple form. I n s t e a d of using t h e full B o l t z m a n n equation for this p r o b l e m it becomes sufficient to use the e q u a t i o n for the second m o m e n t s alone. W e easily obtain, in this way, an explicit form for t h e t h e r m a l c o n d u c t i v i t y of the crystal in t e r m s of molecu.!ar quantifie~ --
2 6 3
--
264
R. B R O U T A N D I. P R I G O G I N E
(cf. §§ 4-5). The phenomenological equation for the change of temperature (Fourier's equation) m a y be justified in the same w a y (cf. § 5). We restrict our study to one-dimensional chains with nearest neighbour interactions in spite of the fact that such systems have a finite thermal conductivity only if four phonon processes are retained 3). The formalism is however so much easier for these systems that we prefer to consider them in detail in the text and to consider the three dimensional extension which presents no special problems in the Appendix. 2. Homogeneous systems. Let us consider the "local" distribution function for n variates Yl . . . . Yn (we use the notation of ref. 1)) /(Yl, Yz . . . . Yn, t). The central result of our previous paper 1) ~.) was to show that under rather general conditions this distribution function tends assymptotically to the multivariate gaussian distribution
/(yl . . . . .
y . , t) -+ (IAII/2/~ "/2) exp (-- Z a,jy, y~)
(2.1)
iA,~I ---- Y~Yj = cofactor of ai~ The time necessary to reach the assypmtotic distribution is of order (1/o) where co is the characteristic crystal frequency. Further this local system will become spatially homogeneous in these times. We have seen that these statements are true if the initial range of correlation of variates is finite, i.e. there exists a number K such that Y~Yn+,,, = 0,
all n,
all m >~ K
(2.2)
. If initially there were no systematic inhomogeneities extending over a macroscopic range, then 1/co is also the order of magnitude of the time necessary to reach homogeneity for the whole system. This is the situation which we shall investigate in this paragraph. The asymptotic distribution (2. l) differs from the canonical equilibrium distribution only b y the existence of the second moments y,y~. In a homogeneous system the moments are independent of i and we m a y use the notation i k = YnYn+k = limN-~oo 1/N E,, YnYn+~
(2.3)
It was shown that the _Ts are invariants for the harmonic systems 1) ,) and that they are Fourier coefficients of the energy of normal modes *) i s = limN_~oo E l E, cos lk
(2.4)
The persistence of the second moments in (2.1) shows therefore that equipartition of energy has not been reached. We m a y note the simple physical interpretation of the first I s (I 0 mean energ3~, -Yl mean heat (or energy) flow . . . . ). *) All f o r m u l a e in the t e x t are giveix ill t e r m s of non d i m e n s i o t t a l v a r i a b l e s . The t i m e used is the ~lormalized t i m e (2k/m)¿/~ t, and the energies of the n o r m a l m o d e s are d i v i d e d b y Eo.
I R R E V E R S I B L E PROCESSES PART VI : T H E R M A L C O N D U C T I V I T Y OF CRYSTALS 2 6 5
The role of anharmonic forces is therefore reduced to the decay of the correlations (2.3) or to the evolution of the energies of normal modes towards equipartition. For small anharmonic forces this is a far slower process. The evolution towards statistical equilibrium m a y therefore be split into two steps which proceed according to quite different time scales. As the anharmonic evolution of the system is reduced to the s t u d y of second moments only, the examination of the full distribution with a Boltzmann equation becomes superfluous. The decay of the second moments can be described b y the Peierls a) ,) equations dEk/dt = o~ Zk,~,, Ckk'k'" Look E~, Ek,, -- ~ok, Ek Ek,, -- cok,, Ek Ek,7
(2.5)
where Ckk,k," contains the anharmonic constant squared and is symmetric in its indices. The sum in (2.5) is over all processes such that k-- (k'+k") -
:K
+
= 0
(2.6)
K is the vector in the reciprocal lattice necessary to bring the third phonon (say'k") into the basic unit cell of the reciprocal lattice when the two others are given. The Peierls equations (2.5) are valid only in a time smoothed sense for times of 0 (1/(anharmonic constant) ~). Near equilibrium we m a y linearize equations (2.5). We set Ek = E0 + ~ and retain only first order terms in ~k. Using (2.6) we get d~k/dt : X k, Bkk,-e k,
(2.7)
with B ~ , - - B~,~ = 2 [X~,, C~,~,, c o ~ , -- 5~,/2 o~ Z~,,~,,, C~,,~,,,]
(2.8)
The equations of decay of the [ Nare obtained from (2.4) with the result d i n / d t = Zm~>l 2 .... i,~
n > 1
(2.9)
The relaxation constants ~tn,,, are defined '~n, :
~,,, :
Zkk" COS kn Bkk, cos k ' m
(2.10)
Again (2.9) is valid only near equilibrium and for homogeneous systems. The decay constants have again a simple physical meaning. For example 1/211 is the decay time of the heat flow for a system in which all the other second moments [,~ (m > 1) are zero. We shall see that the thermal conductivity can precisely be expressed in terms of the 2,,~. We m a y stress that if initially the range of correlation is infinite (and (2.2) is not satisfied) we cannot subdivide the evolution of the system into the two steps we just discussed. This is the situation (cf. 2.4) if initially only a single normal mode is excited. The subsequent evolution depends therefore to a certain extent on the initial preparation of the system.
266
R.' BROUT AND I. PRIGOGINE
3. Entropy. In this section we study the thermodynamic consequences of the mechanical evolution discussed in the previous section. We define the entropy of a local system of n variables as S I"l = -- k f/c"l log/c-I dy: . . . . . .
dy n
(3.1)
where l ~"1 is the local distribution function. This value of s Inl takes on a well defined unambiguous value only for, n>~range of correlation, in which case the "surface" terms due to correlations with the system exterior to n become negligible in comparison to the number of interior terms. We shall always assume a finite range of correlation and define S I'1 accordingly. It is first observed that in a homogeneous system where only harmonic forces are present that the maximum of (3.1) subject to the dynamical constraints Y,Y,,+m = i m (independent of t) leads to the multivariate Gaussian distribution (2.1) towards which the system asymptotically tends. This is the same distribution which G r a d 6) has found in a generalized microcanonical ensemble where there are a number of dynamical constraints in addition to the total energy. Thus the thermodynamic interpretation of the harmonic evolution is that the entropy increases to its highest possible value within the framework of the mechanical constraints imposed b y the harmonic invariants of the system. When one removes the harmonic invariants b y allowing anharmonic forces, the system may further evolve, increasing its entropy to the final equilibrium value. It is to this evolution that we now turn. It is necessary to consider only small deviations of second moments from their equilibrium value. If deviations from equilibrium are retained to the second order we find from (3.1) (S (") - - S(°))/k = { Z,L1 (y~ - - y~O) _ ¼ Z , L I ( ~ - - y~O) _ ½ E~., (y,y,)2
(3.2)
where the superscript 0 implies the average taken with the equilibrium distribution function. The first two terms on the r.h.s, of (3.2) are unsigned and can be shown to give rise to an entropy flow related to the heat flow in the system. Indeed comparison of (3.2) with the macroscopical expression S - - S O = c~f~, ° d T / T = c~(T - - T o ) I T ) - - ½ c~(T - - To)~/T~
(3.3)
identifies the first two terms in (3.2) as the change in entropy due to a reversible change in temperature. We m a y also study the significance of the different terms in (3.2) b y the following method. We m a y split the change of entropy into two contributions: the entropy flow d , S / d t and the entropy production d i S / d t 6) 7) d S / d t =: d, S / d t + d , S / d t (3.4) The second principle requires that d,S/dt > 0
(3.5)
IRREVERSIBLE
PROCESSES PART VI : THERMAL CONDUCTIVITY OF CRYSTALS
267
Let us write the equations of motion including, in a formal way, the an~ harmonicity (cf. t) dy,/dt = ½(Y~+I -- Y,-1) + (dY,,/dt),,h
(3,6)
The evolution of the second moments in time is then d y . y . + m / d t -- ½ [Y.Yn+,~+t -- Y . - t Y . + ~ + { lYe+mY.+1 -- Y-Y~+,~-I] + + (dy,,y,,/dt)a,, h
(3.7)
If (3.7) is substituted in (3.2) we get d , S " / d t = E,.~ y, yj (dyj y,/dt)..~
(3.8)
where we have included all boundary terms in the entropy flow and isolated the entropy production. For a homogeneous system (with n >~ range o f corralation K) using (2.9) we obtain diS'*/dt ~- -- n Ei. ~ ~i~ I i l s
(3.9)
From (2.8) it can be deduced that the eigenvalues of Bk,, are negative, which is physically obvious. Consequently the eigenvalues of the A~jdefined by (2.10) are also negative and (3.9) is a positive definite form as it must be. 4. Stationary state. In discussing the possibilities of a stationary state one fact becomes immediately apparent. The stationary state cannot exist indefinitely. For eventually the thermal gradient will be smoothed out by harmonic forces. If we confine our attention to a subsystem of length N , over which there exists a temperature gradient but no other inhomogeneities this smoothing takes place for t ,-~ N in the middle of the system.Thus for t ~ N the center of the system maintains a constant gradient. At the same time, the time must be sufficiently long so that equs (2.5) are valid i.e. t >~ 1/co. In this intermediate time scale the possibility of a stationary state presents itself. We further suppose that aside from the presence of the systematic temperature gradient all other moments are initially but irregularly inhomogeneous and are rapidly smoothed out in t = 0 (1/o~). In this case we can write down eqs (3.7) in a time smoothed sense using (2.9)
dil/dt=
d i k / d t = Y'~ ~kz iz
(4.1)
½ (Y,+I * - -- Y~) + Z,~>l ~t, 1,
(4.2)
dye/dr = 0
(4.3)
The average value of the energy y~ is not changed by the decay equations (2.9.). In the stationary state eqs (4.1) and (4.2) have their left hand sides equal to zero and it is possible to solve for the I's in terms of grad -?0 = ½(Y~+I--Y~).~
268
R. B R O U T A N D I. P R I G O G I N E
In particular, we have for the heat flow i 1 = ½(A-1)n grad h
(4.4)
which is Fourier's law. (A-1)n is the 11 component of the inverse of the matrix [lyt~][. B y these considerations we have therefore identified the thermal conductivity coefficient n as = ½(A-1)u
(4.5)
As ;tij is proportional to I o, we have n ,-~ 1/T. As a first approximation we m a y expect n = ½ 1/;tn
(4.6)
where 1/,ln is the relaxation time of the heat flow. w e see that the thermal conductivity is directly expressed in terms of the coefficients a~j related to the second moments i~ and not in terms of the relaxation coefficients Ckk,k,, appearing directly in equations (2.5). As can be expected in the high temperature region (classical mechanics) the inverting processes do not play any special role in a. Their existence must however be assumed to insure the ergodicity of the system 8) 4). If the entropy production (3.9) is minimized subject to constant heat flow eqs (4.1) become zero b y minimization and eq. (4.2) is zero b y constraint. We have thus the well known theorem 6) that the entropy production in the presence of a gradient is minimum in the stationary state.
5. Fourier's equation. In order to derive the phenomenological equations of thermal conductivity out of the stationary state, let us reconsider (3.7). We divide our system of length N into a large number of domains each of which is large with respect to microscopic dimensions and small with respect to macroscopic variations. Since the domain contains a large number of atoms we may use our time smoothed and space smoothed equations of anharmonic decay assuming that all moments b u t the energy are homogeneous within a domain. In this way (3.7) m a y be written in a domain D with boundaries n and n + D d/dt (Y~,,DY,Y,+k)= (d/dt (~,d) Y,Y,+k))anh
d/dt (XieD Y,Y,+x)
= ½ (yn+ -3
-
-
-
+
k > 1
(d/dr
(5.1) (5.2)
feD
d/dt (Y',,D
= ½
Y.Y -I)
(5.3)
Now the structure of equations (5.1) and (5.2) differ in character from (5.3) in that (5.3) only involves surface terms of the 'domain, representing the heat flux in and out of the domain. If one adds up (5.3) over all domains one is finally left with a big difference of two surface terms. Thus small inhomogeneities in the heat flow Y~Yn+I become multiplied into macroscopic pro-
I R R E V E R S I B L E PROCESSES PART VI : THERMAL C O N D U C T I V I T Y OF CRYSTALS 2 6 9
portions by the totality of equations (5.3). In (5.1) and (5.2) this successive addition of terms is each time compensated by a local interior term representing anharmonic decay. These local terms then present the possibility of local equilibrium. In (5.2) we see this brought about by the compensation of the harmonic propagation of the gradient against the local decay of heat flow. Note that in the absence of anharmonic terms such a compensation is impossible; the heat flow m a y continue unimpeded. In (5. l) and (5.2) we substitute the equations of anharmonic decay which are valid in the smoothed region D to give (we call T~ the value of I k in D) diD~ d t = X~=z ~k~ _]D k > 2
(5.3)
diD/dt = ½ grad i D + XT=1 21zi~)
(5.4)
diDo/dt = (l D+z -- l D-z)
(5.5)
where we have assumed all moments approximately homogeneous in D except _70. (5.3) and (5.4) contain only terms in ]D and m a y be solved for ?D in terms of (grad i~) D. For asymptotic terms in the intermediate time scale, the solution tends to the state of local equilibrium and we have I-D , -+½ (A-a)1, (grad i0) D and in particular i D --> ½(A-1)zz (grad i0) D which establishes Fourier's heat law in the domain D again yielding ½(A-1)lz for the thermal conductivity •. Combining with (5.5) we have diD/dt = n [(grad i0) z)+z -- (grad
-T0)D-lj
(5.6)
which is the discrete analogue of the thermal conductivity equation ST/at = -, a~T/a~ 2
(5.7)
The essential argument of the establishment of (5.6) was the concept of local equilibrium in microscopic domains sufficiently large for the smoothed equations of decay to be valid. These smoothed equations of decay counteract the harmonic propagation of the gradient and are responsible for the achievement of local equilibrium. It must be remembered that these arguments are valid only within the intermediate time scale previously discussed. For sufficiently large times t ~-~N all gradients break down under the smoothing out action of harmonic forces. This rather rapid breakdown is due to the high efficiency of conductivity of the thermal bath in which the system is immersed. A more anharmonic bath would give rise to a larger intermediate time scale. Acknowledgments. The research reported in this paper has been made possible through the support and sponsorship of the Air Research and Development Command, United States Air Force, through its European Office.
270
R. B R O U T A N D I. P R I G O G I N E
APPENDIX
In this appendix we wish to point out the necessary generalizations that make the arguments in the text applicable to three dimensions. Many of the considerations have immediate generalizations to the general three dimensional case. The considerations leading to equations of the type (4.6), (4.1), (4.2) are necessarily complicated and we treat the simple case of an cubic crystal with nearest neighbour interactions. This latter supports Umklapp processes and so is sufficiently illustrative for our purposes. We first remark that the three dimensional harmonic problem has been • treated in ref. 2) with the same asymptotic results as in one dimension save for the introduction of quasi-invariants I~. These quasi-invariants rapidly become constant. Local distributions tend to multivariate gaussians and the system becomes homogeneous. Unlike the one-dimensional system homogenity does not imply random phases for there m a y still exist correlations between normal modes with the same value of k b u t of different polarizations. Thus the random phase hypothesis goes further than mere homogenization of the crystal. In the case of random phases the energies of the normal modes are related to the invariants I in the homogeneous system b y 9.) s 1 e~(k) e~(k) I~{ = u .a ua.+.. = Z, Zr=
I~,~1" cos km
which is quite similar to the one dimensional formula and allows one to find the decay equation of I~{ which then take on the same form as in one dimension where now the t . . . also contain the inverse vectors of polarization
4k). The most complicated generalization to three dimensions is the establishment of the analogues of (4.1) and (4.2) in the text. This we proceed to do for the cubic crystal with nearest neighbour interactions and two characteristic interaction constants. The important point is to calculate the time derivatives of energy and heat flow of a region of the crystal due to the harmonic equations of motion. From the outset we will consider a crystal homogeneous in all moments save the energy for such times that the quasi-invariance 3) of the second moments, whose time derivatives do not involve the temperature gradients, has already been established. The Hamiltonian of the system we consider is H
3
,, . + l a - -
.j + (U~+lo---U~)2)J
(A. 1)
to give rise to the equation of motion
dt = (2u~--U~+l~--U~-l~)+e [(2Un--Un+a--Un_l~)+ (2U~--U~_IV-- n--l)] (A.2) *) Here again we use non-dimensional variables, with the variables u, u of the same dimensions as y in the text.
I R R E V E R S I B L E PROCESSES PART VI : T H E R M A L C O N D U C T I V I T Y OF CRYSTALS 271
In these equations u,, is the displacement of the nm atom, with components u~, ua,, u~; 1a is the unit vector in the ~th direction. There are three kinds of second moments that must be considered. Of these the type velocity × velocity and position × position do not involve the temperature gradient and hence are quasi-invariants. Of the moments v e l o c i t y × position, it is seen from (A.2) that only moments of the type "a a u , (u,+l --u~,) involve the temperature gradient and hence all others are quasi-invariant. We shall presently establish this moment as the heat flow. In this way the analogues of (2.2) are established for di~/dt k > 2. We now proceed to find the appropriate relations for energy and heat flow. In order to define the energy content of a region containing (2N) z atoms (we take a cube of length 2N) it is necessary to establish a convention of how to split the potential energy of surface atoms. We shall take the natural definition of } the surface potential energy as belonging to a given volume. With this definition
E(2N)S=
½
Z..,,.,,. N lit2, + Y',,a Aoa (u]
-
-
u~,_la) 2]
i~v~'a7,6= " 1/.~z+N n o n l l = - N A., [(u~o,.a,~V+l~--u~.,.a,x)9--(u~.,.a,_.v--UO.o,.a,_~._l) 2] (A.3) a,B#?
where the matrix A oa is given b y A o a = 1 ~ = / 5 ; A ~ a - - e e : # / 5 . The vectors n in the surface terms have had their components expressed in terms of no, n o, %. In order to find the heat flux through the region it is necessary to differentiate (A.3) with respect to time using (A.2).The resulting expression cancels all interior terms and one is left with an analogue of a surface integral. Thus dEc2N)a/dt = Z.=~ a (q~N -- qr--N-1) (A.4) ,,o,,,a=-N
~ ,,°.aN+l~,~
,,o,,aNJ
(U.a.a,V+l~--U.o.aN)]
(A.5)
(A.5) is the analogue of y , (Y,+I + Y,-1) which was identified as a heat flow in one dimension. In (A.5) this becomes a surface integral In order to establish the analogue of (2.2) for 11, we take the time derivative of (A.5). A tedious calculation in which one uses the property of homogeneity of all moments and neglects terms of 0(1) and 0(N) in comparison with 0(N2), yields
dq~/dt=½ Y,.a . {eZ~ = ~AoaE(u~o+L.aN+I--U~o..a,N+I)~--(U~+~,.a.v--U~o,.a.X)~] + 3 + s Z~= 1 A~
~
~
(Unan#+ l , N + l - - Una, n~,N + l)
+ Z~= a 1 A~ r/t ~ ~ q ~ ~.~,.a,N+I -- /~..,-~,N~
2
--
6
(~.o,.~+1,~ ~.o,.~)~] + (A.6)
Each of the three terms in (A.6) expresses a energy difference. As we have assumed our crystal to be homogeneous in every sense but the presence of a systematic gradient, it is appropriate to assume equipartition between the
272
IRREVERSIBLE
PROCESSES PART VI" THERMAL CONDUCTIVITY
OF CRYSTALS
kinetic and potential energies of a single particle in which case (A.6) becomes
dq~,/dt = const X,,,,,a (Ena, no,N+lr -- En,,na,a, )
(A. I)
which establishes the analogue of (2.2) for (dla/dt)ha~o,~i e. It should be pointed out t h a t the model here t r e a t e d supports inverting processes. For given k along a crystal axis say z the three frequencies of polarization are e sin 2 k, e sin 2 k, sin ~ k which is consistent with the selection rules of inverting processes (see P e i e r 1 s 3)). Received 20-12-55.
1,~EFERENCES 1) K l e i n , G. and P r i g o g i n e , I., P h y s i c a 1~ {1953) 1053. Prigogine, I. a n d B i n g e n , R., P h y s i e a :~! (1955) 299. CI. P e i e r 1 s, R., The QualltUm T h e o r y of Solids, Oxford C l a r e n d o n Press 1955. B r o u t, R. a n d P r i g o g i n e , I., P h y s i c a 2 1 . G x a d, H., Comm. pure and app. Sci..5 (1952) 455. 6) P r i g o g i n e, I., E t u d e t h e r m o d y n a m i q u e des Phfinom~nes irr~versibles, Liege, Desoer (1947). 7) D e G r o o t, S. R., T h e r m o d y n a m i c s of irreversible Processes, N o r t h H o l l a n d P u b l i s h i n g Co. (19521.
2) 31 4) 5)
Physica XXII 263-272
Prigogine, I. 1956
STATISTICAL MECHANICS OF IRREVERSIBLE PROCESSES PART VI : THERMAL CONDUCTIVITY OF CRYSTALS b y R. B R O U T a n d I. P R I G O G I N E Facult~ des Sciences de l'Universit~ Librc de Bruxelles, Belgique
R~sum~
Nous avons montr~ ant~rieurement 1)3) que les fonctions de distribution locales tendent, par suite des interactions harmoniques, vers une forme gaussienne et cela moyenant des hypotheses assez g4n6rales sur la distribution initiale. L'~volution vers la distribution d'6quilibre d'un r4seau 16g~rement anharmonique peut d~s lors ~tre de composs6e en deux 4tapes proc6dant suivant des 6chelles de temps diff4rentes : une *chelle de temps courte li6e aux fr6quences harmoniques et une 4chelle de temps longue liSe ~t l'anharmonicit~ au cours de laquelle s'6galent les Snergies des modes normaux. Cette ,,preparation" initiale du syst~me par les seules forces harmoniques est utilis6e pour obtenir une th~orie tr~s 614mentaire de la conductibilit6 thermique dans les r4seaux. En particulier nous obtenons une expression explicite du coefficient de conductibilit6 thermique en termes de grandeurs mol4culaires. I. Introduction. In parts I I I and IV of this series it has been shown t h a t u n d e r r a t h e r general conditions (summarized in § 2 of this paper) local distribution functions in an infinite h a r m o n i c crystal t e n d a s s y m t o t i c a l l y to m u l t i v a r i a t e gaussian distributions. This result has interesting consequences which are discussed in the present paper. T h e evolution t o w a r d s equilibrium of a slightly a n h a r m o n i c crystal can be subdivided into two steps which proceed according to quite different time scales (cf. § 2). T h e first step corresponds to the a p p e a r a n c e of the m u l t i v a r i a t e gaussian distribution and is r e a c h e d in times of order 1/co (co = characteristic f r e q u e n c y of the h a r m o n i c crystal). T h e second step is the disappearance of all second order correlations (except of course autocorrelations) t h r o u g h the effect of a n h a r m o n i c forces. This p a r t corresponds to the evolution t o w a r d s e q u i p a r t i t i o n of e n e r g y b e t w e e n n o r m a l modes a n d corresponds to a " l o n g " time scale related to t h e a n h a r m o n i c i t y of t h e crystal. W e shall show t h a t if we t a k e a c c o u n t of this initial " p r e p a r a t i o n " of the s y s t e m b y h a r m o n i c forces, the molecular t h e o r y of t h e r m a l c o n d u c t i v i t y t a k e s a r e m a r k a b l y simple form. I n s t e a d of using t h e full B o l t z m a n n equation for this p r o b l e m it becomes sufficient to use the e q u a t i o n for the second m o m e n t s alone. W e easily obtain, in this way, an explicit form for t h e t h e r m a l c o n d u c t i v i t y of the crystal in t e r m s of molecu.!ar quantifie~ --
2 6 3
--
264
R. B R O U T A N D I. P R I G O G I N E
(cf. §§ 4-5). The phenomenological equation for the change of temperature (Fourier's equation) m a y be justified in the same w a y (cf. § 5). We restrict our study to one-dimensional chains with nearest neighbour interactions in spite of the fact that such systems have a finite thermal conductivity only if four phonon processes are retained 3). The formalism is however so much easier for these systems that we prefer to consider them in detail in the text and to consider the three dimensional extension which presents no special problems in the Appendix. 2. Homogeneous systems. Let us consider the "local" distribution function for n variates Yl . . . . Yn (we use the notation of ref. 1)) /(Yl, Yz . . . . Yn, t). The central result of our previous paper 1) ~.) was to show that under rather general conditions this distribution function tends assymptotically to the multivariate gaussian distribution
/(yl . . . . .
y . , t) -+ (IAII/2/~ "/2) exp (-- Z a,jy, y~)
(2.1)
iA,~I ---- Y~Yj = cofactor of ai~ The time necessary to reach the assypmtotic distribution is of order (1/o) where co is the characteristic crystal frequency. Further this local system will become spatially homogeneous in these times. We have seen that these statements are true if the initial range of correlation of variates is finite, i.e. there exists a number K such that Y~Yn+,,, = 0,
all n,
all m >~ K
(2.2)
. If initially there were no systematic inhomogeneities extending over a macroscopic range, then 1/co is also the order of magnitude of the time necessary to reach homogeneity for the whole system. This is the situation which we shall investigate in this paragraph. The asymptotic distribution (2. l) differs from the canonical equilibrium distribution only b y the existence of the second moments y,y~. In a homogeneous system the moments are independent of i and we m a y use the notation i k = YnYn+k = limN-~oo 1/N E,, YnYn+~
(2.3)
It was shown that the _Ts are invariants for the harmonic systems 1) ,) and that they are Fourier coefficients of the energy of normal modes *) i s = limN_~oo E l E, cos lk
(2.4)
The persistence of the second moments in (2.1) shows therefore that equipartition of energy has not been reached. We m a y note the simple physical interpretation of the first I s (I 0 mean energ3~, -Yl mean heat (or energy) flow . . . . ). *) All f o r m u l a e in the t e x t are giveix ill t e r m s of non d i m e n s i o t t a l v a r i a b l e s . The t i m e used is the ~lormalized t i m e (2k/m)¿/~ t, and the energies of the n o r m a l m o d e s are d i v i d e d b y Eo.
I R R E V E R S I B L E PROCESSES PART VI : T H E R M A L C O N D U C T I V I T Y OF CRYSTALS 2 6 5
The role of anharmonic forces is therefore reduced to the decay of the correlations (2.3) or to the evolution of the energies of normal modes towards equipartition. For small anharmonic forces this is a far slower process. The evolution towards statistical equilibrium m a y therefore be split into two steps which proceed according to quite different time scales. As the anharmonic evolution of the system is reduced to the s t u d y of second moments only, the examination of the full distribution with a Boltzmann equation becomes superfluous. The decay of the second moments can be described b y the Peierls a) ,) equations dEk/dt = o~ Zk,~,, Ckk'k'" Look E~, Ek,, -- ~ok, Ek Ek,, -- cok,, Ek Ek,7
(2.5)
where Ckk,k," contains the anharmonic constant squared and is symmetric in its indices. The sum in (2.5) is over all processes such that k-- (k'+k") -
:K
+
= 0
(2.6)
K is the vector in the reciprocal lattice necessary to bring the third phonon (say'k") into the basic unit cell of the reciprocal lattice when the two others are given. The Peierls equations (2.5) are valid only in a time smoothed sense for times of 0 (1/(anharmonic constant) ~). Near equilibrium we m a y linearize equations (2.5). We set Ek = E0 + ~ and retain only first order terms in ~k. Using (2.6) we get d~k/dt : X k, Bkk,-e k,
(2.7)
with B ~ , - - B~,~ = 2 [X~,, C~,~,, c o ~ , -- 5~,/2 o~ Z~,,~,,, C~,,~,,,]
(2.8)
The equations of decay of the [ Nare obtained from (2.4) with the result d i n / d t = Zm~>l 2 .... i,~
n > 1
(2.9)
The relaxation constants ~tn,,, are defined '~n, :
~,,, :
Zkk" COS kn Bkk, cos k ' m
(2.10)
Again (2.9) is valid only near equilibrium and for homogeneous systems. The decay constants have again a simple physical meaning. For example 1/211 is the decay time of the heat flow for a system in which all the other second moments [,~ (m > 1) are zero. We shall see that the thermal conductivity can precisely be expressed in terms of the 2,,~. We m a y stress that if initially the range of correlation is infinite (and (2.2) is not satisfied) we cannot subdivide the evolution of the system into the two steps we just discussed. This is the situation (cf. 2.4) if initially only a single normal mode is excited. The subsequent evolution depends therefore to a certain extent on the initial preparation of the system.
266
R.' BROUT AND I. PRIGOGINE
3. Entropy. In this section we study the thermodynamic consequences of the mechanical evolution discussed in the previous section. We define the entropy of a local system of n variables as S I"l = -- k f/c"l log/c-I dy: . . . . . .
dy n
(3.1)
where l ~"1 is the local distribution function. This value of s Inl takes on a well defined unambiguous value only for, n>~range of correlation, in which case the "surface" terms due to correlations with the system exterior to n become negligible in comparison to the number of interior terms. We shall always assume a finite range of correlation and define S I'1 accordingly. It is first observed that in a homogeneous system where only harmonic forces are present that the maximum of (3.1) subject to the dynamical constraints Y,Y,,+m = i m (independent of t) leads to the multivariate Gaussian distribution (2.1) towards which the system asymptotically tends. This is the same distribution which G r a d 6) has found in a generalized microcanonical ensemble where there are a number of dynamical constraints in addition to the total energy. Thus the thermodynamic interpretation of the harmonic evolution is that the entropy increases to its highest possible value within the framework of the mechanical constraints imposed b y the harmonic invariants of the system. When one removes the harmonic invariants b y allowing anharmonic forces, the system may further evolve, increasing its entropy to the final equilibrium value. It is to this evolution that we now turn. It is necessary to consider only small deviations of second moments from their equilibrium value. If deviations from equilibrium are retained to the second order we find from (3.1) (S (") - - S(°))/k = { Z,L1 (y~ - - y~O) _ ¼ Z , L I ( ~ - - y~O) _ ½ E~., (y,y,)2
(3.2)
where the superscript 0 implies the average taken with the equilibrium distribution function. The first two terms on the r.h.s, of (3.2) are unsigned and can be shown to give rise to an entropy flow related to the heat flow in the system. Indeed comparison of (3.2) with the macroscopical expression S - - S O = c~f~, ° d T / T = c~(T - - T o ) I T ) - - ½ c~(T - - To)~/T~
(3.3)
identifies the first two terms in (3.2) as the change in entropy due to a reversible change in temperature. We m a y also study the significance of the different terms in (3.2) b y the following method. We m a y split the change of entropy into two contributions: the entropy flow d , S / d t and the entropy production d i S / d t 6) 7) d S / d t =: d, S / d t + d , S / d t (3.4) The second principle requires that d,S/dt > 0
(3.5)
IRREVERSIBLE
PROCESSES PART VI : THERMAL CONDUCTIVITY OF CRYSTALS
267
Let us write the equations of motion including, in a formal way, the an~ harmonicity (cf. t) dy,/dt = ½(Y~+I -- Y,-1) + (dY,,/dt),,h
(3,6)
The evolution of the second moments in time is then d y . y . + m / d t -- ½ [Y.Yn+,~+t -- Y . - t Y . + ~ + { lYe+mY.+1 -- Y-Y~+,~-I] + + (dy,,y,,/dt)a,, h
(3.7)
If (3.7) is substituted in (3.2) we get d , S " / d t = E,.~ y, yj (dyj y,/dt)..~
(3.8)
where we have included all boundary terms in the entropy flow and isolated the entropy production. For a homogeneous system (with n >~ range o f corralation K) using (2.9) we obtain diS'*/dt ~- -- n Ei. ~ ~i~ I i l s
(3.9)
From (2.8) it can be deduced that the eigenvalues of Bk,, are negative, which is physically obvious. Consequently the eigenvalues of the A~jdefined by (2.10) are also negative and (3.9) is a positive definite form as it must be. 4. Stationary state. In discussing the possibilities of a stationary state one fact becomes immediately apparent. The stationary state cannot exist indefinitely. For eventually the thermal gradient will be smoothed out by harmonic forces. If we confine our attention to a subsystem of length N , over which there exists a temperature gradient but no other inhomogeneities this smoothing takes place for t ,-~ N in the middle of the system.Thus for t ~ N the center of the system maintains a constant gradient. At the same time, the time must be sufficiently long so that equs (2.5) are valid i.e. t >~ 1/co. In this intermediate time scale the possibility of a stationary state presents itself. We further suppose that aside from the presence of the systematic temperature gradient all other moments are initially but irregularly inhomogeneous and are rapidly smoothed out in t = 0 (1/o~). In this case we can write down eqs (3.7) in a time smoothed sense using (2.9)
dil/dt=
d i k / d t = Y'~ ~kz iz
(4.1)
½ (Y,+I * - -- Y~) + Z,~>l ~t, 1,
(4.2)
dye/dr = 0
(4.3)
The average value of the energy y~ is not changed by the decay equations (2.9.). In the stationary state eqs (4.1) and (4.2) have their left hand sides equal to zero and it is possible to solve for the I's in terms of grad -?0 = ½(Y~+I--Y~).~
268
R. B R O U T A N D I. P R I G O G I N E
In particular, we have for the heat flow i 1 = ½(A-1)n grad h
(4.4)
which is Fourier's law. (A-1)n is the 11 component of the inverse of the matrix [lyt~][. B y these considerations we have therefore identified the thermal conductivity coefficient n as = ½(A-1)u
(4.5)
As ;tij is proportional to I o, we have n ,-~ 1/T. As a first approximation we m a y expect n = ½ 1/;tn
(4.6)
where 1/,ln is the relaxation time of the heat flow. w e see that the thermal conductivity is directly expressed in terms of the coefficients a~j related to the second moments i~ and not in terms of the relaxation coefficients Ckk,k,, appearing directly in equations (2.5). As can be expected in the high temperature region (classical mechanics) the inverting processes do not play any special role in a. Their existence must however be assumed to insure the ergodicity of the system 8) 4). If the entropy production (3.9) is minimized subject to constant heat flow eqs (4.1) become zero b y minimization and eq. (4.2) is zero b y constraint. We have thus the well known theorem 6) that the entropy production in the presence of a gradient is minimum in the stationary state.
5. Fourier's equation. In order to derive the phenomenological equations of thermal conductivity out of the stationary state, let us reconsider (3.7). We divide our system of length N into a large number of domains each of which is large with respect to microscopic dimensions and small with respect to macroscopic variations. Since the domain contains a large number of atoms we may use our time smoothed and space smoothed equations of anharmonic decay assuming that all moments b u t the energy are homogeneous within a domain. In this way (3.7) m a y be written in a domain D with boundaries n and n + D d/dt (Y~,,DY,Y,+k)= (d/dt (~,d) Y,Y,+k))anh
d/dt (XieD Y,Y,+x)
= ½ (yn+ -3
-
-
-
+
k > 1
(d/dr
(5.1) (5.2)
feD
d/dt (Y',,D
= ½
Y.Y -I)
(5.3)
Now the structure of equations (5.1) and (5.2) differ in character from (5.3) in that (5.3) only involves surface terms of the 'domain, representing the heat flux in and out of the domain. If one adds up (5.3) over all domains one is finally left with a big difference of two surface terms. Thus small inhomogeneities in the heat flow Y~Yn+I become multiplied into macroscopic pro-
I R R E V E R S I B L E PROCESSES PART VI : THERMAL C O N D U C T I V I T Y OF CRYSTALS 2 6 9
portions by the totality of equations (5.3). In (5.1) and (5.2) this successive addition of terms is each time compensated by a local interior term representing anharmonic decay. These local terms then present the possibility of local equilibrium. In (5.2) we see this brought about by the compensation of the harmonic propagation of the gradient against the local decay of heat flow. Note that in the absence of anharmonic terms such a compensation is impossible; the heat flow m a y continue unimpeded. In (5. l) and (5.2) we substitute the equations of anharmonic decay which are valid in the smoothed region D to give (we call T~ the value of I k in D) diD~ d t = X~=z ~k~ _]D k > 2
(5.3)
diD/dt = ½ grad i D + XT=1 21zi~)
(5.4)
diDo/dt = (l D+z -- l D-z)
(5.5)
where we have assumed all moments approximately homogeneous in D except _70. (5.3) and (5.4) contain only terms in ]D and m a y be solved for ?D in terms of (grad i~) D. For asymptotic terms in the intermediate time scale, the solution tends to the state of local equilibrium and we have I-D , -+½ (A-a)1, (grad i0) D and in particular i D --> ½(A-1)zz (grad i0) D which establishes Fourier's heat law in the domain D again yielding ½(A-1)lz for the thermal conductivity •. Combining with (5.5) we have diD/dt = n [(grad i0) z)+z -- (grad
-T0)D-lj
(5.6)
which is the discrete analogue of the thermal conductivity equation ST/at = -, a~T/a~ 2
(5.7)
The essential argument of the establishment of (5.6) was the concept of local equilibrium in microscopic domains sufficiently large for the smoothed equations of decay to be valid. These smoothed equations of decay counteract the harmonic propagation of the gradient and are responsible for the achievement of local equilibrium. It must be remembered that these arguments are valid only within the intermediate time scale previously discussed. For sufficiently large times t ~-~N all gradients break down under the smoothing out action of harmonic forces. This rather rapid breakdown is due to the high efficiency of conductivity of the thermal bath in which the system is immersed. A more anharmonic bath would give rise to a larger intermediate time scale. Acknowledgments. The research reported in this paper has been made possible through the support and sponsorship of the Air Research and Development Command, United States Air Force, through its European Office.
270
R. B R O U T A N D I. P R I G O G I N E
APPENDIX
In this appendix we wish to point out the necessary generalizations that make the arguments in the text applicable to three dimensions. Many of the considerations have immediate generalizations to the general three dimensional case. The considerations leading to equations of the type (4.6), (4.1), (4.2) are necessarily complicated and we treat the simple case of an cubic crystal with nearest neighbour interactions. This latter supports Umklapp processes and so is sufficiently illustrative for our purposes. We first remark that the three dimensional harmonic problem has been • treated in ref. 2) with the same asymptotic results as in one dimension save for the introduction of quasi-invariants I~. These quasi-invariants rapidly become constant. Local distributions tend to multivariate gaussians and the system becomes homogeneous. Unlike the one-dimensional system homogenity does not imply random phases for there m a y still exist correlations between normal modes with the same value of k b u t of different polarizations. Thus the random phase hypothesis goes further than mere homogenization of the crystal. In the case of random phases the energies of the normal modes are related to the invariants I in the homogeneous system b y 9.) s 1 e~(k) e~(k) I~{ = u .a ua.+.. = Z, Zr=
I~,~1" cos km
which is quite similar to the one dimensional formula and allows one to find the decay equation of I~{ which then take on the same form as in one dimension where now the t . . . also contain the inverse vectors of polarization
4k). The most complicated generalization to three dimensions is the establishment of the analogues of (4.1) and (4.2) in the text. This we proceed to do for the cubic crystal with nearest neighbour interactions and two characteristic interaction constants. The important point is to calculate the time derivatives of energy and heat flow of a region of the crystal due to the harmonic equations of motion. From the outset we will consider a crystal homogeneous in all moments save the energy for such times that the quasi-invariance 3) of the second moments, whose time derivatives do not involve the temperature gradients, has already been established. The Hamiltonian of the system we consider is H
3
,, . + l a - -
.j + (U~+lo---U~)2)J
(A. 1)
to give rise to the equation of motion
dt = (2u~--U~+l~--U~-l~)+e [(2Un--Un+a--Un_l~)+ (2U~--U~_IV-- n--l)] (A.2) *) Here again we use non-dimensional variables, with the variables u, u of the same dimensions as y in the text.
I R R E V E R S I B L E PROCESSES PART VI : T H E R M A L C O N D U C T I V I T Y OF CRYSTALS 271
In these equations u,, is the displacement of the nm atom, with components u~, ua,, u~; 1a is the unit vector in the ~th direction. There are three kinds of second moments that must be considered. Of these the type velocity × velocity and position × position do not involve the temperature gradient and hence are quasi-invariants. Of the moments v e l o c i t y × position, it is seen from (A.2) that only moments of the type "a a u , (u,+l --u~,) involve the temperature gradient and hence all others are quasi-invariant. We shall presently establish this moment as the heat flow. In this way the analogues of (2.2) are established for di~/dt k > 2. We now proceed to find the appropriate relations for energy and heat flow. In order to define the energy content of a region containing (2N) z atoms (we take a cube of length 2N) it is necessary to establish a convention of how to split the potential energy of surface atoms. We shall take the natural definition of } the surface potential energy as belonging to a given volume. With this definition
E(2N)S=
½
Z..,,.,,. N lit2, + Y',,a Aoa (u]
-
-
u~,_la) 2]
i~v~'a7,6= " 1/.~z+N n o n l l = - N A., [(u~o,.a,~V+l~--u~.,.a,x)9--(u~.,.a,_.v--UO.o,.a,_~._l) 2] (A.3) a,B#?
where the matrix A oa is given b y A o a = 1 ~ = / 5 ; A ~ a - - e e : # / 5 . The vectors n in the surface terms have had their components expressed in terms of no, n o, %. In order to find the heat flux through the region it is necessary to differentiate (A.3) with respect to time using (A.2).The resulting expression cancels all interior terms and one is left with an analogue of a surface integral. Thus dEc2N)a/dt = Z.=~ a (q~N -- qr--N-1) (A.4) ,,o,,,a=-N
~ ,,°.aN+l~,~
,,o,,aNJ
(U.a.a,V+l~--U.o.aN)]
(A.5)
(A.5) is the analogue of y , (Y,+I + Y,-1) which was identified as a heat flow in one dimension. In (A.5) this becomes a surface integral In order to establish the analogue of (2.2) for 11, we take the time derivative of (A.5). A tedious calculation in which one uses the property of homogeneity of all moments and neglects terms of 0(1) and 0(N) in comparison with 0(N2), yields
dq~/dt=½ Y,.a . {eZ~ = ~AoaE(u~o+L.aN+I--U~o..a,N+I)~--(U~+~,.a.v--U~o,.a.X)~] + 3 + s Z~= 1 A~
~
~
(Unan#+ l , N + l - - Una, n~,N + l)
+ Z~= a 1 A~ r/t ~ ~ q ~ ~.~,.a,N+I -- /~..,-~,N~
2
--
6
(~.o,.~+1,~ ~.o,.~)~] + (A.6)
Each of the three terms in (A.6) expresses a energy difference. As we have assumed our crystal to be homogeneous in every sense but the presence of a systematic gradient, it is appropriate to assume equipartition between the
272
IRREVERSIBLE
PROCESSES PART VI" THERMAL CONDUCTIVITY
OF CRYSTALS
kinetic and potential energies of a single particle in which case (A.6) becomes
dq~,/dt = const X,,,,,a (Ena, no,N+lr -- En,,na,a, )
(A. I)
which establishes the analogue of (2.2) for (dla/dt)ha~o,~i e. It should be pointed out t h a t the model here t r e a t e d supports inverting processes. For given k along a crystal axis say z the three frequencies of polarization are e sin 2 k, e sin 2 k, sin ~ k which is consistent with the selection rules of inverting processes (see P e i e r 1 s 3)). Received 20-12-55.
1,~EFERENCES 1) K l e i n , G. and P r i g o g i n e , I., P h y s i c a 1~ {1953) 1053. Prigogine, I. a n d B i n g e n , R., P h y s i e a :~! (1955) 299. CI. P e i e r 1 s, R., The QualltUm T h e o r y of Solids, Oxford C l a r e n d o n Press 1955. B r o u t, R. a n d P r i g o g i n e , I., P h y s i c a 2 1 . G x a d, H., Comm. pure and app. Sci..5 (1952) 455. 6) P r i g o g i n e, I., E t u d e t h e r m o d y n a m i q u e des Phfinom~nes irr~versibles, Liege, Desoer (1947). 7) D e G r o o t, S. R., T h e r m o d y n a m i c s of irreversible Processes, N o r t h H o l l a n d P u b l i s h i n g Co. (19521.
2) 31 4) 5)