- Email: [email protected]
Thermal conductivity of amorphous and crystalline thin films by molecular dynamics simulation
ARTICLE IN PRESS Physica B 404 (2009) 1790–1793
Contents lists available at ScienceDirect
Physica B journal homepage: www.elsevier.com/locate/physb
Thermal conductivity of amorphous and crystalline thin films by molecular dynamics simulation Zhengxing Huang, Zhenan Tang , Jun Yu, Suyuan Bai Department of Electronic Engineering, Dalian University of Technology, Dalian, Liaoning 116024, China
a r t i c l e in f o
a b s t r a c t
Article history: Received 6 November 2008 Received in revised form 10 February 2009 Accepted 17 February 2009
Thermal conductivity (TC) of thin films will be influenced by boundary if the thickness is close to the mean free path (MFP). In this paper, we calculate the TC of crystalline and amorphous SiO2 thin films, which are commonly used materials in micro devices and Integrated Circuits, by nonequilibrium molecular dynamics (NEMD) simulations. The calculation temperatures are from 100 to 700 K and the thicknesses are from 2 to 8 nm. For crystalline thin films, thickness is less than MFP, for amorphous thin films, the thickness is larger than MFP. The TC of crystalline thin films reach their peak values at different temperatures for different thicknesses, the smaller thickness, the larger peak value obtained. But for amorphous thin films, the results show that the temperature dependence of thin films is similar to bulk materials. The obtained temperature dependence of the thin films is consistent with some previous measurements and the theory predictions. & 2009 Elsevier B.V. All rights reserved.
PACS: 66.70.+f 02.70.Ns Keywords: Thin film Molecular dynamics Thermal conductivity
1. Introduction Molecular dynamics (MD) is useful to calculate thermal conductivity (TC) of very thin films, such as the thickness is less than 10 nm. Equilibrium molecular dynamics (EMD) and nonequilibrium molecular dynamics (NEMD) methods have been developed to calculate TC of materials. The EMD methods are based on the use of Green–Kubo formulae, TC is obtained from the integral of the heat current autocorrelation function [1–4]. For the NEMD methods, a hot plate and a cold plate are defined and TC is obtained from the Fourier’s law [5–11]. As is well known, for crystalline materials, TCs are proportional to T3 at low temperatures and to T 1 at high temperatures, TCs of most amorphous materials decrease with the decreasing temperatures. For thin films, the dependence may not be always right. Some previous measurements of crystalline thin films show that TCs reach their peak values at different temperatures for different thicknesses [12]. The smaller thickness the larger peak values obtained. But for amorphous thin films, measurements show that the temperature dependence of thin films is the same as bulk materials and little relative to their thicknesses [13–15]. TCs of crystalline and amorphous SiO2 thin films are calculated by NEMD simulations in this work. Temperatures from 100 to 700 K and the
Corresponding author. Tel.: +86 411 84707945; fax: +86 411 84708918.
E-mail addresses: [email protected] (Z. Huang), [email protected] (Z. Tang). 0921-4526/$ - see front matter & 2009 Elsevier B.V. All rights reserved. doi:10.1016/j.physb.2009.02.022
thicknesses from 2 to 8 nm are used. As proved by some other researches, boundary conditions are not important for NEMD simulations [16,17]. Three-dimensional periodic condition is used in all calculations. The long ranged interactions are handled by the Ewald summation method. A parallel computer Lenovo Shenteng 1800 is employed to do the calculations. In all the calculations, a modified LAMMPS (large-scale atomic/molecular massively parallel simulator) code was used. A constant temperature gradient and momentum not conserved (CTGMNC) temperature controlling method was used in all the calculations [5].
2. Physical models The density of amorphous SiO2 made by thermally grown is 2200 kg/m3. In this paper, an amorphous SiO2 structure with the density of 2200 kg/m3 was used and the unit cell is shown in Fig. 1. The cell is a cubic one, the length is 21.4 A˚ and it has 648 atoms. In the simulations, the X and Y directions have only one cell and the Z direction has different cells to build thin films of different thicknesses. Fig. 2 shows a four-cell simulation structure. To compare the TCs of different structures, we used a crystalline SiO2 with the density similar to the amorphous SiO2. The crystalline unit cell is shown in Fig. 3 and its density is 2180 kg/m3. The cell is also a cubic one, the length is 7.16 A˚ and it has 24 atoms. The super-cell is shown in Fig. 4, the X and Y directions have three cells and the Z direction has different cells to build thin films of different thicknesses.
ARTICLE IN PRESS Z. Huang et al. / Physica B 404 (2009) 1790–1793
1791
The potential function for SiO2 is the potential of van Beest et al. [8,18,19]. The short-range force cutoff is set to 8 A˚ and the method used to calculate long-range force is Ewald [20,21]. In all the simulations, the Verlet integration method was used. Equilibrium is determined by the velocity distribution becoming the Maxwell distribution or by the calculated TC becoming constant. The former condition for equilibrium is usually not enough for a NEMD calculation but the latter is used in this study. Fig. 5 shows a calculated TC vs the calculation time. It can be seen that the TC converges to a constant value when the simulation time is longer.
3. Results and discussion
Fig. 1. A unit cell of the amorphous SiO2 simulated.
Fig. 2. The super-cell of amorphous SiO2.
Fig. 3. A unit cell of the crystalline SiO2 simulated.
The calculated TCs of amorphous SiO2 are shown in Fig. 6. The results from Boltzmann transport equation (BTE) [22] are compared with our results in Fig. 6. Since the BKS potential will got a little larger TC, the BTE results are normalized to the bulk TC 2 W/m K [1]. The change trends of results from BTE and MD are the same but a large deviation found in the case of 2.1 and 6.4 nm. Since the phonon concept is questionable for disordered materials, whether the BTE is an appropriate method to interpret the thermal transport still need to be examined. So the deviation from different thermal transport mechanisms is acceptable. The results also tell us that the TCs decrease with the decreasing thicknesses. The variation of the calculated TCs with temperatures is shown in Fig. 7. It can be seen that TCs of the two thinner films decrease with the decreasing temperatures, similar to temperature dependence of bulk amorphous SiO2 thermal conductivities. For two thicker films, temperature dependence of TCs seems not clear. If we accepted a relative larger error of MD simulations, about 15%, the TCs can be seen as temperature independent. The results for crystalline SiO2 are shown in Fig. 8. The results also show that TC decrease with the decreasing thicknesses. The bulk material TC can be deduced from the thin films according to the size effect formula [8,23]. The deduction procedure is shown in Fig. 9 and the obtained bulk TCs are 29.5, 25.6, 10.0 and 8.9 W/m K for 100, 300, 500 and 700 K, respectively. The mean free path (MFP) of bulk crystalline SiO2 can be obtained from the kinetic theory of gases [10] as 66, 35, 10, and 8 nm for different temperatures. The smaller temperatures, the larger MFPs obtained. Thus, all the simulations of crystalline SiO2 are in the situation that thickness is less than or close to MFP. It is obvious that TCs of bulk materials increase with the decreasing temperatures. The variation of the calculated crystalline SiO2 thin films TCs with temperatures is shown in Fig. 10. It can be seen that different temperature dependences for thin films of different thicknesses. The thicker thin films reach their values at the smaller temperature. In the simulated temperature range, TCs of the two thinner crystalline SiO2 films still not reach peak value. From the variance trend, it will appear in a temperature higher than 700 K. By comparing Figs. 7 and 10, a different temperature dependence of TC between amorphous and crystalline thin films can be found. For crystalline thin films, when the thickness is close to MFP, with the decreasing thicknesses, the phonon boundary scattering will dominate over phonon–phonon umklapp scattering which decreases the TC with an increase in temperature. That is to say, in higher temperatures, the umklapp processes dominate and thermal conductivities of different thicknesses are closer. But in lower temperatures, boundary scattering dominate and thermal conductivities will be influenced by thicknesses more remarkable. Thus, thermal conductivities tend to reach different
ARTICLE IN PRESS 1792
Z. Huang et al. / Physica B 404 (2009) 1790–1793
Fig. 4. The super-cell of crystalline SiO2.
Fig. 5. The calculated TC varied with the time used.
Fig. 7. The thermal conductivities of amorphous SiO2 thin films versus temperatures.
Fig. 6. The calculated TCs of amorphous SiO2 and compared with those from BTE.
Fig. 8. Thickness dependence of TC of crystalline SiO2 thin films from MD.
peak values for different thicknesses thin films. A theory calculation of TC of diamond thin films obtained the same result [24]. A similar result has been obtained for silicon nanowires from experimental measurements [25]. For amorphous materials, heat
conduction mechanisms is much complex [22,26,27], the kinetic theory may be over-simplified. If we accepted the kinetic theory here and obtained MFP value from Zeng et al. The thickness of simulated amorphous SiO2 thin films is larger than MFP. From the
ARTICLE IN PRESS Z. Huang et al. / Physica B 404 (2009) 1790–1793
1793
4. Conclusions By MD calculations, TCs of amorphous and crystalline SiO2 thin films were obtained. From the results, TC’s variance with thickness and temperature is MFP dependent. As for the simulated amorphous SiO2 thin films, thickness is larger than MFP, the temperature dependence is similar to the bulk materials. But for the simulated crystalline SiO2 thin films, thickness is less than or close to MFP, the obtained TCs reach their peaks at different temperatures. The thicker films reach their values at the smaller temperatures. This characteristic is important for thermal management and thermal design of micro devices and ICs.
Acknowledgments
Fig. 9. The bulk TCs deduced from TCs of thin films.
Z. Huang thanks John E. Carpenter at the Mayo Clinic College of Medicine for help running LAMMPS and Chantrenne Patrice at the Thermal Center of Lyon for help modifying LAMMPS. The authors thank the financial supports from the National Natural Science Foundation of China (90607003) and 863 High Technology Program (2006AA040106 and 2006AA040102). All the calculations have been done by Lenovo 6800 supercomputer, which is provided by Department of Engineering Mechanics of Dalian University of Technology. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14]
Fig. 10. The temperature dependence of TCs of crystalline SiO2 thin films from MD.
results of this paper and some other measurements [14,15], temperature dependence of TC of amorphous thin films is similar to bulk materials.
[15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27]
A.J.H. McGaughey, M. Kaviany, Int. J. Heat Mass Transfer 47 (2004) 1799. A.J.H. McGaughey, M. Kaviany, Int. J. Heat Mass Transfer 47 (2004) 1783. K.V. Tretiakov, S. Scandolo, J. Chem. Phys. 120 (2004) 3765. S.G. Volz, G. Chen, Phys. Rev. B 61 (2000) 2651. Z.X. Huang, Z.A. Tang, Phys. B 373 (2006) 291. R.D. Mountain, R.A. MacDonald, Phys. Rev. B 28 (1983) 3022. P. Jund, R. Jullien, Phys. Rev. B 59 (1999) 13707. Y.G. Yoon, R. Car, D.J. Srolovitz, et al., Phys. Rev. B 70 (2004) 012302. S. Maruyama, Phys. B 323 (2002) 193. S.H. Choi, S. Maruyama, K.K. Kim, et al., J. Korean Phys. Soc. 43 (2003) 747. F. Muller-Plathe, J. Chem. Phys. 106 (1997) 6082. A.D. McConnell, K.E. Goodson, Annu. Rev. Heat Transfer 14 (2005) 129. D.G. Cahill, M. Katiyar, J.R. Abelson, Phys. Rev. B 50 (1994) 6077. M.B. Kleiner, S.A. Kuhn, W. Weber, IEEE Trans. Electron. Devices 43 (1996) 1602. R.M. Costescu, M.A. Wall, D.G. Cahill, Phys. Rev. B 67 (2003) 054302. J.R. Lukes, D.Y. Li, X.G. Liang, et al., J. Heat Transfer 122 (2000) 536. P. Chantrenne, J.L. Barrat, J. Heat Transfer 126 (2004) 577. B.W.H. van Beest, G.J. Kramer, R.A. van Santen, Phys. Rev. Lett. 64 (1990) 1955. P. Jund, R. Jullien, Philos. Mag. A 79 (1999) 223. A.Y. Toukmaji, J.A. Board, Comput. Phys. Commun. 95 (1996) 73. E.L. Pollock, J. Glosli, Comput. Phys. Commun. 95 (1996) 93. T. Zeng, G. Chen, J. Heat Transfer 123 (2001) 340. P.K. Schelling, S.R. Phillpot, P. Keblinski, Phys. Rev. B 65 (2002) 144306. T. Kalisik, P. Majumdar, IEEE Trans. Compon. Pack. Tech. 30 (2007) 683. D.Y. Li, Y.Y. Wu, P. Kim, et al., Appl. Phys. Lett. 83 (2003) 2934. D.G. Cahill, R.O. Pohl, Phys. Rev. B 35 (1987) 4067. P. Sheng, M. Zhou, Science 253 (1991) 539.
Contents lists available at ScienceDirect
Physica B journal homepage: www.elsevier.com/locate/physb
Thermal conductivity of amorphous and crystalline thin films by molecular dynamics simulation Zhengxing Huang, Zhenan Tang , Jun Yu, Suyuan Bai Department of Electronic Engineering, Dalian University of Technology, Dalian, Liaoning 116024, China
a r t i c l e in f o
a b s t r a c t
Article history: Received 6 November 2008 Received in revised form 10 February 2009 Accepted 17 February 2009
Thermal conductivity (TC) of thin films will be influenced by boundary if the thickness is close to the mean free path (MFP). In this paper, we calculate the TC of crystalline and amorphous SiO2 thin films, which are commonly used materials in micro devices and Integrated Circuits, by nonequilibrium molecular dynamics (NEMD) simulations. The calculation temperatures are from 100 to 700 K and the thicknesses are from 2 to 8 nm. For crystalline thin films, thickness is less than MFP, for amorphous thin films, the thickness is larger than MFP. The TC of crystalline thin films reach their peak values at different temperatures for different thicknesses, the smaller thickness, the larger peak value obtained. But for amorphous thin films, the results show that the temperature dependence of thin films is similar to bulk materials. The obtained temperature dependence of the thin films is consistent with some previous measurements and the theory predictions. & 2009 Elsevier B.V. All rights reserved.
PACS: 66.70.+f 02.70.Ns Keywords: Thin film Molecular dynamics Thermal conductivity
1. Introduction Molecular dynamics (MD) is useful to calculate thermal conductivity (TC) of very thin films, such as the thickness is less than 10 nm. Equilibrium molecular dynamics (EMD) and nonequilibrium molecular dynamics (NEMD) methods have been developed to calculate TC of materials. The EMD methods are based on the use of Green–Kubo formulae, TC is obtained from the integral of the heat current autocorrelation function [1–4]. For the NEMD methods, a hot plate and a cold plate are defined and TC is obtained from the Fourier’s law [5–11]. As is well known, for crystalline materials, TCs are proportional to T3 at low temperatures and to T 1 at high temperatures, TCs of most amorphous materials decrease with the decreasing temperatures. For thin films, the dependence may not be always right. Some previous measurements of crystalline thin films show that TCs reach their peak values at different temperatures for different thicknesses [12]. The smaller thickness the larger peak values obtained. But for amorphous thin films, measurements show that the temperature dependence of thin films is the same as bulk materials and little relative to their thicknesses [13–15]. TCs of crystalline and amorphous SiO2 thin films are calculated by NEMD simulations in this work. Temperatures from 100 to 700 K and the
Corresponding author. Tel.: +86 411 84707945; fax: +86 411 84708918.
E-mail addresses: [email protected] (Z. Huang), [email protected] (Z. Tang). 0921-4526/$ - see front matter & 2009 Elsevier B.V. All rights reserved. doi:10.1016/j.physb.2009.02.022
thicknesses from 2 to 8 nm are used. As proved by some other researches, boundary conditions are not important for NEMD simulations [16,17]. Three-dimensional periodic condition is used in all calculations. The long ranged interactions are handled by the Ewald summation method. A parallel computer Lenovo Shenteng 1800 is employed to do the calculations. In all the calculations, a modified LAMMPS (large-scale atomic/molecular massively parallel simulator) code was used. A constant temperature gradient and momentum not conserved (CTGMNC) temperature controlling method was used in all the calculations [5].
2. Physical models The density of amorphous SiO2 made by thermally grown is 2200 kg/m3. In this paper, an amorphous SiO2 structure with the density of 2200 kg/m3 was used and the unit cell is shown in Fig. 1. The cell is a cubic one, the length is 21.4 A˚ and it has 648 atoms. In the simulations, the X and Y directions have only one cell and the Z direction has different cells to build thin films of different thicknesses. Fig. 2 shows a four-cell simulation structure. To compare the TCs of different structures, we used a crystalline SiO2 with the density similar to the amorphous SiO2. The crystalline unit cell is shown in Fig. 3 and its density is 2180 kg/m3. The cell is also a cubic one, the length is 7.16 A˚ and it has 24 atoms. The super-cell is shown in Fig. 4, the X and Y directions have three cells and the Z direction has different cells to build thin films of different thicknesses.
ARTICLE IN PRESS Z. Huang et al. / Physica B 404 (2009) 1790–1793
1791
The potential function for SiO2 is the potential of van Beest et al. [8,18,19]. The short-range force cutoff is set to 8 A˚ and the method used to calculate long-range force is Ewald [20,21]. In all the simulations, the Verlet integration method was used. Equilibrium is determined by the velocity distribution becoming the Maxwell distribution or by the calculated TC becoming constant. The former condition for equilibrium is usually not enough for a NEMD calculation but the latter is used in this study. Fig. 5 shows a calculated TC vs the calculation time. It can be seen that the TC converges to a constant value when the simulation time is longer.
3. Results and discussion
Fig. 1. A unit cell of the amorphous SiO2 simulated.
Fig. 2. The super-cell of amorphous SiO2.
Fig. 3. A unit cell of the crystalline SiO2 simulated.
The calculated TCs of amorphous SiO2 are shown in Fig. 6. The results from Boltzmann transport equation (BTE) [22] are compared with our results in Fig. 6. Since the BKS potential will got a little larger TC, the BTE results are normalized to the bulk TC 2 W/m K [1]. The change trends of results from BTE and MD are the same but a large deviation found in the case of 2.1 and 6.4 nm. Since the phonon concept is questionable for disordered materials, whether the BTE is an appropriate method to interpret the thermal transport still need to be examined. So the deviation from different thermal transport mechanisms is acceptable. The results also tell us that the TCs decrease with the decreasing thicknesses. The variation of the calculated TCs with temperatures is shown in Fig. 7. It can be seen that TCs of the two thinner films decrease with the decreasing temperatures, similar to temperature dependence of bulk amorphous SiO2 thermal conductivities. For two thicker films, temperature dependence of TCs seems not clear. If we accepted a relative larger error of MD simulations, about 15%, the TCs can be seen as temperature independent. The results for crystalline SiO2 are shown in Fig. 8. The results also show that TC decrease with the decreasing thicknesses. The bulk material TC can be deduced from the thin films according to the size effect formula [8,23]. The deduction procedure is shown in Fig. 9 and the obtained bulk TCs are 29.5, 25.6, 10.0 and 8.9 W/m K for 100, 300, 500 and 700 K, respectively. The mean free path (MFP) of bulk crystalline SiO2 can be obtained from the kinetic theory of gases [10] as 66, 35, 10, and 8 nm for different temperatures. The smaller temperatures, the larger MFPs obtained. Thus, all the simulations of crystalline SiO2 are in the situation that thickness is less than or close to MFP. It is obvious that TCs of bulk materials increase with the decreasing temperatures. The variation of the calculated crystalline SiO2 thin films TCs with temperatures is shown in Fig. 10. It can be seen that different temperature dependences for thin films of different thicknesses. The thicker thin films reach their values at the smaller temperature. In the simulated temperature range, TCs of the two thinner crystalline SiO2 films still not reach peak value. From the variance trend, it will appear in a temperature higher than 700 K. By comparing Figs. 7 and 10, a different temperature dependence of TC between amorphous and crystalline thin films can be found. For crystalline thin films, when the thickness is close to MFP, with the decreasing thicknesses, the phonon boundary scattering will dominate over phonon–phonon umklapp scattering which decreases the TC with an increase in temperature. That is to say, in higher temperatures, the umklapp processes dominate and thermal conductivities of different thicknesses are closer. But in lower temperatures, boundary scattering dominate and thermal conductivities will be influenced by thicknesses more remarkable. Thus, thermal conductivities tend to reach different
ARTICLE IN PRESS 1792
Z. Huang et al. / Physica B 404 (2009) 1790–1793
Fig. 4. The super-cell of crystalline SiO2.
Fig. 5. The calculated TC varied with the time used.
Fig. 7. The thermal conductivities of amorphous SiO2 thin films versus temperatures.
Fig. 6. The calculated TCs of amorphous SiO2 and compared with those from BTE.
Fig. 8. Thickness dependence of TC of crystalline SiO2 thin films from MD.
peak values for different thicknesses thin films. A theory calculation of TC of diamond thin films obtained the same result [24]. A similar result has been obtained for silicon nanowires from experimental measurements [25]. For amorphous materials, heat
conduction mechanisms is much complex [22,26,27], the kinetic theory may be over-simplified. If we accepted the kinetic theory here and obtained MFP value from Zeng et al. The thickness of simulated amorphous SiO2 thin films is larger than MFP. From the
ARTICLE IN PRESS Z. Huang et al. / Physica B 404 (2009) 1790–1793
1793
4. Conclusions By MD calculations, TCs of amorphous and crystalline SiO2 thin films were obtained. From the results, TC’s variance with thickness and temperature is MFP dependent. As for the simulated amorphous SiO2 thin films, thickness is larger than MFP, the temperature dependence is similar to the bulk materials. But for the simulated crystalline SiO2 thin films, thickness is less than or close to MFP, the obtained TCs reach their peaks at different temperatures. The thicker films reach their values at the smaller temperatures. This characteristic is important for thermal management and thermal design of micro devices and ICs.
Acknowledgments
Fig. 9. The bulk TCs deduced from TCs of thin films.
Z. Huang thanks John E. Carpenter at the Mayo Clinic College of Medicine for help running LAMMPS and Chantrenne Patrice at the Thermal Center of Lyon for help modifying LAMMPS. The authors thank the financial supports from the National Natural Science Foundation of China (90607003) and 863 High Technology Program (2006AA040106 and 2006AA040102). All the calculations have been done by Lenovo 6800 supercomputer, which is provided by Department of Engineering Mechanics of Dalian University of Technology. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14]
Fig. 10. The temperature dependence of TCs of crystalline SiO2 thin films from MD.
results of this paper and some other measurements [14,15], temperature dependence of TC of amorphous thin films is similar to bulk materials.
[15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27]
A.J.H. McGaughey, M. Kaviany, Int. J. Heat Mass Transfer 47 (2004) 1799. A.J.H. McGaughey, M. Kaviany, Int. J. Heat Mass Transfer 47 (2004) 1783. K.V. Tretiakov, S. Scandolo, J. Chem. Phys. 120 (2004) 3765. S.G. Volz, G. Chen, Phys. Rev. B 61 (2000) 2651. Z.X. Huang, Z.A. Tang, Phys. B 373 (2006) 291. R.D. Mountain, R.A. MacDonald, Phys. Rev. B 28 (1983) 3022. P. Jund, R. Jullien, Phys. Rev. B 59 (1999) 13707. Y.G. Yoon, R. Car, D.J. Srolovitz, et al., Phys. Rev. B 70 (2004) 012302. S. Maruyama, Phys. B 323 (2002) 193. S.H. Choi, S. Maruyama, K.K. Kim, et al., J. Korean Phys. Soc. 43 (2003) 747. F. Muller-Plathe, J. Chem. Phys. 106 (1997) 6082. A.D. McConnell, K.E. Goodson, Annu. Rev. Heat Transfer 14 (2005) 129. D.G. Cahill, M. Katiyar, J.R. Abelson, Phys. Rev. B 50 (1994) 6077. M.B. Kleiner, S.A. Kuhn, W. Weber, IEEE Trans. Electron. Devices 43 (1996) 1602. R.M. Costescu, M.A. Wall, D.G. Cahill, Phys. Rev. B 67 (2003) 054302. J.R. Lukes, D.Y. Li, X.G. Liang, et al., J. Heat Transfer 122 (2000) 536. P. Chantrenne, J.L. Barrat, J. Heat Transfer 126 (2004) 577. B.W.H. van Beest, G.J. Kramer, R.A. van Santen, Phys. Rev. Lett. 64 (1990) 1955. P. Jund, R. Jullien, Philos. Mag. A 79 (1999) 223. A.Y. Toukmaji, J.A. Board, Comput. Phys. Commun. 95 (1996) 73. E.L. Pollock, J. Glosli, Comput. Phys. Commun. 95 (1996) 93. T. Zeng, G. Chen, J. Heat Transfer 123 (2001) 340. P.K. Schelling, S.R. Phillpot, P. Keblinski, Phys. Rev. B 65 (2002) 144306. T. Kalisik, P. Majumdar, IEEE Trans. Compon. Pack. Tech. 30 (2007) 683. D.Y. Li, Y.Y. Wu, P. Kim, et al., Appl. Phys. Lett. 83 (2003) 2934. D.G. Cahill, R.O. Pohl, Phys. Rev. B 35 (1987) 4067. P. Sheng, M. Zhou, Science 253 (1991) 539.