- Email: [email protected]
Calculation of the thermal conductivity of superlattices by molecular dynamics simulation
Physica B 316–317 (2002) 247–249
Calculation of the thermal conductivity of superlattices by molecular dynamics simulation Brian C. Daly*, Humphrey J. Maris Physics Department, Brown University, P.O. Box 1843, Providence, RI 02912, USA
Abstract We report on molecular dynamics studies of heat flow in superlattices. The computer simulations are performed using classical mechanics with periodic boundary conditions. We present results for the variation of the conductivity with the repeat distance. We discuss the relation of these results to experimental data in the literature. r 2002 Elsevier Science B.V. All rights reserved. Keywords: Thermal conductivity; Superlattices; Molecular dynamics
The thermal conductivity of superlattices is of great current interest due to the use of superlattice structures in semiconductor lasers [1] and thermoelectric devices [2,3]. In addition to these device applications, the study of heat flow in superlattices offers insight into the effect that modifications of the phonon dispersion relation have on lattice conductivity [4]. Measurements on Si/Ge [5] and GaAs/AlAs [6] superlattices have shown that the thermal conductivity in the direction perpendicular to layers is reduced by as much as an order of magnitude compared to the conductivity of the bulk constituents. Part of the decrease in the thermal conductivity can be attributed to the reduction in the group velocity of phonons due to zone-folding [7]. However, quantitative calculations of this effect show that the reduction in thermal conductivity should increase as the thickness of the *Corresponding author. Tel.: +401-863-2642; fax: +401863-2024. E-mail address: brian [email protected] (B.C. Daly).
layers making up the superlattice is increased within the range from one to ten monolayers. Experimental results [5,6] show the opposite effect in this range of layer thickness. However, the extent to which interfacial roughness and other defects in the superlattice affect the experimentally measured thermal conductivity is not yet known. In order to investigate the effects of the different superlattice parameters on thermal conductivity, we have performed molecular dynamics simulations on a simple model of a superlattice and present some of the preliminary results here. The simulations were designed to model heat flow in GaAs/AlAs superlattices. For these simulations we have used a simple FCC lattice of atoms with nearest neighbor harmonic and anharmonic interactions. Calculations based on a model of this type have been presented in papers by Maradudin and co-workers [8], and by Maris and Tamura [9,10]. The parameters that enter for this model are the atomic mass M; the lattice parameter a0 ; and the second and third derivatives of the interatomic potential which we will denote by b and b0 : For
0921-4526/02/$ - see front matter r 2002 Elsevier Science B.V. All rights reserved. PII: S 0 9 2 1 - 4 5 2 6 ( 0 2 ) 0 0 4 7 6 - 3
B.C. Daly, H.J. Maris / Physica B 316–317 (2002) 247–249
GaAs, we select the value of the mass M to be the average of the atomic masses of Ga and As. The parameters a0 and b were then chosen so as to give the correct density and bulk modulus for GaAs. The bulk modulus of AlAs differs from that of GaAs by only a few % and so the same value of b was used for AlAs. Finally, the value of b0 was chosen so as to give a reasonable value for the Gruneisen constant in the two materials. In the simulations performed so far the effect of disorder in the isotopic mass of Ga has not been included. Simulations were performed for lattices of dimensions 3200 20 20, with periodic boundary conditions applied in all three directions. The heat flow was along the long direction (the xdirection) of the lattice. The simulations were performed using personal computers with processor speeds of 600 MHz, and each took between 3 and 5 days to complete. In order to determine the thermal conductivity, we set up an initial temperature distribution of the form Tðx; t ¼ 0Þ ¼ T0 þ DT0 cosð2px=Lx Þ;
ð1Þ
where Lx is the length of the sample in the xdirection, and the magnitude DT0 of the temperature perturbation is 10% of the base temperature T0 : We then monitored the decay of the temperature perturbation as a function of time in order to determine the thermal conductivity. The position and momentum of the atoms are then calculated as time progresses using a simple ‘‘leapfrog’’ algorithm [11]. At selected time intervals, we calculate the energy of each atom in the sample. We then find the average of the energy of the 400 atoms with a given x-coordinate, and from this average, we determine the temperature TðxÞ as a function of x: Finally, we multiply this temperature distribution by the weighting function cosð2px=Lx Þ and sum over all values of x in order to determine the current magnitude DTðtÞ of the temperature perturbation. Results for DTðtÞ=DT0 are shown in Fig. 1. These cooling curves exhibit a number of interesting features. First, we note that DTðtÞ contains a small amplitude oscillation. When the temperature distribution is applied at t ¼ 0; a spatially varying thermal stress is set up. This
1.0 0.8
∆ T (t) /∆ T0
248
0.6
0.4 BULK
1x1
2x2
4x4
0.2 0 0
20
40
60
80
100
TIME (psec)
Fig. 1. Molecular dynamics results at 300 K. Decay curves are plotted for the bulk GaAs lattice, as well as the 1 1, 2 2, and 4 4 superlattices.
excites the sample into a low frequency vibrational mode in which the atoms vibrate in the xdirection. The energy associated with this motion makes a small contribution to the total energy which leads to an oscillatory contribution in DTðtÞ: Neglecting this oscillation, we now try to analyze the cooling curve to determine k: Based on Fourier’s law, we expect that DTðtÞ=DT0 ¼ expðp2 Dt=L2x Þ;
ð2Þ
where D is the thermal diffusivity. If the value of D is chosen so as to give a fit at large t; we find that that the fit based on Eq. (2) gives too high a cooling rate for small t: As we will discuss in a more detailed paper [12], this discrepancy occurs because the Green’s function solution to the one dimensional diffusion equation gives an infinite velocity at t ¼ 0; whereas clearly heat can never travel faster than the phonon velocity v: A modification of Eq. (2) that allows for this is 2 p 4Dv2 t2 DTðtÞ=DT0 ¼ exp 2 : ð3Þ Lx 4D þ v2 t The results for DTðtÞ=DT0 were fit to this form using D and v as adjustable parameters. The value of k was then calculated from D and the specific heat. We have tested the effects of changes in the lateral dimensions of the lattice and have found that the decay rate is unchanged when we increase the lateral dimensions from 20 20 to 40 40. We
THERMAL CONDUCTIVITY (W/cm K)
B.C. Daly, H.J. Maris / Physica B 316–317 (2002) 247–249
300K 400K
BULK 300K
1.0 BULK 400K 0.5
0.2 1x1
2x2
4x4
8x8 16x16 32x32
SUPERLATTICE PERIOD Fig. 2. Thermal conductivity versus superlattice period for 300 and 400 K. The solid lines represent the thermal conductivity calculated from simulations of bulk GaAs at the two temperatures. The closed circles (K) and open circles (J) represent the thermal conductivity calculated from the superlattice simulations at 300 and 400 K respectively.
therefore elected to use 20 20. In tests of the length dependence of the conductivity, we find that even up to a length of 3200 atom spacings, the value of k still has not converged to a stable value. We presume that this is because the length of the lattice in the x-direction is not sufficiently long compared with the phonon mean free path. An estimate of the mean free path for GaAs based on the formula of Stoner and Maris [13] gives a value ( whereas the length of our sample is of B600 A, ( 7168 A. Results for k as a function of superlattice repeat distance are plotted in Fig. 2, along with results for bulk GaAs. The conductivity of GaAs at 300 K is approximately twice the experimental value. The conductivity of the all of the superlattices we have studied is reduced by more than a factor of two compared to the bulk material. For the longer period superlattices the decrease in k relative to bulk GaAs is in reasonable agreement with experimental data. For 8 8, 16 16 and 32 32 superlattices the reduction factor is larger the shorter the period, as is found experimentally. The reduction for the 16 16 at 300 K is a factor of 3.5,
249
compared to a reduction of 4.4 measured for a 17 17 sample [6]. For the short period superlattices, the calculation gives a conductivity that increases as the period becomes shorter, whereas in the experiments the conductivity is found to decrease with decreasing period. This may indicate that for short period superlattices, k is significantly affected by defects at the interfaces. We plan to perform simulations to investigate the possible effects of interfacial roughness in the near future, and to include the effects of disorder in the isotopic mass of Ga. This work was supported by the Air Force Office of Scientific Research MURI grant entitled ‘‘Phonon Enhancement of Electronic and Optoelectronic Devices’’ (Grant No. F4962-00-1-0331). We would like to thank S. Tamura of Hokkaido University for valuable discussions.
References [1] T.E. Sale, Vertical Cavity Surface Emitting Lasers, Research Studies Press, Taunton, Somerset, England, 1995. [2] L.D. Hicks, T.C. Harman, M.S. Dresselhaus, Appl. Phys. Lett. 63 (1993) 3231. [3] P.J. Lin-Chung, T.L. Reinecke, Phys. Rev. B 51 (1995) 13244. [4] S.Y. Ren, J.D. Dow, Phys. Rev. B 25 (1982) 3750. [5] S.M. Lee, D.G. Cahill, R. Venkatasubramanian, Appl. Phys. Lett. 70 (1997) 2957. [6] W.S. Capinski, et al., Phys. Rev. B 59 (1999) 8105. [7] S. Tamura, Y. Tanaka, H.J. Maris, Phys. Rev. B 60 (1999) 2627. [8] A.A. Maradudin, P.A. Flynn, R.A. Coldwell-Horsfall, Ann. Phys. 15 (1961) 360. [9] H.J. Maris, S. Tamura, Phys. Rev. B 47 (1993) 727. [10] H.J. Maris, S. Tamura, Phys. Rev. B 51 (1995) 2857. [11] R.P. Feynman, R.B. Leighton, M. Sands, in: The Feynman Lectures on Physics, Vol. 1, Addison-Wesley, Reading, MA, 1963. [12] B.C. Daly, H.J. Maris, K. Imamura, S. Tamura, Phys. Rev. B, submitted (2002). [13] R.J. Stoner, H.J. Maris, Phys. Rev. B 47 (1993) 11826.
Calculation of the thermal conductivity of superlattices by molecular dynamics simulation Brian C. Daly*, Humphrey J. Maris Physics Department, Brown University, P.O. Box 1843, Providence, RI 02912, USA
Abstract We report on molecular dynamics studies of heat flow in superlattices. The computer simulations are performed using classical mechanics with periodic boundary conditions. We present results for the variation of the conductivity with the repeat distance. We discuss the relation of these results to experimental data in the literature. r 2002 Elsevier Science B.V. All rights reserved. Keywords: Thermal conductivity; Superlattices; Molecular dynamics
The thermal conductivity of superlattices is of great current interest due to the use of superlattice structures in semiconductor lasers [1] and thermoelectric devices [2,3]. In addition to these device applications, the study of heat flow in superlattices offers insight into the effect that modifications of the phonon dispersion relation have on lattice conductivity [4]. Measurements on Si/Ge [5] and GaAs/AlAs [6] superlattices have shown that the thermal conductivity in the direction perpendicular to layers is reduced by as much as an order of magnitude compared to the conductivity of the bulk constituents. Part of the decrease in the thermal conductivity can be attributed to the reduction in the group velocity of phonons due to zone-folding [7]. However, quantitative calculations of this effect show that the reduction in thermal conductivity should increase as the thickness of the *Corresponding author. Tel.: +401-863-2642; fax: +401863-2024. E-mail address: brian [email protected] (B.C. Daly).
layers making up the superlattice is increased within the range from one to ten monolayers. Experimental results [5,6] show the opposite effect in this range of layer thickness. However, the extent to which interfacial roughness and other defects in the superlattice affect the experimentally measured thermal conductivity is not yet known. In order to investigate the effects of the different superlattice parameters on thermal conductivity, we have performed molecular dynamics simulations on a simple model of a superlattice and present some of the preliminary results here. The simulations were designed to model heat flow in GaAs/AlAs superlattices. For these simulations we have used a simple FCC lattice of atoms with nearest neighbor harmonic and anharmonic interactions. Calculations based on a model of this type have been presented in papers by Maradudin and co-workers [8], and by Maris and Tamura [9,10]. The parameters that enter for this model are the atomic mass M; the lattice parameter a0 ; and the second and third derivatives of the interatomic potential which we will denote by b and b0 : For
0921-4526/02/$ - see front matter r 2002 Elsevier Science B.V. All rights reserved. PII: S 0 9 2 1 - 4 5 2 6 ( 0 2 ) 0 0 4 7 6 - 3
B.C. Daly, H.J. Maris / Physica B 316–317 (2002) 247–249
GaAs, we select the value of the mass M to be the average of the atomic masses of Ga and As. The parameters a0 and b were then chosen so as to give the correct density and bulk modulus for GaAs. The bulk modulus of AlAs differs from that of GaAs by only a few % and so the same value of b was used for AlAs. Finally, the value of b0 was chosen so as to give a reasonable value for the Gruneisen constant in the two materials. In the simulations performed so far the effect of disorder in the isotopic mass of Ga has not been included. Simulations were performed for lattices of dimensions 3200 20 20, with periodic boundary conditions applied in all three directions. The heat flow was along the long direction (the xdirection) of the lattice. The simulations were performed using personal computers with processor speeds of 600 MHz, and each took between 3 and 5 days to complete. In order to determine the thermal conductivity, we set up an initial temperature distribution of the form Tðx; t ¼ 0Þ ¼ T0 þ DT0 cosð2px=Lx Þ;
ð1Þ
where Lx is the length of the sample in the xdirection, and the magnitude DT0 of the temperature perturbation is 10% of the base temperature T0 : We then monitored the decay of the temperature perturbation as a function of time in order to determine the thermal conductivity. The position and momentum of the atoms are then calculated as time progresses using a simple ‘‘leapfrog’’ algorithm [11]. At selected time intervals, we calculate the energy of each atom in the sample. We then find the average of the energy of the 400 atoms with a given x-coordinate, and from this average, we determine the temperature TðxÞ as a function of x: Finally, we multiply this temperature distribution by the weighting function cosð2px=Lx Þ and sum over all values of x in order to determine the current magnitude DTðtÞ of the temperature perturbation. Results for DTðtÞ=DT0 are shown in Fig. 1. These cooling curves exhibit a number of interesting features. First, we note that DTðtÞ contains a small amplitude oscillation. When the temperature distribution is applied at t ¼ 0; a spatially varying thermal stress is set up. This
1.0 0.8
∆ T (t) /∆ T0
248
0.6
0.4 BULK
1x1
2x2
4x4
0.2 0 0
20
40
60
80
100
TIME (psec)
Fig. 1. Molecular dynamics results at 300 K. Decay curves are plotted for the bulk GaAs lattice, as well as the 1 1, 2 2, and 4 4 superlattices.
excites the sample into a low frequency vibrational mode in which the atoms vibrate in the xdirection. The energy associated with this motion makes a small contribution to the total energy which leads to an oscillatory contribution in DTðtÞ: Neglecting this oscillation, we now try to analyze the cooling curve to determine k: Based on Fourier’s law, we expect that DTðtÞ=DT0 ¼ expðp2 Dt=L2x Þ;
ð2Þ
where D is the thermal diffusivity. If the value of D is chosen so as to give a fit at large t; we find that that the fit based on Eq. (2) gives too high a cooling rate for small t: As we will discuss in a more detailed paper [12], this discrepancy occurs because the Green’s function solution to the one dimensional diffusion equation gives an infinite velocity at t ¼ 0; whereas clearly heat can never travel faster than the phonon velocity v: A modification of Eq. (2) that allows for this is 2 p 4Dv2 t2 DTðtÞ=DT0 ¼ exp 2 : ð3Þ Lx 4D þ v2 t The results for DTðtÞ=DT0 were fit to this form using D and v as adjustable parameters. The value of k was then calculated from D and the specific heat. We have tested the effects of changes in the lateral dimensions of the lattice and have found that the decay rate is unchanged when we increase the lateral dimensions from 20 20 to 40 40. We
THERMAL CONDUCTIVITY (W/cm K)
B.C. Daly, H.J. Maris / Physica B 316–317 (2002) 247–249
300K 400K
BULK 300K
1.0 BULK 400K 0.5
0.2 1x1
2x2
4x4
8x8 16x16 32x32
SUPERLATTICE PERIOD Fig. 2. Thermal conductivity versus superlattice period for 300 and 400 K. The solid lines represent the thermal conductivity calculated from simulations of bulk GaAs at the two temperatures. The closed circles (K) and open circles (J) represent the thermal conductivity calculated from the superlattice simulations at 300 and 400 K respectively.
therefore elected to use 20 20. In tests of the length dependence of the conductivity, we find that even up to a length of 3200 atom spacings, the value of k still has not converged to a stable value. We presume that this is because the length of the lattice in the x-direction is not sufficiently long compared with the phonon mean free path. An estimate of the mean free path for GaAs based on the formula of Stoner and Maris [13] gives a value ( whereas the length of our sample is of B600 A, ( 7168 A. Results for k as a function of superlattice repeat distance are plotted in Fig. 2, along with results for bulk GaAs. The conductivity of GaAs at 300 K is approximately twice the experimental value. The conductivity of the all of the superlattices we have studied is reduced by more than a factor of two compared to the bulk material. For the longer period superlattices the decrease in k relative to bulk GaAs is in reasonable agreement with experimental data. For 8 8, 16 16 and 32 32 superlattices the reduction factor is larger the shorter the period, as is found experimentally. The reduction for the 16 16 at 300 K is a factor of 3.5,
249
compared to a reduction of 4.4 measured for a 17 17 sample [6]. For the short period superlattices, the calculation gives a conductivity that increases as the period becomes shorter, whereas in the experiments the conductivity is found to decrease with decreasing period. This may indicate that for short period superlattices, k is significantly affected by defects at the interfaces. We plan to perform simulations to investigate the possible effects of interfacial roughness in the near future, and to include the effects of disorder in the isotopic mass of Ga. This work was supported by the Air Force Office of Scientific Research MURI grant entitled ‘‘Phonon Enhancement of Electronic and Optoelectronic Devices’’ (Grant No. F4962-00-1-0331). We would like to thank S. Tamura of Hokkaido University for valuable discussions.
References [1] T.E. Sale, Vertical Cavity Surface Emitting Lasers, Research Studies Press, Taunton, Somerset, England, 1995. [2] L.D. Hicks, T.C. Harman, M.S. Dresselhaus, Appl. Phys. Lett. 63 (1993) 3231. [3] P.J. Lin-Chung, T.L. Reinecke, Phys. Rev. B 51 (1995) 13244. [4] S.Y. Ren, J.D. Dow, Phys. Rev. B 25 (1982) 3750. [5] S.M. Lee, D.G. Cahill, R. Venkatasubramanian, Appl. Phys. Lett. 70 (1997) 2957. [6] W.S. Capinski, et al., Phys. Rev. B 59 (1999) 8105. [7] S. Tamura, Y. Tanaka, H.J. Maris, Phys. Rev. B 60 (1999) 2627. [8] A.A. Maradudin, P.A. Flynn, R.A. Coldwell-Horsfall, Ann. Phys. 15 (1961) 360. [9] H.J. Maris, S. Tamura, Phys. Rev. B 47 (1993) 727. [10] H.J. Maris, S. Tamura, Phys. Rev. B 51 (1995) 2857. [11] R.P. Feynman, R.B. Leighton, M. Sands, in: The Feynman Lectures on Physics, Vol. 1, Addison-Wesley, Reading, MA, 1963. [12] B.C. Daly, H.J. Maris, K. Imamura, S. Tamura, Phys. Rev. B, submitted (2002). [13] R.J. Stoner, H.J. Maris, Phys. Rev. B 47 (1993) 11826.