Molecular Dynamics Simulations and Thermal Transport at the Nano-Scale

Molecular Dynamics is a two-volume compendium of the ever-growing applications of molecular dynamics simulations to solve a wider range of scientific and engineering challenges. The contents illustrate the rapid progress on molecular dynamics simulations in many fields of science and technology, such as nanotechnology, energy research, and biology, due to the advances of new dynamics theories and the extraordinary power of today's computers. This first book begins with a general description of underlying theories of molecular dynamics simulations and provides extensive coverage of molecular dynamics simulations in nanotechnology and energy. Coverage of this book includes: Recent advances of molecular dynamics theory Formation and evolution of nanoparticles of up to 106 atoms Diffusion and dissociation of gas and liquid molecules on silicon, metal, or metal organic frameworks Conductivity of ionic species in solid oxides Ion solvation in liquid mixtures Nuclear structures


Introduction
This chapter presents an overview of the Molecular Dynamics (MD) simulation technique to predict thermal transport properties of nanostructured materials. This covers systems having characteristic lengths of the order of a few nanometers like carbon nanotubes, nanowires and also superlattices, i.e. composite materials made of submicronic thickness of solid layers. The common features of these systems is the small ratio between their characteristic system size and the phonon mean free path, which leads to ballistic heat transport and deviations from the classical Fourier law. Also when the density of interfaces gets large, the energy transport properties of the materials can not longer be described solely by the thermal conductivities of the constituents of the material, but depend also on the thermal boundary resistance which measures the transmission of phonons across an interface. In this context, molecular dynamics was proven to be a very useful technique to study heat transport in nanostructured materials. The main reasons are; the length scale probed by the method is in the nanometer range, and it does not make any assumption on the phonons dynamics except their classical nature.
In this contribution, we present two MD methods, the equilibrium and the non-equilibrium method, which are now commonly used to determine both the thermal conductivity and the thermal boundary resistance of nanostructured materials. We focus on superlattices and discuss how the structural features of the interfaces like height, shape, inter-diffusion phenomena and the layer thickness affect the thermal conductivity of the superlattice. We show how these complex phenomena can be predicted by simple models of Lennard-Jones crystals with a mass ratio corresponding to the acoustic impedance ratio of Si/Ge and GaAs/AlAs superlattices.

Molecular dynamics
The development of molecular simulations began in the early fifties after the considerable development of computer facilities in the United States during World War II. A few years after the first Monte-Carlo simulation, Alder and Wainwright, first introduced in the late 50's the Molecular Dynamics (MD) method (Alder & Wainwright, 1957, 1959. The aim of the first simulations both Monte-Carlo or MD was to probe the different phases of model hard spheres. The necessity to model liquids motivated the development of realistic potentials. Rahman was the first to model Argon using Lennard-Jones potential , which is still considered as a standard potential for MD (Rahman, 1964). This opened the way to consider a broad range of condensed matter systems, ranging from liquid water first modelled by Rahman and Stillinger to silicon (Rahman & Stillinger, 1974). The efforts towards realistic modelling have permitted to applicate MD to characterize the collective excitations in solids (Hensen 1976). The last step came with the implementation of MD thermostats opening the way to probe the phonon dynamics in rare-gas solids (Ladd 1986). More recently, there has been an increasing number of MD studies on nano-scale heat transport motivated by fine measurements of energy transport in nano-materials (Volz 1999). We refer the reader to the classical textbooks on MD (Allain &Tildesley 1987 andFrenkel &Smit 1996) for a thorough introduction.
MD is a simulation method based on the numerical integration of Newton's equation of motion: where: m , i r  and V denote respectively the mass of the particle i, its position vector and the inter-atomic potential. The integration scheme commonly used is the Verlet algorithm (Verlet 1967a(Verlet , 1967b  (3) and the error is larger than the error on the positions, but as soon as the velocities are only used to compute the instantaneous kinetic energy the consequences are minor. Of course there are other integration schemes, but the Verlet algorithm has the advantage to be simple, easy to implement, stable and time-reversible. The typical time step used in the integration algorithms is a fraction of the characteristic atomic period where m is the mass of the particles and  denotes a typical atomic diameter.
Given the power of modern computers, it takes a few hours on a mono-processor machine to follow the trajectory of a set of 10000 particles over a time comparable to 1 ns. Using parallel MD codes or GPUs (graphics processing units) opens the way to model larger www.intechopen.com systems or reach orders of magnitude longer times . The system represented in a MD simulation has microscopic dimensions, and usually one is not so much interested in simulating a nano system with free surfaces, but rather a part of a bulk system. To embed the small system simulated in a bulk-like system, one can use periodic boundary conditions. These conditions assume that the system is repeated periodically in all space directions and thus any particle close to the boundaries of the simulation box interacts with an image particle, and when a particle crosses one of the face of the central simulation box, one of its images will enter the central box through the opposite face. One of the consequences of the use of periodic boundary conditions is the cut-off of long wavelengths fluctuations, i.e. those having a wavelength larger than the simulation box length. In the context of heat transport simulations this will imply that phonons with long wavelengths are not present.
The strength of MD lies on the versatility and the flexibility in the choie of the atomic potential (eq.1). Generally speaking, the total potential can be decomposed in the sum where the first term represents the one body contribution to the potential and is often associated to an external field like e.g. an electric field. The second and third terms represent respectively two-body and three-body potentials of the inter-atomic interaction. The choice of the inter-atomic potential depends on the problem at hand (type of atomic structure and on the properties to be studied), and on the necessity to model a real system or an experiment. One of the simplest potential used in condensed matter physics is the Lennard-Jones (LJ) potential, where the total potential decomposes in a sum of pair potentials: This is the LJ pair potential, which depends only on the distance between neighbouring particles ij ij r= r  . The interaction energy  identifies with the depth well of the potential, and  is the atomic diameter. Although very simple, the LJ potential is often thought to be a good potential to model rare gas and in particular Argon, using the following parameters: . Due to the algebraic decay of the LJ potential, all the atoms of the simulation box interact. From a practical point of view however, it would be computationally costly to estimate 2 /2 N inter-atomic forces at each time step, with 10000 N  . To bring back the computational cost of the force calculation to   ON , the interatomic pair potential is often truncated at a cut-off radius, usually 2,5 c r=  , which represents typically a neighbourhood of 50 atoms per particle. It should be mentioned at this point, that even if the potential at the radius of truncation is small compared to the well depth:   0,016 c Vr   , truncating the potential may have a non negligible effect on the thermodynamics of the system modelled. As an example, let us estimate the effect of truncating the inter-atomic potential respectively on the internal energy per particle and on the pressure. The respective errors are given by: where we have used the previous value of the cut-off radius and we have assumed that beyond r c , the medium is structureless, i.e. the pair correlation function may be considered close to unity. The previous errors are found to be typically 10% the values of the internal energy and pressure for a condensed phase. This is not worrying since the calculation of thermodynamic quantities can be corrected using the previous estimations, but this has to be kept in mind when calculating a phase diagram for instance.
Although the LJ potential is quite simple, some situations require to use more sophisticated semi-empirical potentials whose parameters are chosen to reproduce either microscopic or macroscopic properties of a model system. These semi-empirical potentials are usually no longer pair potentials like the LJ potential but many-body, where the many-body terms describes how does the potential energy of an atom depends on its coordination.
As an illustration, let us mention the embedded atom model (EAM) used in the context of metals (Daw & Baskes 1984), where the potential wr depends on the microscopic and macroscopic properties to be reproduced, usually the lattice constant, the elastic constants and the sublimation energy of the metal (Daw & Baskes 1984).
Another limitation of LJ potentials is the description of solids crystalizing in non-compact structures. Indeed, Lennard Jones atoms form close-packed structures at low temperatures (usually fcc at low pressure), and thus LJ potentials are also not adapted to model materials which crystallize in opened structures such as diamond. Semi-conductors like silicon or germanium are usually modelled using many-body semi-empirical potentials, in which the many body terms account for the local preferential bond ordering of the semi-conductor. Among the most famous potentials for silicon, let us quote the Tersoff potential (Tersoff 1986(Tersoff , 1988a(Tersoff , 1988b and the Stillinger-Weber potential (Stillinger & Weber 1985). The parameters of the former are chosen so as to reproduce the elastic properties of bulk silicon, while the latter describes satisfactorily the structural properties of liquid silicon. Of course the list of potentials is non exhaustive and there exists plethora of empirical potentials to model as disparate systems as charged systems, liquid crystals, polymers, surfactants, granular media,... and mixtures of them! As we have already briefly mentioned, the output of MD simulations can be compared to experiments if we know how to relate the microscopic state of the material under study to the macroscopic observables typically measured in an experiment. Generally speaking, one can distinguish three classes of information that can be extracted from MD. On the one hand, MD allows to compute global quantities which correspond usually to the thermodynamics of the system modelled. For instance, the internal energy can be expressed as: and the total kinetic energy: where the brackets denote an ensemble average. If sufficient time has been left to the system to reach equilibrium, the ensemble average may be performed by averaging over a sufficient long time thanks to ergodicity. On the other hand, in non-equilibrium situations, different initial conditions (usually the positions and the velocities of the atoms) may be used to generate independent trajectories in phase space. The latter definition of the kinetic energy is intimately related to the definition of the local temperature : where the sum runs over particles i in a small volume V(r) centred around r and which contains n(r) particles. Another important quantity in heat transport simulations is the energy flux, defined as where E i is the total energy of particle i, yielding for pair potentials: For simulations in crystals at low temperatures, the previous energy flux displays large oscillations due to optical phonons which carry a negligible amount of heat, and sometimes it can be more suited to work with the equivalent definition of the flux: where the superscript "0" denotes the equilibrium positions of the atoms. From the knowledge of the local temperature and the energy flux, it is possible to measure thermal conductivity and we would come back to the measurement of transport coefficients in MD in the next section, with a particular emphasis on heat transport.
The other type of information that one can extract from MD is structural. Since the output of a MD simulation is the set of positions ij r  , it is possible to compute the pair distribution where  is the mean number density of the system. Physically, g(r) represents the probability to find two particles separated by a distance r, and   2 4 rg r d r   is the average number of particles located at a distance between r and r+dr from a given particle. Note that sometimes, it may be more useful to compute the structure factor: which is simply related to the pair distribution function through a Fourier transform: Finally, MD may be used to compute the vibrational properties of a model system. In particular, the vibrational density of states (DOS) where the sum runs over the eigenmodes of the system can be estimated using the Fourier transform of the velocity autocorrelation function: where here   vt  is the velocity of an atom. One has to keep in mind that the previous equation is only approximate, as it assumes small deviations from harmonicity, and thus is only valid strictly speaking at low temperatures. Note also that the resulting DOS mixes information about the different polarizations of the atomic vibrations, e.g. transversal or longitudinal, and to determine the DOS relative to a particular polarization, a more refined analysis of the atomic displacement should be undertaken. Apart from the DOS, MD is a common technique to characterise the propagation of phonons in a crystal.
The wavelength dependance of the phonon lifetime may then be compared to the theoretical predictions, given either by Callaway 1959or Holland 1963, see e.g. (McGaughey 2004).
As we have already mentioned, the conversion of the microscopic information given by a MD simulation to macroscopic observables requires to evaluate averages over phase space. The phase space is a multidimensional space generated by the positions and the momenta of a classical system, and has dimension of 6N for a system of N particles. A molecular dynamics simulation generates a sequence of points in phase space as a function of time, whose points belong to the same thermodynamic ensemble. The first MD simulations were performed in the micro-canonical ensemble or NVE ensemble, which corresponds to a fixed number of atoms N, a fixed volume V, and a fixed energy E. Although these simulations are very easy to implement, they do not represent a realistic system as the set of molecules studied is completely isolated. The need to model real situations led to the development of thermostats for MD, which allow to simulate the dynamics of a system in the canonical or the grand-canonical ensemble for instance (Frenkel & Smit, 1996).
Although MD can be used to calculate the thermodynamics, structural and vibrational properties of systems at equilibrium, it is more designed to study the non-equilibrium situations where the system under study is driven by an external force; the latter can be of thermodynamic nature. In the language of out of equilibrium statistical physics, MD can help in determining the relation between the forces and the fluxes. Statistical physics predicts that there is proportionality between the forces and the fluxes if the system under study is driven weakly out of equilibrium. The coefficients of proportionality are called transport coefficients and measure the susceptibility of the system to respond to a given thermodynamic force. For instance, the thermal conductivity quantifies the amount of energy flowing in a material submitted to a temperature gradient. MD is a very powerful tool to calculate the transport coefficients of a model system. To this end, two routes can be traditionally followed: either the model system is submitted to an external force or the thermally induced fluctuations of an internal variable are probed at equilibrium. Statistical physics states indeed that the typical time decay of the spontaneous thermal fluctuations is www.intechopen.com 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 nanostructure.

Limitations of molecular dynamics
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 input in a Boltzmann equation.
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 we couple the system with heat reservoirs.
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 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 modes (Dammak2009).
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 accounted for.

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 semiconductors . 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(Evans, , 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).

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: 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 where c i s t he heat cap a c ity per un it 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 www.intechopen.com 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: 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(Plimpton, , 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).  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.

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): 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: is the low pulsation cut-off introduced by the periodic boundary conditions, and the upper bound is 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 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.

Thermal boundary resistance
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: 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 1 0,1MWm K   (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 www.intechopen.com 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): 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): 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.

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 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 www.intechopen.com 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 / Wm K   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): 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. 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 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 , 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).

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 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, Termentzidis et al, 2009,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.

Modelling Si/Ge superlattices with rough interfaces
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 40 0 , 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 12 0 , which is more than the half of the superlattice period for superlattices with period of 20 0 and from 1ML up to 24 0 for superlattices with period of 40 0 (see fig. 4). The maximal interface height is thus more than half the superlattice's period. 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 0 and we considered only one height of interfaces, the 6 0 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 like interfaces are also examined.

Thermal conductivity of Si/Ge like superlattices with rough interfaces
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 www.intechopen.com roughness. Again, two superlattice periods have been considered 20 0 (left) and 40 0 (right). The figures show also some points named "intra-plane" thermal conductivity, which is defined as the thermal conductivity in the direction of 45° both of the in-plane or cross-plane directions. For the case of infinite roughness (with isosceles periodic triangles) the thermal conductivity in both the in-plane and the cross-plane directions is expected to be equal to this of the intra-plane conductivity, which is verified in our simulations.
For smooth interfaces, it is expected that the transmission of phonons across the interface is specular and depends only on the acoustic impedance mismatch between the bulk materials of the superlattice (Swartz andPohl, 1987, 1989). This is actually not the case since it is observed that the thermal conductivity of the layer decreases when the film thickness decreases (at least enough to be of the same order of magnitude as the phonon mean-free path). For rough surfaces with small roughness the transmission of phonons becomes more diffusive and the transmitted phonons are distributed over a wide range of angles, which induces an additional resistance to in-plane transport. This is consistent with the MD results which conclude to a decrease of the in-plane conductivity with the interfacial roughness. For rough surfaces with large roughness there is a combination of specular and diffusive transmission. This last case shows some similarities with the smooth surface case but now specular reflection is accompanied by back scattering.
This back-scattering explains the existence of a minimum in the thermal conductivity observed for free surfaces, as well as for the in-plane conductivity of superlattices. A further increase in the surface roughness leads to higher thermal conductivity. This further increase in the thermal conductivity is related to the fact that on average, phonons experience less reflections on the asperities of the superlattice, and the mean free path length between scattering events increases, leading to an enhanced in plane transport. For the randomlike roughness, the monotonous decrease in the thermal conductivity in www.intechopen.com increasing the roughness's height can be interpreted if we assume that the phonon scattering at the interface remains diffusive, phenomena of back scattering and specular reflection, playing a secondary role here. Hence, we conclude that the variation in the ratio of interface roughness to the superlattice period can tailor the thermal properties of superlattices. Fig. 6. In-plane and cross-plane thermal conductivity as a function of the interface roughness obtained by EMD and NEMD for superlattice with period 20 0 (top) and with 40 0 (bottom).
The cross-plane and in-plane thermal conductivities obtained by NEMD and EMD for superlattices with various shapes of interfaces are plotted in figure 7. Further details about the modelling of these interfaces and the physical explanation of the results are given in Termentzidis et al, 2011b. It is striking to note that for rough interfaces, the anisotropy of the thermal conductivity is drastically reduced. Regarding the anisotropy between the in-plane and cross-plane directions, we can categorize the interfaces in three different groups, first the smooth interfaces which displays the maximum anisotropy, with the in-plane thermal conductivity being more than twice larger than the cross-plane thermal conductivity; A second group contains totally random interfaces for which the two thermal conductivities exhibit their minimum values and a third group with periodic rough interfaces of a specific shape (triangular, square or cosine) where the thermal anisotropy is negligible compared to the uncertainties of the methods.

Thermal conductivity predictions for GaAs/AlAs superlattices, impact of rough interfaces
The GaAs/AlAs superlattices are important materials mainly for their optical properties.
Simulations of the thermal conductivity of superlattices as a function of their period exhibit a minimum for period around 8 to 10 monolayers (Daly et al 2002, Imamura et al 2003, Chen 2005b, but this minimum is not observed in experimental measurements (Capinski et al 1999). The quality of the interfaces might be the reason to explain this discrepancy. For perfectly smooth interfaces the minimum exists, while for rough interfaces with a height as small as one monolayer the minimum disappears and the thermal conductivity increases with the periodicity of the superlattice. The works of Daly et al 2002 andal 2003 are based on rough interfaces of one atomic layer and with a stochastic distribution of the two types of atoms that compose the superlattice. Fig. 7. In-plane and cross-plane thermal conductivities for several shapes of the interfaces obtained by EMD and NEMD for superlattices with period 20a 0 and with constant height of interfaces 6a 0 .

Modelling rough interfaces for the GaAs/AlAs superlattices
A new modeling of realistic interfaces is considered with the present study (Termentzidis et al 2010(Termentzidis et al , 2011a. Interfaces with square formed islands of one monolayers and pyramide like islands of two monolayers are modelled. Furthermore there are two characteristic lengths of square islands depending if GaAs or AlAs is on top. It has been proven that large scale islands are formed when an AlAs layer grows on a GaAs layer while small scale islands are formed when a GaAs layer grows on an AlAs layer (Tanaka & Sakaki 1987, Jusserand et al 1990. Furthermore interfaces with interdiffusion parts are also considered. Figure Figure 10 shows the predicted cross-plane thermal conductivity as a function of the superlattice period, for a variety of interface configurations. In figure 10 at left the thermal conductivity is plotted for smooth interfaces, for rough interfaces with height of one monolayer and three coverage factors (1%, 10% and 50%) and finally for interfaces with www.intechopen.com height of 2 monolayers. The coverage factor of 1% exhibits the minimum in thermal conductivity as the smooth interfaces. For coverage factor of 10% and 50% the minimum dissapears and the results confirm previous theoretical observations (Daly et al 2002, Imamura et al 2003, while the two defect concentrations do not influence the thermal conductivity. For height of roughness of 2MLs, the thermal conductivity is higher than for smooth interfaces and exhibits a maximum for a 8 0 . In the same figure at the right the influence of interdiffusion of the two species is shown, with similar behavior as the interfaces of height of two MLs. These unexpected results are related with the fact that inelastic scattering could enhance the thermal conductivity through interfaces (Lepi et al 2003). The study shows that the thermal conductivity depends strongly on the detailed description of the interfaces, including height, shape of roughness, interdiffusion of species etc. Changing the structure of the interface can favor or deteriorate the thermal transport through interfaces, and we showed that the interface structure is a relevant parameter which controls the thermal properties of superlattices. Fig. 10. Cross-plane thermal conductivity as a function of the superlattice period with NEMD for GaAs/AlAs systems. Left: smooth interfaces, rough interfaces of height of 1ML with three different concentrations and rough interfaces of height of 2ML are presented. Right: rough interfaces with interdiffusion. The smooth and rough interfaces with a concentration of 50% are also shown for comparison.

Conclusion
We hope we have helped in showing the possibilities of the molecular dynamics technique to probe heat transport in solids, and in particular nanomaterials. Molecular dynamics is a relatively simple and flexible method to be used especially today when stable optimized open source codes have become avalaible: LAMMPS, DLPOLY, GROMACS to name a few. The ever increasing number of publications has helped in resolving controversies regarding heat transfer at the nanoscale, and also getting physical insights in classical problems. A physical understanding of energy transport across two solids-a very simple question-still poses a challenge ! MD simulations may help in observing the scattering of phonons at interfaces, which certainly complements experimental investigations. In the context of nanomaterials, MD is well adapted to characterize ballistic heat transport in nano-objects, although care should be taken not to introduce spurious sources of scattering ! Let us hope that further experimental measurements may improve the modeling of the nanomaterials that we have considered. Molecular Dynamics is a two-volume compendium of the ever-growing applications of molecular dynamics simulations to solve a wider range of scientific and engineering challenges. The contents illustrate the rapid progress on molecular dynamics simulations in many fields of science and technology, such as nanotechnology, energy research, and biology, due to the advances of new dynamics theories and the extraordinary power of today's computers. This first book begins with a general description of underlying theories of molecular dynamics simulations and provides extensive coverage of molecular dynamics simulations in nanotechnology and energy. Coverage of this book includes: Recent advances of molecular dynamics theory Formation and evolution of nanoparticles of up to 106 atoms Diffusion and dissociation of gas and liquid molecules on silicon, metal, or metal organic frameworks Conductivity of ionic species in solid oxides Ion solvation in liquid mixtures Nuclear structures