**3. Predicting thermal conductivity with EMD and NEMD**

In general the heat energy is transmitted through a solid by electrons (mainly in metals) or phonons (mainly in insulators) or other excitations as spin waves. In this paragraph, we focus on the lattice thermal conductivity, which is the dominant mechanism in

proportional to the appropriate transport coefficient. This will be illustrated in the context of the heat transfer in the next sextion. Before focusing on energy transfer applications of MD, let us mention that MD can be used also beyond the linear response domain of out of equilibrium systems. This includes situations of large external forcing such as e.g. polymer melts flowing at large shear rates (Vladkov 2006b). But this includes also systems for which intrinsically the linear relationship between forces and fluxes is violated. As an example consider a nano-structured system at low temperatures for which the phonon mean free path is comparable with the characteristic system size. Heat transfer in such systems is no longer described by Fourier law, but is rather described by an effective conductivity which depends on the strength of the thermal flux flowing across the system. Molecular dynamics has been shown to predict the effective conductivity under conditions where the heat carriers travel ballistically in the system and being scattered only by the boundaries of the

As any simulation technique, molecular dynamics suffers from some intrinsic limitations. The most obvious is the limitation in system size. This may be critical in the modeling of real scale devices with disparate length scales ranging from 1nm to microns. To address such situations, MD should be coupled with a more mesoscopic method, such as the Boltzmann Transport Equations, where the microscopic information on the phonon lifetime is used as

We have also mentioned that a finite system size cuts the long wavelength phonon modes. As we will see, this affects only mildly the measurement of the thermal conductivity unless

The second limitation of MD is its classical nature. Each phonon mode is equally populated, and the heat capacity is given to a good approximation by the Dulong-Petit law. This is a good approximation if one ever considers solids above or just below their Debye temperature. However, many materials do not obey this condition at ambient temperature. Quantum effects in MD may be accounted for in different ways. The phonon lifetimes computed from MD may be used in a Boltzmann transport equation which includes quantum statistics in the phonon occupation number (McGaughey 2004). Quantum effects may be also directly incorporated in the course of a MD simulation, using a Langevin thermostat with a colored noise consistent with a Bose-Einstein distribution for the phonon

Finally, the electronic degrees of freedom are not explicitly simulated in MD. Hence, it is not possible to probe thermal transport in electrical conductors, and in its basic version the contribution of electron-phonon scattering to the transport in semi-conductors is not

In general the heat energy is transmitted through a solid by electrons (mainly in metals) or phonons (mainly in insulators) or other excitations as spin waves. In this paragraph, we focus on the lattice thermal conductivity, which is the dominant mechanism in

**3. Predicting thermal conductivity with EMD and NEMD** 

nanostructure.

**2.1 Limitations of molecular dynamics** 

we couple the system with heat reservoirs.

input in a Boltzmann equation.

modes (Dammak2009).

accounted for.

semiconductors (Yang, 2004). There are three principal techniques used to evaluate the thermal conductivity with molecular dynamic simulations (Allen & Tildesley 1987, Chantrenne, 2007); a. the equilibrium approach based on the Green-Kubo formulae (Frenkel & Smit, 1996), b. the non-equilibrium MD, also called the direct method, which uses a heat source and a heat sink, or the so-called direct method, based on the creation of a temperature gradient, and c. the homogeneous non-equilibrium MD, where a heat flux is induced (Evans, 1982, 1986). The comparison between the two methods has been undertaken by many researchers, who have concluded that the two methods give consistent results (Shelling et al, 2002; Mahajan et al, 2007; Termentzidis et al, 2011b).

#### **3.1 Non-equilibrium molecular dynamics method (NEMD)**

Kotake and Wakuri proposed the direct method (Kotake & Wakuri, 1994), which is similar to the hot-plate experiment setup. A temperature gradient is imposed across the structure under study by allowing thermal power exchange between the heat source and sink and measure the resulting heat flux (Chantrenne & Barrat, 2004a, 2004b, Termentzidis et al, 2009). The thermal conductivity is then obtained as the ratio of the heat flux and the temperature gradient. An alternative, but equivalent way consists in inducing a heat flux and to measure the resulting temperature gradient (Muller-Plathe, 1997). In both cases the system is first allowed to reach a steady state, after which prolong simulations are conducted allowing to obtain correct statistical measurements (Stevens et al, 2007). The NEMD method is often the method of choice for studies of nanomaterials (Poetzsch & Botter, 1994) while for bulk thermal conductivity, particularly of high-conductivity materials, equilibrium method is typically preferred due to less severe size effects. The size effects in the determination of the thermal conductivity using NEMD can be understood using the following simple line of thought (Schelling & al 2002). Imagine a phonon travelling in the crystal between the two heat reservoirs. This phonon can experience at least two kinds of scattering events: either it can be scattered by other phonons travelling in the material, or it can be scattered by the reservoirs which are seen by the considered phonon as a different material with an almost infinite thermal conductivity. The mean free path associated with this former mechanism can be roughly approximated by / 2 *reservoir = L* , where L is the distance between the reservoirs, because on average a phonon will travel ballistically a distance L/2 before being scattered by the heat source/sink. Now invoking Mathiessen rule to predict the effect of reservoir scattering on the overall thermal conductivity, the effective mean free path is given by:

$$\mathbf{1} / \Lambda(\mathbf{L}) = \mathbf{2} / \mathbf{L} + \mathbf{1} / \Lambda(\mathbf{L} \to \infty) \tag{26}$$

where the last term *L* describes the phonon-phonon interaction in a bulklike medium. The length dependence of the thermal conductivity can now be estimated using kinetic theory of gas: *L = cv L* /3 where *c* is the heat capacity per unit volume, *v* is the mean group velocity. Of course, this analysis is simplified because we have assumed that the effect of phonon-phonon scattering can be described in terms of a single mean free path, and we should have rather done the analysis mode by mode (Sellan et al 2010).

However, the previous simple kinetic model predicts:

$$1 \,\,\lambda\left(L\right) = 1 \,\,\lambda\left(L \to \infty\right) + 6 \,\,\left(\text{cvL}\right) \tag{27}$$

Molecular Dynamics Simulations and Thermal Transport at the Nano-Scale 83

Fig. 1. Geometric configuration for the determination of the cross-plane (left) and in-plane (right) thermal conductivity simulations with NEMD and periodic boundary conditions are

Fig. 2. To calculate the thermal conductivity (in the example the in-plane TC) of the

superlattices the energy exchange (b1) and the temperature profile (b2) is extracted for each of the structures (at least 4 structures with increasing length see (a)). Then the inverse thermal conductivity is plotted as a function of the inverse of the size of slab (c). The last diagram helps in calculating the thermal conductivity of an infinitive size superlattice.

used in all directions.

and the size dependence is found to be consistent with MD NEMD simulations (Schelling et al 2002, Termentzidis et al 2009). Because of the algebraic decay, size effects in NEMD are severe, and sometimes it is faster to run equilibirum simulations to determine the thermal conductivity.

We now illustrate the determination of the thermal conductivity using the direct method with the exemple of superlattices composed of a regular alternative arrangement of solid layers, having different physical properties. Thermal transport in superlattices is characterized by the cross-plane and the in-plane conductivities, which correspond respectively to heat flowing in the direction perpendicular or parallel to the interfaces.

The geometry used for the determination of the cross-plane and in-plane thermal conductivities using NEMD is given in figure 3.1.1. Periodic boundary conditions are used in all directions, while in all cases the heat flux is imposed in z-direction. The interatomic interactions are described with Lennard-Jones (LJ) potentials, with energy unit ij=1.0 and length unit ij=1.0. The use of LJ potential is justified by the interest in focusing on the main phenomena introduced by the interface roughness of the super-lattices (for all the following results). The two types of materials A and B may represent respectively the Si and Ge, if one considers Si/Ge-type superlattice (paragraph 5) or the GaAs and AlAs in case of GaAs/AlAs-type superlattice (paragraph 6). The two solids have the same lattice constants and they differ only by the mass ratio of the atoms constituting the superlattices layers. For the Si/Ge-type superlattices the mass ratio is taken to be 2.0, while for GaAs/AlAs this ratio is equal to 1.5. This mass ratio is consistent with the acoustic impedance ratio in the case of Si/Ge (1.78) and in the case of GaAs/AlAs (1.2), as the ratio of acoustic impedances is equal to the square root of the ratio of masses (Swartz & Pohl, 1989).

In this chapter the equilibrium (EMD) and the non-equilibrium (NEMD) methods are presented and the results of the thermal conductivity are reported. For both methods the molecular dynamics code LAMMPS is used (Plimpton, 1995, 1997).

In the NEMD method or direct method, a temperature gradient is imposed across the structure under study. Depending on the boundary conditions the geometry of the thermostats can change. For periodic boundary conditions there is one thermostat (cold or hot) at the middle of the slab and a second thermostat (hot or cold) is divided in two parts, which are positioned at the two edges of the slab (fig. 1). This configuration gives the opportunity to increase the statistics (double the results). The thermal conductivity is then calculated using the Fourier's law monitoring the thermal power exchange and the temperature profile of the system (fig. 2). The infinitive system size thermal conductivity then can be extrapolated by plotting the inverse of the thermal conductivity as a function of the inverse of the system size (Schelling et al, 2002).

single mean free path, and we should have rather done the analysis mode by mode (Sellan

and the size dependence is found to be consistent with MD NEMD simulations (Schelling et al 2002, Termentzidis et al 2009). Because of the algebraic decay, size effects in NEMD are severe, and sometimes it is faster to run equilibirum simulations to determine the thermal

We now illustrate the determination of the thermal conductivity using the direct method with the exemple of superlattices composed of a regular alternative arrangement of solid layers, having different physical properties. Thermal transport in superlattices is characterized by the cross-plane and the in-plane conductivities, which correspond respectively to heat flowing in the direction perpendicular or parallel to the interfaces.

The geometry used for the determination of the cross-plane and in-plane thermal conductivities using NEMD is given in figure 3.1.1. Periodic boundary conditions are used in all directions, while in all cases the heat flux is imposed in z-direction. The interatomic interactions are described with Lennard-Jones (LJ) potentials, with energy unit ij=1.0 and length unit ij=1.0. The use of LJ potential is justified by the interest in focusing on the main phenomena introduced by the interface roughness of the super-lattices (for all the following results). The two types of materials A and B may represent respectively the Si and Ge, if one considers Si/Ge-type superlattice (paragraph 5) or the GaAs and AlAs in case of GaAs/AlAs-type superlattice (paragraph 6). The two solids have the same lattice constants and they differ only by the mass ratio of the atoms constituting the superlattices layers. For the Si/Ge-type superlattices the mass ratio is taken to be 2.0, while for GaAs/AlAs this ratio is equal to 1.5. This mass ratio is consistent with the acoustic impedance ratio in the case of Si/Ge (1.78) and in the case of GaAs/AlAs (1.2), as the ratio of acoustic impedances is equal to the square root of the ratio of masses (Swartz & Pohl,

In this chapter the equilibrium (EMD) and the non-equilibrium (NEMD) methods are presented and the results of the thermal conductivity are reported. For both methods the

In the NEMD method or direct method, a temperature gradient is imposed across the structure under study. Depending on the boundary conditions the geometry of the thermostats can change. For periodic boundary conditions there is one thermostat (cold or hot) at the middle of the slab and a second thermostat (hot or cold) is divided in two parts, which are positioned at the two edges of the slab (fig. 1). This configuration gives the opportunity to increase the statistics (double the results). The thermal conductivity is then calculated using the Fourier's law monitoring the thermal power exchange and the temperature profile of the system (fig. 2). The infinitive system size thermal conductivity then can be extrapolated by plotting the inverse of the thermal conductivity as a function of

molecular dynamics code LAMMPS is used (Plimpton, 1995, 1997).

the inverse of the system size (Schelling et al, 2002).

*L = L + cvL* 6 / (27)

 

However, the previous simple kinetic model predicts:

1/ 1/ 

et al 2010).

conductivity.

1989).

Fig. 1. Geometric configuration for the determination of the cross-plane (left) and in-plane (right) thermal conductivity simulations with NEMD and periodic boundary conditions are used in all directions.

Fig. 2. To calculate the thermal conductivity (in the example the in-plane TC) of the superlattices the energy exchange (b1) and the temperature profile (b2) is extracted for each of the structures (at least 4 structures with increasing length see (a)). Then the inverse thermal conductivity is plotted as a function of the inverse of the size of slab (c). The last diagram helps in calculating the thermal conductivity of an infinitive size superlattice.

#### **3.2 Equilibrium molecular dynamics method (EMD)**

Alternately, the thermal conductivity of a model system can be determined using equilibrium simulations (Shelling et al, 2002). The principle relies on the fact that the regression of the thermal fluctuations of an internal variable, in our case the thermal flux, obeys macroscopic laws. Hence, the time decay of the fluctuations of the flux is proportional to the thermal conductivity. This is mathematically expressed by the Green-Kubo formulae which state that the time integral of the heat flux autocorrelation function is proportional to the thermal conductivity tensor (Evans & Morris, 1990):

$$\mathcal{A}\_{\alpha,\beta} = \frac{1}{Vk\_BT^2} \int\_0^{\star x} \{f\_\alpha(t)f\_\beta(0)\}dt\tag{28}$$

Molecular Dynamics Simulations and Thermal Transport at the Nano-Scale 85

and if we assume that for sufficient low wavevectors that the medium behaves like a Debye

 *<sup>c</sup> L aL* / and thus decreases faster as the error in NEMD *δλ Λ NEMD L L* / . This allows obtaining good estimates of the thermal conductivity with rather small systems. The main drawback is intrinsic to the method, as we need to probe small thermal fluctuations around equilibrium over long time scales, which requires also performing

Apart from the thermal conductivity, both NEMD and EMD simulation techniques may be used to calculate the thermal boundary resistance characterizing heat transport across the interface between two media (Barrat, 2003). The thermal boundary resistance (TBR), also known as the Kapitza resistance is defined as the ratio of the temperature jump at the

> *K <sup>T</sup> R = J*

The Kapitza resistance is thus a measure of the ability of phonons to be transmitted by the interface, with common values falling in the range between 0,01 and <sup>1</sup> 0,1*MW mK* (Cahill et al 2003) depending on the dissimilarity between the two media. There have been several models to predict the TBR between two solids (Swartz & Pohl 1987, 1989). Among the most popular, let us mention the acoustic mismatch model (AMM) (Khalatnikov 1952) in which the transmission of phonons is assumed to be specular and depends only the acoustic impedance mismatch between the two media, in a way similar to Snell Descartes's law ruling the transmission of an electromagnetic wave across the interface between two dielectrics with different optical indexes. The other popular model is the so-called diffuse mismatch model (DMM) (Swartz & Pohl 1989), in which it is assumed that the phonons are diffusively scattered by the interface, and the transmission coefficient depends in this limit in the mismatch of the density of states characterizing the two media. Generally speaking, both models fail to predict the thermal boundary resistance of real interfaces (Cahill et al 2003), and a deep physical understanding of interfacial heat transport between two solids is still missing. In this context, MD may be a convenient tool to probe the transmission of phonons across ideal surfaces, and also has the flexibility to introduce defects at the interface and study their effect on interfacial heat transport. Again, two routes may be followed to calculate the interfacial resistance: in NEMD simulations, the interface to be charaterized is placed between two heat reservoirs fixed at a distance a few atomic layers away from the interface. A net heat flux in the direction perpendicular to the interface is created, and the resulting temperature jump at the interface is measured, allowing to derive the interfacial resistance. With NEMD method only the TBR of smooth interfaces can be

( )

 

*<sup>c</sup> L*

 

several statistical averages over different initial conditions.

interface *T* over the heat flux *J* crossing the interface:

solid, <sup>2</sup> *g* 

 

to obey Callaway model: <sup>2</sup>

**3.3 Thermal boundary resistance** 

0

*gcv d*

2

 (32)

, the group velocity is constant and the phonon relaxation time is supposed

, then the error in the conductivity scales as

(33)

where *V* is the volume of the system, and *J* denotes the component of the heat flux vector along the direction . The equilibrium method consists then in computing the corresponding autocorrelation, which requires following the dynamics of a system over time scales a few times larger than the longest relaxation time present in the system. In the case of heat transfer in solid materials, the longest relaxation times correspond to long wavelengths phonons (a few nm) with a lifetime on the order of 100 picoseconds.

The advantage of the equilibrium method is that it allows to compute the full conductivity tensor from one simulation, which may be appreciated in superlattice simulations for instance, which display large thermal anisotropies. Another advantage of the equilibrium method is that it does not suffer from severe finite size effects as NEMD, mainly because in EMD the phonons are not strongly scattered by the boundaries of the simulation box. We can estimate the finite size effects in EMD using the following assumptions: the thermal conductivity measured in a MD simulation is given by:

$$\lambda = \int\_{a\_{\mathbb{L}}(L)}^{a\_{\mathbb{D}}} g(\boldsymbol{\alpha}) \boldsymbol{c}(\boldsymbol{\alpha}) \boldsymbol{v}^2(\boldsymbol{\alpha}) \boldsymbol{\tau}(\boldsymbol{\alpha}) d\boldsymbol{\alpha} \tag{29}$$

where

$$
\alpha \rho\_c \left( L \right) = \alpha \left( k \equiv \pi \;/ \; L \right) \tag{30}
$$

is the low pulsation cut-off introduced by the periodic boundary conditions, and the upper bound is

$$
\alpha o\_{\rm D} = o(k = \pi \;/\; a) \tag{31}
$$

where *a* is the interatomic step, since high frequency phonons are supposed to contribute weakly to the overall thermal conductivity. In the latter expression of the thermal conductivity, *<sup>B</sup> c =k* is the heat capacity which obeys Dulong and Petit law, *v* is the phonon group velocity, and is the relaxation time of phonon having a pulsation . In the following, we have ignored the possible different polarization states of the phonons to simplify the discussion. The error in the determination of the thermal conductivity is thus:

$$\delta \mathcal{S} = \int\_0^{\alpha\_\varepsilon(L)} g(\alpha) c(\alpha) v^2(\alpha) \tau(\alpha) d\alpha \tag{32}$$

and if we assume that for sufficient low wavevectors that the medium behaves like a Debye solid, <sup>2</sup> *g* , the group velocity is constant and the phonon relaxation time is supposed to obey Callaway model: <sup>2</sup> , then the error in the conductivity scales as *<sup>c</sup> L aL* / and thus decreases faster as the error in NEMD *δλ Λ NEMD L L* / .

This allows obtaining good estimates of the thermal conductivity with rather small systems. The main drawback is intrinsic to the method, as we need to probe small thermal fluctuations around equilibrium over long time scales, which requires also performing several statistical averages over different initial conditions.

#### **3.3 Thermal boundary resistance**

84 Molecular Dynamics – Theoretical Developments and Applications in Nanotechnology and Energy

Alternately, the thermal conductivity of a model system can be determined using equilibrium simulations (Shelling et al, 2002). The principle relies on the fact that the regression of the thermal fluctuations of an internal variable, in our case the thermal flux, obeys macroscopic laws. Hence, the time decay of the fluctuations of the flux is proportional to the thermal conductivity. This is mathematically expressed by the Green-Kubo formulae which state that the time integral of the heat flux autocorrelation function is proportional to

0

corresponding autocorrelation, which requires following the dynamics of a system over time scales a few times larger than the longest relaxation time present in the system. In the case of heat transfer in solid materials, the longest relaxation times correspond to long

The advantage of the equilibrium method is that it allows to compute the full conductivity tensor from one simulation, which may be appreciated in superlattice simulations for instance, which display large thermal anisotropies. Another advantage of the equilibrium method is that it does not suffer from severe finite size effects as NEMD, mainly because in EMD the phonons are not strongly scattered by the boundaries of the simulation box. We can estimate the finite size effects in EMD using the following assumptions: the thermal

<sup>2</sup>

*g cv d*

 (29)

<sup>+</sup>

 

. The equilibrium method consists then in computing the

(28)

denotes the component of the heat flux vector

*<sup>c</sup> L = k= L* / (30)

*<sup>D</sup> k a* / (31)

is the relaxation time of phonon having a pulsation

is the

> . In

<sup>1</sup> <sup>0</sup>

*= J t J dt*

**3.2 Equilibrium molecular dynamics method (EMD)** 

the thermal conductivity tensor (Evans & Morris, 1990):

where *V* is the volume of the system, and *J*

conductivity measured in a MD simulation is given by:

along the direction

where

bound is

is thus:

conductivity, *<sup>B</sup> c =k*

phonon group velocity, and

 

, 2

*B*

wavelengths phonons (a few nm) with a lifetime on the order of 100 picoseconds.

( )

 

 

is the low pulsation cut-off introduced by the periodic boundary conditions, and the upper

where *a* is the interatomic step, since high frequency phonons are supposed to contribute weakly to the overall thermal conductivity. In the latter expression of the thermal

the following, we have ignored the possible different polarization states of the phonons to simplify the discussion. The error in the determination of the thermal conductivity

 

is the heat capacity which obeys Dulong and Petit law, *v*

*<sup>c</sup> L*

  *D*

*Vk T*

Apart from the thermal conductivity, both NEMD and EMD simulation techniques may be used to calculate the thermal boundary resistance characterizing heat transport across the interface between two media (Barrat, 2003). The thermal boundary resistance (TBR), also known as the Kapitza resistance is defined as the ratio of the temperature jump at the interface *T* over the heat flux *J* crossing the interface:

$$R\_K = \frac{\Delta T}{J} \tag{33}$$

The Kapitza resistance is thus a measure of the ability of phonons to be transmitted by the interface, with common values falling in the range between 0,01 and <sup>1</sup> 0,1*MW mK* (Cahill et al 2003) depending on the dissimilarity between the two media. There have been several models to predict the TBR between two solids (Swartz & Pohl 1987, 1989). Among the most popular, let us mention the acoustic mismatch model (AMM) (Khalatnikov 1952) in which the transmission of phonons is assumed to be specular and depends only the acoustic impedance mismatch between the two media, in a way similar to Snell Descartes's law ruling the transmission of an electromagnetic wave across the interface between two dielectrics with different optical indexes. The other popular model is the so-called diffuse mismatch model (DMM) (Swartz & Pohl 1989), in which it is assumed that the phonons are diffusively scattered by the interface, and the transmission coefficient depends in this limit in the mismatch of the density of states characterizing the two media. Generally speaking, both models fail to predict the thermal boundary resistance of real interfaces (Cahill et al 2003), and a deep physical understanding of interfacial heat transport between two solids is still missing. In this context, MD may be a convenient tool to probe the transmission of phonons across ideal surfaces, and also has the flexibility to introduce defects at the interface and study their effect on interfacial heat transport. Again, two routes may be followed to calculate the interfacial resistance: in NEMD simulations, the interface to be charaterized is placed between two heat reservoirs fixed at a distance a few atomic layers away from the interface. A net heat flux in the direction perpendicular to the interface is created, and the resulting temperature jump at the interface is measured, allowing to derive the interfacial resistance. With NEMD method only the TBR of smooth interfaces can be calculated, as for rough interfaces there are geometrical issues about the crosssection/effective surface of the interfaces.

Equilibrium simulations rely again on the Green-Kubo formula for the interfacial conductance G=1/R (Puech 1986):

$$G = \frac{1}{Ak\_B T^2} \int\_0^{+\infty} \langle q(t)q(0) \rangle dt \tag{34}$$

Molecular Dynamics Simulations and Thermal Transport at the Nano-Scale 87

dimensional confinement of the phonon modes. This effect shares some similarities with the electron confinement in a quantum well (Stroscio & Dutta 2001). Phonon engineering in nanostructures can be succeeded in tailoring the phonon modes through the designing of

Fig. 3. Several types of nanostructured materials, including nano-wires, nano-objects,

regime offers the possibility to tailor the thermal properties of nanostructures.

Nanoscale heat transfer has indeed an old history, which can be traced back to the Fermi-Pasta-Ulam study of heat transport in a 1D anharmonic chain (Dhar 2008). It is theoretically predicted that the thermal conductivity of model 1D chain assumptotically increases with

conserving process (Dhar 2008). The possibility of a diverging conductivity together with thte fantastic developement of nanotechnologies in the nineties motivated experimental investigation of energy transfer in 1D systems. There has been however considerably fewer

, at least as soon as phonon scattering is ensured by momentum

Nanomaterials are composite materials with at least one characteristic dimension smaller than 0.5µm. Super-lattices, nanofilms, nanowires, nanotubes and nanoparticles are typical examples of such materials (fig. 3). Nanomaterials exhibit transport properties different from their constituents and a fundamental understanding of heat conduction at the nanoscale is absolutely necessary to tailor their properties. Heat conduction on nanoscale differs significantly from heat transfer laws describing macroscopic transport. Among others, nanoscale energy transport can be controlled by phonon confinement (Montroll, 1950), the modification of the density of states induced by the confinement and also the ballistic behaviour of phonons. At low temperatures, phonons may travel ballistically in the bulk of the nanostructured material being scattered only by the interfaces, the latter mechanism providing the scattering mechanism which controls the thermal conductivity of the material. The crucial parameter describing ballistic transport is the ratio of the phonon mean free path (MFP) over the characteristic lengths of the nanostructure (eg. for 3Dstructures: period and roughness of the interfaces of superlattices, for 2D: interfacial roughness in case of thin films or 1D: length of nanowires). Three different regimes may be distinguished, depending on the ratio of the characteristic length of the nanostructures, over the mean free path of the energy carriers. When the characteristic length is larger than the MFP of carriers the transport is diffusive and can be described by Fourier law, while when it is smaller, the heat carriers feel the nanostructured material as a homogeneous medium with a low thermal conductivity. There is an intermediate regime, where the transport becomes ballistic and diffusion scattering becomes predominant. As we will see, this third

the dimensions and the roughness of nanostructured materials.

superlattices etc.

the chain length L *L*

 

where A is the interfacial area and q(t) is the instantaneous value of the flux flowing across the interface.

The latter quantity may be computed in a MD simulation using the power of interfacial forces (Barrat 2003):

$$\eta(t) = \sum\_{i \in \mathbf{1}, j \in \mathbf{2}} \vec{F}\_{ij} \cdot \vec{\boldsymbol{\upsilon}}\_{i} \tag{35}$$

where the indexes 1 and 2 denote the two media separating the interface to be studied, and we have assumed pair potentials to simplify. Equilibrium simulations have been performed for interface between simple Lennard Jones solids, where the dissimilarity between the two solids is introduced by changing the acoustic impedance ratio between the two solids. The results were found to disagree with NEMD determinations of the Kapitza conductance (McGaughey 2006), especially for solids with a weak acoustic mismatch. These discrepancies are not surprising and can be traced back to the formulation of the equations themselves (Pettersson 1990), as the Green-Kubo formula above predicts a finite conductance for the interface between two identical media, while of course in NEMD one measures an infinite conductance in this situation. Maybe for this reason, most of the MD works on Kapitza resistance have considered the direct method.

Among significant work, let us mention Stevens et al. 2007, who showed that the DMM and AMM models underpredicts the thermal boundary conductance obtained by NEMD, which is found to increase with the temperature, in contrast with the theoretical models predicting a constant value, at least above the Debye temperature of the softer solid. The authors showed also that the existence of a large lattice mismatch induces interfacial stress, which deteriorates thermal transport significantly. The existence of grain boundaries has been shown also to increase the interfacial resistance (Schelling 2004).

Further work is needed, in particular using equilibrium molecular dynamics to compare MD with the available theoretical models.

#### **4. Heat transport in nanostructured materials**

Thermal transport in nanostructured materials has attracted an increasing international interest in the last decades. From a theoretical point of view, nanostructured materials are platforms for testing novel phonon and electron transport theories. From the research and development point of view, nanomaterials are promising candidates for nanoscale on chip coolers (Zhang & Li, 2010). A high density of interfaces, which is the source of phonon scattering, appeared in advanced technological devices affecting their reliability and performance. Phonon interactions are modified significantly in nanostructures due to the

calculated, as for rough interfaces there are geometrical issues about the cross-

Equilibrium simulations rely again on the Green-Kubo formula for the interfacial

*G q t q dt*

where A is the interfacial area and q(t) is the instantaneous value of the flux flowing across

The latter quantity may be computed in a MD simulation using the power of interfacial

where the indexes 1 and 2 denote the two media separating the interface to be studied, and we have assumed pair potentials to simplify. Equilibrium simulations have been performed for interface between simple Lennard Jones solids, where the dissimilarity between the two solids is introduced by changing the acoustic impedance ratio between the two solids. The results were found to disagree with NEMD determinations of the Kapitza conductance (McGaughey 2006), especially for solids with a weak acoustic mismatch. These discrepancies are not surprising and can be traced back to the formulation of the equations themselves (Pettersson 1990), as the Green-Kubo formula above predicts a finite conductance for the interface between two identical media, while of course in NEMD one measures an infinite conductance in this situation. Maybe for this reason, most of the MD works on Kapitza

Among significant work, let us mention Stevens et al. 2007, who showed that the DMM and AMM models underpredicts the thermal boundary conductance obtained by NEMD, which is found to increase with the temperature, in contrast with the theoretical models predicting a constant value, at least above the Debye temperature of the softer solid. The authors showed also that the existence of a large lattice mismatch induces interfacial stress, which deteriorates thermal transport significantly. The existence of grain boundaries has been

Further work is needed, in particular using equilibrium molecular dynamics to compare

Thermal transport in nanostructured materials has attracted an increasing international interest in the last decades. From a theoretical point of view, nanostructured materials are platforms for testing novel phonon and electron transport theories. From the research and development point of view, nanomaterials are promising candidates for nanoscale on chip coolers (Zhang & Li, 2010). A high density of interfaces, which is the source of phonon scattering, appeared in advanced technological devices affecting their reliability and performance. Phonon interactions are modified significantly in nanostructures due to the

*B*

*Ak T*

 <sup>2</sup> 0 <sup>1</sup> <sup>0</sup>

(34)

1, 2 *ij <sup>i</sup> i j <sup>q</sup> t= F v* (35)

section/effective surface of the interfaces.

resistance have considered the direct method.

MD with the available theoretical models.

**4. Heat transport in nanostructured materials** 

shown also to increase the interfacial resistance (Schelling 2004).

conductance G=1/R (Puech 1986):

the interface.

forces (Barrat 2003):

dimensional confinement of the phonon modes. This effect shares some similarities with the electron confinement in a quantum well (Stroscio & Dutta 2001). Phonon engineering in nanostructures can be succeeded in tailoring the phonon modes through the designing of the dimensions and the roughness of nanostructured materials.

Fig. 3. Several types of nanostructured materials, including nano-wires, nano-objects, superlattices etc.

Nanomaterials are composite materials with at least one characteristic dimension smaller than 0.5µm. Super-lattices, nanofilms, nanowires, nanotubes and nanoparticles are typical examples of such materials (fig. 3). Nanomaterials exhibit transport properties different from their constituents and a fundamental understanding of heat conduction at the nanoscale is absolutely necessary to tailor their properties. Heat conduction on nanoscale differs significantly from heat transfer laws describing macroscopic transport. Among others, nanoscale energy transport can be controlled by phonon confinement (Montroll, 1950), the modification of the density of states induced by the confinement and also the ballistic behaviour of phonons. At low temperatures, phonons may travel ballistically in the bulk of the nanostructured material being scattered only by the interfaces, the latter mechanism providing the scattering mechanism which controls the thermal conductivity of the material. The crucial parameter describing ballistic transport is the ratio of the phonon mean free path (MFP) over the characteristic lengths of the nanostructure (eg. for 3Dstructures: period and roughness of the interfaces of superlattices, for 2D: interfacial roughness in case of thin films or 1D: length of nanowires). Three different regimes may be distinguished, depending on the ratio of the characteristic length of the nanostructures, over the mean free path of the energy carriers. When the characteristic length is larger than the MFP of carriers the transport is diffusive and can be described by Fourier law, while when it is smaller, the heat carriers feel the nanostructured material as a homogeneous medium with a low thermal conductivity. There is an intermediate regime, where the transport becomes ballistic and diffusion scattering becomes predominant. As we will see, this third regime offers the possibility to tailor the thermal properties of nanostructures.

Nanoscale heat transfer has indeed an old history, which can be traced back to the Fermi-Pasta-Ulam study of heat transport in a 1D anharmonic chain (Dhar 2008). It is theoretically predicted that the thermal conductivity of model 1D chain assumptotically increases with the chain length L *L* , at least as soon as phonon scattering is ensured by momentum conserving process (Dhar 2008). The possibility of a diverging conductivity together with thte fantastic developement of nanotechnologies in the nineties motivated experimental investigation of energy transfer in 1D systems. There has been however considerably fewer thermal measurement in low dimensional systems compared to electrical measurements for instance, because of the difficulty to measure directly the thermal current. Carbon nanotubes is a good paradigm of 1D objects, because all the heat carriers travel in the axial direction. The first experiments probing energy transfer in isolated multi-wall carbon nanotubes (Kim 2001) reported values of the thermal conductivity 3000 / *W mK* even larger than the conductivity of diamond at room temperature. These large values have been confirmed by other experiments on isolated single wall (Yu 2005, Pop 2006).

This experimental effort has been accompanied by several molecular dynamics works (Berber 2000, Maruyama 2002, Zhang 2005, Yao 2005, Donadio 2007), which have concluded that the phonon mean free path may be m and that the thermal conductivity converges for very long nanotubes (Mingo 2005, Donadio 2007). We note besides that Green-Kubo equilibrium simulations are preferred because in the direct method, it is hard to compute an effective conductivity in a material where heat is transported ballistically. Also the coupling between the nanotube and heat reservoirs in NEMD would certainly affect the conductivity measurement.

Thermal transport in nanowires has been also extensively studied during the last decade. The motivations of the first experimental works was to measure the quantum of conductance (Schwab 2000):

$$\mathbf{g}\_0 = \pi^2 k\_B^2 T \text{ / (3h)}\tag{36}$$

Molecular Dynamics Simulations and Thermal Transport at the Nano-Scale 89

field, with the exception of the modeling of self-assembled monolayers (Luo 2010), Henry and Chen used equilibrium simulations to show that the conductivity of polyethylene chains may be orders of magnitude larger than bulk polyethylene (Henry 2008). A recent study concluded that the conductance of molecular chain is also strongly affected by its

Finally, heat transfer in the vicinity of nano-particles has also shown an increasing interest during the last years. This interest has been motivated first by the early measurements of large thermal conductivity in the so-called nanofluids, i.e. suspensions of nano-particles in a liquid solvent (Keblinski 2008). Molecular simulations have helped in solving the controversy and showed that for well dispersed nano-particles, no enhancement is expected with respect to the effective medium theory (Vladkov 2006a). These simulations have also helped in determining the interfacial thermal resistance characterising the nanoparticle/liquid interface (Vladkov 2006a). Heat transfer around strongly heated nanoparticles has also attracted attention after the development of ultrafast optical techniques which allow to selectively heat nanoparticles in suspension (Plech 2004). When the metallic particles are excited at wavelengths close to the maximum of their optical extinction, their temperature can be raised by hundreds of Kelvin, while the liquid environment in the immediate vicinity may remain at ambient temperature, thus creating very large temperature gradients and energy fluxes flowing from the nanoparticles. This raises interesting new questions regarding nanoscale heat transfer, e.g. regarding the validity of Fourier law at very large temperature gradients (Merabia 2009b), and the competition between heat transfer and boiling of the fluid surrounding the nanoparticles (Merabia 2009b). Apart from the academic interest, these questions have important biomedical applications in hyperthermia, as the appearance of vapor bubbles with submicronic radius would concentrate large thermomechanical stresses which may destroy tumors for instance. Further work is needed to understand the conditions of formation of these "nanobubbles". Let us conclude saying that nanoparticles may also served as thermal contacts to measure the conductance of molecular chains as demonstrated using MD (Merabia2011).

environment (Merabia 2011), with a transition between ballistic and Fourier regime.

**5. Thermal conductivity predictions for Si/Ge superlattices, impact of** 

Heat transport in superlattices, which are materials composed of a periodic or a random arrangement of different alternating materials with a submicronic thickness, has attracted a large scientific interest, as they exhibit low thermal conductivity (Cahill et al, 2003), at least in one direction, usually the direction normal to the interfaces. This makes superlattices promising materials for applications in MEMS and NEMS devices such as semiconductor lasers (Sale, 1995), optical data-storage media (Kim et al, 2000), thermoelectric (Hicks et al, 1993; Lin-Chung & Reinecke, 1993) and thermomechanic devices (Ezzahri et al, 2008). For the latter two categories the thermal-conductivity characteristics are very important to ensure the correct function of the device(Daly et al, 2002). Depending on the use, the highest possible thermal conductivity is required for example to remove the Joule heat in electronic devices, or very low thermal conductivity for thermoelectric applications (Mahan 2004). With a clever tailoring of the properties of superlattice, one can succeed in both directions. The phonon thermal conductivity of superlattices has been started to be in the center of scientific interest quite early as superlattice are promising structures for new electronic devices (Ren & Dow 1982). It has been reported that the thermal conductivity of

**rough interfaces** 

at very low temperatures. During the years 2000, the thermal conductivity of Si nanowires at room temperatures has been measured (Li 2003) reporting values two orders of magnitude smaller than bulk Si. These low values of the thermal conductivity may be understood because in a nanowire the free surfaces induce diffuse scattering of the phonons, contrary to carbon nanotubes where phonons can only travel in the axial direction. As a consequence, the transport properties of semi-conductor nanowires depend on the state of the nanowire surface, and in particular its roughness, opening the way to achieve materials combining low thermal conductivity but electrical transport properties comparable to bulk semiconductor, with promising applications in thermoelectric conversion (Hochbaum 2008).

On the computational side, the first MD simulation of heat transport in a nanowire has been performed by Volz and Chen (Volz 1999), who already measured a two orders of magnitude reduction of the conductivity compared to the bulk. Recent simulations have confirmed the reduction, but yielding contradicting results for very small nanowires diameters (Ponomareva 2007, Donadio 2009, Abs-Da-Cruz 2011). Molecular simulations have been also pointed out the role of surface disorder on the conductivity reduction (Donadio 2009). Comparatively, there has been fewer studies on heat transport in molecular junctions, probably because of the difficulty to measure a thermal current flowing across a molecule. Molecular junctions have been however proposed to be good candidates for thermal rectifiers (Chang 2006, Casati 2007). Whang et al. recently used ultrafast thermal to probe ballistic heat transport in alkane-thiol chains supported on a gold substrate (Whang 2007). This opens the way to measure the conductance of a molecular chain as a function of its length. Theoretical studies have first considered simple one-dimensional chains (Dhar 2008, Segal 2003). Realistic models of molecular junctions have been studied recently (Mingo 2006) using Green function technique. There is relatively few molecular dynamics studies in the

thermal measurement in low dimensional systems compared to electrical measurements for instance, because of the difficulty to measure directly the thermal current. Carbon nanotubes is a good paradigm of 1D objects, because all the heat carriers travel in the axial direction. The first experiments probing energy transfer in isolated multi-wall carbon nanotubes (Kim

conductivity of diamond at room temperature. These large values have been confirmed by

This experimental effort has been accompanied by several molecular dynamics works (Berber 2000, Maruyama 2002, Zhang 2005, Yao 2005, Donadio 2007), which have concluded that the phonon mean free path may be m and that the thermal conductivity converges for very long nanotubes (Mingo 2005, Donadio 2007). We note besides that Green-Kubo equilibrium simulations are preferred because in the direct method, it is hard to compute an effective conductivity in a material where heat is transported ballistically. Also the coupling between the nanotube and heat reservoirs in NEMD would certainly affect the conductivity

Thermal transport in nanowires has been also extensively studied during the last decade. The motivations of the first experimental works was to measure the quantum of

at very low temperatures. During the years 2000, the thermal conductivity of Si nanowires at room temperatures has been measured (Li 2003) reporting values two orders of magnitude smaller than bulk Si. These low values of the thermal conductivity may be understood because in a nanowire the free surfaces induce diffuse scattering of the phonons, contrary to carbon nanotubes where phonons can only travel in the axial direction. As a consequence, the transport properties of semi-conductor nanowires depend on the state of the nanowire surface, and in particular its roughness, opening the way to achieve materials combining low thermal conductivity but electrical transport properties comparable to bulk semiconductor, with promising applications in thermoelectric conversion (Hochbaum 2008).

On the computational side, the first MD simulation of heat transport in a nanowire has been performed by Volz and Chen (Volz 1999), who already measured a two orders of magnitude reduction of the conductivity compared to the bulk. Recent simulations have confirmed the reduction, but yielding contradicting results for very small nanowires diameters (Ponomareva 2007, Donadio 2009, Abs-Da-Cruz 2011). Molecular simulations have been also pointed out the role of surface disorder on the conductivity reduction (Donadio 2009). Comparatively, there has been fewer studies on heat transport in molecular junctions, probably because of the difficulty to measure a thermal current flowing across a molecule. Molecular junctions have been however proposed to be good candidates for thermal rectifiers (Chang 2006, Casati 2007). Whang et al. recently used ultrafast thermal to probe ballistic heat transport in alkane-thiol chains supported on a gold substrate (Whang 2007). This opens the way to measure the conductance of a molecular chain as a function of its length. Theoretical studies have first considered simple one-dimensional chains (Dhar 2008, Segal 2003). Realistic models of molecular junctions have been studied recently (Mingo 2006) using Green function technique. There is relatively few molecular dynamics studies in the

 <sup>2</sup> <sup>0</sup> / 3h *<sup>2</sup> <sup>B</sup> g = kT* 

3000 / *W mK* even larger than the

(36)

2001) reported values of the thermal conductivity

measurement.

conductance (Schwab 2000):

other experiments on isolated single wall (Yu 2005, Pop 2006).

field, with the exception of the modeling of self-assembled monolayers (Luo 2010), Henry and Chen used equilibrium simulations to show that the conductivity of polyethylene chains may be orders of magnitude larger than bulk polyethylene (Henry 2008). A recent study concluded that the conductance of molecular chain is also strongly affected by its environment (Merabia 2011), with a transition between ballistic and Fourier regime.

Finally, heat transfer in the vicinity of nano-particles has also shown an increasing interest during the last years. This interest has been motivated first by the early measurements of large thermal conductivity in the so-called nanofluids, i.e. suspensions of nano-particles in a liquid solvent (Keblinski 2008). Molecular simulations have helped in solving the controversy and showed that for well dispersed nano-particles, no enhancement is expected with respect to the effective medium theory (Vladkov 2006a). These simulations have also helped in determining the interfacial thermal resistance characterising the nanoparticle/liquid interface (Vladkov 2006a). Heat transfer around strongly heated nanoparticles has also attracted attention after the development of ultrafast optical techniques which allow to selectively heat nanoparticles in suspension (Plech 2004). When the metallic particles are excited at wavelengths close to the maximum of their optical extinction, their temperature can be raised by hundreds of Kelvin, while the liquid environment in the immediate vicinity may remain at ambient temperature, thus creating very large temperature gradients and energy fluxes flowing from the nanoparticles. This raises interesting new questions regarding nanoscale heat transfer, e.g. regarding the validity of Fourier law at very large temperature gradients (Merabia 2009b), and the competition between heat transfer and boiling of the fluid surrounding the nanoparticles (Merabia 2009b). Apart from the academic interest, these questions have important biomedical applications in hyperthermia, as the appearance of vapor bubbles with submicronic radius would concentrate large thermomechanical stresses which may destroy tumors for instance. Further work is needed to understand the conditions of formation of these "nanobubbles". Let us conclude saying that nanoparticles may also served as thermal contacts to measure the conductance of molecular chains as demonstrated using MD (Merabia2011).

### **5. Thermal conductivity predictions for Si/Ge superlattices, impact of rough interfaces**

Heat transport in superlattices, which are materials composed of a periodic or a random arrangement of different alternating materials with a submicronic thickness, has attracted a large scientific interest, as they exhibit low thermal conductivity (Cahill et al, 2003), at least in one direction, usually the direction normal to the interfaces. This makes superlattices promising materials for applications in MEMS and NEMS devices such as semiconductor lasers (Sale, 1995), optical data-storage media (Kim et al, 2000), thermoelectric (Hicks et al, 1993; Lin-Chung & Reinecke, 1993) and thermomechanic devices (Ezzahri et al, 2008). For the latter two categories the thermal-conductivity characteristics are very important to ensure the correct function of the device(Daly et al, 2002). Depending on the use, the highest possible thermal conductivity is required for example to remove the Joule heat in electronic devices, or very low thermal conductivity for thermoelectric applications (Mahan 2004). With a clever tailoring of the properties of superlattice, one can succeed in both directions.

The phonon thermal conductivity of superlattices has been started to be in the center of scientific interest quite early as superlattice are promising structures for new electronic devices (Ren & Dow 1982). It has been reported that the thermal conductivity of

Molecular Dynamics Simulations and Thermal Transport at the Nano-Scale 91

In this subsection, we study the effect of the superlattice period and the structure of the interface on both the in-plane and the cross-plane superlattice thermal conductivities. Simulations have been held for two periods of superlattices 20 and 400, which are comparable to the phonon meat free path at the temperature we are working at. To understand the role of the interfacial structure, we have considered two types of interfaces: smooth interfaces on the one hand, and rough interfaces with a right-isoscele-triangles shape, as shown in fig. 4. In the case of rough interfaces, the height of the interfaces has been varied between 1 monolayer (ML) which is half of a lattice constant for a fcc crystal, up to 120, which is more than the half of the superlattice period for superlattices with period of 200 and from 1ML up to 240 for superlattices with period of 400 (see fig. 4). The maximal

Fig. 4. Several triangular shaped interfaces with varing height, for superlattice period

**5.2 Thermal conductivity of Si/Ge like superlattices with rough interfaces** 

We examined also the effect of the shape of the interface on the cross-plane and the in-plane thermal conductivities. In this case, the superlattice period is kept constant and equal to 20<sup>0</sup> and we considered only one height of interfaces, the 60 or 12MLs. Fig. 5 shows the smooth interfaces, the periodic triangular isosceles interfaces and the additional 4 other shapes. The structure shown in 5.iii is obtained by superposing the periodic isosceles of small lengths with the periodic triangular interfaces shown in 5.ii. Cosine like, random like and square

The in-plane and cross-plane thermal conductivities have been calculated using the EMD and the NEMD methods. The results are displayed in figure 6 as a function of the interface

of 400.

like interfaces are also examined.

**5.1 Modelling Si/Ge superlattices with rough interfaces** 

interface height is thus more than half the superlattice's period.

superlattices may be dramatically smaller than the corresponding values of the constituent materials in their bulk form. This decrease has been related to the folding of the Brillouin zone and the related mini-umklapp three-phonon scattering process. Tamura & al 1999 analysed the effect on the phonon spectra in superlattices by three major reasons: a. folding of the phonon branches caused by the periodicity of the superlattices, b. the formation of the mini-band and c. the confinement of the acoustic phonons in the different layers due to the mismatch of the spectra. These three reasons impose reduction of the group velocity in the cross-plane direction, leading to the decrease of the cross-plane thermal conductivity. Chen & Yang 2005a claimed that the group velocity reduction is not sufficient to explain the dramatic decrease of the thermal conductivity, and they argued that one should add the diffusive scattering at the interfaces and treat phonons as incoherent particles. Chen & Neagu 1997 solving the Boltzmann Transport Equation for specular and diffuse interfaces showed that depending on the superlattice period, the thermal conductivity might be influenced either by the diffuse interface scattering or by the scattering induced by the dislocations. The literature is rich in this subject and a series of articles appeared with a lot of experimental (Capinski et al 1999, Huxtable et al 2002, Lee et al 1997) or theoretical results using lattice dynamics method or Equilibrium (Volz 2000, Landry 2008, Termentzidis 2011b, Termentzidis 2011c) and Non-Equilibrium Molecular Dynamics method (Liang & Shi 2000, Chen 2004, Termentzidis 2009, Termentzidis 2010).

One very interesting property of superlattices is their thermal anisotropy. In the paragraphs below the distinction between the in-plane (parallel to the interfaces) and cross-plane (perpendicular to the interfaces) thermal conductivity has been underlined. In general one expects a value of the in-plane thermal conductivity close to the bulk conductivity especially for superlattices with smooth interfaces, where phonons are expected to specularly reflect from the interfaces, and in which each layer behaves like a phonon waveguide. For the cross-plane direction the picture is totally different, with thermal conductivity even smaller than a random alloy of the same material (Mahan 2004). A key point for the physical explanation of the phenomena related to the superlattice thermal conductivity is the thermal boundary resistance or the Kapitza resistance, which has been discussed before.

Molecular dynamics simulations have been performed recently to understand the physical mechanisms ruling the transport properties of superlattices (Landry et al, 2008, 2009, Termentzidis et al, 2009,2010,2011a,2011b,2011c). In this contribution, we will discuss the influence of the interface roughness and of the superlattice period on the in-plane and the cross-plane thermal conductivities of Si/Ge superlattices using both the EMD and NEMD methods. This study proves that heat transport in superlattices is controlled by the interfaces. An atomic knowledge (or description) of the interfaces is necessary for the correct prediction of the thermal conductivity. In turn, understanding the link between the interfacial structure and the thermal conductivity will certainly help in tailoring and controlling the phonon behaviour in nanostructures. This can lead to the augmentation of the lifetime and to optimize the working of several nano-devices. The state of the interfaces is crucial for the determination of the behavior of phonons within the nanostructures. When the layer thickness of the superlattice is comparable with the MFP, the thermal conductivity is controlled by the transmission of phonons across the interfaces of the superlattice. In particular, the thermal boundary resistance or the Kapitza resistance, which has been discussed before, will play a key role in the thermal transport in superlattices with thin layers.
