Phonon thermal conductivity of a nanowire with amorphous structure

Computational Materials Science 59 (2012) 152–157

Contents lists available at SciVerse ScienceDirect

Computational Materials Science journal homepage: www.elsevier.com/locate/commatsci

Phonon thermal conductivity of a nanowire with amorphous structure Tian Long See, Rui Xing Feng, Cheuk Yu Lee, Z.H. Stachurski ⇑ Research School of Engineering, College of Engineering and Computer Science, Australian National University, Canberra, Australia

a r t i c l e

i n f o

Article history: Received 27 January 2012 Accepted 28 February 2012 Available online 30 March 2012 Keywords: Amorphous nanowire Molecular dynamics Thermal diffusivity

a b s t r a c t Thermal conductivity of a model nanowire, composed of Zr–Ti–Cu–Ni–Be amorphous alloy, has been studied by computer simulations and theoretical calculations. The results from the molecular dynamics simulations are compared to predictions from Fourier continuum mechanics theory, and with published experimental data. Analysis of the theoretical phonon thermal conductivity follows the previously published incoherent particle model. The novelty of this study is in the employment of amorphous structure, lacking any order or superlattice. The simulated thermal conductivity is significantly lower than that measured by experiments on bulk alloy. It appears that amorphous structure and side-wall scattering reduce thermal diffusivity significantly. Velocity auto correlation time constant increases during heating cycle in proportion to the ratio of atomic weight divided by atomic scattering cross-sectional area. Ó 2012 Elsevier B.V. All rights reserved.

1. Introduction Single or multi-phase nanowires of crystalline atomic structure are grown routinely in many laboratories. No nanowires with amorphous atomic structure have been made yet (however, twodimensional thin films of both crystalline and amorphous structures can be easily produced). In this paper we report on thermal conductivity of an amorphous virtual nanowire simulated by molecular dynamics. Nanowires are a relatively new form of materials, grown on a very small scale [1]. They can be used for lasers and light emitting diodes [2], as well as for the detection of hydrogen gas has been published recently [3], and other uses have been reported. Since the discovery of bulk metallic glasses there have been many experiments and theoretical predictions carried out, especially for the Zr-based Vitreloy [4–6]. It has medical, defence, aerospace and applications [7]. The thermo-physical properties of the material constitute essential information required in designing the processing cycle [8]. Thermal conductivity of the material is an important property and its prediction involves the inverse heat conduction problem where the property is being estimated through the temperature profile taken from experiment [9–11]. Measurements of thermal conductivity of a number of amorphous alloys have been published [12–15]. One-dimensional nanostructures such as carbon nanotubes and semiconductor nanowires have recently received attention. While most of the current research is being focused on electronic and optical properties,

⇑ Corresponding author. Tel.: +61 2 6125 5681; fax: +61 2 6125 0506. E-mail address: [email protected] (Z.H. Stachurski). 0927-0256/$ - see front matter Ó 2012 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.commatsci.2012.02.041

thermal transport is also of interest for basic science as well as for technological applications [16]. In a bulk solid with thermal conductivity, K, the transport of heat energy, expressed as the flux JU, is driven by temperature gradient [17].

J U ¼ K grad T:

ð1Þ

The coefficients of thermal conductivity and diffusivity, D, are related through the heat capacity per unit volume, CV:

D ¼ K=C V ;

ð2Þ

Eq. (1) provides the basis for experimental measurements of thermal properties, and for theoretical calculations based on Fourier’s law of heat conduction. On atomistic scale, the coefficient, K, is related to the characteristic parameters of a solid:

K ¼ C V av ‘ ¼ C V D;

ð3Þ

where ‘ is the mean free path of energy transport without collision, v is the spatial average speed of the particles, and a is a geometrical constant, varying from 1/3 to 1.23 [17]. Eq. (3) provides for two new ways to obtain thermal conductivity using molecular dynamics computations: From the velocity correlation function,

Di ¼

1 3

Z 0

1

h~ v i ðtÞ  ~ v i ð0Þi dt ¼

1 2 h~ v ð0Þi 3 i

Z

1

zi ðtÞ dt;

ð4Þ

0

where ~ v i ðtÞ and ~ v i ð0Þ are the velocities of the ith atom at time t and zero, respectively, and h~ v i ðtÞ  ~ v i ð0Þi=h~ v 2i ð0Þi is the normalised velocity autocorrelation function (VACF). And from Einstein relationship for mean square displacement (MSD) of individual atoms:

T.L. See et al. / Computational Materials Science 59 (2012) 152–157

Di ¼

1 dðMSDi Þ lim ; 6 t!1 dt

MSDi ¼

1P 2 ½~ r i ð0Þ ; r i ðtÞ  ~ N i

ð5Þ

where ~ rðtÞ and ~ rð0Þ are the positions of the ith atom at time t and time zero, respectively, and N is the number of atoms. MSD represents the random walk of the mean free path, ‘. From published thermal conductivities of the pure constituents one can calculate the effective value in series and in parallel arrangements of the pure phases, giving 118 and 37 W/m K, respectively, at 300 K. The experimentally measured value for this alloy is 5 W/m K [6,14], significantly lower than the above values. For other amorphous metallic alloys the value ranges from 4 to 12 W/m K [12,13]. We use molecular dynamics simulations to estimate the thermal conductivity of this amorphous alloy sample in the shape of a nanowire. The atomistic simulations of thermal conductivity predict much lower values than those measured experimentally for the bulk metallic glass. We consider the reduced thermal conductivity in terms of the local character of atomic environments and incoherent scattering mechanisms. The nanowire is the starting point for the molecular dynamics computations, and its analysis follows the method described by Dames and Chen [18]. 2. Starting material and system specification The description of an ideal amorphous structure of the bulk metallic glass has been published elsewhere [19]. The nanowire model with atomic composition: Zr41.2%, Ti13.8%, Cu12.5%, Ni10%, Be22.5%, was extracted from the central portion of the Round Cell, with dimensions of 2 nm in diameter and slightly over 6 nm in length. It contains approximately 1000 atoms; it is shown in Fig. 1. The geometrical model of an ideal amorphous atomic structure is obtained by addition of spheres onto an initial random seed cluster. The result is a topologically random structure. The pair distribution function (PDF) for such an ideal packing is shown in Fig. 2a. It was obtained for the Round Cell (containing 106 atoms) using the following expression [19,20]:

gðr þ DrÞ ¼

nðrþDrÞ 1 ; 4pr 2 Dr q

ð6Þ

where n(r+Dr) is the number of particles (particle centres) within the spherical shell of radius r and thickness Dr, and q is the average numerical density of the solid. The expression allows for different resolution of the PDF depending on the value of Dr. In Fig. 2 Dr = 0.23 nm was used. The first peak (of split appearance) is an envelope of 15 underlying peaks from each of the possible pairs; from the closest separation of Be–Be (at 0.228 nm) to the most distant of Zr–Zr (at 0.318 nm), each of different intensity related to applicable concentrations. Since the structure is completely random, there is no split in the second peak, usually associated with fragments of fcc or hcp arrangements.

Fig. 1. A side view of the cylindrical model of the nanowire. Axes scale in (nm). Atom colours: Zr – green, Ti – red, Cu – blue, Ni – grey, Be – yellow. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

153

The nanowire, initially having the ideal amorphous structure, was subjected to molecular dynamics. During the run the purely random packing of hard spheres altered slightly because of the assigned enthalpic interactions between the atoms, which vary for the different atomic pairs [21]. The probability of contacts with Zr has increased at the expense of contacts between atoms of smaller atomic radius. As a result the PDF changed slightly to that shown in Fig. 2b. The shape of this peak is almost identical to that shown in Fig. 4 in a publication by Hui et al. [22] showing PDF for the same alloy, obtained independently, thus supporting the correctness of our model. There is no change of shape in the second peak, indicating that no short range ordering has taken place. There is good experimental evidence for nearest neighbour coordination in this alloy. Local clusters of atoms have a wide range of neighbour’s numbers as well as variation in atomic species [21,22,19]. Therefore, it can be assumed to a good approximation that the interactions between atoms are non-directional and are reasonably described by the LJ potentials. As a further verification of the atomic structure we have carried out experimental X-ray scattering on the amorphous alloy, and also simulated X-ray Debye scattering on its IAS model. The results, characteristic of scattering from an amorphous material, are shown in Fig. 3. The experimental X-ray patterns were obtained on Siemens diffractometer with CuKa radiation, 30 kV/40 mA, 40 min exposure. For the simulated X-ray scattering pattern, shown in Fig. 3b, we have applied computations using Debye equation on the IAS model of 106 spheres, with appropriate atomic scattering factors. The excellent agreement between the two scattering patterns gives credence to the IAS packing as a good model for the real metallic amorphous alloy. To obtain a reasonably smooth pattern the scattering sample should contain at least 109 scattering centres. This is achieved by starting with the 106 particles model, and for each summation cycle shifting the origin from the central to the next (1 + i) atom 1000 times. (The code to carry out this computation is available from the authors on request.)

3. Methodology of simulation The molecular dynamics approach is based on 12:6 Lennard– Jones potential. The initial equilibration was carried out at a constant temperature of 300 K for 7000 cycles and the coordinates of the final system were then used for the thermal experiment. MD computations for the heating cycle were carried out in steps of 0.12  1015 s (0.12 fs), with a total simulation run of approximately 11  1012 s. Simple calculations based on Fourier theory showed that this interval of time should be sufficient for temperature rise of 200 K over the length of the wire. To run molecular dynamics, atoms are assigned random velocities according to Boltzmann distribution corresponding to the given temperature of the system, and directions of displacement chosen at random. The potential of each atom with regard to other atoms in the system is calculated. The force, position, velocity and acceleration of each atom are calculated by solving Newton’s equations. Velocity scaling was applied to each simulation cycle. For studies of thermal conductivity the cylindrical nanowire was completely isolated from the environment on its outer surfaces. To achieve thermodynamic isolation, atoms positioned in outer layers of the surfaces (within 1.5 radius of the largest atom) were immobilised. This had the effect of creating an impervious boundary; no atoms or heat could escape through this layer, which ensured conservation of mass and energy within the nanowire interior. However, the fixed atoms acted as elastic barrier since the forcefield were active. There were 665 atoms in the fixed layer, 64 atoms in the heat zone (s), and 283 atoms in the interior.

154

T.L. See et al. / Computational Materials Science 59 (2012) 152–157

Fig. 2. Characterisation of the nanowire structure by pair distribution function (PDF). On left, PDF for the IAS model of the bulk metallic glass [19], middle, PDF for the IAS model after subjecting it to molecular dynamics run [19].

Fig. 3. X-ray characterisation of the bulk metallic glass. On the left, X-ray scattering profile of the Zr-based metallic glass [23]. On the right, simulated X-ray Debye scattering profile using the IAS model. CuKa radiation.

Heat baths were applied to both ends, and the temperature change was monitored in the central region of the nanowire, as shown schematically in Fig. 4. The initial condition of the system was set to and equilibrated at a temperature of 300 K. Next, 500 K was applied to both heat baths. The temperature of the heat bath was maintained at 500 K for every cycle of the simulation through re-scaling the velocity of the atoms in that region to fit the Boltzmann distribution for the specific temperature. The temperature in the central region, C, was measured for each cycle as a space average:

T obs ¼ hTitime ¼

tobs 1 P TðCðtÞÞ t obs t¼1

ð7Þ

The temperature of the heat bath is maintained at 500 K for every cycle of the simulation through re-scaling the velocity of atoms in that region to maintain the Boltzmann distribution for the specific temperature. The temperature, like entropy, is only meaningful at equilibrium. The temperature in any chosen region of the nanowire is evaluated (for every cycle) by taking into account the kinetic energy of atoms according to the equipartition principle:

T¼

N 1 P mi v 2i 3NkB i¼1

ð8Þ

4. Results 4.1. Fourier calculations of thermal conductivity Fig. 4. The nanowire with insulation layers, heat bath at either end, and temperature sensing central region.

The law of heat conduction is expressed by Fourier Eq. (1), coupled with that of conservation of energy:

T.L. See et al. / Computational Materials Science 59 (2012) 152–157

q Cp

@Tðx; tÞ @ 2 Tðx; tÞ þ Qðx; tÞ ¼K @t @x2

ð9Þ

where Q is the external heat energy input and Cp is the specific heat of the material at constant pressure. The standard method of solution is to assume that the function describing temperature can be represented as a product of two independent functions, T(x, t) = J(x) G(t), where J(x) is only a function of x and G(t) is only a function of t, both satisfying the initial and boundary conditions of the system [24]: The boundary conditions for the nanowire are set as follows:

@T @T ðL; tÞ ¼ ð0; tÞ ¼ 0; @x @x

ð10Þ

where L is the length of the system in the heat flow direction. Since the boundary conditions are time independent, then a nontrivial solution is.

Tðx; tÞ ¼ A0 þ

1 P n¼1

An cos

npx ðnLpÞ2 qKCp t e ; L

n ¼ 1; 2; 3; . . .

ð11Þ

The initial and final conditions at x = 0 are set to: T(0, 0) = 300 K, and T(0, 1) = 500 K. By substitution of the initial conditions into Eq. (11) the constants A0 and An are evaluated as 500 K and 200 K, respectively. At the centre of the nanowire x = 0, therefore the cosine term is always = 1. To a good approximation it is sufficient to take n = 1, and neglect values n > 1 since the exponential term will decay rapidly for higher values of n. Therefore, the final solution to this particular heat conduction problem is given by Eq. (12) below.

Tð0; tÞ ¼ 500  200eðt=sÞ

ð12Þ

where s = (L/p)2(qCp/K). The time constant for the system is estimated as 1012 s = 103 ps. Fig. 5 shows the variation of temperature in the centre of the nanowire, predicted by Eq. (12) for different values of K varying from 1 to 6 W/m K. It appears that a fit will occur for K between the two values. 4.2. Atomistic simulations of thermal conductivity The temperature profiles of the heat bath and the central region of the wire are shown in Fig. 5. The temperature of the heat bath is being maintained at the required level with relatively small fluctuations compared to those of the central region, achieved through velocity scaling.

155

It is observed that the central region of the wire is not reaching the equilibrium temperature in a monotonic fashion, but rather has periodic fluctuations, as expected of a small system, and the temperature rise seems to be govern by a higher order equation rather than the first order equation solved above. The relatively large fluctuations of temperature are acceptable as the system is very small and the heat is being transferred through phonon-like waveforms simultaneously from both ends. Since the system is enclosed, the reflection of heat waves within the system is expected. Hence, superposition of the heat waves with different wave-vector is a possible cause of the fluctuation of temperature at the centre of the wire. To characterise the dynamics of the system we have calculated the mean square displacement (MSD) of atoms as well as their velocity autocorrelation functions (VACF) for individual elements (i.e. Zr, Ti, Cu, Ni, Be), as well as cumulative averages for the central and hot zones. Velocity autocorrelation function is a time dependent function which reveals the underlying nature of the dynamical process. The relative value of VACF, h~ v i ðtÞ  ~ v i ð0Þi=h~ v i ð0Þi2 , starts at 1 at the onset of molecular dynamics simulation and tends to zero at long time. Some of these results are shown in Figs. 6 and 7. Separately, from data of recorded MSD of all interior atoms, plotted against time, the gradient was used in Eq. (5) to obtain diffusivity, found to be approximately 0.18  106 m2/s. Using the published value for heat capacity at constant pressure, Cp = 0.41  103 J/kg K, and the density of the alloy, q = 6.1  103 kg/ m3, we calculate the value of CV = 2.5  106 J/m3 K. Given this value, and using Eq. (3), thermal conductivity of the nanowire is estimated as K = 0.45 W/m K. This is significantly lower than the experimentally measured value of 5 W/m K [6]. Applying Eq. (4) to results similar to those shown in Figs. 6 and 7 gives the diffusivity as 0.153  106 m2/s, consistent with the result calculated above. From other sources the measured diffusivity of the alloy at 300 K is given as 2  106 m2/s [14]. 5. Discussion In a glassy solid the velocity rapidly de-correlates with time, resulting in oscillations that are decaying due to perturbing forces. We take the time at the first negative peak as measure of a time constant of the process. For example, referring to Figs. 6 and 7 for nickel atoms, the time of reversal (lowest point) is taken as 3  1015 s at constant temperature (300 K), whereas the same point is 7  1015 s during the heating cycle. Therefore, the relative change is (7–3)/3 = 133%. In all cases this time constant was found

Fig. 5. On the left, temperature–time profile of the heat bath. On the right, temperature–time profile of the central region of the nanowire, with theoretical curve plots for different values of K according to Eq. (12).

156

T.L. See et al. / Computational Materials Science 59 (2012) 152–157

Fig. 8. Random walks of selected atoms in the nanowire.

Fig. 6. The auto velocity correlation function (VACF) of Ni atoms in the central region T = 300 K.

ally drifting to other positions. For example, the displacements of the Cu-1, Ti-1 and Ti-2 atoms show the essential supercooled liquid characteristics. After a period of localised vibrations near one position, the atoms begin to drift away without any locus. The Zr-1 atom also behaves in a similar way. The gradual displacement (drift) of atoms in the solid amorphous structure involves changing association of an atom from one cluster to another. This can be occur by the a/b mechanisms as described elsewhere [19,20], and it is driven, in the first place, by inhomogeneities of dynamics which in the supercooled liquid can be orders of magnitude slower or faster than dynamics in other nearby clusters [25], but also influenced by individual interatomic enthalpies, given in general by summing all the local interactions:

V ij ¼

Fig. 7. The auto velocity correlation function (VACF) of Ni atoms in the central region. Temperature rising from 300 to 500 K.

to increase during the heating cycle, as shown in Table 1. The increase in velocity de-correlation time can be understood as a result of the heat flux in the direction of the temperature gradient. There appears to be a correlation with the ratio of atomic weight divided by the atomic cross-section for scattering, as can be seen in Table 1. The appearance of the VACF graphs in Figs. 6 and 7 indicates strong interatomic forces, characteristic of supercooled liquid-like behaviour. Observed atomic displacements in the system, examples of which are seen in Fig. 8, show the movements to be markedly different to that in a crystalline structure where the atom’s positions are centred around the lattice points. Without any defined lattice in the amorphous structure, atoms move about randomly, gradu-

P

PP

i

i j>i

v 1 ðri Þ þ

v 2 ðri ; rj Þ þ

PP P

v 3 ðri ; rj ; rk Þ þ . . . ;

ð13Þ

i j>i k>j>i

where the summation extends over constituent atoms of the cluster. The situation is different, however, for the adjacent Cu-1 and Ni atoms (upper left in the figure). Due to the chemical affinity of the two elements they remain bound together despite their oscillations. Zirconium is the element of highest concentration, and therefore, a noticeable corresponding peak should appear in the vibrational spectrum for this alloy. The vibrational density spectrum for the whole alloy can be obtained by Fourier transform of the VACF, with the result shown in Fig. 9. A first approach of its analysis can be as follows. Let us consider a local cluster with Zr as its

Table 1 Selected and measured properties of the alloying elements.

VAC time relative increase (%) At. wt./pr2 (au/nm2) Atomic weight (au) Atomic radius (nm)

Ni

Cu

Ti

Zr

Be

133 1215 58.7 0.124

128 1234 63.5 0.128

120 1215 47.9 0.147

110 1134 91.2 0.160

47 228 9.01 0.112

Fig. 9. Vibrational density of states in the nanowire calculated from Fourier transform of VACF.

T.L. See et al. / Computational Materials Science 59 (2012) 152–157

central atom. The most common coordination of clusters in this alloy ranges from 11 to 15 [22,19], therefore likely to include a full selection of the alloying elements, thus making the potential (forcefield) for a central Zr atom in such a cluster approximately the same. Suppose the bonds have metallic character with bond stiffness of the order of k = 40 N/m, then pffiffiffiffiffiffiffiffiffi ffi the vibrational frequency of Zr can be calculated: mZr ¼ 21p k=m  280  1012 Hz. A prominent broad peak, centred around 400  1012 Hz, is clearly seen in Fig. 9. It is conjectured that this peak is associated with vibrational modes of Zr, and that the breadth of the peak reflects the variety of arrangements of local clusters around Zr as the central atom. Other peaks, seen between 1000 and 2000  1012 Hz, are tentatively ascribed to atomic vibrations in the other types of local clusters, with a higher bond stiffness (for Ni, Cu and Ti) or lower atomic mass (Be). Atomic oscillations in glasses have anharmonic character that can be usually modelled as double well potentials. Allen et al. [26] have reported density of vibrational states for amorphous silicon with two broad peaks, which they call ‘propagons’, ‘diffusons’ and ‘locons’, and make the point that thermal conductivity is most strongly affected by localisation. Phonons have no meaning in a glass and the diffusion modes show undirected displacement patterns. The mode coupling theory [27] is an approach to the study of complex behaviour in the supercooled liquids. It corroborates the view that heat conductivity is strongly affected by localisation, with incoherent scattering, and therefore theoretical interpretation of K is contentious. Dames and Chen [18] nominate confinement effects and scattering of phonons by defects as the main causes of limited thermal conductivity. From their analysis the wavelength of acoustic phonons is, k = hv/ kBT = 1 nm, for a cylindrical nanowire of small dimensions. It should be said that thermal conductivity, K, consists of two components: the electronic contribution, Ke, and the photon contribution, Kph. In the simulation experiment described here, the electronic component is omitted since the model uses a simple molecular dynamics with incoherent particle model for phonon thermal conductivity of amorphous nanowire. Diffuse scattering and incoherent three-dimensional dispersion are justified. Side-wall scattering decreases thermal conductivity by a factor of 2. Clearly, there are many effects that act together to reduce thermal conductivity of amorphous solids, especially when in the form of a small nanowire. These effects will be cumulative and may even have amplifying effect on each other. It is assumed that these account for the factor of 10 or so between the experimentally measured and computationally simulated values of thermal conductivity of this alloy described in this paper. It is conjectured that there may be another factor which contributes to this effect, and that is temperature scaling routinely used in molecular dynamic simulations [28,29]. Defining temperature at this small length scale requires some care because the usual criterion of local thermodynamic equilibrium is violated [18]. We use the averaged energy density based on the condition that energy is conserved in our computations. However, the smallness of the system makes it subject to fluctuations. There are important implications from the fluctuation theorem (FT) [30]. One is that small systems will spend part of their time actually running in ’reverse’. By reverse it is meant that they function so as to run in a way opposite to that for which they were presumably designed. This applies to our nanowire as it is an isolated system (no periodic boundary condition). The algorithm used to simulate and calculate the temperature is based on the velocity distribution of the atoms in the system that fits to the shape of the Boltzmann distribution curve of a certain temperature. As the number of atoms in the system is small, there will be some fluctuations occurring when calculating the temperate as the distribution is not evenly spread out as it suppose to be. Consequently, temperature scaling is applied to the hot zones, where

157

the temperature is defined as constant (500 K). The temperature profile of the hot zone(s) is well maintained at a constant temperature by this method (Fig. 5a). In the heated zone, the fluctuations are much larger, but the temperature rises more or less as predicted by Fourier calculations. Furthermore, it can be insisted that the atomic oscillatory wanderings, shown in Fig. 8, are a physically realistic, as are the increases in time constant of the velocity correlations. At this stage of the experiment, it is hard to draw a definite conclusion as to the magnitude of this effect, and how to test it. 6. Conclusions  Atomistic dynamics simulations of small systems, such as a nanowire, require special attention to the set up the boundary conditions, since the usual periodic boundary arrangement may not be appropriate.  The computation of thermal conductivity through MD simulation in a metallic glass nanowire with amorphous structure gives results an order of magnitude less than experimental value, thought to be due to a combination of structural scattering, side-wall effects and volume limited.  Velocity scaling should only be applied to zones maintained a constant temperature. The magnitudes of both MSD and AVCF are significantly affected by temperature scaling, resulting in a large uncertainty.

Acknowledgement Both TLS and RXF gratefully acknowledge Summer Research Scholarships provided by the College of Computer Science and Engineering at the ANU. References [1] S. Paiman, Q. Gao, H. Joyce, Y. Kim, H. Tan, C. Jagadish, X. Zhang, Y. Guo Y, J. Zou, J. Phys. D: Appl. Phys. 43/44 (2010) 6. [2] J.B. Schlager, N.A. Sandord, K.A. Bertness, J.M. Barker, A. Roshko, P.T. Blanchard, Appl. Phys. Lett. 88 (2006) 213106. [3] R.M. Penner, MRS Bullet. 35 (2010) 10. [4] F. Szuecs, C.P. Kim, W.L. Johnson, Acta Mater. 49 (2001) 1507–1513. [5] H. Choi-Yim, W.L. Johnson, Appl. Phys. Lett. 71/26 (1997) 3808–3810. [6] A. Peker, W.L. Johnson, Appl. Phys. Lett. 63/17 (1993) 2342–2344. [7] M. Telford, Mater. Today March (2004) 36–43. [8] C.L. Chang, M. Chang, Int. Commun. Heat Mass Trans. 33 (2006) 1013–1020. [9] G.C. Bussolino, J. Spisiak, F. Righini, A. Rosso, P.C. Cresto, R.B. Roberts, Int. J. Thermophys. 14 (1993) 525–539. [10] W.K. Yeung, T.T. Lam, Int. J. Heat Mass Trans. 39/17 (1996) 3685–3693. [11] J.H. Lim, C.K. Chen, Y.T. Yang, J. Thermophys. Heat Trans. 15 (1) (2001) 34–41. [12] A. Jezowski, J. Mucha, G. Pompe, J. Phys. D: Appl. Phys. 20 (1987) 1500–1506. [13] C.L. Choy, K.W. Tong, H.K. Wong, W.P. Leung , J. Appl. Phys. 70 (1991) 4925– 4929. [14] M. Yamasaki, S. Kagao, Y. Kawamura, Appl. Phys. Lett. 84 (2004) 130–138. [15] M. Yamasaki, S. Kagao, Y. Kawamura, Scripta Mater. 53 (2005) 63–67. [16] D.G. Cahill et al., J. Appl. Phys. 93/2 (2003) 793–818. [17] C. Kittell, Thermal Physics, John Wiely & Sons, New York, 1969. [18] C. Dames, G. Chen, J. Appl. Phys. 95 (2004) 682. [19] C.Y. Lee, T.R. Welberry, Z.H. Stachurski, Acta Mater. 58 (2010) 615–625. [20] Z.H. Stachurski, J. Mater. 4 (2011) 1564–1598. [21] I. Martin, T. Ohkubo, M. Ohnuma, B. Deconihout, K. Hono, Acta Mater. 52 (2004) 4427–4435. [22] X. Hui, H.Z. Fang, G.L. Chen, S.L. Shang, Y. Wang, J.Y. Qin, Z.K. Liu, Acta Mater. 57 (2009) 376–391. [23] G. Wang, J. Shen, J.F. Sun, B.D. Zhou, J.D. FitzGerald, D. Llewellyn, Z.H. Stachurski, Scripta Mater. 53 (2005) 641–645. [24] Richard Haberman, Elementary Applied Partial Differential Equations: With Fourier Series and Boundary Value Problems, third ed., Prentice Hall, USA, 1998. [25] M.D. Ediger, Annu. Rev. Phys. Chem. 51 (2000) 99–128. [26] P.B. Allen, J.L. Feldman, J. Fabian, F. Wooten, Philos. Magaz. 70 (1999) 1715– 1731. [27] D.R. Reichman, P.J. Charbonneau, Stat. Mech. 25 (2005) 1–23. [28] L.V. Woodcock, Chem. Phys. Lett. 10 (1971) 257. [29] K. Hermansson, G.C. Lie, E. Clementi, J. Comput. Chem. 9 (1988) 200203. [30] D.J. Evans, D.J. Searles, Phys. Rev. E 50 (1994) 16451648.