LIBR.AR.Y OF THE U N IVERSITY Of ILLI NOIS 510.84 no. Z50^^5(b cop. 2. Digitized by the Internet Archive in 2013 http://archive.org/details/thermalconductiv250jack 5 / ^> ^7^ Report No. 250 ^ ,. . , w ^ THERMAL CONDUCTIVITY OF ONE-DIMENSIONAL LATTICES by Z S^ October 1, 1967 E. ATLEE JACKSON JOHN R. PASTA JOHN F. WATERS ,^ -z ^^ THE UBRARY OF THE AUG 13 i::j UKIVERSIIY 0^ ILUNUIS Report No. 250 THERMAL CONDUCTIVITY OF ONE -DIMENSIONAL LATTICES E. Atlee Jackson John R. Pasta John F. Waters"^ October 1, I967 Department of Computer Science University of Illinois Urhana, Illinois 618OI E. Atlee Jackson, Departments of Physics and Mechanical Engineering, University of Illinois, Urbana, Illinois John R. Pasta, Departments of Computer Science and Physics, University of Illinois, Urbana, Illinois John F. Waters, Department of Computer Science, University of Illinois, Urbana, Illinois. (Present address: Hydrospace Research Cor-poration, Rickville, Maryland. ^ ^ TABLE OF CONTENTS Page List of Tables iii List of Figures iv List of Symbols ...... v Abstract vi 1. Introduction. , 1 2. The Lattice and Reservoir Model .... ...... 2 3. Computational Features. ... ..» ....... . 7 ^. Equilibrium Results ................ 10 5. Nonequilibri-um Results. ................... 22 6. Conclusion. ...,,......, ........ 3^ References. ........................... 35 Appendix A. ....,,....,.......,,....,. . 36 Appendix B. . ..................... 37 Appendix C. .......................... . 37 Footnotes ............................ 38 LIST OF TABLES Table Page I 16 II 25 III 26 111 LIST OF FIGURES Figure Page 1 One -dimensional lattice k- 2 Displacement of particles (times three) in the harmonic lattice 12 3 Displacement of particles (times three) in an anharmonic lattice ......... 13 h Time-averaged kinetic energy of the particles in an equilibrium lattice 15 5 The normalized energy flux autocorrelation function ... l8 6 The temperature of the particles in a nonequilibrium harmonic lattice ,,..... 2k 7 The temperature in an anharmonic lattice with defects . . 28 8 The same conditions as in Figure 7^ except the lattice has no defects 29 9 The temperature in an anharmonic lattice with a smaller applied gradient 30 IV List of Symbols 1^ i.e. Greek mu X i.e. Greek lambda K i.e. Greek kappa A eapital Greek delta T >; i.e. Greek tau eapital Greek sigma VT gradient of T f capital seript J ^. H, \, K tilde over letters Abstract A numerical study is made of the equilibrium properties and thermal conductivity of one -dimensional lattices with various degrees of anharmonicity and coupling defects. The equilibrium energy flux autocorrelation function is determined, and its relation- ship to the coefficient of thermal conductivity is discussed. The coefficients of thermal conductivity, together with other non- equilibriiim properties, are determined for several applied temperature gradients and lattice conditions. VI 1. Introduction Despite the numerous studies of lattice thermal conductivity since the classic investigation of Peierls [l], there remains a n"amber of unresolved conceptual and analytic problems [2,3] • The theoretical clarification of these problems is impeded in the case of real solids by effects resulting from complicated interatomic forces, lattice structures, impurities, defects, and quantum effects--all of which are largely secondary to many of the conceptual difficulties, and even to many of the analytic problems. The purpose of the present study is to report on mamerical calculations of the thermal conductivity and equilibri\im fluctuations for some simple one -dimensional lattices which have a minimum number of complicating factors. Thus our present objective is not to reproduce results which can be better obtained from real experiments on real solids, but to obtain information about systems which are simple enough for future detailed theoretical analysis. It is hoped that by simplifying the system, and obtaining detailed information which is impossible to obtain by physical experiments, the present results will aid in clarifying some of the conceptual' aspects of lattice thermal conductivity and act as a detailed check of the standard analytic approximations which are used in the case of real solids. It should be made clear, however, that the degree to which studies of the present type can aid in clarifying theoretical questions is rather limited. The theoretical idealization of infinite lattices can only be roughly approximated by numerical calculations, which are at present limited to rather small lattices ( 100-1000 particles) due to computation time. While this size limitation does not appear to have any significant effect on the intensive nature of the coefficient of thermal conductivity, it does have a profound effect on the "behavior of the energy fl-ox autocorrelation function. Because of this, it is very difficult, if not impossible, to check the very interesting and important theoretical relationships arising from fluctuation- dissipation theories. On the other hand, the very occurrence of these difficulties indicates certain limitations resulting from idealizations in the standard fluctuation-dissipation theories. The model lattice and reservoir used in the present calcula- tions are described in the following section. Some aspects of the computations are discussed in section 3- In section k a n'omber of equilibrium properties of this model are reported, including the behavior of the energy fl\xx: autocorrelation fimction. Section 5 contains the results of nonequilibrium situations involving various applied temperature gradients and anharmonicities. 2. The Lattice and Reservoir Model 2 The system which is studied is a linear monatomic lattice with the Hamiltonian N N-1 H i=l i=l = Y ^ mq^ + Y \l n(q. ,-q.-£ f - 7 V (q. ,-q.-^ )' Z_i 2 1 /__, 2 ^^1+1 ^10 3 11+1 1 o + r x- (q. n -q.-i ) h 1^^1+1 ^1 o (1) where q. is the location of the i'th particle, and H is the equilibrium separation distance. Introducing the dimensionless time t \l\i/in. , displacements from equilibrium x. = {^./l- ) - (i-l), and anharmonic coefficients A.. = X. i /\i, k. = k.I /u, the dimensionless Hamiltonian, 1 1 o' ' 1 1 o' ' ' H = H/|i£^, is N N~l H = y ^x^+ y [i (x. -x.f - \ X.(x. -x.)5 4- I K.(x. ~x.)^ /_, 2 1 Z_.[2'i+1 i' 3 I'l+l 1 4 I'l+l 1- i=l i=l (2) Typically the values N = 100, A.. ~ 10, and k. ~ 67 were used in the calculations. The relationship "between these nondimensional units and typical physical lattices is discussed in Appendix A. While the present lattice is monatomic, it can exhibit many of effects of "defects" through the variable anharmonic coefficients X. , and K. . This way of introducing "defects" has the analytic advantage of not complicating the frequency spectrum of the normal modes, but it increases the variety of normal mode interactions [h]. It, of course, does not represent the usual defects in real solids. In fact the present method is intended to circumvent this complication while, at the same time, probing the role of the anharmonic terms in irreversible processes. The lattice is located between two thermal reservoirs which are separated by a distance (N-1)j? or less. If the distance is (N-1)£ , as shown in figure 1, then the lattice just fits between the reservoirs with no applied pressure at temperature T = 0, In this case x = and X-j. = corresponds to contact between the left and right particles and their respective reservoirs. To increase the contact between the reservoirs and the lattice, the case in which the reservoirs were moved inward slightly so that x-, > 0. 5 and x < -0. 5 was also considered. It (T - T^ ) for particle 1, and v < (T = T^) for particle N. This corresponds to Maxwellian reservoirs with particles of the same mass as the lattice particles. The factor v in (3) takes into account the fact that a slow reservoir particle is less likely to "be at the wall when the lattice particle arrives there (the same effect as in molecular _2 effusion). The mean value of T was taken to he 10 , which corresponds to a physical temperature of roughly 300 K (see Appendix A). In order to obtain temperature gradients which are observable over the background statistical fluctuations, it was necessary to take 2(T„ - T^ )/(T_ + T^. ) > 0.^. Therefore the temperature differences were of the same order of magnitude as the mean temperature. However, within the statistical accuracy of the present study (~ 10'^), no consistent dependency of the coefficient of heat conductivity on the gradient was observed. The total dimensionless energy flux (instantaneous) in the lattice is given by [5] N-1 r J(t) = - > - (x. ,+x.)mx. -X.) - A..(x. -X.) + K.(x. -X.)' ^ Z_i 2 ^ 1+1 ±' y 1+1 \' 1^ 1+1 ±' 1^ 1+1 1^ i-1 *- (^) This expression was used to evaluate the energy flux autocorrelation function for the equilibrium system. The harmonic approximation of J(t), namely N-1 J^(t) = - > ^ (x. ,-x.)(x. -X.) (5) H ' Z_i 2 ^ 1+1 1''^ 1+1 i' ^ ' 1=1 which is frequently used in theoretical investigations [1,3^5^6], was also examined. In the cases studied, J-rjCt) was not found to be the n dominant part of j(t)--which casts considerable doubt on many theo- retical treatments of the energy flux. The average heat flux, J, through the nonequilibrium lattice was determined directly from the time-averaged energy exchange at the two reservoirs. This was used to obtain the dimensionless coefficient of heat conductivity, K , through Fourier's law J = -K^ ^T (6) The temperature of each particle in the lattice is defined to be twice 3 the time-averaged kinetic energy T. = v^(t) (7) The internal temperature gradient (dimensionless) in (6) is determined by a fit of the values of T. to a straight line. The internal gradient is always considerably less than the applied temperature gradient, (T - T^ )/lOO, because of large temperature differences between the lattice ends and the reservoirs. It should be also emphasized that (t) is simply a definition of T., and does not imply local equilibrium. It would be much better, but rather difficult, to relate the local temperature to energy (or velocity) distributions. What this relation- ship should be is, however, not entirely clear--the concept of local equilibrium may be too naive for a nonequilibrium lattice. 5. Computational Features The problem was programmed in FORTRAN II except for the inner loop integration of the dynamical equations. This Runge-Kutta integra- tion procedure was written in NICAP, the assembly language for the University of Illinois computer, ILLIAC II. It was found that the actual running time on ILLIAC II was almost exactly equal to the dimensionless time of the problem itself. The boundary interactions play an important role in the selection of the method of integration of the 2N equations of motion which derive from the dimensionless Hamiltonian (2); V, = F. -, - F. , X. = V. 1 1+1 1 ' 11 2 5 where F. =Ax.-\.(Ax,) +/<.('^x.) with/!^>c.=x. -x, and the con- 1 111 11 1 1 1-1 dition of free ends requires Ax, = Ax., -, = 0. Each impulsive inter- 1 N+1 ^ action with the reservoir erases the past history of the particle motion and necessitates a restart of the integration- In order to allow high interaction rates with the attendant better statistical behavior, the system is studied under a high pressure configuration with reservoir interactions every 10 to 20 iterations of the difference equations. Rather than mix integration schemes a Runge-Kutta method was adopted for the entire problem with no multistep intermediate calculation. In deciding on the step size and the order of the method to be employed there is the usual compromise between accuracy and running time. A fifth order Runge-Kutta method due to Butcher has been shown [7] to have the best accuracy of a typical set of integration methods for a fixed computing time on this problem and in a suitable error range. A step size of 0.1 time uniiss was chosen. It is by no means clear that this selection is optimal. The strategy adopted is a compromise of conflicting requirements. On the one hand, an important characteristic of the system under study is the possibility of performing an energy check by balancing the reservoir energy fluxes against the calculated energy resident in the lattice. The change in internal energy will be exactly accounted for by the external flux except for errors in the integration which are mainly truncation errors in our computations. Roundoff error is not significant because of the IJ-significant-figure precision of the ILLIAC IE computer used in this work. Thus, by maintaining high accuracy, the energy check can be used to detect possible machine errors. On the other hand, if we relax the requirement of high accuracy in the energy check and use a longer time step with perhaps a less accurate integration scheme the whole effect of truncation error could be compared with a physical experiment with leakages between the system and an external heat bath. As long as these leakages are small compared to the main fluxes in the system the steady state configuration will not be sensibly different from the case of an isolated system . The statistical errors which arise because of the small number of inter- actions with the reservoir per hour of computing time, about 2000 in our 8 ccmputations, can be reduced by increasing the time step. Since the statistics will improve roughly as the square root of the number of interactions or dynamical time a twofold improvement in statistics will require a fourfold increase in the time step. Using a fifth order method this means an error term three orders of magnitude larger for a halving of the statistical error. The advantages of a small error term with a good energy check leads to the step size of 0.1 units we selected. The number of interactions per unit of computing time can be increased by reducing the n\amber of particles. No improvement is realized by this, however, since the lifetime of these pulses is shorter and we have fewer particle parameters to determine temperature gradients and other functionals. It would appear that the information to be obtained is a function of the computing time alone. A rough feeling for the computation involved can be had by observing that for 100 particles we have 200 equations with 5 Runge- Kutta steps per equation. A single time step of 0.1 units takes about 0.1 sec of computing time so each basic integration step takes 100 micro- seconds per particle including all overhead in the problem and in addition to the necessary function evaluation and steps in the integra- tion itself. Again, since we are interested in statistical quantities, dramatic improvement on problems of this type must await the advent of new parallel computers of the ILLIAC IV class. h, Equililprium Results To facilitate the discussion of the results in the remaining sections, it is useful to introduce some terminology and basic numerical values: N = 100; K. =%}?. (Appendix A); T = ^ (T^ + tJ = 0.01 ' 1 3 1 av 2 R L "no defects" means A.. = \ = 10 1 "defects" means that \% of the A., have a value of 1.5^, 15^ have a value of 0.5^., and ^QPJo have the value of 10 (Appendix B) "pressure applied" refers to T = and means that x, > 0.5, =<„<-0.5 "no pressure applied" means that x > and x < "force" (on each reservoir) equals the time-averaged momentum exchange between the reservoir and the end lattice particle CFL and CFR is the collision frequency of the left and right particle with its respective reservoir "Applied VT" = (T^ - T )/l00; Internal VT is obtained from a K li fit of the T. to a straight line. These values will be implied in all the following cases unless other- wise stated. The time averages were over various times, up to ^0,000 seconds (dimensionless) , which represents ^00,000 integration steps. It should be noted that ^0,000 seconds is the time it takes a sound wave to cross the harmonic lattice UOO times and roughly TOO times in the anharmonic lattice. Some of the equilibrium properties of the harmonic and anharmonic lattice were determined, both for comparison with the non- equilibrium cases and because of their intrinsic importance. 10 Typical examples of the equilibrium dynamics of the lattice _2 between the thermal reservoirs (T_ = T^ - 10 ) are shown in Figures 2 and 3. These figures show the displacement (multiplied by three) of each of the one hundred particles as a function of time. Figure 2 illustrates the behavior of the harmonic lattice. Two features are notable in this figure. First, the lattice spends a considerable amount of time not in contact with the reservoirs (i.e., contracted in length), and tends to "slosh" back and forth between the two reservoirs. Secondly, the very energetic disturbances (which trace a line across the figure) have a velocity very near the expected value of 1.0 (one lattice spacing/sec). However, when two of these disturbances interact, there is a delay in the propagation of approximately 10 seconds (indicated by a shift in the heavy lines where they intersect). Presumably this delay is due to the appearance of high Fourier components during the interaction, with a corresponding reduction in the group velocity. The effect of this delay is that the time required for a disturbance to cross a harmonic lattice is not equal to the length of the lattice divided by the sound velocity, but somewhat longer (depending on the niamber of such interactions). The situation in the case of an anharmonic lattice with defects, shown in Figure J), is quite different. The contraction of the lattice is sometimes quite large, but it lasts for a much shorter period of time. The large amplitude disturbances now cross the lattice with a velocity of 1.68 to 1.55. Not only is the velocity larger than in a harmonic lattice, but there is no longer any delay when two large amplitude pulses intersect. Both of these results indicate that there 11 (U o •H -P -P 05 U •H a o CD -P O) -P w o •H -P W o •H -P 05 p- O P> (U o 05 rH P. w -H n OJ pi •H 3INI1 12 UJ OD 3 UJ _J o < Q. a) o •H -p CtJ CIS •H 0) 0) •H -P w 0) H O •H -P o -p •H 3V^I1 13 has teen a considerable change m the frequency vs, wave length relation- ship (to the extent that this relationship is significant in the nonlinear case-'-e.go^ it will he amplitude dependent.). The implication is that ojCk) is mere nearly linear m k^ and has a larger initial slope than in the harmonic case. These conclusions also agree with previous theo- retical studies of these lattices [^j. The situation is not simple, however, for in the case of the anharmonic lattice with no defects the velocity of disturhances is 1.20 to 1.2^, with large interaction delays (~ lO sec,,}, and considerable contraction. The conclusion is that the defects have a large effect on smoothing out the dynamics, and increasing the velocity of disturbances. The magnitude of the statistical fluctuations which occur, when shorter computations are used, is illustrated by the plot of T. for the various particles (averaged over 8,000 seconds), shown in Figure ^, It can be seen that the average temperature of the particles deviates from the applied temperature by only a few per cent. Other significant results, which can be compared to the nonequilibrium values in section 5> SiVe tabulated in Table I. The force on the reservoirs fluctuates very ,L.ittle for different runs, whereas the collision frequencies and energy flux (which should be zero in the present case) exhibit modest fluctuations over snorter run times. These results are consistent with the assumption that the end lattice particle has a Maxwej.lian distribution. If the particle which strikes the reservoir has the distribution 2 f(v) = (|v|/T') e" Ik -v 72T ' o o o o 00 o o CO o in O o rO m Z UJ _i o h- cr < Si •H ra (U iH O •H ■P Fh 0) -P o d) Si 0) o •H . -P (U (1) O G -H •H -P ^:! P> TIJ H bO S U -H > ^ I iH •H a^ Eh 0) O CM -4- 0) I •H in to Od — o cn 00 N ID lo "" ~ o o o o o 001 X 3dniVd3dl^31 15 TABLE I Pressure _ , Run Defects applied Force x 10 J x 10 CFL CFR time yes yes 6. 71 -0.12 0.270 0.270 8,000 yes no 5-15 -0.75 0.210 0.201+ i+,000 16 as it approaches the reservoir, and the form (5) as it leaves the reservoir, then the average momentum exchange per collision is N/n/2 (n/t' + ^/t ). This should equal (force/CF), The values in Table I yield (force/CF) = 0.2^9 for the longer run, and 0.2^+^ to 0.251 for -2 the short runs. For T' = T = 10 , the theoretical value is v2n X 10 = 0.25. Obviously this excellent agreement is only a weak test of the Maxwellian character of the distribution function--nonethe- less, it is of some interest. Because of the theoretical relationship between the energy flux autocorrelation function and the coefficient of heat conductivity [8,9]^ considerable effort was made to determine this function to a reasonable degree of accuracy. The normalized autocorrelation function is defined by t / t o / ° A(t) = / = lim / j(t)j(t+T)dt // J^(t)dt ' t - 00 ^o /^o (8) where j(t) is given by (5)- The numerical approximation used for (8) was 11,080 / 11, 080 A(t) - ^ J(5n)j(5n+T)/ ^ J^(5n) (9) n:=0 / n-0 so the current was sampled every five seconds, for 55^^00 seconds (t < 500). The result obtained for a lattice with defects and applied pressure is shown in Figure 5» The first minimum in Figure 5? at IT o •H -P cd t-J (U o o o 0) •H 1^ 18 roughly 65 seconds, is very different from the corresponding value of 100 seconds in the case of a harmonic lattice. In the latter case this first minimum is related to the time required for a sound wave to cross the system; and for shorter times A(t) has the very simple approximate form A(t) ~ 1 - 1.7 ^ (0 < T < N) (10) In the case of Figure 5^ the first minimum again corresponds roughly to the crossing time in the anharmonic lattice (as determined from Figure 3). Moreover it is clear from Figure 5 that the initial behavior of A(t) is not simply a linear decrease with time, as in (lO), but that it also contains a part with an exponential time behavior. The hope would be that one could separate the boundary effect producing the linear decrease in A(t) from the intrinsic relaxation-time effect. Geometrically this would mean that A(t) should initially behave in an exponential fashion in a coordinate system which is suitably rotated with respect to the one in Figiare 5- This simple expectation proves to be wrong. Instead, it is found that A(t) behaves exponentially only if the coordinates are both rotated and shifted by a suitable amount. The significance of this shift in the coordinates is not clear at present. In the transformed coordinates the autocorrelation curves is given by A'(t') = .55 e"-°^^'^' (0 < T' < 60) (11) if the transformed axis have the same scale as in Figure 5« If the axis had only been rotated (by nearly t^/^) , then the coefficient in (ll) would be 0.68 rather than 0.55« When (ll) is expressed in terms of the 19 original coordinates, it is found that A(t) satisfies the implicit equation a(t) ^ 0.176 - o.oiott + o.82i^ exp {-.00933t + 1.05 (a(t)-i)3 (12) (0 < T < 60) These results do not seem to give any conclusive answer to the value of the relaxation time for the energy flux correlation (if it exists). The relaxation time from (ll), namely (.01^) =^ 70 seconds, may not be significant because of the appreciable shift between the coordinates. However, it represents the only value which can be assigned at present. The energy flux autocorrelation function was also determined using the harmonic approximation for the energy flux, Eq. (5). While the normalized autocorrelation function bears some resemblence to Figure 5^ the normalization constants were very different in the two cases. Specifically it was found that 11,080 11,080 ^ J^(5n) = 160.0 ; ^ j|(5n) = 30.8 (lO) n=0 n=0 The large difference in these values shows that the harmonic current only represents roughly one half of the total energy flux in the anharmonic lattice. It would appear, therefore, that the frequently used [1,3,5^6] harmonic energy flux bears little relationship to the actual energy fliix in real lattices that have large anharmonicities. According to a n-umber of theoretical investigations [8], the coefficient of heat conductivity and the energy flux autocorrelation function should be related by the expression 20 1 K. 00 / dT (11) '^ NT^ "O where the factor N is essentially the diraensionless length of the lattice. While the question of the size of the lattice does not seem to have been given any particular attention, it is clear that (ll) should be interpreted to hold in the "thermodynamic limit," N -* 0°, in order to insure K is an intensive property of the system. While this is rather trival from a theoretical point of view, it is important in the present context, and probably in real solids. The difficulty, in the extreme, is illustrated by considering a finite , but arbitrarily large, harmonic lattice. In this case the integral on the right side of (ll) can be shown to vanish--the energy flux reverses its sign within the time required for a sound wave to cross the lattice. If (ll) applied to this case, then one would have K = 0--which is quite different from the usual statement [l] that /<„ -* °° in the harmonic limit. Note that taking the thermodynamic limit after the integral does not change these results. On the other hand, if one begins with an infinite lattice (e.g., using periodic boundary conditions), then j(t) is rigorously a constant for a harmonic lattice, and the integral in (ll) is infinite. It is clear from this that some care must be taken in applying (ll) to finite lattices (such as in the present case, or in real experiments). As can be seen from Figure 5^ there remains a strong tendency for the integral in (ll) to vanish, even in the case of strong anharmonicity. This tendency is again a result of disturbances (or energy) being reflected at the ends of the lattice. Such a tendency will persist as 21 long as the correlation time is longer than the time required for a sound wave to cross the system. Moreover, the correlation time should presumably "be of the same order as the lifetime of sound waves fphonons). Since sound waves do, in fact, easily survive while crossing real (experimental size) lattices, the correlation time can be expected to be large compared to this crossing time. Thus the results shown in Figure 5 should represent, fairly realistically, the behavior of the autocorrelation function in real finite solids. Since a very large amount of cancellation still occurs in the integral in (ll), it seems very doubtful that the theoretical relationship (ll) can be expected to apply in many real situations. It will be shown later that (ll) apparently does not apply to the present system. 5. None qui lib rium Results When the temperatures of the two reservoirs are made different from one another, a systematic energy exchange between the ends of the lattice and their respective reservoirs is observed. From these results, it is easy to obtain the time-averaged energy flux, J, through the lattice. The more difficult problem is to obtain a measurable internal temperature gradient. The gradient must be larger than the statistical fluctuations (e.g., shown in Figure ^) , which either requires long run times, large applied temperature differences, and/or large anharmonic coefficients with defects. A number of cases involving increasing an- harmonicities were investigated to .deteDrmine the weakest anharmonicity which would yield an observable temperature gradient. 22 In the case of a harmonic lattice, illustrated in Figure 6, the time-averaged kinetic energy of the lattice particles exhibit a rather remarkable behavior. Figure 6 is for a run of ^,000 seconds, with an applied pressure, and ;,T - 'T ) = 6067 x 10 , The applied temperature gradient is indicated by the inclined straight line in the figure. The average flux through the system is J = -0.2^ x 10 , with a force of .01?^+, and CFR = .06k, CFL = .079 (see Tables I, II and III for comparisons). The most remarkable feature of Figure 6 is the high degree of symmetry in the time-averaged kinetic energies of the particles--despite nearly 60O interactions with asymmetric reservoirs. Moreover, the average temperature of the lattice is significantly lower than ;r (T + T ). Other similar computations indicate that these 2 R L features are real and not a statistical accident. At present there appears to be no theory to explain these features of the none qui librium finite harmonic lattice. The obvious dominance of a few modes in Figure 6 may, however, be due to the "sloshing" effect noted in Figure 1. The symmetry is probably due to the fact that pulses, traveling in either direction, are transmitted without alteration. Thus the time average energy of all particles would presumably be the same if it were not for this "sloshing" effect. The only predictable feat'uT-e of Figure 6 is that there is no 'uniform temperature gradient inside the lattice. Obviously Fourier's law does not apply to such a system. When a weak anharmonicity (X. = l) is introduced in this lattice, the plot of T. remained very similar to the one for the harmonic lattice. The principal difference is that the average tempera- ture of the lattice is closer to the average applied temperature than 25 o o o o 00 •H ^ ,i3 o •H H r»- •H § o 05 u) Si a: •H UJ m ^r + {l/h)K r^. Another method for differentiating between the cubic and quartic terms is to allow the system to occupy a different volume. The results in Table III are for a system with defects but no applied 31 pressure (at T = O), and may te compared with the first two entries in Table II. The increased vol-ume of the system effectively reduces the anharmonicity, resulting in an increased value of /<„• Both the force and the collision frequencies are decreases "because of the weaker contact with the reservoirs. The interesting question arises as to whether the above values of K can be shown to be related to the equilibriiom energy flux autocorrelation function by Eq. (ll). If this relationship is valid, then it can also be written in the form 00 ^ NT ^0 where A(t) is the normalized autocorrelation function, (8). Using the value in (lO), we obtain = l6o/ll,080. Using the observed value of K = 27 (first entry in Table II ), (12) would imply that / A(t) dT ~ 19.^ (13) On the other hand, as already noted, the integral of A(t) as determined from Figure 5 is apparently very nearly zero. In any case, over the interval < t < 500, the integral of A(t), determined from Figure 5^ is at least an order of magnitude smaller than required by (13). It seems very unlikely that the difference can be accounted for by the remaining integral of A(t) for t > 500. From these results it appears that the fundamental relationship (ll) does not hold for the present finite lattice. 32 Histograms of the velocities of the end particles, at the time of impact with the reservoir, were also obtained. Because of the large statistical fluctuations, these did not accurately determine the nature of the distributions, particularly for high energies. Thus, while the distributions frequently had a similarity to the Maxwellian p -V /2T ' form, f(v) = (l/T')ve ' , an elementary check of the data in Tables II and III indicates that this is not exactly the case. If the end particles have a distribution of this form, approaching the reservoirs, the average momentum exchange with the reservoir per collision is VJt/2 ( V T ' + vT), whereas the average energy exchange per collision is T' - T. Here, T' refers to the end particle and T to the reservoir. These quantities should respectively equal (force/CF) and (j/CF). Using the values in Tables II and III it is generally found that different values of T' must be used in order to satisfy both of these conditions. This implies that the distribution function for the end particle as it approaches the reservoir is not Maxwellian. Even when a single value of T • is found to suffice, the time-averaged kinetic energy of the end particle is generally not equal to p (T' + T), The values of this average kinetic energy of the end particles are (in the order of entries in Tables II and III) \ .85 ' .87 .80 .87 .9^ .80 .80 .88 ^N 1.19 1.13 1.21 1.13 1.10 1.24 1.19 1.10 No systematic method of correlating these various quantities has yet been found. 33 6. Conclusi on The present study has not only obtained a numlDer of results which may he compared with existing theories, but it has also indicated several limitations of these theories and the present approach. The principal limitation of the present method arises from available compu- tation time (or, what is the same, lattice size). While the present lattice size does not appear to significantly affect the coefficient of thermal conductivity, it does have a strong influence on the form of the energy flux autocorrelation function. The comparison of this function with the coefficient of thermal conductivity must await a theory which explicitly takes into account the finite size of the lattice. Another limitation of the present approach concerns the rather poor statistics, particularly as regards the velocity distribution function of the lattice particles. It would be of considerable interest, for example, to know this distribution for various locations in the lattice. At present there appears to be no adequate theory for such nonequilibrium distributions. Despite these limitations, sufficient information has been obtained for comparison with present theories, and possibly as a guide for as yet undeveloped theories (e.g., for a finite harmonic lattice between two thermal reservoirs). We would like to express our appreciation to Richard Hicks for his assistance in the final computations in this study. 3^ REFERENCES [1] R. E. Peierls, Ann. Physik ^^ IO55 (1929); Quantum Theory of Solids (Oxford Univ. Press, London, 195^) [2] D. K. Co MacDonald, In Transport Processes in Statistical Mechanics; Editor, I. Prigogine (interscience Publishers, New York, 1958), pg. 63 [5] P. Carruthers, Rev. Mod. Phys. 33, 92 (1961) J. M. Ziman, Electrons and Phonons (Oxford Univ. Press, London, 1962), Chapter VIII [k] E. A. Jackson, J. Math. Phys. k,6Q6 (1963) [5] R. Hardy, Phys. Rev. 132, 168 (1963) [6] R, Brout and I. Prigogine, Physica 22, 263 (1956) [7] J. Waters, Comm. A. CM. 9, 293 (1966) [8] See, e.g., R. Kubo, J. Phys. Soc. Japan 12, 57O (195T); L. P. Kadanoff and P. C. Martin, Annals of Physics 2h, kl^ (1963); J. M. Luttinger, Phys, Rev. 135, AI505 (196^) and references therein. [9'J J. A. McLennan, Jr., Phys. Fluids k, 1319 (1961) 35 Appendix A As a representative system take velocity of sound; v = i N/ia/in = 3 x 10 cm/sec. so' ' _o equililDri-uiii separation; i = 3 x 10 cm -23 mass; m = 30 AMU ~ 5 x 10 gm. Then one dimensionless unit corresponds to 2 2 -12 energy: \ill = mv = ^.5 x 10 ergs. time: vm/Li - I N = 10 seconds. ' o' s -k Force: |j..£ = 1. 5 x 10 dynes Temperature: mv /k = 3.25 x 10 K, where k is Boltzmann's constant. The dimensional form of Fourier's law is cm -sec where the tilde represents a dimensional quantity. The energy flux per linear chain (for simple cubic lattice) is 2 ~ ,2 sec (f o To so the corresponding dimensionless energy flux J = j/{\ill. ) vfi/m satisfies J - -K^V T ~ 2 / ~ ~, 2 , 3x where K = K i /kv . Therefore k T=KT{li /mv j xxos i_Los = I X 10-9 ~j (_^ig^) . 6.67 X 10-^° K f (-^^£^). 3 T cm-sec^ T cm-sec For a typical value of K T = 28 watts/cm., this yields kJ£ = .187. 36 A pressure of 1 atmosphere, corresponds to an average force per linear chain of 9 ^ 10 dynes, or a dimensionless force of 6 x 10 Thus a dimensionless force of .06 corresponds to 10 atmospheres press\ire. Appendix B For possible future theoretical investigations, we list the fifteen nonlinear coefficients with X.. = 15 and K. = 5 in the case of defects. They are respectively i = ^,8, 9, 13, 23, 30, 57, ^^,51^57, 63, 6^1,73, 91,9^ and i - 3, 11, 15,20,25,29,3^,36,^+7, 60, 75, 77, 82,8^4,96,98. Appendix C The numerical value of the normalized autocorrelation function A(5n), for increasing integer values of n, is 1.000, O.696, 0.520, 0.390, 0.280, 0.190, .015, -.05^4, -.127, -.188, -.250, -.31^, ~o356, -.338, -.300, -.250, -.191, -.1^9, --.109, --071, -.033, +.019, .062, .097, .1^3, «l80, .,202, .213, .20^, .185 (for t - 150}. 37 Footnotes 1. Present address: Hydrospace Research Corporation, Rickville, Md. 2. A recent numerical study of lattice thermal conductivity, which investigates the effects arising from a disordered anharmonic lattice, has been made by D. N. Payton III, M. Rich, and W, M. Visscher [h]. 3o It might be mentioned that the ratio of the average potential energy to the average kinetic energy is approximately 0.75 when X. = 10 and T = .01, so the time-averaged potential energy differs significantly from T/2. k. Numerical values of A(t) are given in Appendix C. 5. A recent study of a finite nonequilibrium harmonic lattice has been made by Z. Rieder, J. L. Lebowitz, and E. Lieb (Yeshiva University). Their system differs somewhat from the present system, and their results are not similar to Figure 6. 38 / i