potential with period a s . . . .

Acknowledgement

First and foremost, I would like to express my sincerest gratitude to my thesis advisor, Dr. Sebastian Volz, for his guidance and encouragement throughout my PhD studies. He is not only an advisor, but also a role model and a friend to me. Back to the year of 2010, I was a student who barely knows anything about nano-scale heat transfer nor computer simulation. It was Sebastian who brought me to this interesting scientific field and unreservedly shared his expertise and experience with me. He is a very busy researcher, but he always spends time discussing with me every week, replies the emails and answers the questions as early as he can. I also appreciate that he provided me with many opportunities to attend conferences from which I benefited a lot. His guidance and personality will continue to inspire me in my future life.

I am greatly indebted to all my thesis committee members: the reviewers Prof. Bernard Perrin, Yuriy Kosevich, Philippe Dollfus, and Dr. Michel Kazan, Konstantinos Termentzidis, Samy Merabia, and Yann Chalopin, for their insight, valuable suggestions, and time. I would also like to express my appreciations to Dr.Yann Chalopin for his help and insightful discussions. I feel very fortunate to have worked with our team members: Shiyun Xiong, Wassim Kassem, Haoxue Han, Benoit Latour, Jordane Soussi, Laurent Tranchant, Horne Kyle, Thomas Antoni, Alexis Pokropivny and Jose-Ediberto Ordonez Miranda ... I gained a lot from our discussions and will always cherish our friendship. I would also like to thank numerous friends for their company and support. Their friendship is a precious lifetime fortune. Particular thanks go to Da Xu and Zhongsen Li, for being there for me when I most needed help and comfort, and also for their company on the tennis court.

I owe a huge thank-you to my wife, Wanni, for her selfless support in the past nine years. Without her encouragement and care, this thesis would not have been possible. Finally, my deepest appreciation goes to my parents and grandmother, for their unconditional support throughout my life, particularly in the past three years of studies abroad. My grandfather helped me make my mind to study abroad and gave me most encouragement, but he left me forever in the first year of my PhD study. This thesis is dedicated to him.

Abstract

This thesis is devoted to the calculation of thermal contact resistance in various molecular systems based on carbon nanotubes (CNTs) and few layer graphene (FLG). This work has been performed through equilibrium molecular dynamics (EMD) simulations. We adopted the temperature difference fluctuations method in our EMD calculations. This method only needs the input of the temperatures of the subsystems whereas the heat flux, that is involved in all the other approaches, remains more difficult to compute in terms of simulation time and algorithm. Firstly, three cases were studied to validate this method, namely: (i) Si/Ge superlattices; (ii) diameter modulated SiC nanowires; and (iii) few-layer graphenes. The validity of the temperature difference fluctuations method is proved by equilibrium and non-equilibrium MD simulations. Then, by using this method, we show that an azide-functionalized polymer (HLK5) has a lower contact resistance with CNT than the one between CNT and PEMA, because HLK5 could form covalent bonds (C-N bonds) with CNT through its tail group azide, while only weak Van der Waals interactions exist in the case of CNT-PEMA contact. The data from our EMD simulations match with the results from experiments in a reasonable range. We then report the thermal contact resistance between FLG and a SiO 2 substrate, which could be tuned with the layer number. Taking advantage of the resistive interface, we show that a SiO 2 /FLG superlattices have a thermal conductivity as low as 0.30 W/mK, exhibiting a promising prospect in nano-scale thermal insulation. In the last part, we investigated the layer number dependence of the crossplane thermal resistances of suspended and supported FLGs. We show that the existence of a silicon dioxide substrate can significantly decrease the crossplane resistances of FLGs with low layer numbers, and the effective thermal conductivities were increased accordingly. The Frenkel-Kontorova model was introduced to explain the substrate-induced band gaps in FLG dispersion relations and the corresponding thermal energy transfer. The enhanced thermal conduction in the cross-plane direction is ascribed to the phonon radiation that occurs at the FLG-substrate interface, which re-distributes the FLG in-plane propagating energy to the cross-plane direction and to the substrate.

Résumé

Cette thèse se concentre sur le calcul des résistances thermiques de contact dans plusieurs systèmes moléculaires à base de nanotubes de carbone (NTCs) et de quelques couches de graphène (QCG). Ce travail a été réalisé en utilisant la méthode de simulation par dynamique moléculaire à l'équilibre (DME). Nous avons utilisé la méthode basée sur les fluctuations de la différence de température dans nos simulations de DME. Cette méthode ne nécessite que l'entrée des températures des sous-systèmes quand le flux thermique, qui intervient dans toutes les autres approches, reste plus difficile à calculer en terme de durée de simulation et d'algorithme. Premièrement, trois cas ont été étudiés pour valider cette méthode : (i) des super-réseaux Si/Ge ; (ii) des nanofils de SiC de diamètres modulés; et (iii) QCGs. La validité de la méthode par fluctuation de la différence de température est démontrée par des simulations de dynamique moléculaire à l'équilibre et hors-équilibre. Ensuite, avec cette méthode, nous montrons qu'un polymère fonctionnalisé azoture (HLK5) a une plus faible résistance de contact avec un NTC que la résistance entre un NTC et un PEMA, car HLK5 forme des liaisons covalentes (C-N) avec un NTC par le groupement azoture, quand seul de faibles interactions de Van der Waals existent dans le cas d'un contact NTC-PEMA. Les données de nos simulations de DME concordent raisonnablement avec les résultats expérimentaux. Nous rapportons ensuite la résistance thermique de contact entre QCGs et un substrat de SiO 2 , qui peut être contrôlée par le nombre de couches de graphène. Avec l'avantage d'une interface résistante, nous montrons que des super-réseaux SiO 2 /QCGs ont une conductivité thermique descendant sous 0.30 W/mK, ce qui est une perspective prometteuse pour l'isolation thermique à l'échelle du nanomètre. Dans la dernière partie, nous recherchons la dépendance de la résistance thermique inter-plan avec le nombre de couches de graphène pour QCGs suspendues ou supportés. Nous montrons que la présence d'un substrat de dioxyde de silicium peut significativement réduire les résistances inter-plan de QCGs possédant peu de couches de graphène, et la conductivité thermique effective est augmentée. Le modèle de Frenkel-Kontorova a été introduit pour expliquer les bandes interdites induites par le substrat dans la relation de dispersion des QCGs et le transfert d'énergie thermique correspondant. L'augmentation de la conduction thermique inter-plan est attribuée au rayonnement de phonons à l'interface QCGs-substrat, qui redistribue l'énergie se propageant dans le plan du QCGs Fig. 3.10). The black and red lines correspond to the groups of atoms far from the surface for small (8.81 nm 2 ) and large (15.67 nm 2 ) cross sections. The green line refers to the DOS in the small section and the first layer, which is the closest to the interface, while the blue one is the DOS in the large section and in the first layer, which is closest to the interface. . . . . . . . . . . . . . . . . . . 3.10 The four considered groups for the calculation of the partial density of states. 1: Group of atoms far from the interface for the small cross-section nanowire (8.81 nm 2 ); 2: group of atoms for the first bilayer closest to the interface; 3, 4: groups of atoms for the large (15.67 nm 2 ) cross-section nanowire far from and the first bilayer closest to the interface. . . . . . . . . . . . . . . . . 3.11 Density of states for four distinct groups (see Fig. 3.10) for frequencies (a) 0-10 THz and (b) 28-40 THz. . . . . . . . . . . . 3.12 Silicon-on-insulator (SOI) circuits with graphene lateral heat spreders attached to the side heat sinks. Microscopic and macroscopic irregularities are present in all practical solid surfaces. Two solid surfaces apparently in contact, therefore, touch each other only at a few individual spots. Since the heat flux lines are constrained to flow through the sparsely spaced actual contact spots, there exists an additional resistance, i.e. an interface resistance to heat flux at a joint. This resistance manifests itself as a sudden temperature drop at the interface [START_REF] Madhusudana | Thermal Contact Conductance[END_REF]. The knowledge of the interface resistance has been more and more important for the successful design of heat dissipation materials and energy conversion devices.

Thermal Interface Materials

With the development of more sophisticated electronic components, including those capable of increasing processing speeds and higher frequencies, having smaller size and more complicated power requirements, relatively extreme temperatures can be generated. However, microprocessors, integrated circuits and other sophisticated electronic components typically operate efficiently only under a certain range of threshold temperatures, typically below 100 • C. The excessive heat generated during operation of these components can not only harm their own performance, but can also degrade the performance and reliability of the overall system and can cause system failure. With the increased need for heat dissipation from microelectronic devices caused by these conditions, thermal management becomes an increasing important element of the design of electronic products. The discontinuity of heat flow at the interface of different components in the system makes the thermal management even more complicated. To reduce the thermal interface resistance, highly thermal conducting mediums are used to eliminate these interstitial air gaps from the interface by conforming to the rough and uneven mating surfaces. These conducting media are typically referred as thermal interface materials (TIMs).

Because the TIM has a greater thermal conductivity than the air it replaces, the heat flow lines are allowed to pass through the gaps and the resistance across the joint decreases, as is shown in Fig. 1.1.

A typical use for a thermal interface material is to thermally connect a computer chip to a cooling module (i.e., a heat sink) to overcome the contact resistance and lack of surface conformity between the heat sink, or the cooling module and the chip or other heat sources. A variety of material types have been developed in response to the changing needs of the electronic packaging market, such as thermal greases, phase change materials and tapes.

Thermoelectric Materials

As for energy conversion devices, one of the types of energy conversion technologies that has received renewed attention is thermoelectric energy conversion. Materials which are able to generate power using the Seebeck effect or refrigerate using the Peltier effect are known as thermoelectric materials [START_REF] Minnich | Bulk nanostructured thermoelectric materials: current research and future prospects[END_REF].

There are three well-known major effects involved in the thermoelectric phenomenon: the Seebeck, Peltier, and Thomson effects. In 1821, Thomas Johann Seebeck discovered that a conductor generates a voltage when subjected to a where V is thermoelectric voltage, ∆T is temperature difference, and α is the Seebeck coefficient (as shown in Fig. 1.2). The Peltier effect refers to the temperature difference induced by a voltage difference. The Thomson effect relates the reversible thermal gradient and electric field in a homogeneous conductor [START_REF] Zheng | Recent advances on thermoelectric materials[END_REF].

The Peltier effect is the basis for many modern-day thermoelectric refrigeration devices, and the Seebeck effect is the basis for thermoelectric power generation devices. The versatility of thermoelectric materials is illustrated in Fig. 1.3, which shows a thermoelectric couple composed of an n-type (negative thermopower and electron carriers) and a p-type (positive thermopower and hole carriers) semiconductor material connected through metallic electrical contact pads. Both refrigeration and power generation may be accomplished using the same module, as shown in Fig. 1.3 (Tritt and Subramanian 2006). The potential of a material for thermoelectric applications is determined in large part to a measure of the material's dimensionless figure of merit,

ZT = α 2 σT k , (1.2)
where α is the Seebeck coefficient, σ is the electrical conductivity, T is the temperature, and k is the thermal conductivity. For a material to be an efficient thermoelectric material, ZT is expected to be as high as possible. Unfortunately, nature does not provide many materials with these properties. Metals have very high electrical conductivity but also very high thermal conductivity. Glasses are the opposite, having very low thermal conductivity but also very low electrical conductivity. But fortunately, the thermoelectrics field has seen remarkable progress in improving thermoelectric properties due to the nanostructured materials such as superlattices and nanocomposites. These materials contain large amount of interfaces, introducing significant interface resistances which decrease the total thermal conductivity of the material while maintaining its electrical properties.

From the examples of thermal interface materials and thermoelectric materials, we can see that thermal interface resistance is playing an important role in the thermal properties of heat dissipation materials and energy conversion materials. Calculating the thermal interface resistance and understanding the corresponding mechanisms thus becomes essential in the future design of micro and nano scale electrical devices. At any finite temperature, the atoms (in an element solid) vibrate about their equilibrium positions. These lattice vibrations may be decomposed into plane and monochromatic waves. Those lattice vibrational waves are characterized by a repetitive and systematic sequence of atomic displacements (longitudinal, transverse, or a combination of the two), which may be characterized by a propagation velocity (v), a wavelength (λ), a wave-vector ( q) (with | q|=2π/λ), a linear frequency (ν) or angular frequency ω=2πν=vq. The energy of lattice vibrational (or elastic) waves is quantized and this quantum is called a phonon.

Transmission of a displacement wave in a solid may be regarded as the coherent motion of a set of phonons, each carrying energy ω(=hν) and momentum q (where is the reduced Planck constant). Thus, thermal vibrations are thermally excited phonons and the thermal conduction in non-metallic solids can be regarded as governed by the creation or annihilation of a phonon. The energy of phonons is of the order of 0.1 eV [START_REF] Keer | Principles of the Solid State[END_REF].

As phonon states comply to the Boltzmann definition of entropy and as those states exist within the canonical ensemble, phonons are governed by the Bose-Einstein statistics. The probability of a phonon with energy ω is given by:

f BE = 1 exp ω k B T -1 , (1.3)
where is the reduced Plank constant, ω is frequency, k B is the Boltzmann constant, and T is temperature. Note that phonons do not have a chemical potential as they are quasi-particles, not real ones.

Phonon dispersion relation

We now consider a linear chain of similar atoms of mass M spaced with lattice constant a as shown in Fig. 1.4. Let u n denote the displacement of the nth atom from its equilibrium position. Considering only the nearest neighbor interactions, the force on the nth atom F n is 1.4) the first term on the right is the increase in length of the bond between atoms n and n + 1, the second term is the increase in length of the bond between atoms n and n -1. β refers to the force constant. We now examine the propagation of waves along the chain of particles. The equation of motion of the nth atom is, using Eq. 1.4, 1.5) we now decompose the displacements into a set of normal modes having the following plane and monochromatic wave form:

F n = β(u n+1 -u n ) -β(u n -u n-1 ), (
M ün = β(u n+1 + u n-1 -2u n ), (
u n = A 0 e i(qna-ωt) , (1.6)
A 0 is the amplitude of the wave, a being the nearest neighbor distance and n the running index. On substituting this trial function in the equation of motion, we can find the solution [START_REF] Kittel | Introduction to Solid State Physics[END_REF]):

-ω 2 M = β(e iqa + e -iqa -2).

(1.7)

Now

e iqa + e -iqa -2 = (e iqa/2 -e -iqa/2 ) 2 = -4sin 2 (qa/2).

(1.8) so that we have a solution if

ω = ±(4β/M ) 1/2 sin(qa/2), (1.9) ω (4β/m) 1/2 = sin 1 2
qa .

(1.10)

The dependence of ω on q is shown in Fig. 1.5, and is called dispersion relation.

Normal mode density of states and energy

The normal modes are equivalent to harmonic oscillators, each of which has an associated wave vector, frequency, and polarization, as defined by Eq. 1.6. Normal mode frequencies and wave vectors are thus defined by the above mentioned dispersion curves. They are completely non-localized spatially. The specification of the normal modes (which is based on the crystal structure and system size) is known as the first quantization. It has nothing to do with quantum mechanics, but indicates that the frequencies and wave numbers available to a crystal are discrete and limited [START_REF] Mcgaughey | Phonon transport in molecular dynamics simulations: formulation and thermal conductivity prediction[END_REF]. This idea can be understood by considering a one-dimensional arrangement of four atoms in a periodic system, as shown in Fig. 1.6. The atoms marked with black dots are equivalent through the application of periodic boundary conditions, which any allowed vibrational mode must satisfy. Two of such waves are shown, with wavelengths of 4a and 2a. Solutions also exist for wavelengths of a and 8a/3. One important distinction between the lattice dynamics problem and the solution of a continuum system (e.g., elastic waves) is that the smallest allowed wavelength is restricted by the spacing between the atoms. A mode with a wavelength of a/2 will be indistinguishable from the mode with wavelength a.

The minimum wavelength (i.e., the maximum wave number) defines the extent of the first Brillouin zone (BZ), the reciprocal space volume accessed by the system. In a mode continuum, there is no lower limit to the wavelength. In both cases, the longest allowed wavelength is determined by the size of the system.

SCIENTIFIC BACKGROUND

For a three dimensional system consists of N atoms, there are 3N independent normal modes with their specific amplitudes and polarizations. To one specific atom, the actual vibration is composed of all normal modes coupled together. The total energy of the lattice vibration is in the form of

E = 3N i=1 ε i , (1.11)
ε i being the ith normal mode energy. According to quantum mechanics, the energy of a normal mode with angular frequency ω is

ε = n + 1 2 ω, (1.12)
where n is the quantum number defining the number of phonons in the mode (positive integers including zero), thus Eq. 1.11 could be written as

E = 3N i=1 n i + 1 2 ω i , (1.13)
ω i is called the phonon energy of the ith mode. The number of phonon modes is frequency dependent, governed by phonon density of states (DOS). The phonon density of states g(ω) of a crystalline solid is defined as follows. The function g(ω)dω corresponds to the fraction of the total number of phonon states in the frequency interval (ω, ω + dω) if g(ω)dω is normalized to unity. Phonon DOS can also be regarded as a projection of the different normal modes of the first Brillouin zone onto the frequency domain.

In the three-dimensional bulk for instance, it is thus proportional to a volume element k 2 dk of the reciprocal space. If now the vibrations only occurs in two dimensional system like films, the density of states becomes proportional to the surface element, kdk. Finally, in the case of a nanowire, where the vibrations can only propagate in one direction, the density of states is proportional to a length element, dk. To sum up, it can be shown that D(k) = k d-1 dk, where d is the dimension of the structure [START_REF] Volz | Thermal Nanosystems and Nanomaterials[END_REF]. The knowledge of g(ω) is essential for understanding the thermodynamic properties.

1.2.2 Size effects in nanostructured materials

1.2.2.

Confinement

The phonon propagation in nanostructures can be largely dependent on the system size. The size effects in nanostructured (nanocrystalline, nanophase or nanocomposite) materials is of great importance both for fundamental considerations and modern technology [START_REF] Uskokovic | The tenth annual conference "YUCOMAT 2008[END_REF]. As previously shown, the amplitude u n of the lattice vibrational wave can be decomposed into plane monochromatic waves, that is, complex exponentials which arguments contain the wave vector q and a time dependence associated with the frequency ω, as given by Eq. 1.6: -ωt) .

u n = A 0 e i(qx
(1.14)

If the wave amplitude is annihilated at the boundary, the wave incident at the stopping point will be reflected with reversal of its phase, as is shown in Fig. 1.7. The superposition of the incident and reflected waves could be modeled as a sum of two exponentials [START_REF] Volz | Thermal Nanosystems and Nanomaterials[END_REF]:

u n ∼ expi(qx -ωt) + expi(-qx -ωt) = cos(qx)e -iωt .
(1.15) The zeros of the cosine function do not depend on time and define the nodes of stationary waves. These stationary waves have zero propagation speed and the phonons are confined rather than propagating. Phonon spatial confinement leads to modification of the phonon dispersion and phonon group velocity, and hence to a change of the phonon scattering rates and transport.

Phonon boundary scattering

The role of the system size effects on thermal properties of nanostructred materials is not only reflected in phonon spatial confinement, but also in phonon boundary scattering. The scattering of phonons on the crystal boundaries has been known since the work of Casimir [START_REF] Casimir | Note on the conduction of heat in crystals[END_REF]), but for many years it was thought to be essentially a low-temperature phenomenon. However, it is now thought that boundary scattering may occur at larger grain sizes and higher temperatures than was previously thought possible. In the case of a nanowire for example, the diameter is much smaller than the phonon mean free path, the scattering at the boundary of the system will enter into the calculation of the thermal conductance. For fully diffusive rough surfaces, the mean free path will be given by the diameter D of a cylindrical nanowire, i.e., Λ = D [START_REF] Casimir | Note on the conduction of heat in crystals[END_REF]). When the mean free path of the phonons Λ is much smaller than the characteristic length scale D of the system, i.e., in the diffusive regime (Λ ≪ D), phonons scatter predominantly with each others and can be treated as a gas at thermal equilibrium. This allows us to define a local temperature T and calculate the effective thermal conductivity from the Fourier law. (see Fig. 1.8(a)) However, when the density of phonons decreases, the distance traveled by a phonon between two collisions can exceed the characteristic length of the structure. Phonons will then enter into more collisions with the walls of the system than with their counterparts within the system. This regime is no longer Brownian, but ballistic, because the phonons will be propagating unhindered in crystalline solids [START_REF] Poole | Encyclopedic Dictionary of Condensed Matter Physics[END_REF]). The ballistic regime can be achieved in a region with dimensions less than the mean free path of phonons. In the pure ballistic regime, the sample size then act as the effective mean free path (see Fig. 1.8(b)).

Ballistic region and diffusive region

Kapitza resistance

Beside the system size, the interface between two phases is another important factor that affects the thermal properties of the whole system. Even if two phases are in perfect contact, the temperature may no longer be continuous across the interface. Kapitza [START_REF] Kapitza | [END_REF] suggested that the temperature at the interface undergoes a change proportional to the normal component of

T out -T in = -R k × Q. (1.16)
Here, R k is Kapitza resistance (or thermal contact resistance) and Q is the heat flux through the surface. The Kapitza resistance appeared first in the context of liquid helium physics, but it exists at any interface [START_REF] Minkowycz | Nanoparticle Heat Transfer and Fluid Flow[END_REF].

The thickness of a liquid or solid layer equivalent to the interface from the thermal point of view is called the Kapitza length:

l K = λ ef f R K , (1.17)
where λ ef f is the effective thermal conductivity. When the Kapitza length is large compared to the problem characteristic dimension, the interfacial contact resistance will play an important role in thermal transport [START_REF] Minkowycz | Nanoparticle Heat Transfer and Fluid Flow[END_REF]).

The Acoustic Mismatch Model and the Diffuse Mismatch Model

There are two main models for evaluating the interface resistance, describing the two limits of the heat exchange process between the media in perfect contact: the acoustic mismatch model (AMM) and the diffuse mismatch model (DMM) [START_REF] Swartz | Thermal boundary resistance[END_REF]).

The Acoustic Mismatch Model

In the acoustic mismatch model, phonons behave like plane waves so that the transmission and reflection coefficients can be calculated from plane-waves dynamics. Moreover, the interface is assumed perfect and no scattering takes place. The reflection is caused by the difference of mass density and speed of sound on the two sides of the interface, and is analogous to the propagation of light across the interface between two optically different materials. For phonons with wavelength much greater than typical interatomic spacings, this continuum approximation might be expected to be accurate. Given this assumption, there are only a few results possible when a phonon is incident on the interface: the phonon can specularly reflect, reflect and mode convert, refract, or refract and mode convert. (That is, the final state is a superposition of these possibilities. The transmission probability is the total fraction of the energy transmitted across the interface.) Fig. 1.9 shows these possibilities [START_REF] Xu | Carbon Nanotube Array Thermal Interfaces[END_REF]. Using the AMM theory, the transmission coefficient for phonon transport between materials A and B may be found to be:

T AB = 4Z A Z B (Z A + Z B ) 2 .
(1.18)

Where Z = ρv is the acoustic impedance, ρ is mass density and v is speed of sound, respectively.

The Diffuse Mismatch Model

In 1987, [START_REF] Swartz | Thermal resistance at interfaces[END_REF] proposed the Diffuse Mismatch Model (DMM) to account for phonon diffusion at the interface. This model assumes that the incident phonons encounter elastic diffuse scattering at the interface and emit to both sides. The probability of being scattered to either side of the interface is proportional to the phonon density of states of that region. From the definition of diffuse scattering, the transmissivity for a phonon from side 1 to side 2, α 1→2 , should be independent of incident angle, wave vector and phonon mode. Hence,

α 1→2 = 1 -α 2→1 , (1.19)
which means that because phonons have no memory of their history, the reflectivity from one side must equal the tranmissivity from the other. With the consideration of incident angle (θ) effects, Swartz and Pohl suggested that the averaged transmissivity should be defined as, 1.20) where p is the phonon mode, and this equation is universal for both AMM and DMM models. Then with the detailed balance principle, the Debye approximation and generalization of different longitudinal and transverse phonon velocities, the averaged transmissivity can be calculated as

Γ 1→2,p = ω 0 α 1→2 (θ, p) • cos(θ) • sin(θ)dω, (
Γ 1→2,p = 1 2 p (v 2→1,p ) -2 p (v 1→2,p ) -2 + p (v 2→1,p ) -2 .
(1.21)

The elastic scattering assumption does not allow phonon frequency change, i.e. phonons cannot convert from mode to mode. Generally this assumption is not true. For instance, if two materials A and B have a large mismatch in their phonon dispersion relations, a certain range of phonons cannot be transmitted across the interface unless some mode conversion occurs at the interface. Thus, some inelastic scattering has to happen to change the phonon frequency. This has been proven by MD simulations (Chalopin and Volz 2013a;[START_REF] Termentzidis | Modulated sic nanowires: Molecular dynamics study of their thermal properties[END_REF].

The main difference between these models is that the AMM assumes a completely specular interface and the DMM assumes a completely diffuse interface.

The AMM works well at low temperatures (<1K), when the phonon wavelengths are much longer than the interface roughness. Reversely, the DMM works better for higher temperatures when the phonon wavelengths are comparable to the interface roughness. It should be noted that neither model accurately fits all experimental data. This is probably because of the fact that they do not consider roughness, polycrystallinity, contaminants and defects that occur at real interfaces. Despite this, these models are important because they provide insight into the conduction mechanisms across interfaces [START_REF] Wang | Energy Transport and Conversion in Nanostructured Materials[END_REF]).

Ambiguity in defining thermal resistance

Thermal resistance is an ambiguous concept relating equilibrium and nonequilibrium quantities: consider two systems with different temperatures T 1 and T 2 (see Fig. 1.10), their temperatures being well defined since they are related to a thermal bath respectively. But if the two systems are exchanging a given amount of heat flux, it is very clear that this heat flux should not be too strong to disturb the thermal equilibrium in the two systems. Otherwise we could not define properly the equilibrium temperatures. This challenge always remains when defining the thermal resistance, because there are always some perturbations at the connection between the two systems, especially when the thermal coupling between the two is very strong as is the case with films in perfect contact. 1.2.6.1 Atomistic Green's Function AGF is best suitable for the cases where harmonic scattering is predominant, such as in the case where the size of a system is comparable or shorter than the bulk mean free path, where phonons travel nearly ballistically and anharmonic effects are less important. It offers the flexibility to incorporate complicated geometries and material compositions, as well as atomic-scale representation.

The motion of a generic harmonic system is completely determined by its force constant matrix K. This matrix corresponds to the second derivatives of the total energy of the system, E, with respect to any pair of degrees of freedom, u i and u j : .22) For the Newtonian momentum conservation condition to be fulfilled, the matrix K satisfies various sum rules, essential when checking the correctness of the force constant matrix [START_REF] Mingo | Phonon transmission through defects in carbon nanotubes from first principles[END_REF]). The vibrational normal modes of the system, Ψ n e -iωnt , satisfy the equation

K ij = ∂ 2 E ∂u i ∂u j . ( 1
M i ω 2 n Ψ (n) i - j K ij Ψ (n) j = 0. (1.23)
Here, M i is the mass of the atom to which degree of freedom i belongs. Let a bar on top of a symbol denote that the symbol is a rank 1 tensor, symbols like K, I, etc. denote rank 2 tensors (matrices). The above equation can be transformed into the matrix dynamical equation:

(ω 2 n I -K)Ψ n = 0, (1.24)
where we have defined the harmonic matrix:

K ij = (M i M j ) -1/2 ∂ 2 E/∂u i ∂u j .
The resolvent Green's function matrix for the system is defined as

G(ω 2 ) = (ω 2 + iδ)I -K -1 , (1.25)
with δ → 0 + . The resolvent is a very useful tool to deal with the effect of perturbations to the system. One can compute how the eigenfunctions are modified when the force constant matrix is modified on a certain region.

The major limitation of this method is the underlying harmonic assumption.

A publication by Mingo [START_REF] Mingo | Anharmonic phonon flow through molecular-sized junctions[END_REF]) introduces a theoretical framework to include anharmonic effects. Another limitation of the method is the fact that the lattice coordination is fixed from the beginning, so that no atomic reorganization is allowed. Thus, the approach can not be applied to systems where atomic diffusion or convection processes take place [START_REF] Zhang | Simulation of phonon transport across a non-polar nanowire junction using an atomistic green function method[END_REF]).

Boltzmann Transport Equation

If the wave nature of the particles is not important in the energy transport, the transport phenomena can be modelled using particle transport theory. The governing equation for the particle transport is the Boltzmann Transport Equation (BTE). The BTE describes the evolution of a particle probability distribution over time, space, and wave vector, as particles experience a series of scattering events and external forces. In principle the BTE is applicable to all energy carriers, including electrons, phonons, or photons as long as the wavelengths of the carriers are small compared to the characteristic length of the object of interest [START_REF] Wong | Thermal Transport for Applications in Micro/Nanomachining[END_REF]. The BTE can be shown to yield the Fourier conduction equation at large scales. The BTE can be written in a general form as:

∂N ∂t + v • ∇N + a • ∂N ∂ v = ∂N ∂t scatter , (1.26)
where v is the particle velocity, a is the particle acceleration rate and N (t, r, v) is the distribution function. At equilibrium, N would be a Maxwell-Boltzmann distribution for gas molecules and a Bose-Einstein distribution for phonons, but in the absence of equilibrium, it may depart significantly from these distributions. The left hand side of Eq. 1.26 contains the transient and drift terms respectively, and the right hand side is the scattering term due to phonon interaction processes.

Molecular Dynamics

Molecular dynamics (MD) is another important method to tackle the thermal properties of materials at finite temperature and for bridging of real and phonon space analysis techniques. A brief introduction to MD method will be given in Chapter 2. MD simulations allow for the natural inclusion of anharmonic effects and for atomic-level observations that are not possible in experiments.

As large-scale computing capabilities continue to grow and component sizes to decrease, the dimensions of systems accessible with MD approach those of real devices (McGaughey et al. 2006). A diagram of the tools for thermal property calculation with corresponding spatial and temporal characteristic dimension is shown in Fig. 1.11. MD is the tool that is used for calculating thermal properties of interesting systems throughout this report.

Organization

We begin Chapter 2 with a brief introduction to MD method. Non-equilibrium and equilibrium MD methods for calculating thermal interface resistance are then presented and compared.

ORGANIZATION

In Chapter 3, three cases were studied to validate the equilibrium MD simulation method based on the temperature difference fluctuations. Firstly, the interface resistance in Si/Ge supterlattice was calculated from this method, and the results were compared with the ones obtained from the exchanged power fluctuations method. The derived thermal conductivities were validated by the values calculated from Green-Kubo formula. Secondly, we defined the thermal resistance of diameter modulated SiC nanowires as a constriction resistance and showed that it is worth twice the internal thermal resistance of the constant diameter nanowire. The constriction resistances calculated from equilibrium MD simulation are in accordance with Non-equilibrium MD simulation data in term of tendency. Thirdly, we used this method to calculate the thermal resistance between few layer graphenes (FLGs) with different layer numbers.

The obtained results were validated by non-equilibrium molecular dynamics.

We then take advantage of the equilibrium technique providing the contact resistance to tackle molecular systems especially carbon nanotubes and graphene.

In Chapter 4, the data of interface resistances between a carbon nanotube (CNT) and polymers (PEMA and azide polymer) calculated from MD are shown. We found that compared with Van der Waals interactions, covalent bonds between CNT and polymer could reduce the resistance. Phonon density of states was calculated to better understand the mechanism of the resistance reduction. The simulated results were further validated by experiments.

Chapter 5 addresses the thermal resistance related to graphene. The contact resistances between FLG and SiO 2 as a function of FLG layer number were reported and compared with experimental data. A varying and large contact resistance between graphene and SiO 2 was found. Taking advantage of such a resistive interface, MD simulations show that SiO 2 /FLG superlattices have thermal conductivity as low as 0.30 W/mK, exhibiting a promising prospect in nano-scale thermal insulation.

Finally, in Chapter 6, we report the layer number dependence of the cross-plane thermal resistances (R c ) of suspended and supported few layer graphene (FLG) by equilibrium MD simulations. The concerned layer number ranges from 2 to 30. The value of R c of suspended FLG was found to be decreasing against layer number, and approaches the graphite limit. The existence of a silicon dioxide substrate significantly decreased the R c of FLG, compared to the suspended ones. But this effect is eliminated when the layer number reaches seven. For FLG layer number larger than seven, R c tends to be consistent in both suspended and supported FLG configurations. The effects of the substrate on the FLG cross-plane thermal properties were discussed. Phonon dispersion relations and density of states were shown to help understanding the mechanisms involved in these phenomena. The Frenkel-Kontorova model and sound lines were introduced to explain the substrate-induced band gaps in FLG dispersion relation and the corresponding effects on FLG cross-plane thermal properties.

Chapter 2

Thermal contact resistance calculation from Molecular Dynamics simulation

Introduction to molecular dynamics simulations

Overview

Molecular dynamics (MD) simulations consist in the computation of the motions of individual molecules in order to model the behaviors of solids, liquids and gases. It provides a technique to track the trajectories and velocities of atoms and molecules in a nanoscale system. MD method can be used effectively to investigate a variety of nanoscale transport processes. The two main components in MD modeling are the atomic structure, and the interatomic potential function (u). The functional form of u can be fitted to match experimental data, or quantum mechanics calculations. Ideally, a potential energy function must combine simplicity, accuracy, and transferability. Generally, the total potential energy (U ) of a system of MD particles can be expressed as a sum of terms involving single, pairs, triplets, and so forth, of atoms: (Allen and Tildesley 1989)

U = i u 1 (r i ) + i j>i u 2 (r i , r j ) + i j>i k>j>i u 3 (r i , r j , r k )+ ••• (2.1)
The summation prevents the inclusion of the interaction potential between a set of atoms more than once. The first term in the potential energy funciton, u 1 (r i ), can be used to incorporate the effect of an external field acting on the system. All the remaining terms represent particle-particle interactions.

The motion of each atom in MD is governed by Newton's second law, which is given by:

d 2 r i dt 2 = F i m i (2.2)
where t is time, m i is the mass of atom i, r i is its position vector, and F i is the total force vector acting on the atom. In the case of the simple pair potential, F i is calculated as follows:

F i = i =j F ij = - i =j ∂u(r ij ) ∂r ij = - i =j ∂u(r ij ) ∂r ij ∂r ij ∂r ij = - i =j ∂u(r ij ) ∂r ij rij , (2.3)
where F ij is the force applied on atom i by atom j, r ij = r i -r j is the separation distance vector, with rij as its unit vector. A depiction of these vectors is shown in Fig. 2.1 (Issa 2012).

Equations of motion for MD systems

Eq. 2.2 must be solved to advance the trajectory of the system in time. One of the most commonly used methods in solving the atomic equations of motion is the finite-difference Verlet algorithm [START_REF] Verlet | Computer "experiments" on classical fluids. i. thermodynamical properties of lennard-jones molecules[END_REF]:

r t+δt = 2r t -r t-δt + δt 2 a t (2.4) v t+δt = 1 2 δt[r t+δt + r t-δt ] (2.5)
where δt is the time step, v, and a = F/m, are the velocity and acceleration vectors, respectively. In order to include velocity in Eq. 2.4 so as to facilitate the adding and removing kinetic energy from a MD system, the velocity-Verlet algorithm was proposed [START_REF] Swope | A computer simulation method for the calculation of equilibrium constants for the formation of physical clusters of molecules: Application to small water clusters[END_REF]:

r t+δt = r t + δt v t + 1 2 δt 2 a t (2.6) v t+ 1 2 δt = v t + 1 2 δt a t (2.7) v t+δt = v t+ 1 2 δt + 1 2 δt a t+δt (2.8)
The position around time t could be obtained from Taylor's expansion:

r t+δt = r t + ṙt δt + 1 2 rt δt 2 + 1 3! ... r t δt 3 + • • • (2.9) r t-δt = r t -ṙt δt + 1 2 rt δt 2 - 1 3! ... r t δt 3 + • • • (2.10)
Adding Eq. 2.9 to Eq. 2.10 and removing the higher order terms, one obtains the expression of the position of one atom in terms of its acceleration:

r t+δt ≈ 2r t -r t-δt + rt δt 2 (2.11)
Subtract Eq. 2.9 by Eq. 2.10, one gets the expression of velocity:

ṙt ≈ r t+δt -r t-δt 2δt (2.
12)

The velocity-Verlet algorithm is mathematically equivalent to the original Verlet algorithm. The choice of the time step should be limited to ensure that the system dynamics with the highest frequencies are sufficiently captured and that the total energy of the system is conserved during a run (Issa 2012).

Initial state

In order for MD to serve a useful purpose, it must be capable of sampling a representative region of the total phase space of the system. An obvious corollary of this requirement is that the results of a simulation of adequate duration are insensitive to the initial state, so that any convenient initial state is allowed. A particularly simple choice is to start with the atoms at the sites of the desired lattice to give the desired density.

Once initial coordinates are specified, it remains to set the initial velocities. The initial velocities are assigned random directions and a fixed magnitude based on temperature; they are also adjusted to ensure that the center of mass of the system is at rest, thereby eliminating any overall flow. The speed of equilibration to a state in which there is no memory of this arbitrarily selected initial configuration is normally quite rapid, so that more careful attempts at constructing a 'typical' state are of little benefit. The assignment of initial velocities is generally done by "sampling" the velocities from a Maxwell-Boltzmann distribution, taking care to ensure that the sampled velocities are consistent with any constraints imposed on the system. The Maxwell-Boltzmann distribution for the velocity v of a particles system of mass m at temperature T is

f (v) = m 2πk B T 1/2 e -mv 2 /2k B T .
(2.13)

The distribution f (v) is an example of a Gaussian probability distribution. More generally, if x is a Gaussian random variable with zero mean, its probability distribution is

f (x) = 1 2πσ 2 1/2 e -x 2 /2σ 2 , (2.14)
where σ is the width of the Gaussian (see Fig. 2.2). Here, f (x)dx is the probability that a given value of the variable, x, will lie in an interval between x and x + dx. (Tuckerman 2008) The size of systems accessible to conventional MD simulations are limited to a few hundred nanometers. This limitation can be overcome by employing periodic boundary conditions as depicted in Fig. 2.3. As one particle moves through a boundary, its images will move across their corresponding boundaries. The number density in the central box is conserved. The use of periodic boundary conditions inhibits the occurrence of long-wavelength fluctuations. For a cube of side L, the periodicity will suppress any density waves with a wavelength greater than L, otherwise there will be a substantial interaction between a particle and its own images in neighboring boxes. The methods used to cope with long-range potentials are discussed in [START_REF] Rapaport | The Art of Molecular Dynamics Simulation[END_REF]). Problems can also arise if there are strong correlations between atoms separated by distances approaching the region size, because periodic wraparound can then lead to spurious effects. One example is the vibration of an atom producing what are essentially sound waves; the disturbance, if not sufficiently attenuated, can propagate around the system and eventually return to affect the atom itself.

Periodic boundary conditions

Even with periodic boundaries, finite-size effects are still present, so how big does the system have to be before they can be neglected? The answer depends on the kind of system and the properties of interest. As a minimal requirement, the size should exceed the range of any significant correlations, but there may be more subtle effects even in larger systems.

Thermodynamic properties

In MD simulations, thermodynamic properties are calculated from the dynamics of atoms in the system, in what could be viewed as a "bottom up" approach. The total energy of the system, E, is the sum of the total potential energy U , and kinetic energy K:

E = U + K = N -1 i N j>i u(r ij ) + 1 2 N i m i v 2 i , (2.15)
where the brackets denote ensemble average, N is the total number of atoms in the system, and the pair potential assumption is implied in Eq. 2.1. In heat transfer studies, two important thermodynamic parameters are the temperature (T), and pressure (P), of the system. For an atomic system these are calculated using:

T = 2 K 3N k B = N i m i v 2 i 3N k B (2.16) P = N k B T V + 1 3V N -1 i N j>i r ij • F ij , (2.17)
where k B is the Boltzmann constant. A quick check of Eq. 2.16 shows that K = (3N ) 1 2 k B T , i.e. each degree of freedom contributes on average an equal amount of 1 2 k B T , in accordance with the equipartition principle. The two terms in the pressure expression Eq. 4.2 represent contributions from the system kinetic and potential energies, respectively. MD simulation suffers from several important and fundamental challenges, the most important of which is that the simulations are classical, that is, the trajectories of the atoms are integrated according to Newtonian mechanics without any quantum-mechanical phenomena included. This is usually justified by stating that for most elements at room temperature the atoms are sufficiently heavy that their de Broglie wavelength is negligible and the atoms can be treated classically. Nevertheless, the permitted energy of vibrational modes is quantised with occupation E n = ω (n+ 1 2 ), and the quanta of energy (phonons) are bosons and so fill the available states according to the Bose-Einstein distribution, as is stated in Chapt.1. Not all modes have a uniform amount of energy: lowfrequency modes will hold more energy than higher-frequency modes according to the relation E i (T ) = ω/(exp(β ω i ) -1), where β = 1/k B T , is the reciprocal temperature. At temperatures below the Debye temperature, T D , the high-frequency modes will be unoccupied until reaching the zero temperature where they only possess zero point energy. This stands in direct contrast to a classical system in which all the modes have on average the same energy k B T /2. For this reason MD simulations are usually performed at temperatures above T D , where at least the participation of all vibrational modes is a reasonable assumption [START_REF] Bichoutskaia | Computational Nanoscience[END_REF]. Alternatively, simulations may be performed at low temperatures if the phenomena of interest only involve low-frequency modes, such as the cross-plane interactions of graphene sheets.

Methods for calculating thermal contact resis-

tance from Molecular Dynamics simulations

Contact resistance from Non-equilibrium molecular dynamics simulations

Typically molecular dynamics (MD) simulations are conducted under equilibrium (or near equilibrium) conditions of constant energy (E) or temperature (T), constant atomic number (N) or chemical potential (µ), and constant volume (V) or pressure (P). The most widely adopted statistical sampling ensembles are the microcanonical ensemble (NVE), the canonical ensemble (NVT) and isobaric-isothermal ensemble (NPT) for solids. However, there is increasing interest in conducting MD simulations which do not fall within the classification of these classical ensembles. These broadly fall into the category of non-equilibrium molecular dynamics (NEMD). The temperature, pressure, volume and energy of the system vary spatially and temporally during the simulation in a complex way (Dumitrica 2010).

METHODS FOR CALCULATING THERMAL CONTACT RESISTANCE FROM MOLECULAR DYNAMICS SIMULATIONS

Steady-state NEMD simulation

The direct method to predict the thermal resistance is a non-equilibrium, steady-state approach. By imposing a one-dimensional temperature gradient on a simulation cell and measuring the resulting heat flux (or, imposing a heat flux and measuring the resulting temperature gradient), one can predict the interfacial thermal conductance, G, by using:

q = G△T, (2.18)
where q is the heat flux, which is defined as the amount of thermal energy transfer rate per unit area perpendicular to the direction of the heat flux and △T is the temperature difference at the interface. The thermal contact resistance is the inverse of conductance.

In a NEMD simulation for thermal resistance calculation, either a constant temperature difference or an applied heat flux is imposed by altering the atomic dynamics in the localized heat sink and source regions. In the constant temperature difference method, the temperatures of the sink and source regions are controlled at the two different values by scaling the atomic velocities in these two regions. The energy differences before and after velocity scaling processes in both the sink and source regions yield the heat flux when divided by the time step. In the heat flux method, the velocities of the atoms in the sink and source regions are scaled so that the same amount of energy is added to the sources and taken out from the sinks. The amount of energy difference before and after the scaling is determined by the energy difference between the fastest atom in the sink region and the slowest atom in the source region. The direct method has a physical intuitiveness and is analogous to an experimental technique.

Transient NEMD simulation

The direct method provides the most efficient means of obtaining the non-linear response of the system, but the major limitation is that the signal to noise ratio goes to zero in the weak field limit. Transient Time Correlation Function (TTCF) formalism was developed to bridge the gap between the Green-Kubo formalism and the direct method [START_REF] Kanellopoulos | Recent Advances in Gas Separation by Microporous Ceramic Membranes[END_REF].

We consider an external force field F being switched at time t=0 onto a system at equilibrium. The response of any time-dependent property of the system, B(t), is given by

B(t) = B(0) + 1 k B T t 0 dt ′ B(t ′ ) Ḣ(0) . (2.19)
The time derivative of the Hamiltonian is viewed as "dissipative heat" in the phenomenological form of a flux times a perturbing thermodynamic force, i.e. Ḣ = -J Q • F where J Q is the dissipative flux. In the case of thermal conductivity, the fictitious force or "heat field" F replaces minus the temperature gradient divided by the mean temperature but still induces an additional heat flux △J Q (t) in the system. Setting B = J Q and substituting for Ḣ in Eq. 2.19 gives the so-called transient time correlation function formula for the thermal conductivity at arbitrary applied field. We then have

λ(F) = lim t→+∞ △J Q (t) F = V 3k B T 2 t 0 dt ′ △J Q (t ′ ) • J Q (0) .
(2.20)

TTCF has been extended to account for systems subject to time-periodic external fields [START_REF] Petravic | Nonlinear response for time-dependent external fields[END_REF] and has been used successfully [START_REF] Todd | Application of transient-time correlation functions to nonequilibrium molecular-dynamics simulations of elongational flow[END_REF] to calculate the non-linear response of a system of simple spheres to elongational flow.

Contact resistance from equilibrium molecular dynamics simulations

For NEMD simulations, the boundary conditions would ideally be applied at a considerable distance from the region of interest to minimize their effect on its dynamics. However, the size of an atomistic system that can be considered is limited by finite computational resources due to their large numbers of degrees of freedom. Therefore most NEMD simulations are affected by their boundaries during the course of a simulation [START_REF] Dumitrica | Trends in Computational Nanomechanics: Transcending Length and Time Scales[END_REF]. Constant rescaling of the atomic velocities in the thermostats introduces artificial scattering and, consequently, the characteristic time of the thermostats is a critical parameter. It is possible to move away from the non-equilibrium approaches and instead to simulate a system at equilibrium in the microcanonical ensemble (NVE), and make use of the fluctuation-dissipation theorem that relates the linear response properties of a system out of equilibrium to the dissipation of thermal fluctuations within the system at equilibrium. Equilibrium molecular dynamics simulation does not suffer from the limitations that restraining NEMD, since it relies on the natural fluctuations of the physical quantities rather than perturbing the system. Three different methods for calculating contact resistance from equilibrium will be presented in this section.

Contact resistance from exchanged power fluctuations method

A microscopic approach was proposed recently by Chalopin et al. (Chalopin et al. 2012) for calculating thermal conductance at interfaces based on the tracking of equilibrium fluctuations:

G = 1 k B T 2 dH A (t) dt H B (t) , (2.21)
where H A,B denotes the Hamiltonian of two systems, A and B, in thermal interaction at the same equilibrium temperature. A more familiar looking formula could be derived by considering

H B (t) = t -∞ ḢB (t ′ )dt ′ .
Using the change of variable t ′ = t -τ , we can then re-express the thermal conductance as the time integral over the dissipated power correlator:

G = 1 k B T 2 t -∞ ḢA (t) ḢB (t ′ ) dt ′ = - 1 k B T 2 ∞ 0 ḢA (t) ḢB (t -τ ) dτ = - 1 k B T 2 ∞ 0 ḢA (0) ḢB (-τ ) dτ = - 1 k B T 2 ∞ 0 ḢA (τ ) ḢB (0) dτ = - 1 k B T 2 ∞ 0 ḢA (τ ) ḢA (0) dτ, (2.22)
which looks like the usual Green-Kubo formula. In the above derivation, we have assumed time-translational invariance, energy conservation, ḢB = -ḢA , and used ḢB (t ′ ) = -ḢB (t -τ ), where the time derivative in the left term is with respect to t ′ and in the right term it is with respect to τ . It seems that the use of Eq. 2.21 would be computationally more advantageous, and the conductance can be calculated within a canonical ensemble, requiring just an ensemble average over uncorrelated equilibrium configurations generated even with a thermostat.

The time derivative of H A can be obtained by considering the Hamiltonian equations of motion:

dH A (t) dt = i∈A∪B ∂H A ∂p i ṗi + ∂H A ∂q i qi , (2.23)
where

ṗi ≡ - ∂H ∂q i = F A i + F B i , (2.24) qi ≡ - ∂H ∂p i = p i m i = v i , (2.25
)

∂H A ∂q i = -F A i + 1 2 F B i f or i ∈ A, (2.26) ∂H A ∂q i = - 1 2 F A i f or i ∈ B, (2.27)
F S i being the force from the side S on atom i. Then we have : which is the net power (or heat flux) absorbed by system A. This method is general and can be applied to any potential, as long as the total system potential energy can be written in terms of individual atom potential energy contributions. The thermal interface conductance in a Si/Ge superlattice with different period thickness was calculated, and the effective thermal conductivity λ e deduced from the interfacial contributions can be expressed as λ e = G • L P , here G denotes the thermal conductance of the superlattice interface of period L P predicted from MD. The derived thermal conductivity was validated by that calculated from Green-Kubo method, as shown in Fig. 2.5 (Chalopin et al. 2012).

dH A (t) dt = 1 2   i∈A v i F B i - j∈B F A j v j   , (2.28)

Contact resistance from temperature difference fluctuations method

Rajabpour and Volz developed an EMD simulation method for calculating thermal resistance based on the temperature difference fluctuations in 2010 (Rajabpour and Volz 2010). The thermal conductance G between two heat baths 1 and 2 at equilibrium temperatures T 1 and T 2 respectively and exchanging a net heat flux Q is classically defined by a linear relation including the temperature difference ∆T = T 1 -T 2 as follows:

Q = G∆T , (2.29)
A conventional expression used for measuring and computing the thermal conductance is based on the time dependent energy conservation equation including the exchanged heat flux Q. The substraction between the energy conservation equations of both bodies:

∂T 1 ∂t = 1 k B N 1 G∆T , (2.30)
and

∂T 2 ∂t = 1 k B N 2 G∆T , (2.31)
yields the equation of the temperature difference ∆T :

∂∆T ∂t = 1 k B 1 N 1 + 1 N 2 G∆T , (2.32)
where k B is the Boltzmann constant and N i is the number of degrees of freedom of the kinetic energy in the heat bath i (N i = 3 × N a where N a is the number of atoms). The solution for ∆T is a decaying exponential function including a relaxation time τ , which defines the reverse of the thermal conductance, i.e. the resistance R, as:

Rk B = τ 1 N 1 + 1 N 2 , (2.33)
Chapter 2 -Thermal contact resistance calculation from Molecular Dynamics simulation

33

This expression is typically used in transient (NEMD) [START_REF] Xue | Two regimes of thermal resistance at a liquid-solid interface[END_REF]) experiments where the time τ is estimated from the exponential decay of the temperature difference after releasing the thermal constraints between both bodies initially set at two different temperatures. However, the linear response theory establishes the heat flux as a convolution product between the conductance and the temperature difference. Eq. 2.32 corrected with a Langevin flux F(t) can be written as:

d∆T dt = - 1 k B 1 N 1 + 1 N 2 t 0 G(t -t ′ )∆T (t ′ )dt ′ + F (t), (2.34)
Note that the Langevin term, which has a white noise time evolution, mimics the fluctuations of equilibrium quantities as proposed in Fokker-Planck equations. Multiplying this equation by ∆T (0) and performing phase average allows for removing the F(t) dependent term in Eq. 2.34 because the Langevin flux and ∆T have an infinitely short correlation time: ∆T (0)F (t) = 0. The resulting equation is:

d ∆T (0)∆T (t) dt = - 1 k B 1 N 1 + 1 N 2 t 0 G(t -t ′ ) ∆T (0)∆T (t ′ ) dt ′ , (2.35)
The brackets denote the phase average or the well-known autocorrelation function (ACF). This function here applies to the equilibrium fluctuations of the temperature difference ∆T . To simplify the notations, we propose to define

C(t) = ∆T (0)∆T (t) and C v = k B / 1 N 1 + 1 N 2 to rewrite Eq. 2.35 as: C v dC(t) dt = -G(t) * C(t), (2.36)
where the star * is the convolution product. Taking the Laplace transform of Eq. 2.36 yields:

sC(s) -C(t = 0) = -G(s)C(s) C v , (2.37)
At the static limit (s=0), the reverse of the spectral conductance, R(s), is yielded as: Two comments need to be drawn at this stage. Firstly, Eq. 2.38 does reveal the exact expression of the relaxation time:

Rk B = C(s = 0) C(t = 0) 1 N 1 + 1 N 2 = ∞ 0 ∆T (0)∆T (t) dt ∆T 2 (0) 1 N 1 + 1 N 2 , (2.38)
τ = ∞ 0 ∆T (0)∆T (t) / ∆T 2 (0) dt, (2.39)
Chapter 2 -Thermal contact resistance calculation from Molecular Dynamics simulation

35

as the time integral of the normalized temperature difference autocorrelation. Secondly, to describe transient evolutions, Eq. 2.30-2.32 appear to be valid only in the range where τ is much smaller than the characteristic time of the perturbation ∆T (or of Q), so that the convolution product in Eq. 2.35 can be reduced to a simple product. Also note that Eq. 2.39 is particularly useful in MD simulations that directly process velocities and temperature whereas the heat flux, that is involved in all the other approaches, remains more difficult to compute in terms of simulation time and algorithm.

The thermal resistance between two Lennard-Jones crystals by equilibrium MD was calculated in Ref. [START_REF] Rajabpour | Thermal boundary resistance from mode energy relaxation times: Case study of argon-like crystals by molecular dynamics[END_REF]. The system is composed of an argon crystal wire with a mass mismatch, as is illustrated in Fig. 2.6. Fig. 2.7 reports the normalized auto-correlation functions (ACFs) of temperature difference fluctuations.

It is important to note that the atoms in one subsystem used in computing the mean temperature T i (i=1,2) should be the ones having direct interactions with the other subsystem across the interface. These atoms can be determined by the potential cut-off radius. The heat flux Q obtained from the hamiltonian equation of motion of two interacting phonon systems whose hamiltonians are

H = H A + H B with Q = -ḢA = [H, H A ] (Chalopin and Volz 2013b), Q(t) = 1/2 i∈A,j∈B α,δ∈{x,y,z} k α,δ i,j ( uα i (t ′ )u δ j (t) -u α i (t) uδ j (t ′ ) ), (2.40)
u i refers to the atomic displacement and ui the instantaneous velocity of atom i. The α and δ exponents refer to the x, y or z component. k α,δ i,j is the interatomic force constant between atoms i and j. According to this equation, Q only depends on the interacting atoms. For the non-interacting atoms, k α,δ i,j is zero and those atoms do not contribute to the heat flux.

The frequency dependent conductance G is only depending on the thermodynamic quantities Q and the equilibrium temperature T 0 [START_REF] Rajabpour | Thermal boundary resistance from mode energy relaxation times: Case study of argon-like crystals by molecular dynamics[END_REF]:

G(ω) = 1 k B T 2 0 ∞ 0 Q(0)Q(t) e -iωt dt.
(2.41) T 0 does not include any microscopic information, thus G depends on the same microscopic quantities as Q, and therefore, G (or resistance R) is only dependent on the atoms that are interacting across the interface.

The temperature difference fluctuation method is the approach that is used to calculate thermal contact resistance throughout this report.

Contact resistance generated by an interstitial body from Green-Kubo formula

The Fourier law states that one-dimensional, steady state heat flux rate between two isothermal surfaces is proportional to the temperature gradient causing the heat flow and the area normal to the direction of heat flux [START_REF] Thirumaleshwar | Fundamentals of Heat and Mass Transfer[END_REF]):

Q = -λA dT dx , (2.42)
where Q is the heat flux rate in x-direction, A is the area normal to the direction of heat flux, dT /dx is the temperature gradient, and λ is the thermal conductivity.

The Green-Kubo formula relates the thermal conductivity to the autocorrelation function of the fluctuations in the heat flux q, at equilibrium [START_REF] Kubo | Statistical Physics II[END_REF]:

λ = V k B T 2 +∞ 0 q(0) • q(t) dt.
(2.43)

The brackets denote the heat flux autocorrelation function, k B is Boltzmann constant, and T is the equilibrium temperature.

The thermal conductance G between two heat baths 1 and 2 at equilibrium temperatures T 1 and T 2 , respectively, and exchanging a net heat flux Q is classically defined by a linear relation including the temperature difference ∆T = T 1 -T 2 as follows:

Q = G∆T = 1 R ∆T . (2.44)
From Eq. 2.42 to Eq. 2.44, we can define the thermal resistance as:

1 R = Q ∆T = λA ∆x = A 2 k B T 2 +∞ 0 q(0) • q(t) dt.
(2.45)

For solids, the heat flux expression is as follows [START_REF] Volz | Molecular-dynamics simulation of thermal conductivity of silicon crystals[END_REF]:

q(t) = 1 V N i=1 N j=1,j =i 1 2 r ij • (F ij • v i ) + 1 6 k=1,k =i,j (r ij + r ik ) • (F ijk • V i ) , (2.46)
the first term is the two-body contribution while the second term is the threebody contribution.

Conclusion

NEMD consists in applying a temperature gradient, and measuring the corresponding flux (or, imposing a heat flux and measuring the resulting temperature gradient). This method was proven to be fairly robust in predicting thermal conductivity and even interfacial resistance based on the reproduction of the experimental conditions. However, artefacts will be generated from the artificial changes made on the system dynamics. Also, the definition of local temperature is ambiguous because temperature is an equilibrium quantity while the whole system studied is under global non-equilibrium conditions, as is discussed previously in Chapter 1. Furthermore, the thermal bath imposed to the system introduces size effects of the thermal dynamics, making it difficult to compute molecular or nano-sized systems. Modal quantities such as spectral phonon relaxation times, etc., also remain out of reach because they cannot be specifically excited.

EMD relies on the "natural" fluctuations of the quantities computed in a system that freely evolves with given interaction potentials and boundary conditions. Since there are no artificial force applied on the system, the size effect could be much reduced compared with NEMD simulation. Thus it allows for the treatment of the averaged microscopic quantities that cannot be excited macroscopically. In the following chapter, we will be focused on the validation of the equilibrium temperature difference fluctuations method with several case studies.

Chapter 3

Validation of the temperature difference fluctuations method

In this chapter, the temperature difference fluctuations method, introduced in the previous chapter, is used for calculating thermal contact resistances in three cases: (i) Si/Ge superlattices; (ii) diameter modulated SiC nanowires; (iii) few layer graphenes. The obtained results are validated by EMD and NEMD simulations.

3.1 Interface thermal conductance in Si/Ge superlattices

Introduction

Superlattices are structures composed of alternating layers of two materials that have nearly the same lattice parameter. The most studied superlattices are those composed of two materials with a certain degree of periodicity. Semiconductor superlattices have optical, electronic, and thermal properties that vary significantly from those of the bulk constituent materials [START_REF] Daly | Molecular dy-References 123 namics calculation of the thermal conductivity of superlattices[END_REF].

Their low thermal conductivity makes superlattices promising materials for applications in devices such as semiconductor lasers, optical data-storage media, thermoelectric, and thermomechanic devices, where the thermal conductivity characteristics is important for the device operation [START_REF] Termentzidis | Nonequilibrium molecular dynamics simulation of the in-plane thermal conductivity of superlattices with rough interfaces[END_REF].

SiGe superlattices offer the possibility of efficient thermoelectric generation from materials which are compatible with existing Si-based information and communication technologies. The reduced thermal conductivity from alloy scattering leads to the possibility of improved thermoelectric properties with respect to bulk elemental or compound semiconductors. The thermal conductivity and the thermal resistance across the Si-Ge interface have been studied previously, by both experiments [START_REF] Lee | Thermal conductivity of si/ge superlattices[END_REF][START_REF] Borca-Tasciuc | Thermal conductivity of symmetrically strained si/ge superlattices[END_REF]; Aksamija and

INTERFACE THERMAL CONDUCTANCE IN SI/GE SUPERLATTICES

Knezevic 2011) and MD simulations [START_REF] Termentzidis | Nonequilibrium molecular dynamics simulation of the in-plane thermal conductivity of superlattices with rough interfaces[END_REF]Chalopin et al. 2012). These results could provide a relevant benchmark to validate the thermal resistances that are calculated from the temperature difference fluctuations method.

Molecular dynamics setup

The three-body Stillinger-Weber interatomic potential [START_REF] Stillinger | Computer simulation of local order in condensed phases of silicon[END_REF] was implemented to describe the Si-Si and Ge-Ge covalent interactions. Potential parameters have also been combined according to the mixing rules described by Ref. [START_REF] Zi | Modification of stillinger-weber potentials for si and ge[END_REF][START_REF] Ethier | Epitaxial growth of si1-xgex on si(100)2 ×1: A molecular-dynamics study[END_REF] to model the interfacial Si-Ge interactions. The time step was set to 1 fs. The structure was relaxed using a steepest descent energy minimization algorithm (sd). The simulation first ran in the isothermal-isobaric (NPT) ensemble with a temperature of 300 K and a pressure of 0 bar for 400 ps, then in the canonical ensemble (NVT) for another 200 ps to equilibrate the system. Atomic trajectories were then calculated in the microcanonical ensemble (NVE). Thermal equilibrium has been achieved within 500 ps. Thermal resistance simulations have been operated on Si/Ge superlattices having a geometrical configuration similar to that presented in Fig. 3.2. Heat conduction mechanisms were studied in the Z direction. Periodic boundary conditions were applied in the three directions, which yields a Si/Ge interface of infinite cross section in the plane perpendicular to direction Z and involves a periodicity L P = N (aGe + aSi) in the same direction, where aSi/Ge stands for the lattice parameters and N is the number of unit cells. We included at least two contacts in the periodic supercell. We assume here a perfect interface between Si and Ge in which there is a perfect lattice matching. This is often achieved through the introduction of coherency strain in one or both of the phases. The lattice parameter difference of 4 % between Si (0.543 nm) and Ge (0.565 nm) has been matched to 0.554 nm in the direction perpendicular to the interface. This relaxation procedure has already been presented in various publications [START_REF] Landry | Thermal boundary resistance predictions from molecular dynamics simulations and theoretical calculations[END_REF]. The resulting lattice parameters in the direction parallel to the contact correspond, respectively, to 0.531 and 0.573 nm for Si and Ge phases.

Interfacial thermal conductance of Si/Ge superlattice

The cut-off radius for Si-Ge interaction was set as 3.85 Å, and only the interacting atoms across the interface were included to calculate the temperature, as is shown in Fig. 3.3. The simulation temperature was set to 300 K. The autocorrelation fuction (ACF) of temperature difference fluctuations and the corresponding integration time (integral of the ACFs vs time) of a 3×3× 6 nm Si/Ge superlattice at 300K is presented in Fig. 3.4. The Kapitza resistance, which is proportional to the integration time, could thus be calculated. The interface thermal conductance is the inverse of Kapitza resistance. The thermal conductances of Si/Ge superlattices with period thicknesses ranging from 6 to 60 nm are presented in Fig. 3.5. The results show a satisfying agreement between the temperature difference fluctuation method and the exchanged power fluctuations method. Fig. 3.6 shows the thermal interface conductance of superlattices over a tem-perature ranging from 300 K to 900 K for a period of L P = 10 nm. It could be seen that the temperature does not affect the conductance of the interface and clearly saturates around the value of 6.6×10 8 W/m 2 K -1 .

3.1.4 Cross-plane thermal conductivity of Si/Ge superlattice

In Ref. (Chalopin et al. 2012), an effective cross-plane thermal conductivity based on interface conductance was compared with Green-Kubo approach and obtained consistent values. In this section, to further validate the interface conductance G obtained from the temperature difference fluctuations method, the effective cross-plane thermal conductivity λ was evaluated from λ = G × L p and compared with the data calculated from Green-Kubo formula:

λ = V 3k B T 2 +∞ 0 q(0) • q(t) dt, (3.1)
The brackets denote the heat flux autocorrelation function, k B is the Boltzmann constant, V is the system volume and T is the equilibrium temperature. We have implemented Green-Kubo thermal conductivity calculations for period thickness ranging from 6 to 60 nm using the same equilibrium temperature as the thermal conductance calculations (300 K). Twenty ensemble averages have been considered for each thermal conductivity point. The results are shown in Fig. 3.7. Thermal conductivity of Si/Ge superlattices is fully driven by the interfaces' contribution in the concerned period length range, indicating that the internal phonon-phonon scattering is negligible. The agreement in the thermal conductivity values calculated from Green-Kubo formula and G × L p validates the relevance of the temperature difference fluctuations method.

Conclusion

Interface thermal conductances of Si/Ge superlattices with different period thicknesses were calculated from the temperature difference fluctuations method.

We observed that the thermal conductance is not dependent on the studied period thickness except for the ones with periods thickness lower than 10 nm.

Further confrontations with Green-Kubo cross-plane thermal conductivity calculations demonstrated that the thermal conductivity of Si/Ge superlattices can be deduced from interface conductances in the concerned period thickness range. Thus, thermal conductivity increases with increasing lattice periods as long as the dominant phonon mean free path involved remains larger than the superlattice period. This indicates that interfacial scattering in perfect superlattices remains the predominant source of scattering in very large phonon mean free path materials. The calculated thermal conductance is not dependent on the temperature. These data are consistent with the results obtained from the exchanged power fluctuations method.

3.2 Thermal constriction resistance of diameter modulated SiC nanowires

Introduction

Silicon carbide (SiC) is an important semiconductor with a wide electronic band gap, and has attracted early interest because of its high mechanical strength, high chemical stability, and high thermal conductivity, which make it an interesting candidate for microelectronic applications. It is used in nanoelectronics and its outstanding hardness makes SiC suitable for devices working in extreme environments. Novel hydrogen storage devices and a series of applications in catalysis, biomedical, and optics are some of the uses for this material [START_REF] Wu | Tuning the morphologies of sic nanowires via the control of growth temperature, and their photoluminescence properties[END_REF].

SiC sustains a variety of stable nanostructures, like nanowires, nanotubes or nanocages, nanobelts, nanorods, etc. [START_REF] Lyver | Lattice thermal conductivity in sic nanotubes, nanowires and nanofilaments: A molecular dynamics study[END_REF][START_REF] Wang | Periodically twinned sic nanowires[END_REF][START_REF] Zhou | Simple approach to beta-sic nanowires: Synthesis, optical, and electrical properties[END_REF]. SiC also exhibits a large number of lattice structure and stacking fault polytypes. This pronounced polytypism counts more than 200 crystalline modifications [START_REF] Zywietz | Influence of polytypism on thermal properties of silicon carbide[END_REF]. The synthesis of nanowire superlattices (NWSLs) made of a variety of dissimilar or the superposition of different polytypes of the same material [START_REF] Wu | Tuning the morphologies of sic nanowires via the control of growth temperature, and their photoluminescence properties[END_REF][START_REF] Wang | Periodically twinned sic nanowires[END_REF][START_REF] Choi | Growth and modulation of silicon carbide nanowires[END_REF]) have both been suggested in literature in this field. Chemical vapor deposition, arc discharge, carbon nanotubes in vapor-solid or solid-solid reactions have been used to fabricate SiC nanowires and nanostructured materials [START_REF] Shen | Silicon carbide hollow nanospheres, nanowires and coaxial nanowires[END_REF][START_REF] Zhang | Silicon carbide nanosprings[END_REF]. It has been shown that such nanostructures exhibit different electronic, optical, and transport properties from those displayed by plain nanowires [START_REF] Mizuno | Acoustic phonon modes and dispersion relations of nanowire superlattices[END_REF].

Phonon confinement observed in semiconductor nanowires modifies the phonon velocities, while surface and interfacial scattering increases the thermal resistance [START_REF] Dames | Theoretical phonon thermal conductivity of si/ge superlattice nanowires[END_REF]. Measuring temperatures and heat fluxes with a nanometer scale resolution is difficult, and such measurements are most often carried out using atomic simulations, which may help in understanding the microscopic phenomena at play [START_REF] Volz | Clamped nanowire thermal conductivity based on phonon transport equation[END_REF]. Molecular dynamics (MD) simulations allows the prediction of the thermal transport properties of semiconductors and insulators because in these materials, energy is mainly carried by phonons. Among the MD simulation approaches to calculate the thermal conductivity, there are two main methods: the nonequilibrium (NEMD) method based on forcing a temperature gradient on the system and the equilibrium (EMD) method within the Green-Kubo approach [START_REF] Merabia | Thermal conductance at the in-References 127 terface between crystals using equilibrium and nonequilibrium molecular dynamics[END_REF].

In this section, we shall study the homogeneous nanowire superlattices of SiC with diameter modulation. We show that the band-structure engineering based on the significant geometrical modulation of the nanowires is responsible for the reduction of thermal constriction resistance. The constriction resistance is related to the confinement of low-frequency modes, as shown by phonon density of states calculations.

Molecular dynamics setup

The LAMMPS code is used for MD simulations and the interatomic forces are calculated thanks to the Tersoff potential [START_REF] Tersoff | New empirical model for the structural properties of silicon[END_REF]). Free boundary conditions are used in the x and y directions, while periodic boundary condition is applied in the z direction. The diameter modulated 2H nanowire considered in our study is depicted in Fig. 3.8.

We consider the diameter modulation of the 2H SiC nanowire over just one period (and therefore just one interface) for the sake of computational time.

The thermal resistance of a diameter modulated nanowire is the sum of the internal thermal resistance of each constant cross section of the nanowire and the constriction resistance near the interface. The constriction resistance differs from the Kapitza resistance, which is the thermal resistance at the interface between two different solids. It is calculated from the temperature difference fluctuations method. The final autocorrelation functions were derived from the average of twenty trajectories with different initial velocities.

Thermal constriction resistance of diameter modulated SiC nanowires

The constriction resistance for diameter modulated nanowires is calculated with the EMD method and found to be equal to 2.37×10 -10 m 2 K/W, while the internal thermal resistance for the small and large cross sections are equal to 1.21 and 1.30×10 -10 m 2 K/W, respectively. The NEMD method is based on the calculation of temperature difference at the two sides of the interfaces and the knowledge of the heat flux q in such a way that R = ∆T /q. With the NEMD method, the average internal thermal resistance for 2H SiC nanowires is equal to 0.6×10 -10 m 2 K/W. The diameter modulated nanowires show a constriction resistance of 0.9×10 -10 m 2 K/W. One possible explanation for the discrepancies between the EMD and NEMD values in this study might be as follows: only one period of the 2H SiC nanowire was considered in the NEMD simulation, and the imposed thermal bath will introduce size effects of the thermal dynamics to such a small-sized system. On the other hand, the EMD simulation allows for the treatment of systems with small sizes, as is discussed in Chapter 2. The two methods even if they do not

THERMAL CONSTRICTION RESISTANCE OF DIAMETER MODULATED SIC NANOWIRES

give the same values, provide the same trends.

The diameter modulaiton decreases the thermal conductivity of the SiC nanowire below the one of a nanowire having the smallest diameter [START_REF] Termentzidis | Modulated sic nanowires: Molecular dynamics study of their thermal properties[END_REF]. This reduction of the thermal conductivity can be attributed to the reduction of the group velocity and the correlation between the wire dimensions and the phonon mean free path, which is around 2 nm here.

Phonon density of states.

The phonon density of states (DOS) for the case of the diameter modulated 2H nanowire is studied with cross sections varying from 8.81 to 15.67 nm 2 . The phonon density of states of the particle is numerically obtained by decomposing the time correlation function of the atomic velocities into the Fourier space, as is expressed in Eq. 3.2.

P (ω) = 1 k B T i ∞ 0 m i v i (0) • v i (t) e -iωt dt, (3.2)
Where ω denotes to the angular frequency, m i is the atomic mass. The quantity in brackets represents the velocity autocorrelation function. Velocities of the atoms were calculated from the integration of the second Newton's law at 300 K from 100 ps. The final DOSs were plotted from the average of three trajectories. The derivation of the DOS expression could be found in Appendix A.

Figure 3.9: Density of states for four distinct groups (see Fig. 3.10). The black and red lines correspond to the groups of atoms far from the surface for small (8.81 nm 2 ) and large (15.67 nm 2 ) cross sections. The green line refers to the DOS in the small section and the first layer, which is the closest to the interface, while the blue one is the DOS in the large section and in the first layer, which is closest to the interface.

In Fig. 3.9, the partial density of states is calculated for the small and large cross sections far from the interfaces and for the first adjacent layers on both sides of the interfaces (Fig. 3.10). We focused on three areas: (a) For intermediate frequencies between 25 and 30 THz where there is a decrease in the density of states of the adjacent layers in comparison to the ones of the constant section nanowires; note that this main peak is generated by the SiC bond; (b) for small frequencies between 0 and 10 THz [Fig. 3.11(a)] where confined modes appear in the first layer of the large cross-section and as peaks related to the stationary phonons; (c) for frequencies ranging from 30 to 37 THz where the additional local maximum of the DOS is attributed to the interfacial and surface phonon modes (Fig. 3.11(b)). The DOS at the first layer in the large cross-section exhibits additional lowfrequency peaks [see Fig. 3.11(b) for more details]. These low-frequency peaks are only shown at the cross-section layers, indicating phonon confinement. Similar phenomena are also observed in glasses [START_REF] Zorn | Boson peak in confined disordered systems[END_REF][START_REF] Jain | Influence of confinement on the vibrational density of states and the boson peak in a polymer glass[END_REF]. We identified these frequencies as being very close to frequencies of the blocked transverse and longitudinal phonons. These phonons cannot pass from the small cross section to the larger or vice versa, as they correspond to unavailable frequencies. The phonon interactions might be altered due to this confinement and as a result a higher Kapitza resistance was found for diameter modulated nanowires.

Conclusions

We have defined the thermal resistance of diameter modulated SiC nanowires as a constriction resistance and have shown that it is worth twice the internal thermal resistance of the constant diameter nanowire. The density of states study shows that phonon modification is related to the confinement of the modes and additional surface modes appear with the diameter modulation. We

Chapter 3 -Validation of the temperature difference fluctuations method 51 believe that a certain number of phonons are trapped in the large cross sections and this might be a helpful way of achieving the objective of further reducing the thermal conductivity of nanowires for application with thermoelectric materials.

A combined modulation of diameter and lattice structure might even reduce further the thermal conductivity of nanowires.

Thermal resistance between few layer graphenes

Introduction

Graphene and few layer graphene (FLG) were proposed as materials for heat removal owing to their extremely high in-plane thermal conductivity [START_REF] Balandin | Superior thermal conductivity of single-layer graphene[END_REF][START_REF] Geim | The rise of graphene[END_REF][START_REF] Saito | Ballistic thermal conductance of a graphene sheet[END_REF]). They were used to incorporate chips to remove heat from localized hot spots, such as in silicon-oninsulator (SOI) integrated circuits [START_REF] Subrina | Heat removal in silicon-on-insulator integrated circuits with graphene lateral heat spreaders[END_REF]) and high-power density electronic devices [START_REF] Yan | Graphene quilts for thermal management of high-power gan transistors[END_REF]. Compared with single layer graphene, FLG is more promising from the practical applications' point of view, because its thermal conductivity is less subject to the deterioration due to extrinsic effects such as substrate impurities (Koh et al. 2010a). However, despite progress, the production of large area and high-quality FLG with a strictly targeted number of layers still remains a challenge. A mixture of FLG with different layer numbers is usually obtained during its synthesis on a substrate. The heat removal efficiency of such FLG structure could then be largely affected by the variation of the FLG thickness, and the corresponding variation of the inter-FLG resistance.

The inter-layer resistance in FLG is an unknown intrinsic property that is likely to affect the performance of FLG-based heat removal structures. For example, more in-plane heat flux could be allowed through cross-section from thinner to thicker FLGs, and a lower inter-layer resistance would be more favorable for heat flux to pass across different layers. However, similar to the weak intertube thermal conduction of carbon nanotubes [START_REF] Bozlar | Carbon nanotube microarchitectures for enhanced thermal conduction References at ultra low mass fraction in polymer composites[END_REF][START_REF] Gordiz | Thermal rectification in multi-walled carbon nanotubes: A molecular dynamics study[END_REF], the inter-FLG thermal conduction of FLG is much weaker than that in the basal plane direction. Understanding and controlling the inter-FLG resistance is, therefore, crucial to the development and performance of FLG-based films in thermal management. In this section, a series of atomistic simulations, aimed at assessing the geometry effects, are performed to systematically study the inter-FLG resistances with a total layer number equals to 2, 3, 4, 5, 6, 10 and 20, respectively. We then report the resistances as a function of different temperatures. The obtained results are compared with that from NEMD simulations and Debye model calculations.

Molecular dynamics setup

Equilibrium molecular dynamics (EMD) simulation was used to investigate the thermal resistance between FLGs. Graphene sheets of 4× 4 nm with periodic boundary conditions applied in the x and y directions were chosen in this study.

The FLG-FLG systems are illustrated in Fig. 3.13. Taking 2,3-layer graphene as an example, the two layers on the left were considered as one FLG subsystem and the remaining three layers were considered as another FLG subsystem. All simulations were carried out by using the LAMMPS package. Adaptive intermolecular reactive empirical bond order (AIREBO) potential [START_REF] Stuart | A reactive potential for hydrocarbons with intermolecular interactions[END_REF] was used to simulate the C-C interactions. AIREBO potential describes the bonded C-C interactions (r < 2 Å, where r is the distance between two atoms) with the reactive empirical bond order (REBO) potential of Brenner [START_REF] Brenner | Empirical potential for hydrocarbons for use in simulating the chemical vapor-deposition of diamond films[END_REF]) and uses Lennard-Jones (LJ) potential for the non-bonded Van der Waals C-C interactions (2 Å< r < cutoff radius). In the current study, the

53

LJ potential is expressed by:

Φ(r) = 4ε σ r 12 - σ r 6 , (3.3)
where ε refers to the depth of the potential well and σ to the equilibrium distance at which the inter-particle potential is zero. The parameters ε=2.39 meV and σ=3.414 Å were implemented [START_REF] Girifalco | Carbon nanotubes, buckyballs, ropes, and a universal graphitic potential[END_REF]. The simulations were conducted in the microcanonical ensemble (NVE) with a time step of 1 fs. The desired temperatures of the system were set by explicitly rescaling the velocities of the atoms for 100 ps. Then the thermal bath was removed and the system was relaxed during one million time steps. Another one million time steps were considered thereafter to compute the temperatures of both subsystems.

Inter-FLG thermal resistance as a function of layer number

The ACFs and the integration time (integral of the ACFs vs time) of the 3,3layer graphene model at 300K are presented in Fig. 3.14(a). The inter-FLG 54 3.3. THERMAL RESISTANCE BETWEEN FEW LAYER GRAPHENES thermal resistance, which is proportional to the integration time, could thus be calculated. equilibrium MD (NEMD) simulations [START_REF] Wei | Interfacial thermal resistance in multilayer graphene structures[END_REF]) and Debye model calculations [START_REF] Sun | Graphite c-axis thermal conductivity[END_REF]). Note that the data from these references are the averaged cross-plane thermal resistance of multi-layer graphene, which are derived from the corresponding thermal conductivities k c , according to the following equation:

(N -1) × R = N × L k c , (3.4)
where L is the distance between neighboring layers (i.e. 0.34 nm) and N is the layer number.

As demonstrated in Fig. 3.14(b), EMD simulations give the value of R(1,1) as high as 1.1×10 -8 m 2 KW -1 , which corresponds to an ultra low effective thermal conductivity of 0.062 W/mK in 1,1-layer graphene. The calculated resistances from different methods agree in that the interfacial resistance decreases as the number of FLG increases. The interfacial resistance finally reaches the graphite limit [START_REF] Taylor | The thermal conductivity of pyrolytic graphite[END_REF] at large layer numbers. The drastic jump in resistance between FLGs with the total layer number from two to four is mainly due to the change in the number of interacting layers. The inter-layer Van der Waals force has a cutoff radius of 10.2 Å in our simulations, and one layer in FLG could interact directly with its third nearest neighbor layer. Within the cutoff radius, a larger layer number brings more interacting layers, which increases inter-layer phonon coupling. That may explain why the interfacial resistance is so sensitive to layer numbers in this layer number range. This trend could also be proved in Ref. [START_REF] Ghosh | Dimensional crossover of thermal transport in few-layer graphene[END_REF] where the in-plane thermal conductivity changes from 2800 to 1300 W/mK as the number of atomic planes in FLG increases from two to four. The decrease of the in-plane thermal conductivity is due to the increased cross-plane coupling of the low-energy phonons. It would also be interesting to note that the in-plane thermal conductivity of bilayer graphene could be largely decreased if the inter-layer sp 3 bonds are introduced (Rajabpour and Vaez Allaei 2012), in which case the inter-layer resistance might be lowered due to the stronger inter-layter interactions. For FLGs with a total layer number larger than four, the interfacial resistance decrease is mainly due to the size effect along the cross-plane direction. Although the phonon mean free path of the multilayer graphene along the cross-plane direction is much smaller than its in-plane counterpart, size effects along the cross-plane direction should not be neglected if this structure is only composed of several layers of graphene [START_REF] Wei | Interfacial thermal resistance in multilayer graphene structures[END_REF]). This has also been proved in Fig. 3.14(b) where the interfacial resistance converges to graphite limit, indicating the size effect in the cross-plane direction is eliminated with a sufficient number of layers. The data from EMD are compared with that from NEMD simulation [START_REF] Wei | Interfacial thermal resistance in multilayer graphene structures[END_REF] and Debye model calculation [START_REF] Sun | Graphite c-axis thermal conductivity[END_REF], and a satisfying agreement is obtained.

3.3. THERMAL RESISTANCE BETWEEN FEW LAYER GRAPHENES

The inter-layer resistance is weaker if the FLG is thicker, and it approaches the value of the interfacial resistance between neighboring layers of graphite. This phenomenon indicates that FLG based heat spreaders, if with uneven layer number, should have larger thickness to achieve the lower inter-FLG resistance.

Thicker FLG indeed is expected to have a better heat spreading efficiency than the thinner one, not only due to its lower interfacial resistance, but also because its larger cross-section allows for more in-plane heat flux to pass through [START_REF] Yan | Graphene quilts for thermal management of high-power gan transistors[END_REF]). Note that it is favorable for improving the heat removal efficiency when heat passes from thinner FLG to thicker one due to the enlarged crosssection, whereas this efficiency will be affected if the heat flux comes from the opposite side. To solve this problem will pose significant fabrication challenges.

Discussions

The average inter-FLG distance L a during the simulation is a critical factor that could affect the calculated value of the inter-FLG thermal resistance. Table 3.2 lists L a of the FLG-FLG systems with different total layer numbers.

Table 3.2: The average inter-FLG distance of the FLG-FLG systems with different total layer numbers at 300 K.

Total layer number of FLG-FLG system Judging from Table 3.2, the average inter-FLG distance L a for 1,1-layer graphene (total layer number equals 2) is slightly larger than the other systems. This is because it is the distance between two surface layers, whose atoms have a smaller coordination number than the bulk atoms and have a higher degree of freedom that would allow them to move towards a energetically more favored position. This may affect the calculated resistance to some extent since the interaction force will be modified due to the change of the equilibrium distance. However, this effect should not be predominant given that the largest variation of L a towards the initial equilibrium inter-layer distance is less than 1.4%. With larger total layer number, L a is decreased and approaches 3.40 Å. In our EMD simulation, no additional conditions were applied to the system to preserve L a , because any external constraints will disturb the thermal equilibrium of the system.

The static roughness of graphene, i.e. the out-of-plane deviation of the carbon atoms is another important factor that needs to be considered. The rippled configuration of graphene is observed by both transmission electron microscopy [START_REF] Meyer | On the roughness of single and bi-layer graphene membranes[END_REF]) and MD simulation (Neek-Amal and Peeters 2011). The ripples in FLG are more likely to be formed in FLG with lower layer number Chapter 3 -Validation of the temperature difference fluctuations method 57 and at higher temperatures. In order to check the effects of the ripples on the calculated resistances, we performed the simulations at temperatures ranging from 300 K to 600 K, and the results are shown in Fig. 3.15. It can be seen that the temperature does not affect the inter-FLG resistance.

The temperature dependence of the thermal boundary resistance at interfaces between two solids is expected to mainly follow the temperature dependence of heat capacity [START_REF] Swartz | Thermal boundary resistance[END_REF]. Our calculations are performed in the classical approximation where the heat capacity is a constant and equals to 3N k B [START_REF] Keer | Principles of the Solid State[END_REF], which explains the absence of temperature dependence of the inter-FLG resistance.

In our simulations, it seems that these ripples do not have obvious impact on the inter-FLG resistance because the density of ripples increases with temperature (Neek-Amal and Peeters 2011), while Fig. 3.15 shows that the inter-layer resistance remained constant at different temperatures. One possible reason might be that the graphene sheets are moving in phase during the simulations. Atoms in one layer are moving in response to the movement of the atoms in the neighboring layers due to the inter-layer interaction, and all the layers exhibit similar bending deformation without significantly changing the inter-layer

CONCLUSIONS

distance, which is consistent with our MD visualization of the time-dependent atom trajectories. Nevertheless, systematic studies are needed to fully understand the effects of static roughness on FLG thermal properties.

Conclusion

The inter-FLG resistances with different layer numbers were calculated from equilibrium molecular dynamics simulations. We found that the inter-FLG resistance is smaller if the FLG is thicker, and it approaches the value of the interfacial resistance between neighboring layers of graphite. Thicker FLG is suggested in thermal devices in order to achieve higher heat spreading efficiency.

In the cases where the layer number of FLG could not be strictly controlled, a FLG mixture with larger thickness might be advantageous owing to the lower inter-layer resistance. The observed properties do not depend on the temperature, guaranteeing stable performances of FLG based structures under different working temperatures. These results provide useful information for future designs of FLG based thermal devices.

Conclusions

The temperature difference fluctuations method is validated by calculating the thermal resistances of Si/Ge superlattices, few-layer graphenes and diameter modulated SiC nanowires. The obtained results were found to be in satisfying agreement with the data from equilibrium MD and NEMD simulations.

(1) The interface thermal conductances of Si/Ge superlattices with different period thicknesses are calculated and the results are consistent with the reported MD values obtained from exchanged power fluctuations method. The cross-plane thermal conductivities of Si/Ge superlattice functions period thickness were also calculated from Green-Kubo formula, which agree with the ones derived from λ = G • L p .

(2) We have defined the thermal resistance of diameter modulated SiC nanowires as a constriction resistance and have shown that it is worth twice the internal thermal resistance of the constant diameter nanowire. The density of states study shows that phonon modification is related to the confinement of the modes and additional surface modes appear with the diameter modulation.

There is a factor of two discrepancy between the data calculated from EMD and NEMD, which may be due to the artifacts generated from NEMD at very small system size.

(3) The inter-FLG resistances with different layer numbers are calculated and are found to be decreasing against the total layer number, before it approaches the graphite limit with sufficient layer number. The observed properties do not These results proved the validity of the temperature difference fluctuations method in calculating thermal conductance (resistance).

Chapter 4

Thermal contact resistance between carbon nanotubes and functionalized polymers

For carbon nanotubes (CNTs) based thermal interface materials, one issue is to ensure a low interface resistance at the CNT-superstrate contact. This resistance is indeed identified as a major contribution to the total resistance of such an interface (Cola et al. 2007a). In this chapter, we report on an azidefunctionalized polymer (HLK5) as efficient thermal glue, which could form covalent bonds (C-N bonds) with carbon nanotube (CNT) through its tail group azide. By using temperature difference fluctuations method, MD simulations predict that the thermal contact resistance (TCR) between CNT and polymer should be decreased by a factor of 6 to 3, when they are connected by covalent bond (with HLK5) rather than by Van der Waals force (with PEMA). In accordance with MD prediction, experimental data show that thermal interface resistance (TIR) with PEMA bonding is 1.50 mm 2 KW -1 , and it is 4.36 times as large as TIR with HLK5 bonding. The mechanism behind this high thermal conductance is discussed from the aspect of lattice vibrations. C-N covalent bond is thus suggested to be an effective bonding in reducing the thermal contact resistance between polymer and carbon based materials. These results provide useful information for future designs of thermal glue for carbon based materials with good thermal conduction and structurally stable interfaces.

Introduction

Carbon nanotubes are molecular-scale tubes of graphitic carbon with outstanding properties. They are among the stiffest and strongest fibres known, and have remarkable electronic properties and many other unique characteristics. The bonding in carbon nanotubes is sp 2 , with each atom joined to three neighbors, as in graphite. The tubes can therefore be considered as rolled-up graphene 4.1(a) shows the formation of a (4, 2) nanotube. By connecting the sites O and A, and B and B', a single wall carbon nanotube will be constructed out of graphene sheet. The vectors OA and OB define the chiral vector C h and the translation vector T of a carbon nanotube, respectively. The rectangle OAB'B gives the unit cell of the one dimensional structure. The chiral vector uniquely determines the nanotube structure and can be expressed on the basis of unit vector of graphene lattice: The thermal conductivity of isolated single wall (SW) and multiwall (MW) carbon nanotubes is about 3000 W/mK or higher at room temperature, which is higher than the of diamond [START_REF] Prasher | Turning carbon nanotubes from exceptional heat conductors into insulators[END_REF]. Therefore CNTs are of both fundamental and practical interest in applying in areas that require high Chapter 4 -Thermal contact resistance between carbon nanotubes and functionalized polymers 63 heat transfer efficiency, such as nano-composites [START_REF] Yu | Enhanced thermal conductivity in a hybrid graphite nanoplatelet, a carbon nanotube filler for epoxy composites[END_REF][START_REF] Huang | Aligned carbon nanotube composite films for thermal management[END_REF][START_REF] Wei | Thermal expansion and diffusion coefficients of carbon nanotube-polymer composites[END_REF][START_REF] Winey | Improving electrical conductivity and thermal properties of polymers by the addition of carbon nanotubes as fillers[END_REF]), chip cooling [START_REF] Shenoy | Minichannels with carbon nanotube structured surfaces for cooling applications[END_REF][START_REF] Kordas | Chip cooling with integrated carbon nanotube microfin architectures[END_REF]), nano-fluids [START_REF] Ding | Heat transfer of aqueous suspensions of carbon nanotubes (cnt nanofluids)[END_REF][START_REF] Xie | Nanofluids containing multiwalled carbon nanotubes and their enhanced thermal conductivities[END_REF][START_REF] Baby | Investigation of thermal and electrical conductivity of graphene based nanofluids[END_REF], as well as thermal interface materials (TIMs) (Xu and Fisher 2006b;[START_REF] Qu | Gecko-foot-mimetic aligned single-walled carbon nanotube dry adhesives with unique electrical and thermal properties[END_REF][START_REF] Tong | Dense vertically aligned multiwalled carbon nanotube arrays as thermal interface materials[END_REF][START_REF] Balandin | Thermal properties of graphene and nanostructured carbon materials[END_REF]. In spite of their high longitudinal thermal conductivity, CNTs could not yet be used to produce bulk materials with equivalently high thermal properties, mainly because thermal contact resistances are predominant between CNTs once the percolation fraction is exceeded but also between CNT and substrate (Xu and Fisher 2006c). Designing a thermal glue able to efficiently connect CNTs would lead revisiting and improving several application fields where nanotubes were inadequately involved to enhance heat fluxes.

-→ C h = n - → a 1 + m - → a 2 , ( 4
In this chapter, we report on a polymer as efficient thermal glue, which has significant low thermal contact resistance with CNTs. To highlight this efficiency, we have measured the thermal contact resistance between CNT and polymer by using a macroscopic configuration where vertically aligned CNTs (VACNTs) which are grown on a silicon substrate and glued to a copper superstrate. The measured thermal resistance between substrate and superstrate directly yields the CNT-glue contact resistance by extrapolation.

This configuration was extensively studied to design thermal interface materials (TIMs) that have attracted much research interests especially in the field of electronic cooling. Ivanov et al. [START_REF] Ivanov | Fast and highly anisotropic thermal transport through vertically aligned carbon nanotube arrays[END_REF]) synthesized VACNTs on Si substrates, and they reported that VACNTs exhibit higher thermal diffusivity than any other macroscopic nanotube assembly. Lin et al. [START_REF] Lin | Molecular phonon couplers at carbon nanotube/substrate interface to enhance interfacial thermal transport[END_REF] synthesized VACNTs on copper substrate and chemically bonded them to a silicon surface. This VACNT/mating substrate interface improved the thermal conductivity of the TIM by almost two orders of magnitude with regard to the bare loose contact between VACNT and the substrate. Xu and Fisher (Xu and Fisher 2006c) showed that the thermal contact resistance (TCR) between silicon wafers and copper with a CNT array interface is about 19.8 mm 2 KW -1 , while the combination of a CNT array and a phase change material produces a minimum resistance of 5.2 mm 2 KW -1 .

Contacts between VACNTs and the superstrate

Dry VACNT-superstrate contact with pressure

The application of pressure on the interface contributes to enhance the thermal transfer, see Fig. 4.2 [START_REF] Cola | Contact mechanics and thermal conductance of carbon nanotube array interfaces[END_REF]). The dense arrays with high mechanical compliance could reduce the thermal contact resistance of CNT array interfaces by increasing the real contact area in the interface. For example, by increasing the pressure from 150 to 450 kPa, the thermal contact resistance of Si/Ti-CNT-Cu decreased from 37 to 20 mm 2 K/W (Xu and Fisher 2006a). Chapter 4 -Thermal contact resistance between carbon nanotubes and functionalized polymers

Metal contact

65

Another method is to use a metal as an interface material between the CNT and the superstrate, see Fig. 4.3. It is a solid-state transfer process of thermocompression bonding VACNTs to metalized semiconductor surfaces, such as those that might be found in either thermal or electrical contact applications. This process allows for the direct transfer at 150 • C of separately grown VACNTs to other silicon based metalized substrates, a temperature within the acceptable thermal budget of MEMS or semiconductor processing [START_REF] Johnson | Thermocompression bonding of vertically aligned carbon nanotube turfs to metalized substrates[END_REF]). It is assumed that the larger CNT-metal TBRs arise from the smaller phononmode overlapping between the CNT and the metals at low frequencies, and the low phonon transmission coefficient at the metal-CNT interface in the intermediate and high frequency range.

Polymer contact

Hung et al. proposed a thin film of polymer as a thermal interface material at the VACNT/superstrate interface, as is shown in Fig. 4.5 (Le Khanh The lowest TIR was measured to be 1.8 mm 2 KW -1 with VACNTs of 10 µm length. The CNTs length dependent thermal resistance measurements suggest that the polymer thin film allows both a reduction in the CNT/superstrate contact resistance and an enhancement of the VACNT intrinsic thermal conductivity. The addition of a polymer layer can increase the effective contact surface. The application of heat and pressure can allow the polymer to flow, and the CNTs to enter the polymer layer. In this manner, the polymer replaces the air film for the shorter CNTs, allowing a reduction of the total thermal resistance by the increase of the number of CNTs connected to the copper. The effective contact surface increase seems more beneficial to the thermal properties than the addition of a low conductivity layer. This phenomenon allows the polymer bonded samples to reach the same resistance values than the metal thermocompression bonded samples.

MD calculation of the thermal contact resistance between a carbon nanotube and polymers

It is reported that the TIR at CNT-superstrate contact is a major contribution to the total resistance of such an interface (Cola et al. 2007a). Giving that the PEMA layer is quite thin (only 600 nm) and it is fully contacted with Cu, the resistance originating from PEMA and that between PEMA and Cu are negligible. As a result, the TIR between CNT and PEMA should be the limiting factor to the Si/CNT/PEMA/Cu TIM thermal conduction. However, due to the chemical inertness of nanotube walls, typically only weak dispersion forces exist between the tube and matrix [START_REF] Shenogin | Role of thermal boundary resistance on the heat flow in carbon-nanotube composites[END_REF]). Further more, the tail group of PEMA, methyl, is an inert group. Thus PEMA could not have strong interaction with CNT, and this significantly baffles the further decrease of TIR.

To further improve the thermal contact, we now propose the molecule HLK5, which is known to form covalent bonds (C-N bonds) with CNT through its tail group azide. The thermal contact resistance between CNT and PEMA and that between CNT and HLK5 are calculated from MD simulations.

Molecular dynamics setup

Equilibrium molecular dynamics (EMD) simulations were conducted first to predict the TCR of CNT/PEMA and CNT/HLK5, respectively. The structure of the models are shown in Fig. 4.6. The thermal resistance R between two systems with temperature difference ∆T was calculated from the temperature difference fluctuation method. EMD simulations were conducted using the LAMMPS software package. The models of CNTs were built by the web-accessible nanotube structure generator "TubeGen Online" (version 3.4). The chirality of the CNT was [START_REF] Ni | Substrate Induced Cross-Plane Thermal Propagative Modes in Few Layer Graphene[END_REF][START_REF] Ni | Substrate Induced Cross-Plane Thermal Propagative Modes in Few Layer Graphene[END_REF] and its length was set to 4.87 nm. The models of polymer chains were established by the free graphics program "MacMolPl" [START_REF] Bode | Macmolplt: A graphical user interface for gamess[END_REF]. Adaptive intermolecular reactive empirical bond order (AIREBO) potential [START_REF] Stuart | A reactive potential for hydrocarbons with intermolecular interactions[END_REF] was used to describe the C-C and C-H interactions. This potential is not only widely used in MD simulations of CNTs [START_REF] Ong | Molecular dynamics simulation of thermal boundary conductance between carbon nanotubes and sio 2[END_REF][START_REF] Grujicic | Atomic-scale computations of References the lattice contribution to thermal conductivity of single-walled carbon nanotubes[END_REF]), but also used to simulate polymers such as polyethylene and polystyrene [START_REF] Henry | High thermal conductivity of single polyethylene chains using molecular dynamics simulations[END_REF][START_REF] Delcorte | Sputtering kilodalton fragments from polymers[END_REF]. Dreiding potential [START_REF] Mayo | Dreiding -a generic force-field for molecular simulations[END_REF]) was chosen to simulate the interactions between C-N, C-O, and C=O. The force between CNT and PEMA was described by Lenard-Jones potential. All the atom pairs were interacting within a neighbor cutoff distance of 10.2 angstroms. Periodic boundary condition was applied in the z direction. The time step was set to 0.5 fs. The structure was relaxed using a steepest descent energy minimization algorithm (sd). The simulation first ran in the isothermal-isobaric (NPT) ensemble with a temperature of 300 K and a pressure of 0 bar for 400 ps, then in the canonical ensemble (NVT) for another 200 ps to equilibrate the system. Atomic trajectories were then calculated in the microcanonical ensemble (NVE) at 300 K. The CNT and polymer temperatures could be computed with MD simulation by using the time dependent kinetic energies of those latter. The final ACFs were derived from the average of 20 trajectories calculated in the same NVE ensemble. As the integration time is proportional to TCR, the thermal resistance and thermal conductance could be calculated, as shown in Table 4.1.

MD CALCULATION OF THE THERMAL CONTACT RESISTANCE BETWEEN A CARBON NANOTUBE AND POLYMERS

HLK5 bonding between VACNTs and Cu

The VACNTs are provided by Smart High Tech AB, Sweden. They were grown by thermal CVD process onto wafers with a 10 nm Al 2 O 3 diffusion barrier layer and iron catalyst. Four series of samples with 10, 30, 50 and 100 µm heights were studied. The samples were then diced into 1cm 2 pieces. The HLK5 thin films (600 nm thick) were spin-coated at 3000 min -1 onto the surface of the Cu square from filtered 1,1,2-trichloroethane solutions. The two surfaces : Si-CNTs and Cu-HLK5 were then brought into contact. The samples were submitted to 66 kPa and heated at 120 • C during 90 minutes. The thickness of HLK5 film was fixed at 0.6 µm and the lengths of the CNTs varied.

Thermal resistance characteration

The total thermal resistance of the multi-layered samples were mesured with a thermal impedancementer (QuickLine T M -50 Thermal Interface analyser, Anter Corporation, Pittsburgh) at room temperature. It allows the measurement of the complete thermal resistance of the sample excluding the resistances from the substrate (Si) and the superstrate (Cu).

Experimental results

Experiments and measurements of the TIR were then conducted using HLK5 bonding between VACNTs and Cu, and the results are compared with PEMA bonding. The measured resistance includes the intrinsic resistance of the TIM and the contact resistance at substrate/material interface. By varying the thickness of the VACNTs, a linear relation between measured resistance and VACNTs length could be obtained. Then the extrapolated resistance at zero CNT length is not affected by the CNT length and yields the CNT-copper resistance. ) is also presented as a comparison. When PEMA layer is substituted by HLK5, the total resistance of the sample is decreased. More than 25% decrease of the total resistance was found in all the considered CNT lengths, indicating that optimized thermal paths are generated with CNT-polymer covalent bonding. The contact resistance values at zero CNT length with PEMA bonding and HLK5 bonding were found to be 1.1 mm 2 KW -1 and 0.34 mm 2 KW -1 , respectively. This contact resistance contains contributions of resistances from the Si-CNT, CNT-polymer, polymer-Cu contact, as well as the intrinsic resistance from polymer. As the polymer layer is quite thin and it fully contacts the copper surface, the resistance originating from polymer and that between polymer and Cu can be considered as negligible. As a result, the thermal interface resistance between CNT and polymer should be the dominating factor to the thermal contact resistance. The TIR with PEMA bonding is 1.50 mm

Comparison between MD simulations and experimental data

In order to compare the simulation data to the experimental ones, we now derive an estimation of the interfacial contact resistance between the CNT and the polymer film. In the following discussion, R M H (R M P ) denotes the TCR between CNT and one monomer of HLK5 (or PEMA); R H (R P ) denotes the total TCR between CNT and HLK5 (or PEMA) in the TIM system; and N H (N P ) denotes the total number of CNT-HLK5 (or CNT-PEMA) contacts in the system. These definitions are illustrated in Chapter 4 -Thermal contact resistance between carbon nanotubes and functionalized polymers

73

According to Table 4.1, R U P is three times as large as R U H . Assuming that U H occupies four carbons in one unit cell of CNT, and U P occupies two (Fig. 4.10), then N H in the TIM is twice as U P . However, this is based on the assumption that all the contacts between CNT and U H are covalent bonds. If only a part of HLK5 chains were connected with CNT by covalent bonds, and the rest of them were interacting with CNT via Van der Waals force, then N H /N P should be smaller than 2. But N H should at least equal to N P . We thus propose a range of approximation for the ratio between the number of polymer molecules and the number of HLK5 ones as follows:

N P N H ≈ 1 → 1 2 . ( 4.3)
The total TIR between CNT and polymer equals to the resistance between CNT and one unit polymer chain divided by the number of CNT-polymer contacts, which yields equation 4.4:

R P R H = R U P /N P R U H /N H = 3 → 6. (4.4)
Thus the resistance between CNT and polymer should consequently be decreased by a factor of 6 to 3 when using HLK5 instead of PEMA, according to the EMD simulation results. This range covers the ratio of R P /R H =4.36

given by experiments (Eq. 4.2). Thus the experimental data matches our EMD prediction in a reasonable range of accuracy.

Phonon density of states of carbon nanotube and polymers

The physical mechanism at play in the thermal contact relies here on the covalent bondings, through which heat is transferred as atomic vibrations. It was reported that the heat transfer from substrate to alkane chains is the ratelimiting step [START_REF] Wang | Ultrafast flash thermal conductance of molecular chains[END_REF]. In this study, the rate is controlled by the coupling strength between phonons in CNT and alkane vibrations. The normalized vibrational density of states (VDOS) of CNT, PEMA and HLK5 are plotted to better understand the TIR decrease from the aspect of lattice vibrations, as shown in Fig. 4.11. According to the VDOS spectrum, there are more modes in HLK5 than in PEMA in the intermediate frequency ranges. The modes with intermediate frequency are most important for thermal transport as they are more abundant and have large group velocities. Similar to Hu's work [START_REF] Hu | Interfacial thermal conductance between silicon and a vertical carbon nanotube[END_REF] where the covalently bonded CNT-Cu exhibited smaller TIR, we attribute the difference between TIR of the CNT-polymer interface with and without covalent bonds is mainly due to the increase in transmission coefficient of intermediate frequency modes due to the presence of C-N covalent bonds. The covalent bonds have a simpler mechanical effect on the thermal contact : the enhancement of the adhesion between CNTs and the superstrate. Indeed, when submitting the assembly to mechanical traction, CNTs can be separated from the growth substrate and fully transferred to the copper substrate : The CNT/polymer adhesion is stronger than the CNT/silicon adhesion. As such transfer cannot be performed on the samples bonded with PEMA, it can be concluded that the covalent bonds play a role in this strong interaction.

Discussions

In order to verify the formation of covalent bonding during the interface closing with HLK5, an Auger Electron Spectrocopy (AES) study was conducted. If there is a chemical reaction between CNT and HLK5, H and N atoms should be found on CNTs. However, because of the pollution of H from environment and the quantity of N is so limited in HLK5, it is difficult to detect the existence of HLK5 on CNT directly. Therefore, another polymer was synthesised : HLK9 (Fig. 4.12). This polymer not only has the same functional group azide (-N 3 ) as HLK5, but also has a lot of Si atoms that can be used as AES probes on the surface of the CNTs. After mixing CNTs and HLK9 at 100 • C for 90 min, the CNTs were taken out and washed extensively with dichloromethane (CH 2 Cl 2 ).

It is assumed that the repeated washing of the CNT allows the solubilization of the physically absorbed HLK9. Therefore, the remaining silicon can be attributed to HLK9 chain covalently bonded to the carbon nanotubes. Samples were washed in dichloromethane extensively until no HLK9 molecules were detected in dichloromethane by infrared spectroscopy. The AES surface survey on CNT and on the substrate (GaAs wafer) are presented in Fig. 4.13. The presence of Si and S from the AES survey of CNTs in Fig. 4.13 indicated that HLK9 chains were connected to CNT through covalent bonding, thus HLK5 chains should have covalent bonding with CNT via the same mechanism. Another key property of HLK5 yielding highly efficient thermal contact consists in its low glass transition temperature (8 • C). This glass transition temperature ensures the VACNTs to fully penetrate into it and form C-N covalent bonding, producing a good VACNTs-HLK5 contact. A good surface-to-surface contact at the interface is indeed known to improve the efficiency of TIM in transferring heat across the interface, and thereby lowering the TIR [START_REF] Shaikh | The effect of a cnt interface on the thermal resistance of contacting surfaces[END_REF]). The local TCR at contact can be approximated as reversely proportional to the atomic vibrational frequency (f) as follows [START_REF] Rajabpour | Thermal boundary resistance from mode energy relaxation times: Case study of argon-like crystals by molecular dynamics[END_REF]:

RK B = 1 f × 1 N 1 + 1 N 2 . (4.5)
Thus by knowing the highest vibrational frequency of the bond, we can estimate the upper limit of thermal conductance across the bond. The highest vibrational frequency for C-N is about 40.5 THz [START_REF] Stuart | Infrared Spectroscopy: Fundamentals and Applications[END_REF], which corresponds to a thermal conductance of 0.28 nWK -1 . This data is quite close to our EMD simulation data where the thermal conductance between CNT and HLK5 is 0.20 nWK -1 , as is shown in Table 4.1.

Nitrogen seems to be the relevant candidate to decrease the TCR because heavier atoms yield lower vibrating frequencies, according to the harmonic oscillator model. Atoms that are heavier than nitrogen, such as oxygen and silicon, are thus expected to have higher thermal resistances with carbon, than the one between C and N.

In the following, the C-N covalent bond in HLK5 was substituted by C-O and C-C (Fig. 4.14) and the corresponding phonon transmission at the CNT-Polymer interfaces were calculated. Dreiding potential [START_REF] Mayo | Dreiding -a generic force-field for molecular simulations[END_REF]) was employed to describe these covalent bonds. We used the following expression for phonon transmission:

T (ω) = β i2ω i∈A,j∈B α,δ∈{x,y,z} k α,δ i,j uα i (ω) u * δ j (ω), (4.6)
where β =1/k B T and the T is the system mean temperature, u i and k α,δ i,j refer to the atom i displacement, and interatomic force constant tensor. The derivation of Eq. 4.6 is shown in Appendix B.

Eq. 4.6 can be estimated by using MD simulations in the microcanocical ensemble (NVE). At each time steps (1 fs), each interacting pair of atoms i, j is identified to properly sample the corresponding velocity fluctuations. The phonon transmission of CNT-Polymer interfaces with C-C, C-N and C-O covalent bonds are shown in Fig. 4.15. The contributions of the functional groups in the polymers to the phonon transmission are identified in the figure .   There is not much difference in the phonon transmission spectra of Polymer-CNT contacts with C-C, C-N and C-O covalent bonds because the atomic mass of carbon, nitrogen and oxygen are very close, and the rest of the polymer structure is the same. But the carbon-nitrogen double bond can provide more thermal paths than the single bond. Besides, for atoms that are lighter than nitrogen, few could be covalently bonded to CNT and form a stable structure. As a result, C-N bond might be a better choice to bond with CNT and decrease the thermal contact resistance.

Conclusion

In summary, TIR between CNT and polymer could be simply and effectively calculated from EMD by equilibrium fluctuations of temperature difference. In agreement to the EMD prediction, smaller TIR between CNT and polymer was achieved for CNT-HLK5 surface than for CNT-PEMA surface. This is because the strong covalent bonds between CNT and HLK5 resulted in a lower barrier for heat carrier. The lowest TIR of Si/CNT/HLK5/Cu TIM was measured as low as 1.40 mm 2 KW -1 with VACNTs of 10 µm length. Finally, we report that using materials with functional groups containing nitrogen (such as azide) to bond with carbon based materials (such as CNTs, graphene, carbon nanowires, et al.) by C-N covalent bond may be the most efficient way to reduce the thermal contact resistance between them.

Chapter 5

Thermal conduction in few layer graphene/SiO 2 systems This chapter includes two sections. In the first section, we report the thermal contact resistance between few layer graphene (FLG) and a SiO 2 substrate. While graphene and FLG are considered as having the highest thermal conductivity in their in-plane directions, our molecular dynamics (MD) simulations however show that those systems are also characterized by a superior thermal contact resistance, which could be largely tuned with the layer number when in contact with a silica substrate. Taking advantage of the resistive interface, in the second section, we show that SiO 2 /FLG superlattices have a thermal conductivity as low as 0.30 W/mK, exhibiting a promising prospect in nanoscale thermal insulation. These findings pave the way for an improved thermal management of nanoscale systems such as thermal barrier coatings and phase change memory materials with atomic-scale super-insulators.

5.1 Modulation of the few layer graphene/ SiO 2 thermal contact resistance

Introduction

Few layer graphene (FLG) was proposed as a material for heat removal owing to its extremely high in-plane thermal conductivity [START_REF] Balandin | Superior thermal conductivity of single-layer graphene[END_REF][START_REF] Geim | The rise of graphene[END_REF]Koh et al. 2010c;[START_REF] Saito | Ballistic thermal conductance of a graphene sheet[END_REF][START_REF] Shahil | Graphene, multilayer graphene nanocomposites as highly efficient thermal interface materials[END_REF]. FLG was for instance envisioned to serve as a heat duct from active nano-devices to heat sink regions [START_REF] Chen | Thermal contact resistance between graphene and silicon dioxide[END_REF]). However, the instability of a freestanding graphene (it has an intrinsic three-dimensional structure or ripples) leads graphene to be used on devices by laying it on a substrate [START_REF] Ao | Density functional theory calculations on graphene/α-sio 2 (0001) interface[END_REF]. For example, the insulating substrates made from silicon dioxide (SiO 2 ) are widely used as a dielectric medium in electronic devices [START_REF] Jang | Thicknessdependent thermal conductivity of encased graphene and ultrathin graphite[END_REF]. As is the case in most molecular systems, the excellent properties of the isolated graphene was found to fade when embedded in a matrix or encased on a substrate. Atomic systems are indeed extremely sensitive to their environment since the mechanisms related to the low number of internal atomic interactions could be easily dominated by the ones due to the interactions with external bodies. The intrinsic thermal properties of atomic-scale compounds therefore become less relevant than the ones generated by their interfaces. Obviously, this issue is critical to the wide electronic-device applications, such as field effect transistors. The change of thermal contact resistance between FLG and SiO 2 against graphene layer number thus deserves to be analyzed. In present, two surfaces (SiC and SiO 2 ) have drawn a great attention [START_REF] Fan | Interaction between graphene and the surface of sio 2[END_REF]. Graphene can be grown epitaxially on SiC substrate by high temperature annealing. Therefore, the interaction between graphene and SiC surface have been largely investigated [START_REF] Mattausch | Ab Initio study of graphene on sic[END_REF][START_REF] Ni | Raman spectroscopy of epitaxial graphene on a sic substrate[END_REF][START_REF] Varchon | Electronic structure of epitaxial graphene layers on sic: Effect of the substrate[END_REF]). The other one is SiO 2 surface which is the thin oxide layer of Si substrate and is widely used as an insulating medium for Si-based device design [START_REF] Zhang | Experimental observation of the quantum hall effect and berry's phase in graphene[END_REF][START_REF] Novoselov | Two-dimensional gas of massless dirac fermions in graphene[END_REF][START_REF] Oostinga | Gate-induced insulating state in bilayer graphene devices[END_REF].

After the graphene is discovered on a SiO 2 surface by optical measurement, various experiments about graphene sheets and nanoribbons [START_REF] Han | Energy bandgap engineering of graphene nanoribbons[END_REF]), including the device fabrication and the measurement of fundamental properties, have been taken on SiO 2 substrate [START_REF] Zhang | Experimental observation of the quantum hall effect and berry's phase in graphene[END_REF][START_REF] Novoselov | Two-dimensional gas of massless dirac fermions in graphene[END_REF][START_REF] Oostinga | Gate-induced insulating state in bilayer graphene devices[END_REF]. The properties of graphene absorbed on SiO 2 are possible to be modulated and this phenomenon has been reported by many experiments [START_REF] Ishigami | Atomic structure of graphene on sio 2[END_REF][START_REF] Berciaud | Probing the intrinsic properties of exfoliated graphene: Raman spectroscopy of free-standing monolayers[END_REF]). The stable gate-controlled conduction has proven that the few-layer graphene can be stable on SiO 2 surface [START_REF] Novoselov | Electric field effect in atomically thin carbon films[END_REF].

The thermal transport between FLG and SiO 2 is of interest for both fundamental reasons and for applications. From the fundamental perspective, thermal transport provides insight into the properties of phonons and their interactions at the FLG-SiO 2 interface. For applications, power dissipation often limits the performance of the devices involving such structures. Of particular importance is the issue of interfacial thermal transport between graphene layers and the substrate which plays a critical role in current saturation of graphene devices.

The thermal contact between graphene and the substrate is known to be impeded by the weak Van der Waals interactions at play. Due to the high interfacial resistance produced by those interactions, a system of carbon sheets appears as a unique means to design heat paths at atomic scale. The contact resistance between graphene and SiO 2 was previously studied by different measurement methods, including Raman and electrical [START_REF] Freitag | Energy dissipation in graphene field-effect transistors[END_REF]), Pump-Probe [START_REF] Mak | Measurement of the thermal conductance of the graphene/sio[sub 2] interface[END_REF], as well as 3-omega methods [START_REF] Chen | Thermal contact resistance between graphene and silicon dioxide[END_REF]. A broad interval of contact resistance values ranging from 5.6×10 -9 m 2 KW -1 to 4×10 -8 m 2 KW -1 were reported. Despite a large number of studies on the graphene and FLG-SiO 2 contact resistances, those latter have always been considered as independent to the layer number [START_REF] Mak | Measurement of the thermal conductance of the graphene/sio[sub 2] interface[END_REF] since no correlation between resistance and layer number could be observed due to experimental inaccuracy. In this section, we provide atomic scale simulation demonstrations of the unequally high thermal contact resistance generated by supported graphene and FLG systems but also of the possibility to tune this resistance with the number of layers.

and σ O-C =3.001 Å. The cutoff distances R c of the LJ potential for the Si-C and O-C interactions were set to 2.5 σ [START_REF] Ong | Molecular dynamics simulation of thermal boundary conductance between carbon nanotubes and sio 2[END_REF]. The sensitivity to cutoff radius was checked with R c equals 4 σ and 5 σ, and no substantial variations in the thermal resistance were observed.

Periodic boundary conditions were applied in the x and y directions. The structure was relaxed using a steepest descent energy minimization algorithm (sd). The simulation first ran in the isothermal-isobaric (NPT) ensemble with a temperature of 300 K and a pressure of 0 bar for 400 ps, then in the canonical ensemble (NVT) for another 200 ps to equilibrate the system. Atomic trajectories were then calculated in the microcanonical ensemble (NVE) at 300 K.

The thermal contact resistances between SiO 2 and FLGs with different layer numbers were calculated from MD simulation with the temperature difference fluctuation method. The system was relaxed during one million time steps and another one million time steps were considered thereafter to compute the temperatures in both subsystems. The final ACFs were derived from the average of more than 20 trajectories.

5. 1.3 Thermal contact resistance between FLG and SiO 2 versus layer number The contact resistance between FLG and SiO 2 as a function of graphene layer number is reported in Fig. 5.3. The resistance between one-layer graphene and silica was calculated as 3.34×10 -8 m 2 KW -1 . The value of the contact resistance was found to decrease against layer number, and converges to 1.70× 10 -8 m 2 KW -1 for the six-layer graphene. Chapter 5 -Thermal conduction in few layer graphene/SiO2 systems 85 than three, N i remained unchanged while the integration time has decreased (see blue dashed line in Fig. 5.4). As a result, the calculated contact resistance is decreasing against layer number.

FLG layers have interactions with their neighboring layers within a cutoff radius of 3 × σ C-C . For the third FLG layer, it could interact directly with its third neighbor, i.e. the sixth layer. The method that we used for calculating resistance involves a temperature difference ∆T = T 1 -T 2 , the 3 layers that interact with SiO 2 are in one subsystem with a mean temperature T 1 and the Si, O atoms that interact with C are in another subsystem with mean temperature T 2 . The sixth layer interacts directly with the third one, so it has a direct effect on T 1 , thus affecting the calculated results. Other graphene sheets except these six layers do not have direct effect on T 1 , and their contribution to the contact resistance is negligible. This may be the reason that the contact resistance saturates when layer number reaches six.

Discussions

Koh et al. (Koh et al. 2010a) measured the thermal conductance of Au/Ti/FLG/ SiO 2 interfaces, and found that the total conductance is independent on the FLG layer number, which seems contradictory to our findings in Fig. 5.3 at first glance. However, considering the low cross-plane thermal conductivity of FLG calculated from MD simulations [START_REF] Wei | Interfacial thermal resistance in multilayer graphene structures[END_REF]) and Debye model [START_REF] Sun | Graphite c-axis thermal conductivity[END_REF]) (less than 0.1 W/mK for 6-layer graphene), the total resistance of FLG will increase with layer number, which is likely to cancel the FLG/SiO 2 contact resistance reduction, thus the measured total conductance is not affected by the graphene layer number.

The contact resistance between FLG and SiO 2 reveals the interface phonon scattering, mostly the scattering between FLG flexural modes and the substrate. It is widely reported that the strong interfacial scattering of flexural modes leads to a decrease of the in-plane thermal conductivity in supported graphene [START_REF] Seol | Two-dimensional phonon transport in supported graphene[END_REF][START_REF] Wang | Thermal transport in suspended and supported few-layer graphene[END_REF][START_REF] Ong | Effect of substrate modes on thermal transport in supported graphene[END_REF]. Moreover, the increasing in-plane thermal conductivity of FLG supported on SiO 2 with layer number measured by Jang et al. [START_REF] Jang | Thicknessdependent thermal conductivity of encased graphene and ultrathin graphite[END_REF]) indicates decreased interface scatterings, which corresponds to a decreasing contact resistance between FLG and SiO 2 . These findings are in accordance with our simulation data. It would also be interesting to note that the substrate-induced disruptions of the outermost graphene layers penetrate a finite distance into the core of the flake, and this characteristic distance is reported to be 2.5 nm (about 7 layers) [START_REF] Jang | Thicknessdependent thermal conductivity of encased graphene and ultrathin graphite[END_REF], which is fairly close to the distance at which the contact resistance converges in this study.

For the Lennard-Jones system, a cutoff radius is often set to R c = 2.5 σ in the literature [START_REF] Hurst | Transport Coefficient Computation Based on Input/output Reduced Order Models[END_REF][START_REF] Thomin | The Impact of Non-equilibrium Microstructure on the Mechanical Response of Polymer Nanocomposites[END_REF][START_REF] Mastny | Melting line of the lennard jones system, infinite size, and full potential[END_REF][START_REF] Nagel | High Performance Computing in Science and Engineering[END_REF][START_REF] Wedekind | Nano-Droplets at Birth Computer Experiments on Gas Phase Nucleation[END_REF], and reasonable results could be obtained. This cut-off radius is not set arbitrarily: at 2.5 σ, the potential energy is merely 1.6% of its value at the minimum (-ǫ) [START_REF] Toxvaerd | Communication: Shifted forces in molecular dynamics[END_REF]. The sensitivity to LJ potential cut-off radius between CNT and SiO 2 was also checked by [START_REF] Ong | Molecular dynamics simulation of thermal boundary conductance between carbon nanotubes and sio 2[END_REF], and they found no significant conductance change if R c increases from 2.5 σ to 4 σ. Nevertheless, the accuracy of the simulated data at certain R c may depend on the system that is studied, thus it is worthwhile to check if longer R c gives consistent values. We investigated the thermal resistance between SiO 2 and single layer graphene using 4 σ and 5 σ, as shown below: It seems that longer cut-off radius than 2.5 σ does not obviously affect the calculated results. With R c = 5 σ, we further calculated 1, 3, 5, 7 layer graphene supported on SiO 2 : The results between R c = 2.5 σ and R c = 5 σ are within error bars.

Conclusion

By using molecular dynamics simulations, the contact resistance between FLG and SiO 2 was found to be decreasing with layer number and converged for six graphene layers. The contact resistance decrease is ascribed to the stronger coupling between FLG and substrate when there are more graphene layers interacting with SiO 2 . The contact resistance between SLG and SiO 2 was found to be 3.34×10 -8 m 2 KW -1 , which is more than one order of magnitude larger than the inter-layer resistance in graphite. The calculated results are within the typical range of the measured contact resistance between graphene and SiO 2 .

5.2 Few layer graphene based superlattices as efficient thermal insulators

Introduction

The weak Van der Vaals interactions between FLG and SiO 2 generate a relatively large contact resistance, compared with that between many solids, as is shown in Fig. 5.7. to the other solid-solid interfaces, SiO 2 -Graphene could provide a large contact resistance without acquiring a very large mass mismatch. The graphene encased by SiO 2 could provide twice the contact resistance due to the doubled interface, and this value is as high as 6.64×10 -8 m 2 KW -1 if the largest contact resistance obtained in our MD simulation is considered. Taking the thermal conductivity of SiO 2 as 1.38 Wm -1 K -1 [START_REF] Kleiner | Thermal conductivity measurements of thin silicon dioxide films in integrated circuits[END_REF], the thickness of SiO 2 that provides the same amount of resistance is 91.6 nm. This means that an encased graphene layer with an effective thickness less than 0.7 nm, i.e. about twice the value of σ Si-C , provides the same thermal resistance as a 91.6 nm thick SiO 2 layer, as is shown in Fig. 5.8. Given that SiO 2 is already a good thermal insulator due to its low thermal conductivity [START_REF] Zhang | Fabrication of thick silicon dioxide layers for thermal isolation[END_REF], graphene systems encased by SiO 2 should behave more like efficient thermal insulators than conductors in the direction perpendicular to the interface. Materials with structures such as the one of SiO 2 /graphene superlattices are therefore proposed to be used as planar thermal insulator while requiring a very reduced volume.

Molecular dynamics setup

FEW LAYER GRAPHENE BASED SUPERLATTICES AS EFFICIENT THERMAL INSULATORS

To verify the projection of SiO 2 /graphene thermal insulator, we calculated the thermal conductivity of SiO 2 /Single layer graphene (SLG) and SiO 2 /FLG superlattice (Fig. 5.9) from the well-known Green-Kubo formula [START_REF] Kubo | Statistical Physics II[END_REF]:

k = V k B T 2 +∞ 0 q(0) • q(t) dt,
(5.1)

the brackets denote the heat flux autocorrelation function, V is the system volume and T is the temperature. Both crystal SiO 2 (c-SiO 2 ) and amorphous SiO 2 (a-SiO 2 ) were studied, and the thermal conductivities of the corresponding bulk SiO 2 were also calculated as benchmarks. For superlattices, the thermal conductivities were only calculated along the direction perpendicular to the interfaces, by using the cross-plane component of the heat flux.

It seems that the Tersoff potential is not sufficient for the thermal conductivity calculation since the long-range Coulombic interaction contribution to the thermal conductivity should not be neglected. As a result, BKS potential ( [START_REF] Vanbeest | Force-fields for silicas and aluminophosphates based on abinitio calculations[END_REF]) was used to simulate SiO 2 to include the columbic force in the thermal conductivity calculation. Amorphous SiO 2 was produced by annealing the crystalline structure at 6000 K [START_REF] Ong | Molecular dynamics simulation of thermal boundary conductance between carbon nanotubes and sio 2[END_REF] and then cooling down rapidly to room temperature. Periodic boundary conditions were applied to x, y and z directions to generate the superlattice structure. The period thickness L and the layer number of FLG were varied to attain the lowest thermal conductivity.

Thermal conductivities of FLG/SiO 2 superlattices

As is shown in Fig. 5.10, the thermal conductivities of both c-SiO 2 /SLG and a-SiO 2 /FLG superlattice increase with increasing period thickness and approach the corresponding bulk silica thermal conductivity. The thermal conductivity increase with period thickness is due to the increased silica proportion or reduced interface density. Amorphous SiO 2 /SLG superlattice has generally lower thermal conductivity than crystal SiO 2 /SLG, due to the lack of periodicity of the crystalline structure and the normal modes are no longer plane-waves in such a way that phonons are damped as in anharmonic crystals [START_REF] Hunklinger | Phonons in amorphous materials[END_REF]. Unlike strong interactions between two solids such as covalent bonds (Ni et al. 2012) and nanojunctions [START_REF] Jalabert | Ballistic thermal conductance of a lab-in-a-tem made si nanojunction[END_REF]) which increase the thermal conductance, the weak Van der Waals interactions between graphene sheet and SiO 2 generate large thermal resistance and reduce the thermal conductivity of the corresponding superlattice. The lowest thermal conductivity appeared at the shortest period thickness, where SiO 2 has only one unit cell with a length of 0.6 nm. From Fig. 5.3, the contact resistance between SLG and SiO 2 is about 3.34×10 -8 m 2 K/W. If this is repeated every 1 nm with two such interfaces, then the corresponding thermal conductivity would be roughly approximated as 0.017 W/mK. The corresponding value reported in Fig. 5.10 is around 0.35 W/mK, larger by a factor of 20 and suggesting that the behavior of these samples is more complicated than simple addition of series resistances. This phenomenon might be due to the fact that the period thicknesses are within the characteristic length in which the contributions of the evanescent modes arise. Here, we define an evanescent mode as a nonpropagative surface mode that exists at the interface and that is locally confined over a typical thickness. Similar results were observed in Si/Ge superlattices (Chalopin et al. 2012), where higher thermal conductivities were found at very small period lengths.

Thermal conductivity was then calculated from models with a fixed a-SiO 2 length (0.6 nm) but with different FLG layer numbers. According to Fig. 5.11, the thermal conductivity of a-SiO 2 /FLG superlattice slightly decreases with the layer number, and saturates to 0.30 W/mK at four graphene layers. We also investigated the a-SiO 2 /5-layer graphene superlattice at different temperatures up to 1000 K, and found that the thermal conductivity was not temperature dependent, as shown in Fig. 5.12. The lowest thermal conduc-tivity of 0.30 W/mK found in the a-SiO 2 /FLG superlattice is close to that of a suspended 10-layer graphene in the cross-plane direction [START_REF] Wei | Interfacial thermal resistance in multilayer graphene structures[END_REF][START_REF] Sun | Graphite c-axis thermal conductivity[END_REF]. The low thermal conductivity of a-SiO 2 /FLG superlattice is however conserved whatever its total thickness, while the one of thicker FLG will significantly increase [START_REF] Wei | Interfacial thermal resistance in multilayer graphene structures[END_REF][START_REF] Sun | Graphite c-axis thermal conductivity[END_REF]. Another advantage of a-SiO 2 /FLG superlattice is that it could endure high temperatures of more than 1800 K, while most polymers with equivalent thermal conductivities melt below 400 K, which is the typical working temperature of today's thermal barrier coatings and phase change memory materials.

Conclusion

To summarize, the thermal conductivity of SiO 2 /SLG and SiO 2 /FLG superlattices were calculated from the Green-Kubo formula. The results show that SiO 2 /FLG superlattices have a low thermal conductivity of around 0.30 W/mK, indicating a large potential as nano-scale insulators. The a-SiO 2 /FLG superlattice could endure high temperatures of more than 1800 K, while most insulators such as polymers with equivalent thermal conductivities will fail below 400 K, which is the typical working temperature of today's thermal barrier coatings and phase change memory materials. The low thermal conductivity of a-SiO 2 /FLG superlattice is independent on temperature up to 1000 K, guaranteeing stable performances of a-SiO 2 /FLG superlattice based structures under different working temperatures.

Conclusion

As a conclusion, by using molecular dynamics simulations, the contact resistance between FLG and SiO 2 was found to be decreasing with layer number and converged for six graphene layers. The contact resistance decrease is ascribed to the stronger coupling between FLG and substrate when there are more graphene layers interacting with SiO 2 . The contact resistance between SLG and SiO 2 was found to be 3.34× 10 -8 m 2 KW -1 , which is more than one order of magnitude larger than the inter-layer resistance in graphite. Owing to such a resistive interface, MD simulations show that SiO 2 /FLG superlattices have a low thermal conductivity of 0.30 W/mK, exhibiting a promising prospect in nano-scale thermal insulation. These findings are promising in thermal barrier coatings and phase change memory materials and provide useful information in the thermal design of small-scale devices.

Chapter 6

Cross-plane thermal properties in suspended and supported few layer graphenes

In this chapter, we report the layer number dependence of the cross-plane thermal resistances (R c ) of suspended and supported few layer graphenes (FLGs) by equilibrium molecular dynamics (EMD) simulations. We show that the existence of a silicon dioxide substrate can significantly decrease the cross-plane resistances of FLGs with low layer numbers, and the effective thermal conductivities were increased accordingly. The Frenkel-Kontorova model was introduced to explain the substrate-induced band gaps in FLG dispersion relation and the corresponding thermal energy transfer. The enhanced thermal conduction in the cross-plane direction is ascribed to the phonon radiation that occurs at the FLG-substrate interface, which re-distributes the FLG in-plane propagating energy to the cross-plane direction and to the substrate. These results provide a better insight in the field of heat transfer across FLG layers and yield useful information to help designing FLG based thermal materials.

Introduction

Since its discovery, graphene has garnered great interest due to its many remarkable physical properties [START_REF] Morozov | Giant intrinsic carrier mobilities in graphene and its bilayer[END_REF][START_REF] Ghosh | Extremely high thermal conductivity of graphene: Prospects for thermal management applications in nanoelectronic circuits[END_REF][START_REF] Nika | Lattice thermal conductivity of graphene flakes: Comparison with bulk graphite[END_REF][START_REF] Dorgan | Mobility and saturation velocity in graphene on sio 2[END_REF], among which its high in-plane thermal conductivity has been suggested as a key advantage for applications in microelectronics and thermal management [START_REF] Jang | Thicknessdependent thermal conductivity of encased graphene and ultrathin graphite[END_REF]). The exotic properties of monolayer graphene, however, are not always beneficial for the realization of graphene devices or interconnects, and in many cases, few layer graphenes (FLGs) have more favorable properties for practical applications. For example, transverse electric fields in bi-layer and few layer graphene can open band gaps up to 0.2 eV, which is crucial for the operation of fieldeffect transistors. In addition, FLG provides a better transparent conductive electrode due to the lower sheet resistance, and is less susceptible to the effects of substrate impurities due to interlayer screening. The recent progress in fabricating large area and high quality graphene sheets (Li et al. 2009;[START_REF] Li | Highly conducting graphene sheets and langmuir-blodgett films[END_REF][START_REF] Reina | Large area, few-layer graphene films on arbitrary substrates by chemical vapor deposition[END_REF]) lead to a more promising FLG application prospect. Fig. 6.1 shows one of its applications in field effect transistors as heat spreaders. The knowledge of how thermal properties evolve from monolayer graphene to graphite will facilitate the development of graphene devices (Koh et al. 2010b).

In contrary to its in-plane counterpart, however, the cross-plane thermal conductivity of FLG is supposed to be much lower due to the weak Van der Waals interaction between layers. This strongly limits the heat dissipation efficiency of FLG based thermal materials. The basal plane conduction property is also notably depending on the weak Van der Waals cross-plane interactions that reduce the inter-layer scattering [START_REF] Nika | Theoretical description of thermal transport in graphene: The issues of phonon cut-off frequencies and polarization branches[END_REF]. As a result, the knowledge of the cross-plane thermal properties of FLG becomes essential in understanding the thermal mechanims at play in FLGs. Generally, there are two configurations of FLG, either supported on a dielectric substrate such as silicon dioxide, or suspended over a hole in the substrate. Several works have studied the cross-plane thermal properties of suspended FLG, including a non-equilibrium molecular dynamics (NEMD) study [START_REF] Wei | Interfacial thermal resistance in multilayer graphene structures[END_REF]) and a Debye model calculation [START_REF] Sun | Graphite c-axis thermal conductivity[END_REF]. They found decreased cross-layer thermal resistances of FLGs against the layer number. For supported FLGs, despite a large number of studies on the in-plane thermal properties [START_REF] Lindsay | Flexural phonons and thermal transport in multilayer graphene and graphite[END_REF][START_REF] Wang | Thermal transport in suspended and supported few-layer graphene[END_REF][START_REF] Ong | Effect of substrate modes on thermal transport in supported graphene[END_REF], how the heat propagates across the layers is rarely been reported.

In the second section of this Chapter, we present the calculation of the crossplane thermal resistances R c of FLGs in both supported and suspended configurations by equilibrium MD simulations. We show that R c decreases as the FLG layer number increases for the suspended configuration. The existence of a SiO 2 substrate significantly reduces R c of FLGs with small layer numbers. The effective cross-plane thermal conductivity of FLG, approximated from the product of thermal conductance and FLG thickness, is increased due to the substrate effect.

In the third and fourth sections, discussions will be focused on how the substrate affects the FLG cross-plane thermal properties. Phonon density of states and dispersion relations are calculated to better understand the physical mechanisms behind these phenomena. The Frenkel-Kontorova model and sound line are introduced to explain the substrate effects on the flexural modes and the corresponding changes in the energy transfer. In the MD simulations, suspended armchair graphene layers (3×3.7nm) and that supported on a crystalline SiO 2 (3×3.7 × 2 nm) with periodic boundary conditions in x and y directions were chosen in this study (see Fig. 6.3). The MD simulations were conducted using the LAMMPS software package (Plimpton 1995). Aadaptive intermolecular reactive empirical bond order (AIREBO) potential [START_REF] Stuart | A reactive potential for hydrocarbons with intermolecular interactions[END_REF] was used to simulate FLG, which has been widely used to describe the carbon-carbon interactions in graphene and carbon nanotubes [START_REF] Shenoy | Edge-stress-induced warping of graphene sheets and nanoribbons[END_REF]Ni et al. 2012;[START_REF] Zhao | Size and chirality dependent elastic properties of graphene nanoribbons under uniaxial tension[END_REF]). A modified Tersoff potential was used to describe the Si-Si, Si-O, and O-O interactions in silica [START_REF] Munetoh | Interatomic potential for si-o systems using tersoff parameterization[END_REF]). Lennard-Jones (LJ) potential was adopted for interlayer interactions in FLG and the interfacial interactions between FLG and silica, with ǫ The cross-plane thermal resistance R c of FLG based on the Green-Kubo formula was derived in Chapter 2:

1 R = A 2 k B T 2 +∞ 0 q(0) • q(t) dt. (6.1)
This approach has been used for calculating the Kapitza resistance such as for solid/liquid interfaces [START_REF] Barrat | Kapitza resistance at the liquid-solid interface[END_REF] and superlattices (Chalopin et al. 2012;[START_REF] Mcgaughey | Molecular dynamics prediction of the thermal resistance of solid-solid interfaces in superlattices[END_REF], by considering the fluctuating heat fluxes across the interface. In the same manner, it is possible to use the fluctuationdissipation theorem to express the thermal resistance of a thick system as proposed in Eq. 6.1. Thus the averaged inter-layer resistance of a FLG system can be calculated if the total cross-plane heat fluxes q c are taken into account.

We can then estimate the effective cross-plane thermal conductivity λ e as follows:

λ e = L R c , (6.2)
where L is the inter-layer thickness of the FLG, i.e. 0.34 nm. For suspended FLG, there is a drastic R c jump when the layer number increases from two to four, which may be due to the change in the number of interacting layers. The inter-layer Van der Waals force has a cutoff radius of 10.2 Å in our simulations, and one layer in FLG could interact directly with its third nearest neighboring layer. Within the cutoff radius, a larger layer number brings more interacting layers, which increases inter-layer phonon coupling. For FLGs with layer number larger than four, R c decrease is mainly due to the size effect, which could be eliminated with sufficient layers, and R c approaches the graphite limit [the black dashed line [START_REF] Taylor | The thermal conductivity of pyrolytic graphite[END_REF]] , as is shown in Fig. 6.4(b). Unlike in suspended FLG, no significant change is found for R c in supported FLG, indicating that R c is less sensitive to the layer number due to the existence of the substrate. At small layer numbers, the R c in supported FLG is much

Results and discussion

Chapter 6 -Cross-plane thermal properties in suspended and supported few layer graphenes 101 lower than that in suspended configuration. These results suggest that one can reduce R c in suspended FLG either by increasing graphene layers or by adding a substrate. When FLG layer number reaches seven, R c is not affected by the substrate anymore, and both configurations tend to have the same cross-plane resistance. Jang et al. [START_REF] Jang | Thicknessdependent thermal conductivity of encased graphene and ultrathin graphite[END_REF]) reported that the substrate could affect the in-plane thermal conductivity of FLG up to a characteristic distance of 2.5 nm (7 layers), due to the weak Van der Waals coupling between adjacent graphene layers. We found in this work that the substrate-induced disruptions to cross-plane thermal resistance of FLG also vanish at the same distance.

The effective cross-plane thermal conductivity, λ e , of suspended and supported FLG is shown in Fig. 6.5. The λ e of suspended FLG is increasing with layer number, and an ultra low λ e of 0.067 W/mK was obtained in 2-layer graphene.

The effective thermal conductivity λ e of supported FLG was found to be higher than that of the suspended configuration by a factor of two, until the layer number reaches seven. For FLGs with layer number larger than seven, λ e is not affected by the substrate anymore, and both configurations tend to give the same value. The λ e of suspended FLG matches well with the data from a non-equilibrium molecular dynamics (NEMD) simulation by Wei et al [START_REF] Wei | Interfacial thermal resistance in multilayer graphene structures[END_REF]).

Phonon density of states of FLGs

The phonon density of states (DOS) of suspended and supported FLGs are calculated from the time correlation function of the atomic velocities in the Fourier space [START_REF] Domingues | Thermal relaxation times and heat conduction in beta-cristobalite and alpha-quartz silica structures[END_REF]). Fig. 6.6 shows the DOS of suspended FLGs with different layer numbers. Judging from Fig. 6.6, the change of layer number does not cause obvious differences in the DOS of suspended FLGs. But the insert shows that at very low frequencies, the modes in the low-layer number FLGs are thermally populated with a larger probability than in the high-layer number FLGs. The small-layer number FLGs are mainly consist of the surface layers, and the enhancement of the acoustic modes at lower frequency can be directly attributed to the presence of the weakly bonded surface atoms, characterized by a lower coordination (Şopu et al. 2011).

PHONON DISPERSION RELATIONS OF FLGS

According to Fig. 6.7, for different layer numbers, both FLG configurations showed almost the same DOS in the frequency range below 50 THz. However, there is a red-shift of the peak of the optical transverse modes (TO) from supported to suspended FLG, which could be attributed to the strain effects due to the surface stresses. In the case of Van der Waals interaction, the surface stresses are tensile. Removing the substrate will expose the surface layers of FLG and allow a slight lattice expansion of FLG corresponding to an increase in the inter-atomic bond length which indicates a softening of inter-atomic force constant (Şopu et al. 2011). Since the optical modes provide only a minor contribution to specific heat, the DOS could not provide direct information in explaining the different cross-plane thermal behaviors in suspended and supported FLGs.

Phonon dispersion relations of FLGs

To better understand the mechanisms behind the different behaviors of the cross-plane thermal properties exhibited in the different FLG configurations, the FLG phonon dispersion relations are calculated from MD simulations in this section. The sound line and the Frenkel-Kontorova model are then introduced to explain the substrate effects on the dispersion relation of the flexural acoustic modes and the corresponding energy transfer.

Phonon dispersion relation calculation method

The phonon dispersion relation could be calculated from the mode kinetic energy in the reciprocal space. Twenty unit cells of graphene were modeled along the wave vector direction, as shown in Fig. 6.8: For a given polarization, p=x,y,z, the reciprocal space representation of the velocity for atom i at time t is [START_REF] Heino | Dispersion and thermal resistivity in silicon nanofilms by molecular dynamics[END_REF]:

v p i,k (t) = v p i (t)e -ik•r i (t) , (6.3)
the per-atom kinetic energy in the reciprocal space is

E p i,c (k, t) = mv p i,k (t) • v p * i,k (t) 2 . (6.4)
By Fourier transforming the atomic kinetic energy, the eigen-frequencies for a given wave vector k and the corresponding amplitudes of the mode kinetic energy could be found, as is shown in Fig. 6.9.

E p c (k, ω) = m i +∞ -∞ v p i (t)e ikr(t) e iωt dt 2 2
. (6.5)

Figure 6.9: The mode kinetic energy as a function of frequency and wave vector for longitudinal atomic motions of suspended 3-layer graphene.

The transverse and longitudinal branches of the dispersion relation can be obtained by changing the polarization of the velocities, then they are combined together to form the total dispersion relation, as is shown in Fig. 6.10. The spectrum of phonon density of states (DOS) that calculated from all modes (all velocity components, black line) and only from the flexural modes (cross-plane direction velocity, red line) are also shown. The observed dispersion relation of ZA modes in supported FLG is consistent with the model proposed by [START_REF] Kosevich | Long wavelength surface oscillations of a crystal with an adsorbed monolayer[END_REF], where the dispersion relations of a monolayer either strongly or weakly absorbed on a deformable crystal were theoretically studied. In their model, the adsorbed layer, strongly bonded with the substrate, is characterized by surface waves with a Rayleigh-wave polarization, i.e. the surface modes of the absorbed layer moving in phase with the Rayleigh-wave of the substrate. A linear (acoustic) dispersion was found in the small wave vector domain. Similar results were obtained in Ref. [START_REF] Ong | Effect of substrate modes on thermal transport in supported graphene[END_REF], where the strong coupling between the flexural modes of graphene with surface waves of the substrate has been shown

Chapter 6 -Cross-plane thermal properties in suspended and supported few layer graphenes 109 to lead to a hybridized linear dispersion in the small wave vector region. For the monolayer which is weakly coupled to the substrate, the dispersion relations for the modes with Rayleigh polarization and for the gap modes were predicted in the small wave vector region [START_REF] Kosevich | Long wavelength surface oscillations of a crystal with an adsorbed monolayer[END_REF], which is in exact agreement with the present predictions for supported FLGs.

The concept of a sound line can be used to explain the effects of the lowfrequency band gap in the dispersion relation of FLG on the reduced R c , as is shown in Fig. 6.13. below the transverse sound line is white. In the grayed and patterned regions, the gap modes of FLG are propagating modes and radiate phonon energy into the substrate [START_REF] Benchabane | Evidence for complete surface wave band gap in a piezoelectric phononic crystal[END_REF]. The formation of the gaps at zero wave vector brings more ZA modes to the grayed and patterned regions above the sound lines. Thus the in-plane propagating energy of FLG will be partly transferred to the cross-plane direction and to the substrate and will result in an enhanced thermal conduction in the cross-plane direction.

The Frenkel-Kontorova Model

The emergence of the low group velocity branches of the ZA modes appearing in the dispersion relation of supported FLG can also be explained by the Frenkel-Kontorova model [START_REF] Braun | The Frenkel-Kontorova Model, Concepts, Methods, and Applications[END_REF], which describes the dynamics of a linear chain of particles interacting with their nearest neighbors in the presence of an external periodic potential. The potential energy U of the chain consists of two parts, U = U sub + U int , (6.6)

The first term U sub characterizes the interaction of the chain with an external periodic on-site potential, defined in the simplest form,

U sub = ε s 2 n 1 -cos 2πx n a s , (6.7)
where ε s is the potential amplitude, x n is the coordinate of the n-th particle in the chain, and a s is its period. The second term U int in Eq. 6.6 takes into account a linear coupling with the nearest neighbors of the chain, U int = g 2 n (x n+1 -x n -a 0 ) 2 , (6.8)

Chapter 6 -Cross-plane thermal properties in suspended and supported few layer graphenes 111 and it is characterized by the elastic constant g and the lattice constant a 0 . Introducing the dimensionless variables, the Hamiltonian can be written as:

H = n 1 2
dx n dt 2 + (1 -cos x n ) + g 2 (x n+1 -x n -a 0 ) 2 , (6.9)

where a 0 → (2π/a s )a 0 , x n → (2π/a s )x n , t → (2π/a s )(ε s /2m a ) 1/2 t, and g → (a s /2π) 2 g(ε s /2) -1 . The first term is the kinetic energy. From the Hamiltonian (Eq. 6.9), we obtain the corresponding equation of motion of a discrete chain, d 2 x n dt 2 + sin x n -g(x n+1 + x n-1 -2x n ) = 0. (6.10)

We consider an infinite chain with a 0 = a s when the ground state of the chain (i.e., its stationary state that corresponds to the absolute minimum of the potential energy) describes a commensurate array of particles. This means that each minimum of the substrate potential is occupied by one particle, so that the "coverage" parameter θ (defined for a finite chain as the ratio of the number of particles to the number of minima of the substrate potential) is θ = 1. In this case it is convenient to introduce new variables for the particle displacements, u n , defined by the relation (6.12)

The frequency spectrum of the Frenkel-Kontorova chain is hence characterized by a minimum frequency ω min = (2π/a s )(ε s /2m a ) 1/2 .

Coming back to the case of 3-layer graphene supported on SiO 2 , if we take the substrate lattice constant a s as 4.98 Å [START_REF] Brice | The lattice constants of a-quartz[END_REF], and the value of the effective ε s as 35.4 meV, the ω min can be approximated to 0.61 THz. Four carbon atoms are considered to be interacting with the substrate within the length scale of a s . ε s is evaluated by moving the carbon atom in the harmonic region and calculating the corresponding potential energy depth arising from the substrate. The higher frequency at zero wave vector is around 2 THz as shown in Fig. 6.13, which is close to the the value of 55 cm -1 (1.65 THz) for perpendicular graphene-SiO 2 vibration frequency reported by Persson and Ueba (Persson and Ueba 2010), but is three times as large as the calculated ω min from the Frenkel-Kontorova model. Note that there are a series of simplifying assumptions taken in the Frenkel-Kontorova model which may cause this 6.5. CONCLUSION discrepancy, such as (i) the one-dimensional model of the particle motion only along longitudinal direction; (ii) only the harmonic interactions between the chain and the substrate are taken into account; (iii) only the nearest neighbors interactions are involved; (iv) the substrate is assumed to be a rigid (nondeformable) body, in which phonon propagation is impossible. However, our FLG system is a three-dimensional model on a deformable substrate with more degrees of freedom than implemented in the Frenkel-Kontorova model, and the anharmonicity and long-range Van der Waals interactions are also included.

Conclusion

In summary, by using molecular dynamics simulations, we show that the crossplane thermal resistance in suspended FLG decreases as the layer number increases. The corresponding effective thermal conductivity increased when it is supported on a SiO 2 substrate. Introducing a substrate or more graphene layers to the suspended FLG will limit the free vibrations of ZA modes and pose resistance on in-plane phonon transport, but transfer the energy to the cross-plane. Phonon density of states of suspended and supported FLGs are calculated, and a red-shift of the peak of the optical transverse modes (TO) from supported to suspended FLG is found, which could be attributed to the strain effects due to the surface stresses. Phonon dispersion relations show that part of the ZA modes of supported FLG were confined and standing waves were formed. The Frenkel-Kontorova model and sound lines were introduced to explain the substrate-induced band gaps in the FLG dispersion relation and the corresponding effects on FLG cross-plane thermal properties. The enhanced thermal conduction in the cross-plane direction is ascribed to the additional cross-plane heat fluxes generated by FLG-substrate phonon radiation. The obtained results are important for FLG applications in microelectronics, interconnects, and thermal management structures.

Chapter 7

Conclusion and future works

Conclusion

In this manuscript, we have validated the MD method to calculate thermal resistance from the temperature difference fluctuations at equilibrium. With this method, the thermal contact resistances in various nanostructures have been studied. The main conclusions of this thesis are the followings:

(1) The temperature difference fluctuations method is validated by calculating the thermal resistances of Si/Ge superlattices, few-layer graphenes and diameter modulated SiC nanowires. The obtained results were found to be in satisfying agreement with the data from equilibrium MD and NEMD simulations.

• The interface thermal conductions of Si/Ge superlattices with different period thicknesses are calculated and the results are consistent with the reported MD values obtained from exchanged power fluctuations method. The cross-plane thermal conductivities of Si/Ge superlattice functions period thickness were also calculated from Green-Kubo formula, which agree with the ones derived from λ = G • L p .

• We have defined the thermal resistance of diameter modulated SiC nanowires as a constriction resistance and have shown that it is worth twice the internal thermal resistance of the constant diameter nanowire.

The density of states study shows that phonon modification is related to the confinement of the modes and additional surface modes appear with the diameter modulation. There is a factor of two discrepancy between the data calculated from EMD and NEMD, which may be due to the artefacts generated from NEMD at very small system size.

• The inter-FLG resistances with different layer numbers are calculated and are found to be decreasing against the total layer number, before reaching the graphite limit with sufficient layer number. The observed results comply well with the data calculated from NEMD simulations and Debye model.

These results proved the validity of the temperature difference fluctuations method in calculating thermal conductance (resistance).

(2) The thermal contact resistance between CNT and polymer was calculated from MD simulations by equilibrium fluctuations of temperature difference.

In agreement to the MD prediction, smaller contact resistance between CNT and polymer was achieved for CNT-HLK5 surface than for CNT-PEMA surface. This is because the strong covalent bonds between CNT and HLK5 resulted in a lower barrier for heat carrier. The lowest TIR of Si/CNT/HLK5/Cu TIM was measured as low as 1.40 mm 2 KW -1 with VACNTs of 10 µm length. Finally, we report that using materials with functional groups containing nitrogen (such as azide) to bond with carbon based materials (such as CNTs, graphene, carbon nanowires, et al.) by C-N covalent bond may be the most efficient way to reduce the thermal contact resistance between them.

(3) The thermal contact resistance between few layer graphene and SiO 2 was found to be decreasing with layer number and converged for six graphene layers. The contact resistance decrease is ascribed to the stronger coupling between FLG and substrate when there are more graphene layers interacting with SiO 2 . The contact resistance between SLG and SiO 2 was found to be 3.34×10 -8 m 2 KW -1 , which is more than one order of magnitude larger than the inter-layer resistance in graphite. Owing to such a resistive interface, MD simulations show that SiO 2 /FLG superlattices have a low thermal conductivity of 0.30 W/mK, exhibiting a promising prospect in nano-scale thermal insulation. These findings are promising in thermal barrier coatings and phase change memory materials and provide useful information in the thermal design of small-scale devices.

(4) We have shown that the cross-plane thermal resistance in suspended FLG decreases as the layer number increase. The corresponding effective thermal conductivity increased when it is supported on a SiO 2 substrate. Introducing a substrate or more graphene layers to the suspended FLG will limit the free vibrations of ZA modes and pose resistance on in-plane phonon transport, but transfer the energy to the cross-plane. Phonon density of states of suspended and supported FLGs are calculated, and a red-shift of the peak of the optical transverse modes (TO) from supported to suspended FLG is found, which could be attributed to the strain effects due to the surface stresses. Phonon dispersion relations show that part of the ZA modes

Future works

There are four extensions of the work presented here that should be pursued:

(1) The method based on temperature difference fluctuations was designed to calculate the thermal resistance between two sub-systems. It is necessary to extend this method to calculate the thermal contact resistance of three sub-systems coupling together.

(2) In Chapter 4, we calculated the thermal resistance between CNT and one monomer of polymer. But experimental measurements usually measure the total resistance between CNT and bulk polymers, which should be calculated directly from MD simulations rather than making assumptions.

(3) Finally, it would be interesting to investigate the relation between in-plane and cross-plane thermal conductivity of FLG. Due to the classical feature of MD simulations, it could not calculate the in-plane thermal conductivity of FLG, owning to the high Debye temperature of graphene. However, this could be done with quantum corrections.

frequency dependent transmission:

T (ω) = β/i2ω
i∈A,j∈B α,δ∈{x,y,z} k α,δ i,j uα i (ω) u * δ j (ω), (B.5)

where β = 1/k B T and the T is the system mean temperature. Eq. B.2 highlights that the atoms involved are those in interaction across the interface. The phonon transmission thus corresponds to the velocity cross spectral density of these contact atoms weighted by the interatomic cohesive force constant that maintains the two systems in interaction. From this consideration, an additional characteristic interface length l c arises. It corresponds to the thickness that ensures the cohesive energy of the contact. In MD simulations, l c is artificially prescribed by the cut-off distance of the selected interatomic potential.

1

 1 Modulation of the few layer graphene/ SiO 2 thermal contact resistance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.1.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . 5.1.2 Molecular dynamics setup . . . . . . . . . . . . . . . . . 5.1.3 Thermal contact resistance between FLG and SiO 2 versus layer number . . . . . . . . . . . . . . . . . . . . . . . . 5.1.4 Discussions . . . . . . . . . . . . . . . . . . . . . . . . . 5.1.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . .

  Figure reprinted from Jin-Cheng Zheng, Front Phys. China, 2008, 3(3): 269-279. . . . . . . . . . . . . . . . . . 3 1.3 Illustration of thermoelectric modules. (a) Cooling module. (b) power generation module. . . . . . . . . . . . . . . . . . . . . . 4 1.4 Coordinates describing the deformation of a linear monatomic lattice of lattice constant a. . . . . . . . . . . . . . . . . . . . . 6 1.5 Dispersion relation of one dimensional lattice vibration. q defines the phonon wave-vector in the above dispersion curve. . . 6 1.6 One-dimensional example of how allowed wave vectors are determined. Figure reprinted from Alan J. H. McGaughey, Phonon transport in molecular dynamics simulations: formulation and thermal conductivity prediction, Advances in Heat Transfer, volume 39, 2006. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 1.7 An incident wave (black line) reflects (black line) on the surface with phase reversal. The superposition of the incident and reflected waves produces a stationary wave of twice the amplitude (red line). Figure reprinted from S. Volz, Thermal nanosystems and nanomaterials, Springer 2009. . . . . . . . . . . . . . . . . 10 1.8 (a) When phonon mean free path Λ ≪ system size D, phonon propagation is in the diffusive region, the effective thermal conductivity λ ef f is calculated from Fourier law; (b) When Λ ≫ system size D, phonon propagation is in the ballistic region, interfacial scattering predominates over internal scattering, the effective thermal conductivity is governed by the system size and the interface resistance R. . . . . . . . . . . . . . . . . . . . . . 11 xvi List of Figures

  Figure reprinted from Subrina et al., IEEE electron device letters 30, 1281-1283, (2009) 3.13 Spatial configuration of of FLG systems: the black solid rectangles symbolize the graphene sheets in one subsystem with a mean temperature T 1 ; the open rectangles symbolize the graphene sheets in the other subsystem with a mean temperature T 2 . Red dotted lines indicate the interfaces of the corresponding FLG-FLGs. The inter-layer resistance is marked as R(m,n), where m and n refer to the number of layers in the two FLG subsystems, respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xviii List of Figures 3.14 (a) Normalized ACF of the temperature difference ∆T for calculating R(3,3). The integral of the ACF vs time is also plotted in the ordinate axis on the right hand side. (b) The calculated inter-FLG thermal resistance as a function of layer number at 300 K (black triangles). The averaged cross-plane thermal resistance from NEMD study (red rhombs), Debye model calculations (blue squares) and resistance between neighboring layers in graphite (dotted line) are also shown. . . . . . . . . . . . . . . . . . . . 3.15 Inter-FLG thermal resistance of FLGs with different total layer numbers as a function of temperature. . . . . . . . . . . . . . . 4.1 Nanotube crystal structure. (a) A (4, 2) nanotube is constructed from graphene sheet by connecting sites O and A, and B and B'. OA and OB define the chiral vector C h and translation vector T, respectively. The rectangle OAB'B defines the unit cell of the nanotube. (b) The three examples for armchair, zigzag and chiral carbon nanotube from up to down respectively. Figure from M. Huang, Studies of Mechanically Deformed Single Wall Carbon Nanotubes and Graphene by Optical Spectroscopy, PhD thesis 2009. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Illustrations and SEM images of a substrate-supported CNT array at various degree of compression. Figure reprinted form Cola, B.A. et al. Contact mechanics and thermal conductance of carbon nanotube array interfaces. Journal of heat and Mass Transfer 52, 1-4. (2009) . . . . . . . . . . . . . . . . . . . . . . . . . 4.3 Schematic of sputtered gold on CNT turf. Figure reprinted from Johnson et al., Nanotechnology 20 (2009) 065703. . . . . . . . . 4.4 Relative thermal boundary resistances of carbon nanotubes in contact with metals and polymers. Figure reprinted from Q. Li et al., Nano Lett., 9 (2009) 3805. . . . . . . . . . . . . . . . . . 4.5 (a) Schematic of direct CNT-Cu contact and polymer bonding; (b) SEM image of polymer bonding. Figure from Thales Research and Technology, Campus Polytechnique. . . . . . . . . . 4.6 Schematic of the CNT-polymer models used in the simulation. . 4.7 Normalized ACF of the temperature difference for CNT-PEMA and CNT-HLK5. The integral of the ACF vs time is also plotted in the ordinate axis on the right hand side. . . . . . . . . . . . 4.8 Measured thermal resistances of the samples for PEMA bonding and HLK5 bonding as a function of the CNT lengths. . . . . . 4.9 Schematic diagram of R U H : TCR between CNT and one repeating unit of HLK5; R H : total TCR between CNT and HLK5 in the TIM system; N H : number of CNT-HLK5 contact. For CNT/PEMA, the similar definitions are defined for R U P , R P and N P . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . bonding model between CNT and PEMA, HLK5. . . . 4.11 The normalized vibrational density of states of CNT, PEMA, and HLK5. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.12 The structure of HLK9. . . . . . . . . . . . . . . . . . . . . . . 4.13 Auger electron spectroscopy on washed CNTs and GaAs wafer. 4.14 CNT-Polymer interface with C-N, C-O and C-C covalent bonds for the phonon transmission calculation. . . . . . . . . . . . . . 4.15 Phonon transmission of CNT-Polymer interface: a comparison between C-C, C-N and C-O covalent bonds. . . . . . . . . . . . 5.1 (a) Optical image of a graphene flake consisting of single and multilayer graphene, (b) AFM image, 20×20 mm

Figure 1 . 1 :

 11 Figure 1.1:A thermal interface material fills the air pockets between the two surfaces, forming a continuous layer to conduct heat away. Figure from internet ( http://www.bit-tech.net/hardware/2009/02/16/all-about-tim/)

Figure 1 . 2 :

 12 Figure 1.2: Simplified diagram of the Seebeck effect. Material A is cooled at one end (in blue color) with a low temperature T and heated at the other end (in red color) with a high temperature T + ∆T . This results in a voltage difference as a function of temperature difference (∆T ). Figure reprinted from Jin-Cheng Zheng, Front Phys. China, 2008, 3(3): 269-279.

Figure 1 . 3 :

 13 Figure 1.3: Illustration of thermoelectric modules. (a) Cooling module. (b) power generation module.

Figure 1 . 4 :

 14 Figure 1.4: Coordinates describing the deformation of a linear monatomic lattice of lattice constant a.

Figure 1 . 5 :

 15 Figure 1.5: Dispersion relation of one dimensional lattice vibration. q defines the phonon wave-vector in the above dispersion curve.

Figure 1 . 6 :

 16 Figure 1.6: One-dimensional example of how allowed wave vectors are determined. Figure reprinted from Alan J. H. McGaughey, Phonon transport in molecular dynamics simulations: formulation and thermal conductivity prediction, Advances in Heat Transfer, volume 39, 2006.

Chapter 1 -

 1 Introduction 9

Figure 1 . 7 :

 17 Figure 1.7: An incident wave (black line) reflects (black line) on the surface with phase reversal. The superposition of the incident and reflected waves produces a stationary wave of twice the amplitude (red line). Figure reprinted from S. Volz, Thermal nanosystems and nanomaterials, Springer 2009.

Figure 1 .

 1 Figure 1.8: (a) When phonon mean free path Λ ≪ system size D, phonon propagation is in the diffusive region, the effective thermal conductivity λ ef f is calculated from Fourier law; (b) When Λ ≫ system size D, phonon propagation is in the ballistic region, interfacial scattering predominates over internal scattering, the effective thermal conductivity is governed by the system size and the interface resistance R.

Figure 1 . 9 :

 19 Figure 1.9: Schematic of the many possibilities within the framework of the acoustic mismatch model for phonons incident on an interface.

Figure 1 .

 1 Figure 1.10: Two systems related to thermal baths with temperatures T 1 and T 2

Figure 1 . 11 :

 111 Figure 1.11: Tools for thermal property calculation with corresponding spatial and temporal characteristic dimension. λ denotes the phonon wavelength, Λ is the phonon mean free path, L the system size, and τ r is phonon relaxation time.

22 2 . 1 .

 21 INTRODUCTION TO MOLECULAR DYNAMICS SIMULATIONS

Figure 2 . 1 :

 21 Figure 2.1: Separation distance and force vectors between atoms i and j.

Chapter 2 - 25 Figure 2 . 2 :

 22522 Figure 2.2: Gaussian distribution given in Eq. 2.14.

Figure 2 . 3 :

 23 Figure 2.3: The meaning of periodic boundary conditions in a two-dimensional case. Figure reprinted from Rapaport. D. C., The Art of Molecular Dynamics Simulation, Cambridge University Press.

Chapter 2 -

 2 Thermal contact resistance calculation from Molecular Dynamics simulation 27 2.1.6 Challenges for MD simulations

Figure 2 . 4 :

 24 Figure 2.4: Temperature jump across the interface between two media crossed by an interfacial flux q. The temperature profile in each medium is schematically represented by red solid lines. Figure reprinted from Merabia and Termentzidis, Physical Review B 86, 094303 (2012).

Chapter 2 - 29 Fig. 2 .

 2292 Fig. 2.4 shows a typical model studied by NEMD simulation and the corresponding temperature profile. The direct method has a physical intuitiveness and is analogous to an experimental technique.

Figure 2 . 5 :

 25 Figure 2.5: Computed thermal conductivity as a function of the layer thickness using the Green-Kubo method and that derived from the contact resistance contribution. Figure reprinted fromChalopin et al, Physical Review B 85, 195302 (2012).

Figure 2 . 6 :

 26 Figure 2.6: Schematic of the MD system simulated at equilibrium from temperature difference fluctuations method. The central segment consists in a wire of conventional argon with length L/2 surrounded by two wires of length L/4 made of argon but with a modified mass. Figure reprinted from Rajabpour and Volz, Journal of applied physics, 108, 094324 (2010).

Figure 2 . 7 :

 27 Figure 2.7: Normalized ACF of the temperature difference ∆T (black line) and its exponential fit (gray line). The integral of the ACF having units displayed in the ordinate axis on the right hand side is also plotted vs the integration time. The mass m 2 of materials 2 is equal to four times the argon atomic mass m 1 . Figure reprinted from Rajabpour and Volz, Journal of applied physics, 108, 094324 (2010).

Chapter 2 -

 2 Thermal contact resistance calculation from Molecular Dynamics simulation 37 Note that the unit of resistance calculated from Eq. 2.45 is K/W . This formula has been used for calculating the Kapitza resistance such as for solid/liquid interface[START_REF] Barrat | Kapitza resistance at the liquid-solid interface[END_REF] and superlattices,(McGaughey and Li 2006) by considering the fluctuating heat power across the interface. In the same manner, it is possible to use the fluctuation-dissipation theorem to express the thermal resistance of a thick system as proposed in Eq. 2.45.In order to have a finite thermal conductance, the autocorrelation function should reach zero. Fig.2.8 shows an example of the heat flux autocorrelation functions for liquid and solid argon at 90 K and 60 K.

Figure 2 . 8 :

 28 Figure 2.8: Heat flux autocorrelation function for argon liquid and solid in the Cartesian coordinates. Figure reprinted from S. Yip, Handbook of Materials Modeling: Models. Springer Science(2005).

Figure 3 . 1 :

 31 Figure 3.1: Cross-sectional transmission electron micrograph of Si/Ge superlattice grown on a Si substrate. Figure reprinted from Steenbergen et al., Applied Physics Letters 99, 071111 (2011).

Figure 3 . 2 :

 32 Figure 3.2: Spatial configuration of a 3× 3× 10 nm Si/Ge superlattice, with a period length L equals 5 nm.

Figure 3 . 3 :

 33 Figure 3.3: Temperatures T 1 and T 2 are calculated from the interacting atoms across the interface.

Figure 3 . 4 :

 34 Figure 3.4: Normalized ACF of the temperature difference ∆T for calculating the Kapitza resistance of a 3 × 3× 6 nm Si/Ge superlattice. The integral of the ACF vs time is also plotted in the ordinate axis on the right hand side.

Figure 3 . 5 :

 35 Figure 3.5: Interface thermal conductance as a function of the period thickness, calculated from the temperature difference fluctuation method. A comparison with data calculated from the exchanged power fluctuations method.

Figure 3 . 6 :

 36 Figure 3.6: EMD thermal conductance predictions of a L p = 20 nm Si/Ge superlattice over 300 K to 900 K.

Figure 3 . 7 :

 37 Figure 3.7: Computed thermal conductivity as a function of the layer thickness using the Green-Kubo method and that derived from the contact resistance contribution.

Figure 3 . 8 :

 38 Figure 3.8: Diameter modulated 2H SiC nanowire.

Figure 3 .

 3 Figure3.10: The four considered groups for the calculation of the partial density of states. 1: Group of atoms far from the interface for the small cross-section nanowire (8.81 nm 2 ); 2: group of atoms for the first bilayer closest to the interface; 3, 4: groups of atoms for the large (15.67 nm 2 ) cross-section nanowire far from and the first bilayer closest to the interface.

Figure 3 . 11 :

 311 Figure 3.11: Density of states for four distinct groups (see Fig. 3.10) for frequencies (a) 0-10 THz and (b) 28-40 THz.

Figure 3 .

 3 Figure 3.12: Silicon-on-insulator (SOI) circuits with graphene lateral heat spreders attached to the side heat sinks. Figure reprinted from Subrina et al., IEEE electron device letters 30, 1281-1283, (2009)

Chapter 3 -

 3 Validation of the temperature difference fluctuations method

Figure 3 .

 3 Figure 3.13: Spatial configuration of of FLG systems: the black solid rectangles symbolize the graphene sheets in one subsystem with a mean temperature T 1 ; the open rectangles symbolize the graphene sheets in the other subsystem with a mean temperature T 2 . Red dotted lines indicate the interfaces of the corresponding FLG-FLGs. The inter-layer resistance is marked as R(m,n), where m and n refer to the number of layers in the two FLG subsystems, respectively.

Figure 3 .

 3 Figure 3.14: (a) Normalized ACF of the temperature difference ∆T for calculating R(3,3). The integral of the ACF vs time is also plotted in the ordinate axis on the right hand side. (b) The calculated inter-FLG thermal resistance as a function of layer number at 300 K (black triangles). The averaged cross-plane thermal resistance from NEMD study (red rhombs), Debye model calculations (blue squares) and resistance between neighboring layers in graphite (dotted line) are also shown.

Fig. 3 .Chapter 3 -

 33 Fig. 3.14(b) shows the calculated inter-FLG thermal resistance as a function of layer numbers at 300K, and the values are compared with that from non-

Figure 3 .

 3 Figure 3.15: Inter-FLG thermal resistance of FLGs with different total layer numbers as a function of temperature.

Chapter 3 -

 3 Validation of the temperature difference fluctuations method 59 depend on the temperature, guaranteeing stable performances of FLG based structures under different working temperatures. The observed results comply well with the data calculated from NEMD simulations and Debye model.

  is an individual graphite layer), as is shown in Fig.4.1(a).

Figure 4 . 1 :

 41 Figure 4.1: Nanotube crystal structure. (a) A (4, 2) nanotube is constructed from graphene sheet by connecting sites O and A, and B and B'. OA and OB define the chiral vector C h and translation vector T, respectively. The rectangle OAB'B defines the unit cell of the nanotube. (b) The three examples for armchair, zigzag and chiral carbon nanotube from up to down respectively. Figure from M. Huang, Studies of Mechanically Deformed Single Wall Carbon Nanotubes and Graphene by Optical Spectroscopy, PhD thesis 2009.

Figure

  Figure 4.1(a) shows the formation of a[START_REF] Ni | Substrate Induced Cross-Plane Thermal Propagative Modes in Few Layer Graphene[END_REF][START_REF] Ni | Significant thickness dependence of the thermal resistance between few-layer graphenes[END_REF] nanotube. By connecting the sites O and A, and B and B', a single wall carbon nanotube will be constructed out of graphene sheet. The vectors OA and OB define the chiral vector C h and the translation vector T of a carbon nanotube, respectively. The rectangle OAB'B gives the unit cell of the one dimensional structure. The chiral vector uniquely determines the nanotube structure and can be expressed on the basis of unit vector of graphene lattice:

  are integers and 0 ≤ m ≤ n. As is shown in Fig. 4.1(b), an armchair nanotube corresponds to the case of n = m, and for zigzag nanotube, m = 0. All other indices correspond to chiral nanotubes.

Figure 4 . 2 :

 42 Figure 4.2: Illustrations and SEM images of a substrate-supported CNT array at various degree of compression. Figure reprinted form Cola, B.A. et al. Contact mechanics and thermal conductance of carbon nanotube array interfaces. Journal of heat and Mass Transfer 52, 1-4. (2009)

Figure 4 . 3 :

 43 Figure 4.3: Schematic of sputtered gold on CNT turf. Figure reprinted from Johnson et al., Nanotechnology 20 (2009) 065703.

Figure 4 . 4 :

 44 Figure 4.4: Relative thermal boundary resistances of carbon nanotubes in contact with metals and polymers. Figure reprinted from Q.Li et al., Nano Lett., 9 (2009) 3805.

4. 2 .

 2 CONTACTS BETWEEN VACNTS AND THE SUPERSTRATE et al. 2010). The Si/CNT/PEMA/Cu TIM thermal interface resistance (TIR) was significantly reduced by one order of magnitude compared with the direct CNT/Cu contact.

Figure 4 .

 4 Figure 4.5: (a) Schematic of direct CNT-Cu contact and polymer bonding; (b) SEM image of polymer bonding. Figure from Thales Research and Technology, Campus Polytechnique.

Figure 4 . 6 :

 46 Figure 4.6: Schematic of the CNT-polymer models used in the simulation.

Chapter 4 -

 4 Thermal contact resistance between carbon nanotubes and functionalized polymers 694.3.2 ResultsThe ACFs and the integration time (integral of the ACFs vs the simulation time) of the two different models are presented in Fig.4.7.

Figure 4 . 7 :

 47 Figure 4.7: Normalized ACF of the temperature difference for CNT-PEMA and CNT-HLK5. The integral of the ACF vs time is also plotted in the ordinate axis on the right hand side.

Figure 4 . 8 :

 48 Figure 4.8: Measured thermal resistances of the samples for PEMA bonding and HLK5 bonding as a function of the CNT lengths.

Fig. 4 .

 4 Fig.4.8 shows the measured thermal resistance with HLK5 bonding as a function of VACNT lengths. The same plot with PEMA bonding that has been published in (LeKhanh et al. ) is also presented as a comparison. When PEMA layer is substituted by HLK5, the total resistance of the sample is decreased. More than 25% decrease of the total resistance was found in all the considered CNT lengths, indicating that optimized thermal paths are generated with CNT-polymer covalent bonding. The contact resistance values at zero CNT length with PEMA bonding and HLK5 bonding were found to be 1.1 mm 2 KW -1 and 0.34 mm 2 KW -1 , respectively. This contact resistance contains contributions of resistances from the Si-CNT, CNT-polymer, polymer-Cu contact, as well as the intrinsic resistance from polymer. As the polymer layer is quite thin and it fully contacts the copper surface, the resistance originating from polymer and that between polymer and Cu can be considered as negligible. As a result, the thermal interface resistance between CNT and polymer should be the dominating factor to the thermal contact resistance. The TIR with PEMA bonding is 1.50 mm 2 KW -1 , and it is 4.36 times as large as TIR with HLK5 bonding:

  Fig 4.9.

Figure 4 . 9 :

 49 Figure 4.9: Schematic diagram of R U H : TCR between CNT and one repeating unit of HLK5; R H : total TCR between CNT and HLK5 in the TIM system; N H : number of CNT-HLK5 contact. For CNT/PEMA, the similar definitions are defined for R U P , R P and N P .

Figure 4 .

 4 Figure 4.10: Possible bonding model between CNT and PEMA, HLK5.

Figure 4 . 11 :

 411 Figure 4.11: The normalized vibrational density of states of CNT, PEMA, and HLK5.

Figure 4 .

 4 Figure 4.12: The structure of HLK9.

Figure 4 . 13 :

 413 Figure 4.13: Auger electron spectroscopy on washed CNTs and GaAs wafer.

Figure 4 .

 4 Figure 4.14: CNT-Polymer interface with C-N, C-O and C-C covalent bonds for the phonon transmission calculation.

Figure 4 .

 4 Figure 4.15: Phonon transmission of CNT-Polymer interface: a comparison between C-C, C-N and C-O covalent bonds.

Figure 5 . 1 :

 51 Figure 5.1: (a) Optical image of a graphene flake consisting of single and multilayer graphene, (b) AFM image, 20×20 mm 2 . Region I, monolayer graphene; region II, double-layer graphene; region III, triple-layer graphene; region VI, multilayer graphene; region V, SiO 2 substrate. Figure reprinted from Kim et al. Adv. Mater.2008, 20, 3589.

Figure 5 . 3 :

 53 Figure 5.3: Contact resistance between FLG and SiO 2 substrate versus layer number at 300 K: a comparison between MD simulations and measurements. The red dashed line is the averaged value of the experimental results. The black dashed lines are the typical range of the measured contact resistance between graphene and SiO 2 from Table 5.1.

Figure 5 . 4 :

 54 Figure 5.4: Temperature difference autocorrelation functions (solid lines) and the running time integrals (dashed lines) for 1,3,5-layer graphene-SiO 2 interfaces at 300 K.

Figure 5 . 5 :

 55 Figure 5.5: Thermal resistance between SiO 2 and single layer graphene with different cut-off radii.

Figure 5 . 6 :

 56 Figure 5.6: Contact resistance between FLG and a SiO 2 substrate versus layer number at 300 K: a comparison between R c = 2.5 σ and R c = 5 σ.

Figure 5 . 7 :

 57 Figure 5.7: Compilation of the thermal contact resistance of solid-solid interfaces at 300 K.

Figure 5 . 8 :

 58 Figure 5.8: Thermal contact of a graphene encased by Silicon dioxide: an encased graphene layer provides the same thermal resistance as a 91.6 nm thick SiO 2 layer.

Figure 5 . 9 :

 59 Figure 5.9: Spatial configuration of SiO 2 /SLG and SiO 2 /FLG superlattice.

90

 90

Figure 5. 10 :

 10 Figure 5.10:Period thickness dependence of the thermal conductivities of c-SiO 2 /SLG and a-SiO 2 /SLG superlattices. The thermal conductivities of the corresponding bulk SiO 2 (parallel to c-axis) were also plotted as benchmarks.

Figure 5 . 11 :

 511 Figure 5.11: Thermal conductivity of a-SiO 2 /FLG superlattice as a function of graphene layer number.

Figure 5 . 12 :

 512 Figure 5.12: Thermal conductivity of a-SiO 2 /5-layer graphene superlattice as a function of temperature.

Figure 6 . 1 :

 61 Figure 6.1: FLG as heat spreaders for AlGaN/GaN field effect transistors (FET).Figure reprinted from Z. Yan, et al. Nature Communications 3(827), 2012.

  Figure 6.1: FLG as heat spreaders for AlGaN/GaN field effect transistors (FET).Figure reprinted from Z. Yan, et al. Nature Communications 3(827), 2012.

Figure 6 . 2 :

 62 Figure 6.2: Cross-plane coupling of the low-energy phonons decreases the in-plane thermal conductivity.

  98 6.2. FLG CROSS-PLANE THERMAL RESISTANCE CALCULATIONS FROM MD SIMULATIONS 6.2 FLG cross-plane thermal resistance calculations from MD simulations 6.2.1 Molecular dynamics setup

Figure 6 . 3 :

 63 Figure 6.3: Spatial configuration of (a) 5-layer suspended graphene and (b) 5-layer graphene supported on a SiO 2 substrate.

  C-C =2.39 meV, ǫ Si-C = 8.909 meV, ǫ O-C = 3.442 meV, σ C-C =3.410 Å, σ Si-C =3.326 Å and σ O-C =3.001 Å[START_REF] Ong | Molecular dynamics simulation of thermal boundary conductance between carbon nanotubes and sio 2[END_REF][START_REF] Girifalco | Carbon nanotubes, buckyballs, ropes, and a universal graphitic potential[END_REF]. The cutoff distance of the LJ potential for the Si-C and O-C interactions were set equal to 2.5 σ (Ong and Pop 2010). The structure was relaxed using a steepest descent energy minimization algorithm (sd). The simulation first ran in the isothermal-isobaric (NPT) ensemble with a temperature of 300 K and a pressure of 0 bar for 400 ps, then in the canonical ensemble (NVT) for another 200 ps to equilibrate the system. Atomic trajectories were then calculated in the microcanonical ensemble (NVE) at 300 K. More than twenty Chapter 6 -Cross-plane thermal properties in suspended and supported few layer graphenes 99 independent simulations were performed for each resistance point, until the resistance (averaged from different simulations) converged. The error bars were then defined by looking at the variance of the averaged values.

Fig. 6

 6 Fig.6.4(a) shows the heat flux autocorrelation function of suspended and supported 4-layer graphene. Fig.6.4(b) compares the calculated values of the cross-plane resistance R c of both FLG configurations against layer number. For suspended FLG, there is a drastic R c jump when the layer number increases from two to four, which may be due to the change in the number of interacting layers. The inter-layer Van der Waals force has a cutoff radius of 10.2 Å in our simulations, and one layer in FLG could interact directly with its third nearest neighboring layer. Within the cutoff radius, a larger layer number brings more interacting layers, which increases inter-layer phonon coupling. For FLGs with layer number larger than four, R c decrease is mainly due to the size effect, which could be eliminated with sufficient layers, and R c approaches the graphite limit [the black dashed line[START_REF] Taylor | The thermal conductivity of pyrolytic graphite[END_REF]] , as is shown in Fig.6.4(b).

Figure 6 . 4 :

 64 Figure 6.4: (a) Heat current autocorrelation function of suspended and supported 4-layer graphene. (b) Cross-plane thermal resistance of suspended and supported FLG versus the layer number at 300 K. Dashed line is the resistance between neighboring layers in graphite.

Figure 6 . 5 :

 65 Figure 6.5: Effective cross-plane thermal conductivity of suspended and SiO 2 supported FLG versus the number of graphene layers at 300 K. NEMD data of cross-plane thermal conductivity for suspended FLG is also shown for comparison.

6. 3 .

 3 Figure 6.6: The phonon density of states of suspended FLG with different layer numbers. The insert shows the DOS at low frequency range.

Fig. 6 .

 6 Fig. 6.7 compares the DOS of suspended and supported FLG, with layer numbers of 2, 4, and 10.

Chapter 6 -

 6 Figure 6.7: A comparison of phonon density of states of suspended and supported FLG (a) 2-layer graphene, (b) 4-layer graphene, (c)10-layer graphene.

Figure 6 . 8 :

 68 Figure 6.8: Size of graphene sheets used for calculating dispersion relation.

Figure 6 . 10 :

 610 Figure 6.10: The calculated longitudinal and transverse branches (optical and acoustic) of phonon dispersion relation of supported 3-layer graphene. The spectrums of phonon density of states calculated from all modes (all velocity components, black line) and only from the flexural modes (cross-plane direction velocity, red line) are also shown.

Figure 6 . 11 :

 611 Figure 6.11: Phonon dispersion relations of suspended and supported 1, 3, 5-layer graphene.

Figure 6 .

 6 Figure 6.12: ZA mode branch of 3-layer suspended and supported FLG.

Figure 6 .

 6 Figure 6.13: ZA modes branches of supported 3-layer graphene with longitudinal and transverse sound lines: in the grayed and patterned regions, above the sound lines, graphene modes are coupled with the propagative modes of the substrate. The surface Rayleigh wave dispersion is also shown in red line.

Figure 6 . 14 :

 614 Figure 6.14: Schematic presentation of the Frenkel-Kontorova model: A chain of particles interacting via harmonic springs with elastic coupling g is subjected to the action of an external periodic potential with period a s .

  x n = na s + u n . When the particle displacements u n are small, |u n | ≪ a s , we can linearize the motion equation and obtaind 2 u n dt 2 + u n -g(u n+1 + u n-1 -2u n ) = 0.(6.11)Substituting the plane wave solution u n (t) ∝ exp[iωt -ikn] to Eq. 6.11, the phonon dispersion relation is obtained as:ω 2 (k) = ω 2 min + 2g(1 -cos k).

Chapter 7 -

 7 Conclusion and future works 115 of supported FLG were confined and standing waves were formed. The Frenkel-Kontorova model and sound lines were introduced to explain the substrate-induced band gaps in the FLG dispersion relation and the corresponding effects on FLG cross-plane thermal properties. The enhanced thermal conduction in the cross-plane direction is ascribed to the additional cross-plane heat fluxes generated by FLG-substrate phonon radiation. The obtained results are important for FLG applications in microelectronics, interconnects, and thermal management structures.

  5.2 Few layer graphene based superlattices as efficient thermal insulators . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . 5.2.2 Molecular dynamics setup . . . . . . . . . . . . . . . . . 5.2.3 Thermal conductivities of FLG/SiO 2 superlattices . . . 5.2.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . 5.3 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Phonon density of states of FLGs . . . . . . . . . . . . . . . . . 6.4 Phonon dispersion relations of FLGs . . . . . . . . . . . . . . . .4.3 The Frenkel-Kontorova Model . . . . . . . . . . . . . . . 6.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Future works . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.3 List of publications during PhD study . . . . . . . . . . . . . . Internal thermal resistance and constriction resistance calculated from EMD and NEMD. . . . . . . . . . . . . . . . . . . . . . . 3.2 The average inter-FLG distance of the FLG-FLG systems with different total layer numbers at 300 K. . . . . . . . . . . . . . .List of Figures1.1 A thermal interface material fills the air pockets between the two surfaces, forming a continuous layer to conduct heat away. Figure from internet ( http://www.bit-tech.net/hardware/2009/02/16/allabout-tim/) . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.2 Simplified diagram of the Seebeck effect. Material A is cooled at

	List of Tables
	3.1
	A Derivation of phonon density of states
	B Derivation of phonon transmission
	References

6 Cross-plane thermal properties in suspended and supported few layer graphenes 6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.2 FLG cross-plane thermal resistance calculations from MD simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.2.1 Molecular dynamics setup . . . . . . . . . . . . . . . . . 6.2.2 Results and discussion . . . . . . . . . . . . . . . . . . . 6.3 6.4.1 Phonon dispersion relation calculation method . . . . . 6.4.2 Phonon dispersion relation of suspended and supported FLGs . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 Conclusion and future works 7.1 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.2 4.1 Thermal resistance and thermal conductance of CNT-PEMA bonding and CNT-HLK5 bonding calculated from EMD simulation 5.1 Thermal resistance of the graphene-substrate interfaces. Table reprinted from A. A. Balandin, Proc.SPIE 8101, Carbon Nanotubes, Graphene, and Associated Devices IV . . . . . . . . . . one end (in blue color) with a low temperature T and heated at the other end (in red color) with a high temperature T + ∆T . This results in a voltage difference as a function of temperature difference (∆T ).

  1.9 Schematic of the many possibilities within the framework of the acoustic mismatch model for phonons incident on an interface.1.10 Two systems related to thermal baths with temperatures T 1 and T 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.11 Tools for thermal property calculation with corresponding spatial and temporal characteristic dimension. λ denotes the phonon wavelength, Λ is the phonon mean free path, L the system size, and τ r is phonon relaxation time. . . . . . . . . . . . . . . . . 2.1 Separation distance and force vectors between atoms i and j. . 2.2 Gaussian distribution given in Eq. 2.14. . . . . . . . . . . . . . 2.3 The meaning of periodic boundary conditions in a two-dimensional case. Figure reprinted from Rapaport. D. C., The Art of Molecular Dynamics Simulation, Cambridge University Press. . . . . 2.4 Temperature jump across the interface between two media crossed by an interfacial flux q. The temperature profile in each medium is schematically represented by red solid lines. Figure reprinted from Merabia and Termentzidis, Physical Review B 86, 094303 (2012). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.5 Computed thermal conductivity as a function of the layer thickness using the Green-Kubo method and that derived from the contact resistance contribution. Figure reprinted from Chalopin et al, Physical Review B 85, 195302 (2012). . . . . . . . . . . . 2.6 Schematic of the MD system simulated at equilibrium from temperature difference fluctuations method. The central segment consists in a wire of conventional argon with length L/2 surrounded by two wires of length L/4 made of argon but with a modified mass. Figure reprinted from Rajabpour and Volz, Journal of applied physics, 108, 094324 (2010). . . . . . . . . . 2.7 Normalized ACF of the temperature difference ∆T (black line) and its exponential fit (gray line). The integral of the ACF having units displayed in the ordinate axis on the right hand side is also plotted vs the integration time. The mass m 2 of materials 2 is equal to four times the argon atomic mass m 1 . Figure reprinted from Rajabpour and Volz, Journal of applied physics, 108, 094324 (2010). . . . . . . . . . . . . . . . . . . . . 2.8 Heat flux autocorrelation function for argon liquid and solid in the Cartesian coordinates. Figure reprinted from S. Yip, Handbook of Materials Modeling: Models. Springer Science(2005). . 3.1 Cross-sectional transmission electron micrograph of Si/Ge superlattice grown on a Si substrate. Figure reprinted from Steenbergen et al., Applied Physics Letters 99, 071111 (2011). . . . . . Spatial configuration of a 3 × 3× 10 nm Si/Ge superlattice, with a period length L equals 5 nm. . . . . . . . . . . . . . . . . . . 3.3 Temperatures T 1 and T 2 are calculated from the interacting atoms across the interface. . . . . . . . . . . . . . . . . . . . . . . . . 3.4 Normalized ACF of the temperature difference ∆T for calculating the Kapitza resistance of a 3 × 3× 6 nm Si/Ge superlattice. The integral of the ACF vs time is also plotted in the ordinate axis on the right hand side. . . . . . . . . . . . . . . . . . . . . 3.5 Interface thermal conductance as a function of the period thickness, calculated from the temperature difference fluctuation method. A comparison with data calculated from the exchanged power fluctuations method. . . . . . . . . . . . . . . . . . . . . . . . . 3.6 EMD thermal conductance predictions of a L p = 20 nm Si/Ge superlattice over 300 K to 900 K. . . . . . . . . . . . . . . . . . 3.7 Computed thermal conductivity as a function of the layer thickness using the Green-Kubo method and that derived from the contact resistance contribution. . . . . . . . . . . . . . . . . . . 3.8 Diameter modulated 2H SiC nanowire. . . . . . . . . . . . . . . 3.9 Density of states for four distinct groups (see

	List of Figures	xvii
	3.2

  Compilation of the thermal contact resistance of solid-solid interfaces at 300 K. . . . . . . . . . . . . . . . . . . . . . . . . . . 5.8 Thermal contact of a graphene encased by Silicon dioxide: an encased graphene layer provides the same thermal resistance as a 91.6 nm thick SiO 2 layer. . . . . . . . . . . . . . . . . . . . . 5.9 Spatial configuration of SiO 2 /SLG and SiO 2 /FLG superlattice. 5.10 Period thickness dependence of the thermal conductivities of c-SiO 2 /SLG and a-SiO 2 /SLG superlattices. The thermal conductivities of the corresponding bulk SiO 2 (parallel to c-axis) were also plotted as benchmarks. . . . . . . . . . . . . . . . . . . . 5.11 Thermal conductivity of a-SiO 2 /FLG superlattice as a function List of Figures 5.12 Thermal conductivity of a-SiO 2 /5-layer graphene superlattice as a function of temperature. . . . . . . . . . . . . . . . . . . . . . 6.1 FLG as heat spreaders for AlGaN/GaN field effect transistors (FET). Figure reprinted from Z. Yan, et al. Nature Communications 3(827), 2012. . . . . . . . . . . . . . . . . . . . . . . . . 6.2 Cross-plane coupling of the low-energy phonons decreases the in-plane thermal conductivity. . . . . . . . . . . . . . . . . . .

	Chapter 1
	6.3 Spatial configuration of (a) 5-layer suspended graphene and (b) 5-layer graphene supported on a SiO 2 substrate. . . . . . . . . . 6.4 (a) Heat current autocorrelation function of suspended and sup-Introduction
	1.1 Scientific and Technological Challenges related to
	interfaces
	6.12 ZA mode branch of 3-layer suspended and supported FLG. . .
	6.13 ZA modes branches of supported 3-layer graphene with longitu-
	dinal and transverse sound lines: in the grayed and patterned
	regions, above the sound lines, graphene modes are coupled with
	the propagative modes of the substrate. The surface Rayleigh

2

. Region I, monolayer graphene; region II, double-layer graphene; region III, triple-layer graphene; region VI, multilayer graphene; region V, SiO 2 substrate. Figure reprinted from Kim et al. Adv. Mater. 2008, 20, 3589. . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2 Spatial configuration of a 5-layer graphene system supported by a silicon dioxide substrate. . . . . . . . . . . . . . . . . . . . . . 5.3 Contact resistance between FLG and SiO 2 substrate versus layer number at 300 K: a comparison between MD simulations and measurements. The red dashed line is the averaged value of the experimental results. The black dashed lines are the typical range of the measured contact resistance between graphene and SiO 2 from Table 5.1. . . . . . . . . . . . . . . . . . . . . . . . . 5.4 Temperature difference autocorrelation functions (solid lines) and the running time integrals (dashed lines) for 1,3,5-layer graphene-SiO 2 interfaces at 300 K. . . . . . . . . . . . . . . . . . . . . . . 5.5 Thermal resistance between SiO 2 and single layer graphene with different cut-off radii. . . . . . . . . . . . . . . . . . . . . . . . . 5.6 Contact resistance between FLG and a SiO 2 substrate versus layer number at 300 K: a comparison between R c = 2.5 σ and R c = 5 σ. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.7 of graphene layer number. . . . . . . . . . . . . . . . . . . . . . xx ported 4-layer graphene. (b) Cross-plane thermal resistance of suspended and supported FLG versus the layer number at 300 K. Dashed line is the resistance between neighboring layers in graphite. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.5 Effective cross-plane thermal conductivity of suspended and SiO 2 supported FLG versus the number of graphene layers at 300 K. NEMD data of cross-plane thermal conductivity for suspended FLG is also shown for comparison. . . . . . . . . . . . . . . . . 6.6 The phonon density of states of suspended FLG with different layer numbers. The insert shows the DOS at low frequency range.102 6.7 A comparison of phonon density of states of suspended and supported FLG (a) 2-layer graphene, (b) 4-layer graphene, (c)10layer graphene. . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.8 Size of graphene sheets used for calculating dispersion relation. 6.9 The mode kinetic energy as a function of frequency and wave vector for longitudinal atomic motions of suspended 3-layer graphene.105 6.10 The calculated longitudinal and transverse branches (optical and acoustic) of phonon dispersion relation of supported 3-layer graphene. The spectrums of phonon density of states calculated from all modes (all velocity components, black line) and only from the flexural modes (cross-plane direction velocity, red line) are also shown. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.11 Phonon dispersion relations of suspended and supported 1, 3, 5-layer graphene. . . . . . . . . . . . . . . . . . . . . . . . . . . wave dispersion is also shown in red line. . . . . . . . . . . . .

Table 3 . 1 :

 31 Internal thermal resistance and constriction resistance calculated from EMD and NEMD.

		Internal thermal resistance (R I , 10 -10 m 2 K/W)	Constriction resistance (R C , 10 -10 m 2 K/W)	R C /R I
	EMD	1.25	2.37	1.89
	NEMD	0.6	0.9	1.50

Table 4 .

 4

	1: Thermal resistance and thermal conductance of CNT-PEMA bonding and
	CNT-HLK5 bonding calculated from EMD simulation
		CNT-PEMA CNT-HLK5
	Thermal Resistance (GK/W)	15.0	4.92
	Thermal Conductance (nW/K)	0.07	0.20
	4.4 Comparison between calculated thermal resistances
	from molecular dynamics simulations and exper-
	imental data
	4.4.1 Experimental setup
	4.4.1.1 HLK5 synthesis
	Paratoluene 4-dimethyl amino pyridinium sulfonate (DPTS) was synthesized
	from DMAP (4-dimethyllaminopryridine) and methanol. 6-azido hexanoic acid

  2 KW -1 , and it is 4.36 times as large as TIR with HLK5 bonding:

	R P R H	≈	1.50 0.34	= 4.36.	(4.2)

Table 5 . 1 :

 51 Thermal resistance of the graphene-substrate interfaces. Table reprinted from A. A. Balandin, Proc.SPIE 8101, Carbon Nanotubes, Graphene, and Associated Devices IV

	Interface	Thermal Resistance	Method
	Graphene/SiO 2	∼ 4 × 10 -8 (Km 2 /W )	Raman/Electrical
	Graphene/SiO 2	∼ (0.6 -12) × 10 -8 (Km 2 /W )	Electrical
	Graphene/SiO 2	∼ 2 × 10 -8 (Km 2 /W )	Pump-Probe
	Au/Ti/graphene/SiO 2	∼ 4 × 10 -8 (Km 2 /W )	Raman/Electrical
	Bulk Graphite/Metal	∼ (1 -3) × 10 -8 (Km 2 /W )	Reflectance
	Graphene/a-SiO 2	∼ 4 × 10 -8 (Km 2 /W )	Theory
	Graphene/Oil	∼ (0.4 -4) × 10 -8 (Km 2 /W )	Theory

1.2. SCIENTIFIC BACKGROUND

2.2. METHODS FOR CALCULATING THERMAL CONTACT RESISTANCE FROM MOLECULAR DYNAMICS SIMULATIONS

3.1. INTERFACE THERMAL CONDUCTANCE IN SI/GE SUPERLATTICES

3.2. THERMAL CONSTRICTION RESISTANCE OF DIAMETER MODULATED SIC NANOWIRES

5.1. MODULATION OF THE FEW LAYER GRAPHENE/ SIO2THERMAL CONTACT RESISTANCE

7.1. CONCLUSION

Molecular dynamics setup

The thermal contact resistance between FLG and SiO 2 was studied by MD simulations with LAMMPS software package [START_REF] Plimpton | Fast parallel algorithms for short-range moleculardynamics[END_REF]. Bernal stacking graphene layers (3×3.7 nm) and SiO 2 (3×3.7×2 nm) with periodic boundary conditions in the x and y directions are used, as is shown in Fig. 5.2. We implemented the adaptive intermolecular reactive empirical bond order (AIREBO) potential [START_REF] Stuart | A reactive potential for hydrocarbons with intermolecular interactions[END_REF] to simulate the C-C interactions, which includes a Lennard-Jones (LJ) potential for the non-bonded Van der Waals C-C interactions between layers. The parameters in the LJ potential are ǫ C-C =2.39 meV for the potential well depth and σ C-C =3.41 Å for the distance at which the inter-particle potential is zero [START_REF] Girifalco | Carbon nanotubes, buckyballs, ropes, and a universal graphitic potential[END_REF]. A modified Tersoff potential [START_REF] Munetoh | Interatomic potential for si-o systems using tersoff parameterization[END_REF]) was used to describe the Si-Si, the Si-O, and the O-O interactions in SiO 2 . This potential can reproduce structural and dynamical properties of SiO 2 systems. One of the most important advantages using this potential is the ability to perform large-scale simulations of systems composed of Si, SiO 2 and SiO x . The potential parameters have been determined based on ab initio calculations of small molecules and the experimental data of α-quartz. Although Tersoff potential for SiO 2 does not take into account the long-range Coulombic interaction, this is less likely to matter for the problem at hand as the interfacial thermal transport depends only on the Van der Waals interaction between the graphene and the substrate atoms. [START_REF] Ong | Molecular dynamics simulation of thermal boundary conductance between carbon nanotubes and sio 2[END_REF]) also used the same potential to calculate thermal conductance between carbon nanotubes and SiO 2 and resonable results were obtained. LJ potential was adopted for the interfacial interactions between

Derivation of phonon density of states

For atom i in cell b in direction l, the displacement u could be written in the form of plane wave:

Taking the time derivative of u gives:

Summing over all the atomic contributions after multiplying by the atomic mass m:

(A.4)

Introducing the Dirac delta expression:

When q=-q', l e i(q ′ +q)x(l) b e i(q ′ +q)x(b) =1, otherwise it equqls zero. Then we rewrite Eq. A.4 as

(A.7)

The complex conjugate relation provides:

We have assumed the validity of the classical approximation to make use of the equipartition theorem. The instantaneous kinetic energy associated to each vibrational mode is given by

Taking the Fourier transform of Eq. A.10, we obtain .11) Taking Dirac delta expression of Eq. A.11, q

+∞

-∞ e i(ω-ωq)t dt = q δ(ω -ω q ). (A.12)

Calling g(ω) the density of states, we have

Applying the Wiener-Kintchine theorem to obtain the power spectral density

where v i (ω) is the Fourier Transform of v t over a time-period T s which is the trajectory duration. Finally, g(ω) can be expressed as

Derivation of phonon transmission

To express the frequency and wave-vector dependent phonon transmission, we start from the expression of the exchanged power obtained previously (Chalopin et al. 2012):

u i refers to the atomic displacement of atom i, the α and δ exponents to its x, y or z component and the dot to the time derivative. k α,δ i,j is the interatomic force constant tensor between atoms i and j. Rewriting the correlation as:

introduces a characteristic time τ , which is related to the vibrational mode dynamics whereas t refers to the evolution of the temperature constraint. Applying the Fourier transform on the correlations yields:

Eq. B.3 can be compared to the following Landauer expression: B.4) where ω stands for the energy quantum and T for the phonon transmission. n 0 A,B corresponds to the Bose-Einstein occupation numbers of systems A and B, respectively. Identifying the first term of the difference of Eq. B.3 to the one in Eq. B.4 and considering the classical limit allow extracting the following