M Brice Davier

M David

Lacroix Professeur

M J Ér Ôme Saint-Martin

M Sebastian Volz

Rapporteur M Mathieu

Rapporteur Mme

Jelena Sjakste Charg

M Jay

Theoretical study of thermoelectric properties in nanostructures

List of Tables

Si bulk thermal conductivity at 1000K from EMD and NEMD. . 2

Ge bulk thermal conductivity at 1000K from EMD and NEMD. . [START_REF] Liu | Thermal Conductivity Measurements of Ultra-Thin Single Crystal Silicon Layers[END_REF], [START_REF] Ju | Phonon Scattering in Silicon Films with Thickness of Order 100 Nm[END_REF], [START_REF] Li | Thermal Conductivity of Individual Silicon Nanowires[END_REF] and [START_REF] Hochbaum | Enhanced Thermoelectric Performance of Rough Silicon Nanowires[END_REF]). 33 Thermal conductivity in NW as function of temperature, compared with exp. results from [START_REF] Li | Thermal Conductivity of Individual Silicon Nanowires[END_REF]

Introduction

The thermoelectric conversion is a subject of active research. Its main applications are related to energy harvesting for autonomous power supply system as well as heat management in particular for CPU cooling. Yet, common thermoelectric systems unfortunately rely on rare and toxic compounds such as Bismuth Telluride, Lead Telluride, etc., . Their replacement by Silicon and Germanium which are more abundant and widely used by the microelectronic industry would be appealing. However, their intrinsic thermoelectric properties are poor especially close to room temperature.

In order to improve the thermoelectric conversion of devices based on Si and Ge, nanotechnologies provide new routes of optimization. Indeed, nanostructuring can be used to tune the intrinsic properties of nanostructures by several orders of magnitude.

In this framework, this PhD dissertation investigates the thermal properties of recently fabricated Si and Ge polytype nanowires. They are made of quasiperiodic sequences of cubic and hexagonal phases of Si or Ge. The presence of many polytype interfaces is expected to be beneficial to the thermoelectric properties. Since experimental characterizations of this nanowires are still very challenging, developing in parallel specific numerical tools is of high importance for developing this technology.

During this thesis, an original Monte Carlo simulator dedicated to phonon transport have been developed. It is a full-band and 3D simulator that is able to provide deep insights of thermal transport in complex nanodevices. As the required "semi-empirical" parameters are computed by using ab initio data, a large class of materials and phase can be modeled even the availability of experimental characterization. Besides, studies performed by using Molecular Dynamics are presented in order to parametrize semi-transparent interface modeling in the MC code.

In §1 the context of this thesis is presented. The theoretical models used through this work are presented in §2. The Monte Carlo and Molecular Dynamics simulations methods are detailed in §3. In §4, our methods are validated by studying the properties in bulk materials and thin films. §5 our implementation of rough boundaries in the Monte Carlo simulations is presented as well as the resulting simulations. Finally, in §6, the properties of interfaces between the different phases of Si and Ge are investigated.

Context

Currently, a significant part of the worldwide produced energy is rejected as wasted heat. In 2007, less than half of the produced energy was effectively consumed (cf. [START_REF] Smith | 2007 Estimated International Energy Flows[END_REF]). There is thus a strong motivation to develop technologies that are able to harvest that large amount of lost energy.

Thermoelectricity deals with the coupling between the flux of electrical charges (current) and heat flux, enabling the direct conversion of thermal energy to electrical energy. The thermal energy can be transferred in solid-state systems by several mechanisms: lattice vibrations (related to pseudo particles called phonons), electrons (that is the electronic part of thermal conduction) and electromagnetic fields (i.e. by thermal radiation, cf. G. [START_REF] Chen | Nanoscale Energy Transport and Conversion: A Parallel Treatment of Electrons[END_REF]). In this thesis, we focus on the thermal transport in semiconductors. Only the lattice vibrations are considered in this work, as other mechanisms (electronic conduction and radiation) are negligible in the studied systems. Besides, only the thermal part of the thermoelectric properties are studied.

At the macroscale, the thermal properties of materials can be well characterized by the thermal conductivity that is accurately modeled by using the Fourier's law [START_REF] Fourier | Role of Disorder and Anharmonicity in the Thermal Conductivity of Silicon-Germanium Alloys: A First-Principles Study[END_REF]. However, when the device size approaches the characteristic length scale of microscopic phenomena, the thermal properties are modified, and new theoretical models are required. Due to the high complexity of the nanoscale heat transport, numerical simulations have become very helpful in the study of nanomaterials and devices.

Thermoelectricity

The main application of thermoelectricity that we consider is the harvesting of thermal energy. While thermoelectric generators have a low efficiency with respect to other technologies (cf. [START_REF] Vining | An Inconvenient Truth About Thermoelectrics[END_REF]), they also have several advantages. Their mechanical simplicity makes them more reliable, and independent of gravity or atmospheric pressure. For instance, it makes them attractive to the spatial industry, and they have been used to power several probes and rovers. They also can be easily miniaturized and could be embedded in a wide range of products. By using such generators, temperature differences naturally present in the environment can be exploited to supply low power devices such as autonomous sensors.

Another application of thermoelectrics is the heat management. Current microprocessor technologies are limited by heat dissipation. Indeed, the increase of clock frequencies leads to an increase of the heat production in electrical devices and the limits of passive heat diffusion are reached. Thermoelectric devices could be integrated close to the sources of heat in order to extract it more efficiently [START_REF] Pop | Energy Dissipation and Transport in Nanoscale Devices[END_REF][START_REF] Moore | Emerging Challenges and Materials for Thermal Management of Electronics[END_REF].

The efficiency of a thermoelectric material is characterized by its thermoelectric figure of merit ZT , which is proportional to the Seebeck coeficient S, the electrical conductivity σ and inversely proportional to the thermal conductivity κ.

ZT = S 2 σT κ (1)
In order to improve the conversion efficiency, ZT has to be increased. Consequently, κ should be reduced while S and σ should be increased. However, κ and σ are not independent parameters in common materials, and ZT optimization have remained very limited for decades (cf. [START_REF] Chen | Recent Developments in Thermoelectric Materials[END_REF]). ZT of Bismuth Telluride was then measured around 1. Unfortunately, the best thermoelectric materials happened to be rare and pollutant. This makes difficult their widespread use.

A thermoelectric device based on Si and Ge would be a better candidate, as they are abundant and compatible with the current fabrication technologies. However bulk Si and Ge ZT are too low to allow interesting thermoelectric applications.

Figure 1: Evolution of the maximum ZT over time. Reproduced from [START_REF] Heremans | When Thermoelectrics Reached the Nanoscale[END_REF].

The pioneering work of [START_REF] Hicks | Effect of Quantum-Well Structures on the Thermoelectric Figure of Merit[END_REF] opened new research directions by demonstrating that nanostructures could be used to significantly enhance the thermoelectric efficiency. Following this research trend, this work focuses on the reduction of the thermal conductivity, which is required for improving thermoelectric energy conversion (G. Chen 2005) by using nanostructuring.

Indeed, the characteristic length of electron transport, called the electron mean free path l e mf p , is usually much smaller than those involves in thermal transport referred as l p mf p h. Moreover, boundary scattering mechanisms become significant in nanostructures [START_REF] Ziman | Electrons and Phonons: The Theory of Transport Phenomena in Solids[END_REF]Peierls 1929) in which a reduction of the two previous characteristic lengths can be achieved. Thus, designing a nanostructure with an intermediate length L (l e mf p < L < l p mf p h) ) would have a greater impact on the phonon transport than on the thermal one. Finally, ratios of σ κ and thus ZT higher than its bulk counterpart can be achieved. Experimental measurements in Silicon nanowires of appropriate diameters [START_REF] Qiu | First-Principles Simulation of Electron Mean-Free-Path Spectra and Thermoelectric Properties in Silicon[END_REF] i.e. on the order of 100 nm, have already demonstrated a large reduction of the effective thermal conductivity [START_REF] Li | Thermal Conductivity of Individual Silicon Nanowires[END_REF][START_REF] Boukai | Silicon Nanowires as Efficient Thermoelectric Materials[END_REF]Hochbaum et al. 2008), down to 2 orders of magnitude lower than the bulk conductivity. As the related electrical conductivity reduction is weak, ZT values higher than 1 have been demonstrated (cf. [START_REF] Akiyama | Effects of Polytypism on the Thermoelectric Properties of Group-IV Semiconductor Nanowires: A Combination of Density Functional Theory and Boltzmann Transport Calculations[END_REF]).

Several other nanostructures have been investigated, e.g. nanofilms, superlattices [START_REF] Hu | Si/Ge Superlattice Nanowires with Ultralow Thermal Conductivity[END_REF][START_REF] Mu | Ultra-Low Thermal Conductivity in Si/Ge Hierarchical Superlattice Nanowire[END_REF], core-shell structures [START_REF] Hu | Significant Reduction of Thermal Conductivity in Si/Ge Core-Shell Nanowires[END_REF], phononic crystals [START_REF] Yu | Reduction of Thermal Conductivity in Phononic Nanomesh Structures[END_REF], nanowire networks [START_REF] Verdier | Thermal Transport in 2D and 3D Nanowire Networks[END_REF]. The thermal conductivity can also be reduced by alloying [START_REF] Fourier | Role of Disorder and Anharmonicity in the Thermal Conductivity of Silicon-Germanium Alloys: A First-Principles Study[END_REF] or partial amorphization [START_REF] Donadio | Atomistic Simulations of Heat Transport in Silicon Nanowires[END_REF] of materials.

Recently, [START_REF] Vincent | Novel Heterostructured Ge Nanowires Based on Polytype Transformation[END_REF] developed an original method for fabricating nanowires composed of a quasi-periodic stacking of cubic and hexagonal phases of Si or Ge. Images of a Ge polytype nanowire are shown in Figure 2 (reproduced from [START_REF] Vincent | Novel Heterostructured Ge Nanowires Based on Polytype Transformation[END_REF])), showcasing the two phases and their interface. In addition to the ZT improvement due to the size effect (as in other nanowires), multiple polytype interfaces could affect strongly the thermal conductivity. Meanwhile, the electrical conductivity should be less impacted by the change of phase.

Figure 2: Germanium polytype nanowire [START_REF] Vincent | Novel Heterostructured Ge Nanowires Based on Polytype Transformation[END_REF].

Nanoscale thermal transport

At low dimensions, the distribution of particles can strongly differ from the thermal equilibrium. The classical Fourier's law of conduction then becomes unable to describe the thermal transport, and nanoscale models are required. These models must include the material properties and the effects of nanostructures. They usually consider the heat either from the trajectories of individual atoms, or as collective excitations (i.e. phonons).

DFT (Density Functional Theory) simulations have been used to calculate the thermal properties of bulk materials [START_REF] Seko | Prediction of Low-Thermal-Conductivity Compounds with First-Principles Anharmonic Lattice-Dynamics Calculations and Bayesian Optimization[END_REF] Or small nanostructures [START_REF] Markussen | Heat Conductance Is Strongly Anisotropic for Pristine Silicon Nanowires[END_REF][START_REF] Li | Thermal Conductivity of Diamond Nanowires from First Principles[END_REF]. This ab initio method can be used to study materials or structures that are not previously characterized. However, full quantum simulations such as DFT require a lot of computing resources, and are thus limited to very small systems in the order of hundreds of atoms.

The Atomistic Green's Function (AGF) method can solve the transmission function in a device, but it is generally limited to ballistic transport (N. Mingo and Yang 2003) but can be extended to include all relevant scattering mechanisms [START_REF] Luisier | Atomistic Full-Band Simulations of Silicon Nanowire Transistors: Effects of Electron-Phonon Scattering[END_REF].

Among atomistic simulations, classical Molecular Dynamics (MD) simulation is a very common method to simulate the heat transfer. In this approach, the trajectories of atoms are numerically solved and analyzed. The materials are described by using empirical inter-atomic potentials. This classical approach of transport is accurate only at high temperature, higher than the Debye's temperature, that is 640 K in Si [START_REF] Madelung | Non-Tetrahedrally Bonded Elements and Binary Compounds I[END_REF].

MD simulations have been widely applied for modeling nanowires made of different crystalline structures [START_REF] Termentzidis | Modulated SiC Nanowires: Molecular Dynamics Study of Their Thermal Properties[END_REF] and of diameter up to 20 nm [START_REF] Volz | Molecular Dynamics Simulation of Thermal Conductivity of Silicon Nanowires[END_REF][START_REF] Donadio | Atomistic Simulations of Heat Transport in Silicon Nanowires[END_REF]. The effect of rough interfaces can be modeled directly in real space [START_REF] Liu | Effect of Surface Roughness on Thermal Conductivity of Silicon Nanowires[END_REF][START_REF] Qiu | Lattice Thermal Conductivity Reduction in Bi 2 Te 3 Quantum Wires with Smooth and Rough Surfaces: A Molecular Dynamics Study[END_REF] as well as the resulting reflections of mechanical waves. MD simulations can be performed at equilibrium (EMD) or non-equilibrium (NEMD) conditions. Previous authors have compared them in bulk Si [START_REF] He | Lattice Thermal Conductivity of Semiconducting Bulk Materials: Atomistic Simulations[END_REF][START_REF] Sellan | Size Effects in Molecular Dynamics Thermal Conductivity Predictions[END_REF][START_REF] Schelling | Comparison of Atomic-Level Simulation Methods for Computing Thermal Conductivity[END_REF] and Ge [START_REF] He | Lattice Thermal Conductivity of Semiconducting Bulk Materials: Atomistic Simulations[END_REF]. The conductance of Si/Ge interfaces has also been investigated [START_REF] Merabia | Thermal Conductance at the Interface Between Crystals Using Equilibrium and Nonequilibrium Molecular Dynamics[END_REF] as well as the conductivity of Si/Ge superlattices [START_REF] Termentzidis | Molecular Dynamics Simulations and Thermal Transport at the Nano-Scale[END_REF]. Due their computational consumptions, MD studies are generally limited to devices that are made of ten thousands of atoms.

Methods based on the resolution of the Boltzmann Transport Equation (BTE) for phonons are relevant for larger system and can cover the full temperature range. The linearized BTE can be solved analytically when simple phonon dispersions and scattering terms are considered as done in pioneering works of Callaway [START_REF] Callaway | Model for Lattice Thermal Conductivity at Low Temperatures[END_REF]) and Holland [START_REF] Holland | Analysis of Lattice Thermal Conductivity[END_REF]. When considering simple geometries, it can be solved for instance by a direct approach [START_REF] Nghiêm | Electro-Thermal Simulation Based on Coupled Boltzmann Transport Equations for Electrons and Phonons[END_REF], a discrete ordinate method [START_REF] Yang | Thermal Conductivity Modeling of Periodic Two-Dimensional Nanocomposites[END_REF] or a finite volume method [START_REF] Narumanchi | Comparison of Different Phonon Transport Models for Predicting Heat Conduction in Silicon-on-Insulator Transistors[END_REF].

But for 3D problems, a stochastic particle Monte Carlo method [START_REF] Jacoboni | The Monte Carlo Method for Semiconductor Device Simulation[END_REF]) is much more numerically efficient and complex scattering terms can be implemented [START_REF] Hamzeh | Monte Carlo Study of Phonon Dynamics in III-V Compounds[END_REF]. This versatile approach can solve accurately the BTE much beyond the linear approximation and in complex geometries. It has been used to study bulk materials [START_REF] Mazumder | Monte Carlo Study of Phonon Transport in Solid Thin Films Including Dispersion and Polarization[END_REF][START_REF] Péraud | Efficient Simulation of Multidimensional Phonon Transport Using Energy-Based Variance-Reduced Monte Carlo Formulations[END_REF][START_REF] Klitsner | Phonon Radiative Heat Transfer and Surface Scattering[END_REF][START_REF] Peterson | Direct Simulation of Phonon-Mediated Heat Transfer in a Debye Crystal[END_REF], as well as nanocrystalline materials [START_REF] Yang | Thermal Transport in Nanocrystalline Si and SiGe by Ab Initio Based Monte Carlo Simulation[END_REF], porous nanofilms [START_REF] Wolf | Thermal Conductivity of Silicon Nanomeshes: Effects of Porosity and Roughness[END_REF][START_REF] Madelung | Non-Tetrahedrally Bonded Elements and Binary Compounds I[END_REF][START_REF] Hao | Frequency-Dependent Monte Carlo Simulations of Phonon Transport in Two-Dimensional Porous Silicon with Aligned Pores[END_REF] nanowires (Y. [START_REF] Chen | Monte Carlo Simulation of Silicon Nanowire Thermal Conductivity[END_REF][START_REF] Lacroix | Monte Carlo Transient Phonon Transport in Silicon and Germanium at Nanoscales[END_REF], nanoribbons [START_REF] Mei | Full-Dispersion Monte Carlo Simulation of Phonon Transport in Micron-Sized Graphene Nanoribbons[END_REF], etc. An important improvement of the MC simulations is the introduction of variance reduction methods in [START_REF] Péraud | Efficient Simulation of Multidimensional Phonon Transport Using Energy-Based Variance-Reduced Monte Carlo Formulations[END_REF].

In many MC approaches, simple isotropic phonon dispersions [START_REF] Mazumder | Monte Carlo Study of Phonon Transport in Solid Thin Films Including Dispersion and Polarization[END_REF][START_REF] Lacroix | Monte Carlo Transient Phonon Transport in Silicon and Germanium at Nanoscales[END_REF] are assumed. Few works have used a full-band dispersion, i.e. have considered the phonon modes in the entire 3D Brillouin Zone [START_REF] Kukita | Monte Carlo Simulation of Thermal Conduction in Silicon Nanowires Including Realistic Phonon Dispersion Relation[END_REF][START_REF] Aksamija | Anisotropy and Boundary Scattering in the Lattice Thermal Conductivity of Silicon Nanomembranes[END_REF]. However, significant deviations from isotropic calculations can be obtained [START_REF] Kukita | Monte Carlo Simulation of Thermal Conduction in Silicon Nanowires Including Realistic Phonon Dispersion Relation[END_REF]. The accuracy of semi-empirical methods such as BTE or MD is directly related to the choice of their input parameters (inter-atomic potential, dispersion properties, phonon lifetimes, etc.). They are thus not convenient to study a new material for which such parameters are unknown. A methodology for solving the BTE in thin films without any adjustable parameter has been recently presented in [START_REF] Chaput | Ab Initio Based Calculations of the Thermal Conductivity at the Micron Scale[END_REF][START_REF] Yang | Thermal Transport in Nanocrystalline Si and SiGe by Ab Initio Based Monte Carlo Simulation[END_REF]. In this simulation method, the BTE parameters, i.e. full-band dispersion and full-band phonon-phonon scatterings rates, are preliminarily computed by using DFT calculation.

In MC simulation, the reflection at external rough boundaries is usually implemented by considering a specular reflection. Thus, external boundaries are designed in real space by using a realistic profile, e.g. using a saw-tooth shape [START_REF] Moore | Phonon Backscattering and Thermal Conductivity Suppression in Sawtooth Nanowires[END_REF] or a random surface [START_REF] Ramayya | Thermoelectric Properties of Ultrathin Silicon Nanowires[END_REF]. Other models consider a specific scattering term related to a diffusive reflection at the interface that randomizes the propagation direction of diffused phonons. Casimir has linked the thermal conductivity to the width of structures [START_REF] Casimir | Note on the Conduction of Heat in Crystals[END_REF]. More recently, a probability of specular reflection [START_REF] Berman | The Thermal Conductivity of Diamond at Low Temperatures[END_REF][START_REF] Ziman | Electrons and Phonons: The Theory of Transport Phenomena in Solids[END_REF]) using a diffusive characteristic length for ultra-thin wires [START_REF] Santamore | Effect of Surface Roughness on the Universal Thermal Conductance[END_REF] have been defined. In Soffer's work [START_REF] Soffer | Statistical Model for the Size Effect in Electrical Conduction[END_REF], this probability depends on the incident wave vector of phonons plus two empirical parameters: the surface roughness standard deviation and its tangential correlation. This model provides results close to experimental measurements [START_REF] Maurer | Universal Features of Phonon Transport in Nanowires with Correlated Surface Roughness[END_REF]. The interface transmission is described by two common models: the Acoustic Mismatch Model (AMM) and Diffusive Mismatch Model (DMM) [START_REF] Little | The Transport of Heat Between Dissimilar Solides at Low Temperatures[END_REF]. AMM considers perfect interfaces that conserve the wave coherence, while the DMM [START_REF] Swartz | Thermal Boundary Resistance[END_REF] supposes that the interfaces are completely diffusive. More advanced methods based on MD [START_REF] Merabia | Thermal Conductance at the Interface Between Crystals Using Equilibrium and Nonequilibrium Molecular Dynamics[END_REF], AGF [START_REF] Tian | Enhancing Phonon Transmission Across a Si/Ge Interface by Atomic Roughness: First-Principles Study with the Green's Function Method[END_REF] and ab initio lattice dynamics calculations [START_REF] Alkurdi | Critical Angle for Interfacial Phonon Scattering: Results from Ab Initio Lattice Dynamics Calculations[END_REF] have also been used.

Theoretical models

In this chapter, we detail the theoretical basis of our studies. The concepts of thermal conductivity and interface thermal conductance are presented.

The two first sections of this chapter presents the studied devices and the different kinds of material modeling used in this work. The next sections focus on the heat transfer modeling and present relevant semi-analytical approaches. Numerical methods to solve the heat transfer are presented in details in the next chapter.

Simulated devices

In this thesis, different types of nanostructures, schematized in Figure 3, are investigated: nanofilms, nanowires and heterojunctions. In nanofilms, the external interfaces, separated by a finite distance, are called the top and bottom interfaces. In nanofilms both the in-plane and cross-plane configurations are considered, depending on whether the thermal flux is parallel or perpendicular to the interfaces, respectively.

As schematized in Figure 3, to implement numerically the devices, a cubic mesh is used. A material (of arbitrary crystal orientation) is assigned to each cell. All devices are in between a cold thermostat (blue plane) and a hot one (red plane). A face of a cell can be either transparent (when the adjacent cell is made of the same material), specular (no adjacent cell) or diffusive (green planes).

The reflected angle of a particle colliding with a specular face is equal to the incident one. This reflected angle is of course different when a diffusive face is involved and must be selected carefully, as explained later in this chapter. To mimic infinite dimensions, as needed for instance in cross-plane configuration for all directions except the transport one, specular reflections are implemented at opposite boundaries.

In this study, the heat transport is along the X-axis. The four devices shown in Figure 3 have different types of boundaries. They are: (a) cross-plane nanofilm (CPNF), i.e. oriented along the cross-plane direction, with only specular boundaries;

(b) in-plane nanofilm (IPNF), i.e. oriented along the in-plane direction with a pair of specular (XZ planes) and diffusive (XY planes) boundaries (colored in green in Figure 3);

(c) rough nanowire (NW), with only diffusive boundaries.

(d) heterojunctions (HJ), with only specular boundaries, but composed of two materials separated in the center by an interface (colored in purple in Figure 3).

These devices are parametrized by their length L (distance between thermostats along X axis) and their width W (along the Z axis). In this study, only nanowires with a square cross section are considered (i.e. with the height along Y direction also equal to W ).

Models of materials

To solve the heat transport equation in a solid-state system, the prior knowledge of both the dispersion relation and the scattering rates is required.

The dispersion relation is the relation between the phonon energy and wave vectors, defined in the entire reciprocal space. In a periodic lattice, the reciprocal space is periodic and the dispersion relation can be completely defined in a primitive cell centered on the origin called the Brillouin zone (BZ).

The movement of phonons in a material is perturbed by several scattering mechanisms, which can be characterized by their scattering rates. Usually, the main source of thermal resistance is due to the phonon-phonon scatterings.

It is possible to have an analytical representation of the dispersion relation and the scattering rates by using several approximations. In this thesis, a full band description of the BZ is mainly used, i.e. both the dispersion relations and scattering rates are numerically calculated for a large number of wave vectors in the whole BZ.

Full band models

In order to study the details of thermal transport more accurately than the previous works, all phonon modes are considered, and the anisotropic effects are included in the present work. Thus the energy dispersion, the phonon velocity as well as the phonon scattering rates in all directions are considered.

Indeed, all phonon states belonging to the Brillouin Zone (BZ) are considered in our "full-band" description. To determine all these full-band material parameters, we use two different methods: the Adiabatic Bond Charge Model (ABCM) and the Density Functional Theory (DFT).

ABCM is a semi-empirical model derived from the Force Constant Model (FCM).

In FCM, the force matrices are defined between every neighboring atoms in the primitive unit cell of the lattice, ignoring all other interactions. The energy can thus be calculated for any wave vector from the equation of motion of the phonons.

The ABCM improves its accuracy by accounting for the delocalization of the electron cloud. The resulting negative charge called 'bond charge' accumulates at the center of the interatomic bond. Using the parameters of [START_REF] Weber | Adiabatic Bond Charge Model for the Phonons in Diamond, Si, Ge, and $\alpha$-Sn[END_REF] for Si and Ge, the resulting ABCM dispersions fit experimental measurements with an error of only 2%. The dispersion relations are calculated in this work for both the cubic and hexagonal phases of Si and Ge, as detailed in [START_REF] Larroque | Étude Théorique de L'anisotropie Du Transport Thermique Dans Les Nanostructures à Base de Silicium et de Germanium[END_REF].

DFT is an ab initio method which is relevant to investigate the phonon properties of a large range of materials. It is accurate and does not require any empirical parameter, but consumes heavy computing power. The dispersion relation and scattering rates for the cubic and hexagonal phases of Si (called Si3C and Si2H, respectively) that are used in this thesis have been computed by [START_REF] Chaput | Ab Initio Based Calculations of the Thermal Conductivity at the Micron Scale[END_REF], as explained in [START_REF] Togo | Distributions of Phonon Lifetimes in Brillouin Zones[END_REF].

Since the scattering rates cannot be calculated via the ABCM which considers by principle only the harmonic part of the potentials, the thermal conductivity can only be calculated by using DFT. However, the ABCM dispersion is sometimes used to calculate the interface thermal conductance since we do not have yet a complete set of DFT data for all Ge phases. For both methods, the first BZ is discretized in a set of N wave vectors, with N = 31 × 31 × 31 = 29791 and N = 31 × 31 × 19 = 18259 for cubic and hexagonal phases, respectively. Cubic and hexagonal phases have 6 and 12 phonon modes, respectively. The angular frequency, the group velocity and the phonon-phonon scattering rates are calculated for each discrete state characterized by a couple of a wave vector ⃗ q and a mode m.

The only phonon scattering rates computed here by DFT are those related to the intrinsic phonon-phonon scattering mechanisms that are dominant in bulk materials. These phonon-phonon scattering rates λ are calculated via the finite displacement method detailed in supplementary materials of [START_REF] Chaput | Ab Initio Based Calculations of the Thermal Conductivity at the Micron Scale[END_REF][START_REF] Togo | Distributions of Phonon Lifetimes in Brillouin Zones[END_REF], for 101 temperatures ranging between 0 and 1000 K. In brief, DFT calculations have been performed by using the finite displacement method with PAW pseudopotentials, within the PBE approximation. This approach allows obtaining the forces on atoms, from which the harmonic and anharmonic force constants can be extracted. The phonon frequency and group velocity can then be computed from the dynamic matrix.

During the simulation, the values corresponding to intermediate temperatures are interpolated by using a cubic spline method. We should note that Normal and Umklapp scattering mechanisms are not distinguished in our approach. Thus, our phonon-phonon scattering rates represent actually an average effect due to a combination of the two kind of mechanisms.

Due to the discrete nature of our BZ description, the definition of the iso-energy states that are used for instance to select the final state after any scattering event leads to a fluctuation of energy inside the discrete states due to their finite size. Thus, the strict conservation of angular frequency (or energy) is modulo the frequency discretization step ∆ω. In our simulation ∆ω is defined as ∆ω = ω max /128 where ω max is the maximum phonon frequency in the material (this value of ∆ω leads to negligible average energy loss as shown in §3.1.11).

The resulting iso-energy curves of the first phonon mode from DFT data are plotted in Figure 4. Parts (a) and (b) represent the BZ in face centered cubic and hexagonal lattices (for Si3C and Si2H, respectively) and their high symmetry points (from [START_REF] Setyawan | High-Throughput Electronic Band Structure Calculations: Challenges and Tools[END_REF]). Parts (c), (d), (e), and (f) are cartographies of the angular frequency in the main cutting planes of the BZ.

In the cutting plane ΓXL of cubic silicon in (c), iso-energy curves are far from being circular (spherical in 3D) as in isotropic materials in the whole frequency range. In particular, the anisotropy between the L and U points is strong. In the hexagonal phase in (d), (e), and (f), an isotropic behavior is nearly achieved in all planes (KΓM , M ΓA, and KΓA, respectively) but only at low frequencies.

In contrast, far from the G point the anisotropy becomes strong, in particular between the M and K points.

We can also see that the phonon group velocities, which are the derivatives of the angular frequency with respect to the wave vectors, will not be isotropic. Furthermore, in the optical modes, some regions of the BZ even show a group velocity with a direction opposite to the corresponding wave vectors. Similarly to the previous figure, the scattering rates in the main planes of Si3C and Si2H are represented in Figure 5. It shows an even stronger anisotropy in the dispersion relations. These figures illustrate the necessity of using a full-band representation of the BZ instead of a simple isotropic analytical formula.

Models of rough boundaries

In nanostructures, the contribution of rough boundaries on the thermal resistance needs to be considered. Indeed, this contribution can become stronger than those of phonon-phonon scattering mechanisms.

A first model of the rough boundary contribution in thermal resistance was introduced by [START_REF] Casimir | Note on the Conduction of Heat in Crystals[END_REF]. He considered all rough boundaries as black bodies at a constant temperature gradient. This means that all points in a structure emit as much energy as they absorb, i.e. they are perfectly diffusive boundaries. An average scattering rate within this model can be defined such as:

λ boundary = v W F (2)
where F is a form factor related to the shape of the device independent from the device width W , and v is the phonon group velocity.

This model was improved by the work presented in [START_REF] Berman | The Thermal Conductivity of Diamond at Low Temperatures[END_REF], by adding a probability p specular to have a specular reflection at the rough boundary (thus 1p specular is the probability to have a diffusive reflection). This probability is constant and the rough boundary scattering rates have the form:

λ boundary = v W 1 + p specular 1 -p specular (3)
Soffer further improved the probability of specular reflection in [START_REF] Soffer | Statistical Model for the Size Effect in Electrical Conduction[END_REF], making it dependent on the incident angle of the phonon, the rms of the surface height deviation (called the surface roughness parameter ∆), and the spatial correlation of these deviations.

In the case where the spatial correlation of the surface roughness is negligible, the probability of specular reflection at a rough boundary becomes:

p specular = e -(2cos(θ)∆|⃗ q|) 2 (4)
where ⃗ q is the incident wave-vector, θ its angle from the normal vector of the surface, and ∆ is the surface roughness parameter. [START_REF] Kazan | Thermal Conductivity of Silicon Bulk and Nanowires: Effects of Isotopic Composition, Phonon Confinement, and Surface Roughness[END_REF] were able to fit experimental measurements of the thermal conductivity of nanowires from [START_REF] Hochbaum | Enhanced Thermoelectric Performance of Rough Silicon Nanowires[END_REF] with this model.

Since this model considered isotropic materials, we needed to adapt them to our full-band dispersion relation (see §2.2.1) in which the velocity and wave-vector of a given phonon state are not necessarily collinear. For a particle in a state j, the incident angle θ is changed to:

cos(θ j ) = v ⊥ |⃗ v j | (5)
where v ⊥,j is the component perpendicular to the boundary interface of the velocity ⃗ v j of the incident phonon.

q ∥ being the component of the wave vector tangent to the surface, [START_REF] Chen | Thermal Conductance of Thin Silicon Nanowires[END_REF] set p specular to 1 when W q ∥ > 1 or to 0 otherwise. Also, [START_REF] Ramayya | Thermoelectric Properties of Ultrathin Silicon Nanowires[END_REF] models rough boundaries in real space using only specular reflections.

Heat transfer modeling and thermal conductivity

The thermal conductivity of a device is its main thermal property. In a nanostructure heated by two thermostats, it can be defined as:

κ = Q ∆T L ( 6
)
where Q is the net heat flux density (flux per unit of area) from the hot thermostat to the cold thermostat, ∆T is the temperature difference between the thermostats and L the distance between them.

Isotropic analytical models

In order to calculate bulk thermal conductivity, [START_REF] Callaway | Model for Lattice Thermal Conductivity at Low Temperatures[END_REF] first used an isotopic and linear phonon dispersion relation, with only one artificial mode averaging all actual mode properties. The relaxation times were energy and temperature dependent, fitted from experimental measurements. This simple model managed to fit correctly the experimental thermal conductivities of Si and Ge at low temperatures. [START_REF] Holland | Analysis of Lattice Thermal Conductivity[END_REF] improved this model by separating the transverse and longitudinal acoustic modes. A frequency threshold was also used to better model the properties of higher energy phonons. The resulting bulk Si and Ge conductivities are relevant in the whole temperature range.

For phonon transport, many works assumed such simple dispersion relations [START_REF] Mazumder | Monte Carlo Study of Phonon Transport in Solid Thin Films Including Dispersion and Polarization[END_REF][START_REF] Klitsner | Phonon Radiative Heat Transfer and Surface Scattering[END_REF][START_REF] Peterson | Direct Simulation of Phonon-Mediated Heat Transfer in a Debye Crystal[END_REF]Y. Chen et al. 2005;[START_REF] Lacroix | Monte Carlo Simulation of Phonon Confinement in Silicon Nanostructures: Application to the Determination of the Thermal Conductivity of Silicon Nanowires[END_REF]). However, they require several parameters to fit on experimental measurements, and thus do not have predictive capabilities to investigate materials that have not been previously investigated experimentally.

Besides, these approximations cannot reproduce accurately both the heat capacity and the heat conductivity (cf. the work of Larroque, Dollfus, and Saint-Martin (2017)) using a single normalization parameter in the whole temperature range.

Full-band semi-analytical models

Using the analytical dispersions relation and scattering rates mentionned in §2.2.1, the thermal conductivity can also be calculated with analytical formula.

It is usually derived from the linearized BTE (see §3.1.1).

In this thesis, semi-analytical models (SA) as well as numerical approaches (MC and MD methods presented in the next chapter) are used to estimate the thermal conductivity. Previous full band models have been developed for instance by Mingo et al. (2003), and [START_REF] Aksamija | Anisotropy and Boundary Scattering in the Lattice Thermal Conductivity of Silicon Nanomembranes[END_REF] demonstrated predictive calculations of the thermal conductivity with a full band dispersion relation. One originality of this work is to directly includes both dispersion relation and phonon-phonon scattering rates computed by using DFT calculations.

SA ballistic model

The ballistic SA model estimates the thermal conductivity in infinitely short crossplane nanofilms (CPNF). The resulting thermal conductivity is linearly dependent on the distance L between thermostats. The heat flux densities involved are emitted from the hot Q T hot and cold thermostats Q T cold (at temperatures T hot and T cold respectively). The fluxes are estimated within a ballistic Landauer's formalism, without any diffusive mechanism leading to a heat transmission equal to 1.

The heat flux density injected by a thermostat at the temperature T , normal to the vector ⃗ n is then:

Q T = V s ∑ state j |⃗ vj •⃗ n|>0 ℏω j |⃗ v j • ⃗ n|f BE (ω j , T ) (7)
where ω j and ⃗ v j are the angular frequency and velocity of the phonons in the state j, V s is the reciprocal volume of one state in the Brillouin zone, and f BE is the equilibrium distribution.

The net heat flux density is:

Q ballistic = Q T hot -Q T cold (8)
The resulting thermal conductivity that it is L-dependent has the following form:

κ ballistic = L V s ∆T ⎛ ⎜ ⎜ ⎝ ∑ state j ⃗ vj •⃗ n>0 ℏω j ⃗ v j • ⃗ n f BE (ω j , T hot ) + ∑ state j ⃗ vj •⃗ n<0 ℏω j ⃗ v j • ⃗ n f BE (ω j , T cold ) ⎞ ⎟ ⎟ ⎠ (9)
Usual approximations use the local temperature derivative of the phonon distribution instead of the temperature bias ∆T at the thermostats, leading to:

κ ballistic = L V s ∑ state j ℏω j |⃗ v j • ⃗ n| ∂f BE ∂T (ω j , T eq ) (10)
We found that these two approaches are not equivalent at low L. Similar expressions of the ballistic conductivity have been previously discussed for instance in Mingo et al. (2003) .

SA diffusive model

The SA diffusive model leads to a length-independent conductivity (as in an infinitely long CPNF) and should provide the bulk thermal conductivity. In this model, the conductivity is written as:

k diffusive = V s ∑ statej ℏω j |⃗ v j • ⃗ n| 2 1 λ j ∂f BE ∂T (ω j , T eq ) (11)
where λ j is the phonon-phonon scattering rate for state j.

The diffusive model is a solution of the linearized Boltzmann Transport equation (BTE). The steady-state BTE for each state j has the form:

⃗ v j • ⃗ ∇f j = ∂f j ∂t ⏐ ⏐ ⏐ ⏐ scatt (12)
where f j is the distribution of phonons in the state j.

The wave vector-dependent relaxation time approximation 1/λ j is used to simplify the expression the phonon-phonon scattering rate, as follows:

∂f j ∂t ⏐ ⏐ ⏐ ⏐ scatt = -(f j -f BE (ω j , T ))λ j (T ) (13)
By using this approximation, the system tends "naturally" to recover its equilibrium distribution.

We define f j as the sum of the Bose-Einstein equilibrium distribution plus a small deviation δf j .

f j = f BE (ω j , T ) + δf j (14)
The deviation from the equilibrium distribution is assumed to be much smaller than the equilibirum distrubtion and also independent from the temperature at the first order:

δf j ≪ f BE (ω j , T ) (15) ∂(δf j ) ∂T = 0 (16)
Besides, as the spatial grandient can be expressed as follows:

⃗ ∇f j = ∂f j ∂T ⃗ ∇T (17)
The left hand side of Eq. 12 can be written as:

⃗ v j • ⃗ ∇f j = ⃗ v j • ⃗ ∇T ∂f j ∂T = ⃗ v j • ⃗ ∇T ∂f BE ∂T (ω j , T ) (18)
Eq.12 becomes:

δf j = ⃗ v j • ⃗ ∇T ∂f BE ∂T (ω j , T ) 1 λ j ( 19
)
The heat flux density in the transport direction ⃗ n (i.e. colinear with ⃗ ∇T ) is defined by:

Q = ∑ states j Q j (20) Q j = V s ℏω j ⃗ v j • ⃗ n δf j (21)
We then introduce δf j :

Q j = V s ℏω j ⃗ v j • ⃗ n ( ⃗ v j • ⃗ ∇T ∂f BE ∂T (ω j , T ) 1 λ j ) (22) Q j = ℏω j |⃗ v j • ⃗ n| 2 1 λ j ∂f BE ∂T (ω j , T eq ) ⃗ ∇T • ⃗ n (23)
The conductivity of Eq. 11 is finally obtained from the definition of the thermal conductivity.

SA Matthiessen model

In order to model the thermal conductivities in devices with an intermediate length in which a heat transport regime between ballistic and diffuse limits occurs, we first mix the previous models by using a Matthiessen's rule. It means summing ballistic and diffusive conductivities (of Eq 11 and 9, respectively) as follows:

1

κ Matthiessen = 1 κ ballistic + 1 κ diffusive (24)

Spectral SA models: CPNF, IPNF and NW

More complex models based on a spectral decomposition of the thermal conductivities are usually more accurate. Thus, a specific model for each kind of studied nanostructures have been derived. Then, one model is dedicated to NanoFilms in Cross-Plane configuration (the CPNF model), another to NanoFilms in In-Plane configuration (IPNF) and the last one to nanowires (NW). In each case, the Soffer's model (see §2.3) is used to model the rough boundary scattering mechanism. In this approach, ballistic and phonon-phonon scattering rates are once again summed but their spectral dependences are considered. Indeed, the average distance over which a phonon moves along the transport direction ⃗ n before colliding with a thermostat is:

L = L 2 (25)
with L the distance between thermostats.

Similarly, in IPNF and NW models, the average distance over which a phonon moves in a transverse direction ⃗ n ⊥ before colliding with a rough boundary is

W = W ( 1 1 -p specular - 1 2
) (26)
with W the width of these devices, and by considering the probability of specular reflection p specular from Eq. 4.

For these three models (NFCP, NFIP and NW), the total spectral scattering term for a state j is given by:

λ CPNF,j = λ j + |⃗ v j • ⃗ n| L (27) λ IPNF,j = λ j + |⃗ v j • ⃗ n| L + |⃗ v j • ⃗ n ⊥,1 | W (28) λ NW,j = λ j + |⃗ v j • ⃗ n| L + |⃗ v j • ⃗ n ⊥,1 | W + |⃗ v j • ⃗ n ⊥,2 | W (29)
Finally, the associated thermal conductivities κ NFCP , κ NFIP , and κ NW are calculated from Eq. 11 by replacing λ j by their relevant scattering rates.

Interface thermal conductance modeling

The interface thermal conductance (ITC) is related to the temperature drop at an interface between two materials induced by a non-zero heat flux. Thus, the ITC is defined by:

G = Q inter ∆T inter (30)
where Q inter is the heat flux density flowing though the interface, and ∆T inter is the temperature drop at the interface.

It should be noted that the definition of ∆T inter is not trivial. The temperature at one side of the interface should be calculated (or measured) by considering only the incident phonons and not the total phonon distribution (cf. [START_REF] Swartz | Thermal Boundary Resistance[END_REF]). In the case of devices in (or close to) a diffusive phonon transport regime, these two different temperature definitions are equivalent, but they differ significantly in the case of a ballistic regime.

The two models mainly used for ITC investigation are the Diffusive Mismatch Model (DMM) and the Acoustic Mismatch Model (AMM). The full-band adaptations of AMM and DMM are detailed in [START_REF] Larroque | Phonon Transmission at Si/Ge and Polytypic Ge Interfaces Using Full-Band Mismatch Based Models[END_REF]. In this thesis, we only investigate the DMM, which assumes a bad quality of the interface and leads to a diffusive process for every phonon crossing the interface. Then, each incident phonon can be either transmitted or reflected according to a transmission coefficient t, but in all cases it will undergo a diffusive scattering process inducing a memory loss of the incident phonon properties.

Similarly to the SA ballistic model derivation for the thermal conductivity, the ITC can be theoretically estimated via a Landauer's approach. The net heat flux density at the interface between the materials A and B is determined by using a transmission coefficient t A→B (ω), which is considered to be only dependent on the angular frequency of the particles. The spectral impedance of a material A is defined for a range ∆ω as:

I A (ω) = V s ∆ω ∑ states j in material A ℏω A j ⃗ v A j • ⃗ n δ(|ω -ω A j | < ∆ω 2 ) (31)
The transmission coefficient is then calculated by using the impedance of the two materials (A and B) involved in the interface:

t A→B (ω) = I B (ω) I A (ω) + I B (ω) (32)
The final semi-analytical formula for the ITC can be written as follows:

G DMM (ω, T ) = ∂f BE ∂T (ω, T )I A (ω)t A→B (ω) (33)
where G DMM (ω, T ) is the spectral decomposition of the ITC at temperature T , and f BE is the equilibrium Bose-Einstein distribution.

Monte Carlo and Molecular Dynamics simulations

The main concepts related to the nanoheat transfer are reminded in the previous chapter. This chapter focuses on the numerical methods that have been used during this thesis to investigate the phonons in nanostructures with a particular interest for interfaces. First, an original Monte Carlo simulator partially developed during this thesis is presented, second the Molecular Dynamics simulation methods are reminded.

Monte Carlo method

The particle Monte Carlo method is a stochastic approach to solve the Boltzmann Transport Equation (BTE). In this thesis, we present a new simulator for MC phonon transport, which solves the BTE in 3D for both the real space and reciprocal space (full band approach). The required semi-empirical, i.e. the phonon dispersion and the phonon-phonon scattering rates, are parametrized in the entire 3D Brillouin zone by using ab initio DFT calculations (see §2.2.1).

After introducing the BTE, the model of the simulated particles is presented and the studied devices are described. Then, we detail our temperature definition and our complete algorithm. Our configuration method to adjust simulation parameters is explicated. Finally, we discuss our two phonon treatment of three phonon scattering processes, the assumption of the average elastic scatterings at rough boundaries.

Boltzmann Transport Equation

The Boltzmann Transport Equation (BTE) describes the time evolution of the phonon distribution function in phase space f j (⃗ r, ⃗ q, t), where j is the phonon state, ⃗ r the position in real space, ⃗ q the wave-vector and t the time. Its timedependent expression is given by:

∂f j ∂t + ⃗ v • ⃗ ∇f = ∂f ∂t ⏐ ⏐ ⏐ ⏐ scatt , ( 34
)
where ⃗ v j is the phonon group velocity (⃗ v j = ⃗ ∇ ⃗ q ω j , ω j is the phonon angular frequency) and

∂f ∂t ⏐ ⏐ ⏐ scatt
is the collision term.

Figure 6: Schema of a possible phonon trajectory.

Principles of a particle Monte Carlo simulation

In order to solve the BTE, we use the particle Monte Carlo method for phonon transport [START_REF] Péraud | Efficient Simulation of Multidimensional Phonon Transport Using Energy-Based Variance-Reduced Monte Carlo Formulations[END_REF]. In this stochastic approach, the trajectories of a large number of semi-classical particles are randomly selected. The trajectories for each particle are described as a conituous sequence of free flights (linear trajectories in real space without any change in the reciprocal space) and instantaneous scattering events (no change in real space but with a modification of the state in the reciprocal space). The scattering mechanism ending a free-flight can be either a phonon-phonon scattering or a collision with the device boundary. For each particle, the initial state, the duration of each free flight, the type and the effect of each scattering event are chosen randomly according to all relevant scattering rates. Figure 6 indicates a schematized particle trajectory in a device. The particle is first injected from a hot thermostat on the left contact. Then it undergoes several scattering processes in between its free flights. It finally exits through the cold thermostat on the right contact. At the end of a MC simulation, the phonon distribution is computed by summing over all the particles k belonging to the mode m at a given time:

f m (⃗ r, ⃗ q, t) = ∑ particles k δ(⃗ q -⃗ q k (t))δ(⃗ r -⃗ r k (t)) (35)
It should be mentioned that only the trajectories during the steady state regime are considered here. However, the transient response could be also investigated by using our MC simulator.

Model of pseudo particles

To reduce the particle number and thus the computational resources, in our model each simulated particle represents a packet of N ω phonons with a frequency ω. Considering phonon bundle with a constant number of phonons N ω = N (and thus particles of different energies E p = N ℏω according to ω) would make difficult the treatment of elastic scattering mechanisms [START_REF] Klitsner | Phonon Radiative Heat Transfer and Surface Scattering[END_REF]). Thus, the number of phonons in a bundle N ω at the frequency ω is tuned to keep constant the total energy E p (= N ω ℏω) carried by each simulated pseudo particle (i.e. a phonon bundle). This input parameter E p defines the energy resolution of the simulation.

Variance reduction and reference temperature

Only the deviation of phonon distributions from their equilibrium state are relevant to investigate heat transfer. We would like to simulate only the excess or lack of phonons with respect to their equilibrium distribution in order to limit strongly the number of studied phonons and thus the number of simulated particles. Hence, a reference temperature T 0 close to and usually below the actual temperature in the device is defined. Only the phonon bundles (the simulated particles) differing from the Bose-Enstein distribution at T 0 are simulated. Besides, each simulated particle gets a new parameter s equal to +1 (or -1)

to represent an excess (or a missing particle, respectively) with respect to the reference distribution at T 0 .

For instance, the density of simulated particles in a state j in equilibrium at temperature T is:

n eq j (T ) = ℏω j E p (f BE (ω j , T ) -f BE (ω j , T 0 )) V s (2π) 3 (36)
where f BE is the Bose-Einstein distribution, and V s the volume of a state in the reciprocal space.

This technique reduces both the simulation time and the numerical noise (see §3.1.12). All simulations presented in this work are based on this approach called "energy-based variance-reduced method" and have been developed by Péraud and co-workers [START_REF] Péraud | Efficient Simulation of Multidimensional Phonon Transport Using Energy-Based Variance-Reduced Monte Carlo Formulations[END_REF].

Effective temperature

The temperature in each device cell must be updated during the simulation, as the phonon-phonon scattering rates are temperature-dependent. The energy density-temperature relation E V (T ) at equilibrium is computed by using this formula:

E V (T ) = ∑ state j ℏω j ( f BE (ω j , T ) -f BE (ω j , T 0 ) ) V s (2π) 3 (37)
where V s is the reciprocal volume of the state, and T 0 is the reference temperature. This function is plotted in Figure 7. A linear behavior is observed above 200K (this study is at 300K). The effective local temperature in a device cell is obtained by computing the phonon energy density. Since every pseudo particle (phonon bundle) has a constant energy, the local energy density in the device cell c is proportional to the particle density and it is given by:

E V,c = E p V c ∑ particles k s k (38)
Where V c being the volume of the cell c, and s k the sign of the particle k.

Then, the energy-temperature relationship of Eq. 37 is numerically inverted. Finally, the local temperature is

T c = E -1 V (E V,c ).

Monte Carlo Algorithm

The main Monte Carlo algorithm is presented in Figure 8. First, the initial phonon distributions of the particles are generated ( §3.1.7). Two nested loops are then performed, one over time and the other one over particles. At every time step, the temperature is updated in each cell ( §3.1.5), and particles are injected from the thermostats ( §3.1.8.2). Then, the displacement of each particle is computed and their coordinates are updated ( §3.1.6.1). During a time step, the events that can interrupt a free flight are either a collision with a cell boundary ( §3.1.8.2) or a phonon-phonon scattering event ( §3.1.8). At the end of the simulation, some post-processing calculations are performed to compute the final results ( §3.1.8.3).

Time of free flight and transport

The free flight corresponds to the movement of a particle between two scattering events. The scattering events can be separated in standard (bulk material) scattering mechanisms such as phonon-phonon scattering, and phonon-boundary scattering mechanisms.

Assuming Poissonian processes (cf. [START_REF] Jacoboni | The Monte Carlo Method for Semiconductor Device Simulation[END_REF]), at the beginning of the free flight, the time before the next standard scattering event t scattering is randomly selected according to the following formula:

t scattering = - ln(n random ) λ j (T c ) (39)
where n random is a uniform random number in ]0;1] and λ j (T c ) is the total scattering rate from a phonon in a state j at a local temperature T c . These scattering rates λ j (T c ) are assumed to be constant during a time step. In the general case, λ j (T c ) is the sum of the scattering rates corresponding to all relevant scattering mechanisms that are assumed to be independent. In this work, only the phonon-phonon scattering rates contribute to λ j (T c ). They are calculated by using an ab initio approach as previously mentioned. This method was first developed for electron transport modeling [START_REF] Jacoboni | The Monte Carlo Method for Semiconductor Device Simulation[END_REF].

Besides, the time before the next boundary collision depends on the distance d between the particle and the boundaries along its transport direction. As the phonon velocity is constant during a free flight, we have

t boundary = min ( d x v x j , d y v y j , d z v z j ) (40)
The actual free flight duration t ff for the particle is thus limited by the first event that occurs, i.e.

t ff = min (t remaining , t boundary , t scattering ) ( 41
)
where t remaining is the remaining time before the end of the i-th time step δt for the particle. The interruption of a free flight by the end of a time step has no impact on the other scattering rates since they are Poissonian processes.

t remaining = (t i + δt) -t current,k (42)
At the end of the free flight, the position in real space ⃗ r k of the particle k with a velocity ⃗ v j is updated according to:

⃗ r k (t current,k + t ff ) = ⃗ r k (t current,k ) + ⃗ v j (t current,k )t ff (43)
It should be mentioned that the wave vector is not modified during a free flight (only scattering mechanisms following the free flight have an impact on the wave vectors). Finally, the time counter of the particle is updated. New sequences of free flights and scattering events are selected until the end of the timestep.

t current,k = t current,k + t ff (44)

Initial conditions

In a cell of volume V c , the initial number of particles N_{init,m,c} in a mode m is numerically calculated by using the equilibrium density n eq (ω j , T c ) defined in Eq. 36. Summing over all (m-mode) states, it yields:

N init,m,c = V c ∑ n eq j (ω j , T c ) ( 45
)
where E p and T c are the energy of a simulated particle and the local temperature, respectively.

For each initial particle, the angular frequency is selected according to a equilibrium distribution proportional to the volume of each iso-energy state. Then, its wave-vector (reciprocal space) is randomly and uniformly selected among the iso-energy states. Using this method ensures that the repartition of energy between modes is as exact as possible. The position in the cell (real space) is also randomly selected according to a uniform distribution. Finally, the parameter s (sign) of the particle is positive if the local temperature T c is higher than the reference temperature T 0 , otherwise the parameter s is negative.

Scattering mechanisms

The two main phonon scattering mechanisms in Si and Ge are phonon-phonon scattering and impurity scattering. In this thesis, we ignore the impurity scatterings: the simulated materials are considered isotopically pure. However this mechanism could be easily included in our code by using for instance the relaxation time defined in [START_REF] Holland | Analysis of Lattice Thermal Conductivity[END_REF].

Phonon-phonon scattering mechanism

While the phonon-phonon scattering is a three-phonon process (higher orders being ignored), we approximate them as a two phonon process as proposed in the work of [START_REF] Lacroix | Monte Carlo Transient Phonon Transport in Silicon and Germanium at Nanoscales[END_REF]. This greatly simplifies the algorithm, and since every particle has a fixed energy the conservation of energy is exact. This approximation can be considered as an average behavior of large number of phonons having a large number of collisions. A particle could undergo a phonon-phonon scattering event at the end of a free flight (i.e. t ff = t scattering ). Then, the memory of the initial phonon state is lost (cf. [START_REF] Mazumder | Monte Carlo Study of Phonon Transport in Solid Thin Films Including Dispersion and Polarization[END_REF]), the scattered pseudo particle (phonon bundle) is destroyed and the new one is randomly selected from an equilibrium distribution.

The probability of selecting a new state j is proportional to the equilibrium density of particles weighted by the interaction rate of that state (cf. [START_REF] Lacroix | Monte Carlo Transient Phonon Transport in Silicon and Germanium at Nanoscales[END_REF]):

p j ∝ λ j (T c )n eq j (T c ) (46)

Boundary conditions

The faces of each cell in a device can be a thermostat, a specular boundary, a rough boundary, or a diffusive interface (ignoring the transparent interfaces between adjacent cells of the same material).

Thermostats

Thermostats inject a constant flux of particles that depends on their temperature. The number of particles injected during a time step δt through a surface is:

N thermostat (T ) = ∑ state j |⃗ vj •⃗ n ⊥ |>0 ⃗ v j • ⃗ n ⊥ A dt n eq j (T ) (47)
where ⃗ n ⊥ is the unit vector normal to the thermostated face and A is its area.

As each particle keeps track of its own simulated time, we initialize all of them at the beginning of the time step and make them behave as if they were injected continuously. In this calculation, it should be emphasized that the displacement of phonons during the (short) time step are assumed to be smaller than the mean free path of phonons, i.e. the transport of phonons in the thermostat is assumed to be ballistic.

Another possible method to simulate thermostats is to add a cell behind the thermostated face. Every face of that "black body" cell is specular, so that the phonon distribution inside the cell of the thermostat is constant and is equal to the equilibrium distribution at the thermostat temperature. In this "black body" cell, when a particle collides the thermostated face connected to the device, a standard specular phonon reflection occurs, and its phonon distribution remains unchanged. In parallel, a duplicated particle with the same properties as the incident phonon (position and wave vector) is created and injected (transmitted) inside the device as if the thermostated face were transparent. This last method is more accurate as it simulates completely the transport of phonons inside the thermostat but is more computationally consuming. We have implemented both methods. They all produce the same heat flux using our common simulation parameters. The direct injection at surfaces was chosen, as it is less computationally intensive.

Specular boundary

At smooth boundaries, the particle reflection is always specular, i.e. the wavevector component normal to the surface boundary of the reflected particle (q ′ ⊥0 ) is the reverse of the incident one (q ′ ⊥ = -q ⊥ ). As such specular reflections have no impact on the thermal flux parallel to the interface, they are used in our simulations to emulate semi-infinite boundaries. It should be mentioned that the implementation of specular reflection is obvious within a full-band description only if the orientation of the specular boundary corresponds to a high symmetry plane of the crystal.

Reflection at rough interface

In the case of a collision with a rough boundary, the particle has a given probability to undergo a diffusive reflection that randomizes the final wave vector instead of a specular reflection. In this work, we use the probability of specular reflections derived from the work of [START_REF] Soffer | Statistical Model for the Size Effect in Electrical Conduction[END_REF]. Its adaptation to our full-band model in detailed in §2.3 If a diffusive reflection occurs, the final state is selected among the states belonging to the same iso-energy states and with a relevant orientation of the final velocity in the case of a reflection (i.e. ⃗ v ⊥ • ⃗ v ′ ⊥ < 0 where ⃗ v and ⃗ v ′ are the incident and reflected velocities).

Finally, the probability of each reflected iso-energy state is weighted by two factors: the complementary probability of a specular reflection for a reflected state j ′ with an angle θ j ′ (from Soffer's model) and the related normal component of group velocity (as in the Lambert's law) as follows:

p diffusive,j ′ ∝ (1 -p specular,j ′ )|⃗ v j ′ |cos(θ j ′ ) (48)
This angular distribution of the final state allows the angular distribution conservation of the heat flux and leads to a net flux equal to zero at the steady state along the direction normal to the interfaces. This prevents an unphysical phonon accumulation near the rough boundary.

Transmission at semi-transparent interface

During a collision with a semi-transparent interface, the particle can be either reflected or transmitted. In our code, the phonon transmission probability t A→B is implemented by using the DMM model (see §2.5) assuming fully diffusive interfaces. Thus reflected phonons as well as transmitted ones undergo an elastic diffusion (the memory of the incident phonon is lost). Our implementation of semi-transparent interface is very similar to those of rough boundaries: a new state is selected among the iso-energy states with a valid orientation of the velocity. In order to ensure self-consistency, the transmission probabilities are calculated for each direction (A → B and B → A), with the same angular frequency step ∆ω for the iso-energy states. The probability for selecting a state is weighted by the probability of the particle being reflected, and the component of its velocity normal to the interface.

p DMM,j ′ ∝ (1 -t A→B (ω j ′ )) |⃗ v j ′ |cos(θ j ′ ) (49)

Post-processings

During the simulation, at each time step and in each cell, the local temperature T c is calculated as well as the local thermal heat flux density ⃗ J c [W m -2 ] by summing over the contribution of all particles. The thermal conductivity of the simulated device is then calculated from the average heat flux density ⃗ J c along the transport direction ⃗ n, as defined in Eq. 6.

The confidence interval at 95% calculated for all simulations was found smaller than 1 W m -1 K -1 except in long devices (L = 100µm) for which the precision is reduced to achieve reasonable simulation times.

Criteria for simulation parameters

There are several simulation parameters that must be selected correctly to ensure convergence and accuracy of the simulation. An automatic process is used to define these parameters. First, we need to estimate the thermal relaxation time of the simulated device. The temperature in a cell approximately relaxes to its equilibrium value according to a negative exponential with a constant thermal relaxation time τ :

∆T c (t) ∝ ∆T c (0)exp(- t τ ) (50) ∆T c (t) = T c (t) -T c (∞) (51)
Since there are two thermostats, the furthest cell is at L 2 . This cell is related to the relaxation time of the device according to:

τ = ∂E V ∂T (T ) ( L 2 ) 2 κ (52)
where E V is the energy density, L is the length of the device and κ is an estimation of the thermal conductivity.

In Figure 9, the thermal relaxation time given by MC simulations of a CPNF is plotted as well as the results of equation 52 using the estimated thermal conductivities obtained via the ballistic, diffusive and CPNF SA formula (from §2.4).

In order to extract the relaxation time from simulation, an exponential fit is performed on the temperature of every cell as a function of time. These relaxation times are then extrapolated to the middle of the device, giving the results shown in the figure.

We thus confirm that the thermal relaxation time τ calculated from the SA CPNF thermal conductivity is a good approximation of the behavior of the MC simulations.

10 -9 10 -8 10 -7 10 -6 10 -5 10 -15 Using this approximation, the time step duration δT is set to τ 20 , and the temperature and heat flux are averaged every 5τ . This allows the temperature of the cells to relax progressively towards their steady-state values.

Finally, to choose a relevant particle energy E p , the temperature fluctuation due to one particle displacement in a volume V c is chosen to remain below δT (typically equal to 0.01K) leading to:

E p = δT ∂E V ∂T V c ( 53
)
This requirement on the precision of the temperature in each cell of the device results in an approximately constant number of simulated particles per cell, independent from their dimensions. As we modeled every device with 20 cells, all our simulations contained an average number of 20 000 coexisting particles.

For a typical device of size 1 µm × 100 nm × 100 nm in Si3C, with a 4 K-temperature difference between thermostats and a reference temperature T 0 = 295 K, our selected parameters were E p ≈ 4 10 -18 J and δt ≈ 1 ns, and the simulation lasted 20 minutes on a single thread. Besides, for devices of length ranging between 1 nm and 10 µm, the timestep δt was tuned in the interval between about 0.1 ps and 100 ns.

Validation of two phonon approximation

During the MC simulation, the phonon-phonon scattering mechanism which is actually a three-phonon process is approximated by a two-phonon process. This approximation must conserve, at least on average, the frequency distribution of the phonon gas. By following the frequency change before and after each phononphonon scattering event for each simulated particle, we plot in Figure 10 the time evolution of these cumulated changes in a nanowire in the quasi-ballistic regime (L = 1µm, W = 100nm, ∆ = 0.5nm). The two phonon approximation appears to be consistent as the cumulated change fluctuates and remains negligible with respect to the total energy in the system (2.6 10 -18 J vs. 8.1 10 -14 J, respectively).

Validation of rough boundary scattering

Due to the discretization of the BZ, the elastic scattering due to rough boundaries cannot exactly conserve the angular frequency of simulated phonons. While the energy E p of the simulated particles is exactly conserved in our algorithm, a frequency change equal to the frequency discretization step ∆ω is possible at each scattering event. If it is not balanced, this change could impact the final frequency distribution of phonons.

To check this point, we summed these frequency changes caused by rough boundary scatterings over a complete simulation of a nanowire working in the quasi-ballistic regime (i.e. L = 1µm, W = 100nm, ∆ = 0.5nm). The time evolution of this cumulative change,plotted in Figure 11, shows a constant increasing trend. However, the total energy density corresponding to this drift during a 300 ns long simulation corresponds to an equivalent heat flux of approximately 666 W m -2 , which is negligible as it is five orders of magnitude lower than the energy transfer carried by the actual heat flux of 200 M W m -2 .

Validation of deviational formulation

For all MC simulations in this study, only phonons that are outside the distribution at a reference temperature T 0 are simulated in order to optimize simulation time and accuracy. In order to validate this method, we simulated a CPNF (L = 1µm, W = 100nm, ∆T = 4K) with different T 0 . The resulting temperature profiles are plotted in Figure 12. The temperature profile is less noisy when T 0 increases. Moreover, the simulation time typically decreases from 3 hours for T 0 = 0K to 10 minutes for T 0 = 300K.

T 0 = 0K T 0 = 200K T 0 = 290K T 0 = 300K

Molecular Dynamics method

The Molecular Dynamics methods (MD) involve a classical simulation of atom trajectories. As each atomic position is given at the initialization, it is possible to model accurately the interfaces between materials by this technique.

In a MD simulation, each atom has a defined mass, position, and velocity. A time dependent Newton's second equation is numerically solved to get the atomic trajectories. The instantaneous forces of their interactions are calculated from semi-empirical potentials (depending on atomic positions). Finally, all relevant thermodynamic properties can be obtained from a statistical analysis of the system, e.g. the temperature is derived via:

∑ atoms i 1 2 m i |⃗ v i | 2 = 3 2 N k B T ( 54
)
where m i and ⃗ v i are the mass and velocitiy of the atom i, N is the number of atoms and k B is the Boltzmann constant.

This technique is quite versatile but has some limitations. Firstly, these methods require a lot of computing power and then only small systems of thousands of atoms are manageable. In order to avoid diffusing phonons at the boundaries, periodic conditions are generally applied but the value of the phonon mean free path with respect of the size of the system or the unit cell has to be carefully considered. Secondly, the number of time steps, i.e. the actual simulation time, is commensurate and the time scale, or the frequency scale, is restricted by an empirical cut-off parameter. Finally, we should keep in mind that it is a classical simulation that cannot account for the energy quantization of phonons. MD simulations are thus limited to high temperature regimes above the Debye temperature of the material that can be relatively high in the case of standard semiconductors.

In this thesis, the MD simulation software used was LAMMPS [START_REF] Plimpton | Fast Parallel Algorithms for Short-Range Molecular Dynamics[END_REF].

For interatomic potential for both Si and Ge, the Stillinger-Weber [START_REF] Stillinger | Computer Simulation of Local Order in Condensed Phases of Silicon[END_REF] and Tersoff [START_REF] Tersoff | Modeling Solid-State Chemistry: Interatomic Potentials for Multicomponent Systems[END_REF]) set of parameters, that includes phonon anharmonicity and 3-body interactions, have been used.

Several MD methods exist to study the thermal conductivity of materials. We detail both EMD and NEMD approaches and their advantages and drawbacks will be discussed in following chapters. The calculation of both the thermal conductivity of bulk materials, and the interface thermal conductance (ITC) is detailed.

In every method, periodic boundaries are applied to the system along all dimensions. For simulations with interfaces, this means that two interfaces are present in the system.

Structure and energy relaxation

At the beginning of the simulation, each atom is at its position of minimal potential energy and the distribution of their velocity is Gaussian. Therefore, before any calculation, the system must be relaxed to its equilibrium state.

For both EMD and NEMD, two steps are followed. A first run is performed in the isothermal-isobaric (NPT) ensemble at zero pressure to relax the structure of the system. Since the thermal conductivity is calculated in the micro-canonical (NVE) ensemble, a second run is then performed in this ensemble to ensure that the energy distribution is at its equilibrium.

The Si/Ge interfaces studied in §6 are perfect interfaces, and have a constant strain in each material. The lattice parameter in the direction parallel to the interface is the average of the bulk Si and Ge lattice constants. In the normal direction, the lattice parameter relaxes in both materials to maintain their density [START_REF] Landry | Thermal Boundary Resistance Predictions from Molecular Dynamics Simulations and Theoretical Calculations[END_REF]. The polytype Ge interfaces are not strained, as there is no change of lattice parameter.

Bulk thermal conductivity with EMD

In EMD [START_REF] Sellan | Size Effects in Molecular Dynamics Thermal Conductivity Predictions[END_REF][START_REF] Schelling | Comparison of Atomic-Level Simulation Methods for Computing Thermal Conductivity[END_REF]Volz and Chen 2000), no external perturbations are applied during the simulation, and the thermal conductivity is calculated from a statistical study of atom properties with the Green-Kubo formula invoking the fluctuation-dissipation theorem.

During the simulation, the heat flux is calculated from the energy, stress, and velocity of each atom.

⃗ Q(t) = 1 V ( ∑ i E i ⃗ v i - ∑ S i ⃗ v i ) ( 55
)
where E i is the potential and kinetic energy of atom i, ⃗ v i its velocity, S i its stress tensor and V the volume of the system.

The method to obtain the stress tensor of each atom is detailed in [START_REF] Thompson | General Formulation of Pressure and Stress Tensor for Arbitrary Many-Body Interaction Potentials Under Periodic Boundary Conditions[END_REF].

The average of the autocorrelation of this heat flux is integrated according to the Green-Kubo formula [START_REF] Volz | Molecular Dynamics Simulation of Thermal Conductivity of Silicon Nanowires[END_REF], giving a thermal conductivity expressed by:

κ(t corr ) = V k B T 2 ∫ tcorr 0 ⟨ ⃗ Q(τ ) • ⃗ Q(0)⟩ 3 dτ ( 56
)
where T is the equilibrium temperature.

To explore the whole phase space of the system, several independent simulations are performed [START_REF] Chen | How to Improve the Accuracy of Equilibrium Molecular Dynamics for Computation of Thermal Conductivity?[END_REF]. For each simulation, a black line in Figure 13 shows the evolution of the thermal conductivity as a function of the upper bound of the integral t corr . Finally, an average (red line) is performed over all simulations. The error bars represent the confidence interval at 95%, calculated with:

IC 95% = [ x -2 σ √ n ; x + 2 σ √ n ] ( 57
)
where x is average, σ the standard deviation and n the number of samples.

The value of the thermal conductivity is the average at t corr = t cutoff . At high correlation times, the integrals diverge due to the accumulation of numerical errors (as shown in Figure 13). t cutoff is chosen after the first convergence (plateau of the red line) [START_REF] Li | Atomistic Modeling of Finite-Temperature Properties of Crystalline $\beta$-SiC[END_REF] and before the increase of the standard deviation.

Interface thermal conductance with EMD

The interface thermal conductance (ITC), instead of the previous bulk thermal conductivity, can be calculated with a similar method [START_REF] Chalopin | Thermal Interface Conductance in Si/Ge Superlattices by Equilibrium Molecular Dynamics[END_REF].

The heat flux across the interface between two volumes X and Y is:

Q inter (t) = 1 2A ⎛ ⎝ ∑ i∈X ⃗ v i • ⃗ F Y i - ∑ j∈Y ⃗ v j • ⃗ F X j ⎞ ⎠ (58)
where ⃗ v i is the velocity of atom i velocity, ⃗ F X i the interaction force on atom i from atoms in area X, and A the cross-section area of the interface.

This formula requires that the interaction potential can be calculated separately for each atom pair, and summed accordingly in order to obtain the ⃗ F X i . While this is trivial with the Stillinger-Weber potential which is defined by a two-body and a three-body components, this is not the case for the Tersoff potential.

As explained in §3.2.1, due to the periodic boundaries of the system, there are two interfaces. The area A used for the normalization must account for them.

The thermal conductance is then calculated from the Green-Kubo formula, with a normalization slightly different from Eq. 56:

G(t corr ) = A k B T 2 ∫ tcorr 0 ⟨Q inter (τ ) • Q inter (0)⟩dτ (59)
The results are then processed in the same way as the bulk method: the result is the average of several simulations, and the error is estimated from the confidence interval at 95%.

Bulk thermal conductivity with NEMD

The second kind of method that we study is NEMD. In opposition to EMD methods, their procedures are similar to experimental measurements. Perturbations are applied to the system, in our case the temperature in some regions of the system is controlled by rescaling the kinetic energy of atoms.

These two "external" thermostats induce a temperature gradient in the system and an internal thermal flux that can be calculated. The thermal conductivity is then directly obtained from its definition in Eq. 6.

Since the size of the system is small compared to the mean free path of phonons, the transport is partially ballistic between the thermostats. This means that there is a direct dependence between the calculated conductivity κ and the distance between thermostats L. In Figure 14, the relation between their

inverses 1 κ = f ( 1 L ) is shown in blue.
The error bars are the confidence intervals at 95%. A linear extrapolation (red line) at 1 L = 0 gives the result for bulk conductivity [START_REF] Schelling | Comparison of Atomic-Level Simulation Methods for Computing Thermal Conductivity[END_REF]. The green lines are linear extrapolations of the maximum and minimum errors.

0 5 • 10 -2 0.1 0.15 0.2 2 • 10 -2 4 • 10 -2 6 • 10 -2 8 • 10 -2 0.1 0.12 0.14 Inverse of distance between thermostats [nm -1 ] Inverse of thermal conductivity κ [mKW -1 ]
Error extrapolations Average extrapolation NEMD simulations Figure 14: Inverse of thermal conductivity as a function of the inverse of the distance between thermostats, for Si at 1000K. Results from NEMD simulations (blue line). The red and green lines are the extrapolations at 1 L = 0, for the average and the maximum and minimum errors, respectively.

Interface thermal conductance with NEMD

To calculate the ITC, each thermostat is placed in a different material as represented in Figure 15. White and grey volumes describe two different materials, for instance Silicon and Germanium. The red and blue hashed volume are respectively the hot and cold thermostats. In addition to the thermal gradient inside the materials, there is a temperature drop ∆T at the interfaces due to the Kapitza resistance. In order to estimate it, we calculate the average temperature in slices of the simulation box, along the transport direction. The resulting temperature profile is shown in Figure 16 where the white and grey backgrounds represent the different materials, and the red and blue hashed areas represent the hot and cold thermostats. A linear interpolation of the temperatures far from the interfaces and thermostat is used to extrapolate the temperature at both sides of the interface, giving a precise estimation of ∆T . The ITC is then calculated from its definition in Eq. 30.

Only the property of the interface is taken into account by this temperature drop ∆T inter . In contrast to the bulk conductivity calculations in §3.2.4, the ballistic transport between the thermostats and the interface is not related to length dependence of the ITC. Thus, the calculated value does not need to be extrapolated from several sizes of simulation boxes.

An issue with the NEMD methods is that the rescaling of atom kinetic energy by the thermostats induces an out of equilibrium phonon distribution in the thermostats. Since at our simulation scales the transport is partly ballistic, the equilibrium distribution is never recovered between the thermostats, or at the sides of the interfaces. The definition of temperature is then unclear.

Bulk materials and thin films

This chapter focusses on the thermal properties in both bulk materials and thin films. Molecular Dynamics (MD) simulations have been performed in order to test the different methods and benchmark them against previous results. Meanwhile, our Full Band Monte Carlo (MC) simulator based on ab initio parameters were used to study the anisotropy of Si and Ge in several crystalline phases. By using these advanced numerical simulations, several semi-analytical (SA) models of the thermal conductivity in a quasi-ballistic regime are discussed.

Bulk thermal conductivities by using MD

This section is dedicated to MD simulations, and bulk thermal conductivities of Si and Ge are calculated. For Silicon both the cubic and the hexagonal phases are investigated.

The thermal conductivities are calculated by using both Equilibrium and Non Equilibrium (EMD and NEMD) methods. Stillinger-Weber (SW) and Tersoff (Ter) potentials are used and their related results are discussed. The time step is 1 fs. For each simulation, two preliminary relaxation steps (see §3.2.1) during 40 ps (in simulation time) are performed. The heat flux is calculated by averaging over the next simulated 120 ps. The final thermal conductivity is an average over 20 independent runs. For EMD, we found that a simulation of 8-unit cells along each dimension is the minimum size to compute size independent conductivities.

The results from Si and Ge at T = 1000K are presented in Table 1 andTable 2. As T = 1000K is higher than the Debye temperatures of these materials, that are equal to 640K and 374K, respectively, a MD approach is relevant. Besides, in this temperature regime the phonon mean free path is lower than the system length and the computed thermal conductances are size independent.

For the sake of comparison, previous MD results extracted from [START_REF] Sellan | Size Effects in Molecular Dynamics Thermal Conductivity Predictions[END_REF] and [START_REF] Schelling | Comparison of Atomic-Level Simulation Methods for Computing Thermal Conductivity[END_REF] are also shown. Besides, experimental values from [START_REF] Capinski | Thermal Conductivity of Isotopically Enriched Si[END_REF]) and [START_REF] Ruf | Thermal Conductivity of Isotopically Enriched Silicon[END_REF] for isotopically pure Si (as considered in our simulations) as well as [START_REF] Ozhogin | Isotope Effect in the Thermal Conductivity of Germanium Single Crystals[END_REF] for pure Ge are indicated. These values given at 1000K are derived from an extrapolation from measurements at lower temperatures. Both the hexagonal phases of Si and Ge show a slight reduction in thermal conductivity with respect to their cubic counterpart. Although the hexagonal phases are anisotropic (as we will see later), it should be noted that in EMD, a pseudo isotropic thermal conductivity is calculated because the method averages the contribution of all directions.

κ W m -1 K -1 Exp.
The thermal conductivities are calculated in both cubic and hexagonal Si (resp. named Si3C and Si2H) with the SA diffusive formula (see Eq. 11), at T = 1000K like in MD and T = 300K, in At 1000K, we find that the Si3C results with SA underestimate both the experimental thermal conductivities and the MD ones. We confirm that the hexagonal phases have a reduced conductivity, even more than in MD. A strong anisotropy between Si2Hx and Si2Hz is also observed. One should keep in mind that these values are derived without considering any phonon-impurity scattering mechanism (with just phonon-phonon scattering) and thus this value underestimates the experimental measurements of isotopically pure Silicon [START_REF] Ruf | Thermal Conductivity of Isotopically Enriched Silicon[END_REF]. However, it is a good estimation of the thermal conductivity of natural Silicon at ambient temperature. Hence, no other scattering mechanism is taken into account in our study of bulk Si.

Ballistic to diffusive transition in thin films

By using Monte Carlo simulations, this section studies the transition of thermal properties between bulk material in which the thermal transport is diffusive and cross-plane nanofilms (CPNF) in which the thermal transport is quasiballistic. Several semi-analytical (SA) models are compared to Monte Carlo (MC) simulations. Both the cubic and hexagonal phases of Si are investigated. The angular frequency distribution of both the energy and heat flux are carefully studied. Finally, in order to explore the anisotropy of thermal transport, the angular distribution of the heat flux is checked.

All the following calculations are performed at an average temperature of T = 300K, with a temperature difference between the thermostats of ∆T = 4K, and a reference temperature of T 0 = 295K for the MC algorithm (see §3.1.4).

Comparison between models in cubic Si

In this section and the following, we investigate the influence of the length L on the thermal conductivities κ of CPNF calculated with several methods:

• MC simulations,

• SA ballistic model (cf. Eq. 9). It is linearly dependent on L and corresponds to the short film asymptotic behavior, • SA diffusive model (cf. Eq. 11). It is this independent of L and corresponds to the long film asymptotic behavior, • SA Matthiessen model (cf. Eq. 24). It is a kind of average between the 2 previous models that is widely used to investigate CPNF, • SA CPNF model (cf. Eq. 27). This model is a frequency dependent evolution of the previous model.

Figure 17 shows for cubic Silicon (Si3C) at 300K the thermal conductivities along the [100] lattice orientation. At long L, all the models (but the ballistic one) converge to the diffusive value of thermal conductivity:

κ diffusive = 138 W m -1 K -1
, which corresponds to the bulk conductivity. Likewise, except the diffusive model, they asymptotically converge at short L to the conductance of the ballistic model. This confirms that the models are equivalent at these limits.

We should note that for device length L shorter than 10nm, i.e. at the atomic scale, the considered phonon dispersion relation is not relevant and the indicated MC results are just a guide for the eyes illustrating the asymptotical behavior.

There is a gradual transition between the two asymptotic behaviors. This transition is often modeled in the literature by using the Matthiessen model. However, we find that this modeling exhibits a strong discrepancy with MC results at intermediate length. For instance, at L = 200nm, the SA model overestimates the conductivity by 60%. Meanwhile, the SA CPNF model is much more accurate leading to an error lower than 4% for all film lengths. Thus, the SA CPNF requiring much less computational resources than the complex MC simulations reproduces the MC results with a very good accuracy. This shows that a spectral treatment of the heat flux is necessary to correctly model the transition between ballistic and diffusive regimes.

10 -9 10 -8 10 -7 10 -6 10 -5 10 -4 10 0

10 1 10 2 Device length L [m] Thermal conductivity κ [W m -1 K -1 ] SA Ballistic SA Diffusive SA Matthiessen SA CPNF MC

Comparison between cubic and hexagonal phases of Si

The evolution of the Si thermal conductivity κ in the cubic phase (Si3C, in blue) is plotted in Figure 18 similarly as in the previous section and this evolution is compared to those in the hexagonal phase. Two lattice orientations of the hexagonal phase are considered: [10-10] (Si2Hx, in green) and [0001] (Si2Hz, in red). While all the curves have the same behaviors, Si2Hx and Si2Hz always exhibit lower κ than cubic Si. For instance, at the diffusive limit, i.e. in long films, κ diffusive = 138, 100 and 74W m -1 K -1 for Si3C, Si2Hx and Si2Hz respectively. Besides, these results show an anisotropy between the two hexagonal phase orientations. The anisotropy leads to a κ diffusive 26% higher in the Si2Hx.

We also confirm that the SA CPNF provides impressively accurate results for all these materials. 10 -9 10 -8 10 -7 10 -6 10 -5 10 -4 10 0

10 1 10 2 Device length L [m] Thermal conductivity κ [W m -1 K -1 ] Si3C Si2Hx Si2Hz
Figure 18: Thermal conductivity in CPNF as function of device length L, in Si3C (blue), Si2Hx (green) and Si2Hz (red).

Spectral distributions of energy and heat flux

The previous sections show that the SA CPNF model provides good estimates of the thermal conductivity due to its accurate spectral distribution of the thermal conductivity and thus of the heat flux. In order to study the spectral thermal properties, the average spectral energy and heat flux in a MC simulation are investigated in this section.

In a MC simulation at the steady state, the spectral energy density E MC (ω) of phonons in the nanofilm is computed. In a device cell (i.e. between ⃗ r and ⃗ r + ⃗ dr), the time spent by each particle in all energy states (i.e. between each ω and ω + δω) are summed. This average occupation distribution of each state f MC,j gives the average spectral energy density E MC (ω) as:

E MC (ω) = E p V c ∑ states j f MC,j δ ( |ω -ω j | < ∆ω 2 ) ( 60
)
where V c is the volume of the considered device cell.

At equilibrium, the theoretical spectral energy density E th (ω) is given by the Bose-Einstein distribution. We should note that the distribution at the reference temperature T 0 = 295K is removed as in the MC simulations. It yields:

E th (ω) = V s (2π) 3 ∑ states j ℏω j ( f BE (ω j , T ) -f BE (ω j , T 0 ) ) δ ( |ω -ω j | < ∆ω 2 ) (61)
where V s is the reciprocal volume of a state in the BZ.

Similarly, the spectral heat flux density is calculated by summing for each state the product of their occupation and their related velocity in the transport direction ⃗ n.

Q MC (ω) = E p V c ∑ states j f MC,j ⃗ v j • ⃗ n δ ( |ω -ω j | < ∆ω 2 ) (62)
For the sake of comparisons, in the SA CPNF model, the spectral heat flux density is given by:

Q CPNF (ω) = ∆T L V s (2π) 3 ∑ states j ℏω j (⃗ v j • ⃗ n) 2 1 λ CPNF,j ∂f BE ∂T (ω j , T )δ ( |ω -ω j | < ∆ω 2 ) (63)
We show in Figure 19 the angular frequency distribution of the phonon energy in CPNFs with an intermediate length (at L = 1µm) in Si3C and Si2H. The spectral energy is shown separately for each phonon mode. Crosses represent the MC values and lines are for the theoretical values from the SA CPNF model. We can see that there is no discrepancy between Si2Hx and Si2Hz spectra as the energy is not dependent on the lattice orientation.

We observe that the MC and theoretical approaches give the same result in terms of energy. This means that the particles in the MC simulation are not strongly out of equilibrium, as it is expected for these simulation parameters (length and temperature bias).

In Figure 20, the integral of the angular frequency distribution of the heat flux is plotted for 1 µm long CPNFs made of Si3C, Si2Hx and Si2Hz. The contributions of each mode are indicated with different colors. The integral of the distribution is used here instead of the spectral distribution of the heat flux as it is less noisy and also indicates directly the total contribution of each mode (at the highest frequency).

The evolutions of all spectral distributions are similar to the SA CPNF model, especially in the case of Si3C. However, the total contributions of each mode differ. Moreover, for some modes at high angular frequency (i.e. in the optical modes), a negative contribution to heat flux can be observed in MC results that are not present in the SA model. For instance, in the case of Si2Hz, the 6th mode has a clear negative contribution. We explain this phenomenon by the fact that the diffusive approximation used in SA CNPF model only account for positive velocities, as only their absolute values are used. These 'negative' velocities are particularly present in optical modes at higher angular frequencies and corresponds to the observed negative contributions. Thus, total heat fluxes computed by the SA model are higher than these computed by using MC. This is consistent with the lower SA thermal conductivities previously discussed. The total contributions of each mode to the heat flux are also indicated in Table 5, (a) for Si3C, (b) for Si2Hx and (c) for Si2Hz. The modes are numbered by ascending energy values. As expected, the heat flux is mainly carried by acoustic modes: 96% for Si3C, 61.1% for Si2Hx (with 33.4% in 4th to 6th modes), 66.4% for Si2Hz (with 32.0% in 4th to 6th modes). In Si3C, the 5th and 6th (TO) modes have a negative net contribution. However, the absolute value of the contribution is in the order of magnitude of the margin error of the MC simulation. In Si2H, modes higher than the 8th mode carry a negligible amount of heat.

(a)

Mode

Heat

Angular distribution of heat flux

Our full band description of the phonon dispersion in each material allows to account for the anisotropy of the thermal conductivity. In order to estimate it, the angular distribution of the heat flux is plotted in Figure 22.

In Figure 22 (a), (c), and (d), we represented the 'angular density of state' (aDOS) of Si3C, Si2Hx and Si2Hz, respectively. In our convention, the heat transport direction in a CPNF is always set towards the positive X axis. It is important to note that this angular analysis is performed in real space (not in the wave vector space). We consider a unit sphere, and discretize it along its spherical coordinates, i.e. the polar and azymuthal angles θ and ϕ. The definition of these angles is reminded in Figure 21. For each solid angle (real space direction), we sum the contribution of all the discrete states in the BZ (reciprocal space) of the material that have this phonon velocity orientation i.e. toward that solid angle, around the corresponding θ and ϕ. To summarize, the aDOS is simply the number of states (per unit volume) in each solid angle, normalized by their area S ∆θ,∆ϕ .

aDOS(θ, ϕ) = V s (2π) 3 S ∆θ,∆ϕ ∑ states j δ ( |θ -θ j | < ∆θ 2 ) δ ( |ϕ -ϕ j | < ∆ϕ 2 ) (64)
We observe a high aDOS in the main lattice directions of each phase:

• for Si3C, <110> are the highest, followed by <100> and <111>.

• for Si2H, <10-10> and <0001> have high aDOS, followed by <21-30>.

The aDOS in the hexagonal phase is higher along its main orientations than in its cubic counterpart as there are only 14 main directions in Si2H (6 <10-10>, 6 <21-30> and 2 <0001> directions) lower than the 22 main directions in Si3C (6 <100>, 8 <110> and 8 <111> directions).

In Figure 22 (b), (d), and (f), the angular distribution of the heat flux density is shown in Si3C, Si2Hx and Si2Hz, respectively. The average occupation of states during MC simulations f MC,s is used to calculate the heat flux. Differently from the process presented for the aDOS, to compute the angular heat flux, the (negative or positive) contribution of the state will depend on the sign of the X component of its velocity.

Q MC (θ, ϕ) = E p V s S ∆θ ∆ϕ ∑ states j f MC,j ⃗ v j • ⃗ n δ ( |θ -θ j | < ∆θ 2 ) δ ( |ϕ -ϕ j | < ∆ϕ 2 ) ( 65
)
The figures showing the angular heat fluxes are relatively noisy especially along the directions in which the aDOS is low, as only few MC particles are present. In the three cases, as expected, the main contribution is oriented directly towards the transport direction (θ = 0 and ϕ = 0). However, in Si3C, the heat flux along the <111> directions is higher than along the <110> directions even though their aDOS is lower. Similarly, in Si2Hx, the <21-30> directions contribute more than the <10-10> ones. In an isotropic system, the angular distribution of the flux is expected to be a smooth cosine function (as in the Lambert's distribution). But in both cubic and hexagonal phases, the angular distributions exhibit peaks revealing that the heat flux is mostly transported along specific orientations. Each peak and their related orientations are labeled in the figures.

Angular distribution of heat flux in an isotropic material

In our full band model, state discretization in the reciprocal space that is not very fine (N = 31 × 31 × 31 for Si3C), in particular near the Γ point. In order to check the impact of our mesh on the computed angular heat flux, we have created an artificial isotropic material having the same dispersion along all directions. For this material, the angular frequencies for the wave vectors in the main direction [100] of Si3C is used as a reference for all the 6 modes. Then for each state in the BZ, its angular frequency is computed according to its distance from the Γ point (i.e. the norm of their wave vector). The norm of the velocities is calculated similarly and the direction of the velocity is parallel to the wave vector. Likewise, the isotropic phonon-phonon scattering rates are interpolated along the [100] direction.

The angular DOS and the angular distribution of heat flux computed in a 1µm long CPNF in this isotropic material is shown in Figure 24. The aDOS in Figure 24 (a) has a slight bias in the main direction [100] and in the planes between them, but it is clearly negligible compared to the features visible in the aDOS in previous materials in Figure 22 (a), (c) and (e). In Figure 24 (b), the angular heat flux is very diffuse in the whole hemisphere.

We thus confirm that our discretization of the BZ is sufficient to capture the anistropy of the heat transport in the studied materials.

Rough nanostructures

While the previous chapter details the thermal properties of thin films in the cross-plane configuration, here we consider thin films in the in-plane configuration (IPNF) with 2 rough interfaces and nanowires (NW) with 4 rough interfaces (the surrounding external boundaries), as defined in §2.1. In this chapter, we study the effect of rough boundaries in nanostructures by using our MC method including the phonon-rough interface scattering mechanism (cf. §3.1.8.2). As done in §4 for CPNFs, the evolutions of thermal conductivity, and the spectral and angular distributions of the heat flux are studied. The accuracy of the SA model to model rough IPNFs and NWs is discussed. By using a fully diffusive model for rough boundaries in the MC code, the thermal transport along several lattice orientations is investigated and compared with experimental measurements.

Effect of rough boundaries in IPNFs and NWs

Device geometry

In IPNFs and NWs, the device geometry is mainly characterized by 2 parameters:

• their length L, or distance between thermostats along the X axis, • their width W , or distance between the rough boundaries (cf. Figure 3). This section shows the dependence of the thermal conductivity on these parameters.

In Figure 25, the evolution of the thermal conductivity κ with the device length L is shown for IPNFs and NWs of width W = 100nm and roughness parameter ∆ = 0.5nm. For the sake of clarity, it should be mentioned that the figure is a semilog plot that differs from log-log axes used for similar figures presented in the previous chapter. For L = 100µm, κ that is equal to 138 in CPNFs is reduced down to 87 and 69 W m -1 K -1 in IPNFs and NWs respectively. As expected, the degradation of the thermal conductivity is directly related to the number of rough boundaries in long devices in which the heat transport is diffusive. However, in ultra-short devices in which the heat transport is ballistic, the thermal conductivity evolutions are the same whatever the number of rough interfaces (2 in INPF and 4 in NW). More interestingly for devices having rough boundaries, the semi-analytical models (SA), indicated in solid lines, do not fit MC results for IPNFs and NWs as well as in the case of CPNFs. For L = 10µm, the CPNF, IPNF and NW SA results are 2%, 5% and 12% lower than the corresponding MC result, respectively. And for L = 1µm, the CPNF, IPNF and NW SA results are 3%, 8% and 14% lower than the corresponding MC result, respectively.

10 -9 10 -8 10 -7 10 -6 10 -5 10 -4 For the same 3 different devices with a length L = 1µm and ∆ = 0.5nm, Figure 26 shows the evolution of κ as a function of the device width W . As all devices have a length of L = 1µm, an intermediate heat transport regime occurs mixing ballistic and diffusive transports. For W = 1nm, κ is 16 and 6 W m -1 K -1 for IPNF and NW, resp. The increase of κ with the increase of W in IPNF and NW exhibit a similar behavior. For intermediate width of 1µm, κ reaches 95% and 90% of the CPNF conductivity in IPNF and NW, resp. For devices wider than 1µm in which W > L, the impact of rough boundaries is weak, and the cross-plane conductivity is recovered with a difference lower than 5%. For thinner devices i.e. for W in the range of 10 nm to 200 nm, the width of the device impacts on the thermal conductance by a factor of about 30 W m -1 K -1 per decade (in W ) for both IPNFs and NWs.

Effect of the surface roughness parameter

Figure 27 shows the evolution of κ as a function of the surface roughness empirical parameter ∆ for a device of length L = 1µm and width W = 100nm. Two plateaus can be observed. For ultra-small values of ∆ lower than 0.1 nm, the diffusive reflections are negligible and then the cross-plane thermal conductivity using specular interface reflections for phonons is recovered. For a value of ∆ higher than 1 nm, the conductivities of NWs and IPNFs reach their minimum. This minimum is related to a fully diffusive regime in which all phonon reflections at the external boundaries are diffusive (with a randomized reflected angle). It should be noticed that the overall trends of these Monte Carlo results are reproduced by the relevant SA models. Nevertheless, these SA models systematically underestimate the MC results.

Impact on the spectral and angular distributions of the heat flux

We have performed the same analyses as in §4.2.3 on NWs. The spectral distribution of the thermal energy in a NW follows the equilibrium distribution as shown in Figure 19 for CPNFs.

Figure 28 compares the spectral heat flux in NW and CPNF in Si3C devices with L = 1µm, W = 100nm and ∆ = 0.5nm. The contribution of each mode is plotted with a specific color, solid lines show the results in CPNF while dashed lines stand for NW. We observe a strong reduction of the heat carried via the acoustic modes of the NW with respect to the case of CPNF in accordance with the reduction of the thermal conductance. However, the heat contribution of phonons belonging to the TO mode is higher in NW than in CPNF even if the absolute value remains weak. With respect to the CPNF distribution shown in Table 5 (only differing by the rough boundaries), we observe a heat flux reduction of 49%, 39% and 36% in NW made of Si3C, Si2Hx and Si2Hz, respectively. Indeed, we can explain this difference, as the heat flux in CPNF of Si2H is anisotropic and more intense towards the transport X direction. Besides, the contribution of each mode to the heat flux is indicated in Table 6 (a), (b) and (c) for nanowires made of Si3C, Si2Hx and Si2Hz, respectively. The heat flux reduction ratio is not similar for each phonon mode and depends on the material. In Si3C, the most reduced contribution is related to that of the 2nd and 3rd mode and in Si2Hx the contribution of the 1st, 2nd and 4th ones.

(a)

Mode

Fully diffusive rough boundaries

Previously in this chapter (cf. Figure 27), we have seen that if the phonon boundary scattering mechanism (Soffer's model) is used with a roughness parameter ∆ higher than 1 nm, a plateau of the thermal conductivity is reached. This limit is called the fully diffusive limit of the phonon boundary scattering, i.e. for this high value of ∆ almost all reflections at the interface are diffusive, i.e. all the reflected angles of reflected phonon are randomized.

In order to mimic such very rough external boundaries, a model called "fully diffusive" (FD) has been implemented. Thus, the FD model is equivalent to the previously used Soffer's model but with a zero probability of having a specular reflection (p specular ). To benchmark Soffer's model with ∆ = 0.5nm model and the FD one, the two resulting conductivities in Si3C nanowires are plotted in Figure 30 as a function of the length L (for W = 100nm). A significant discrepancy between the two models can be observed only for L higher than 500 nm. This indicates that in the FD, phonons have an average mean free path around 500nm. The FD model is thus relevant in short NWs.

Rough boundaries and crystalline orientations

Our full-band approach for phonon dispersion allows a priori the study of arbitrary crystalline orientations. Nevertheless, only one crystalline orientation has been investigated for Si3C (and two for Si2H) in the previous sections. Indeed, our implementation of a specular reflection (although it is the standard one) requires the existence of a final state having a wave vector with a negative perpendicular component with respect of the wave vector of the incident phonon (q ⊥0 = -q ′ ⊥0 ). This condition is always fulfilled only if the orientation of the specular boundary is oriented along a high symmetry plane of the crystal. This limitation is not present in the fully diffusive (FD) model since the final state is randomly chosen among all the available final iso-energy states. We are thus able to study the thermal conductivities of Si3C nanowires for several other crystal orientations. In Figure 31, the thermal conductivity as a function of device length is plotted for [100], [110] and [111] lattice orientations by using the FD model. While the conductivities are similar in all directions for devices smaller than a few µm, at higher lengths different values are achieved, revealing some anisotropy in the heat transport. The diffusive limits along the [111] and [110] directions are 6% and 15% lower than the limit for the [100] direction, respectively.

Once again, the semi-analytical models appear quantitatively disappointing in long devices when the fully diffusive approximation is used, nevertheless they are able to capture the qualitatively the trends of the orientation effects.

10 -9 10 -8 10 -7 10 -6 10 -5 10 -4

0 20 40 Device length L [m] Thermal conductivity κ [W m -1 K -1 ] SA Si3C [100] SA Si3C [110] SA Si3C [111] MC Si3C [100] MC Si3C [110] MC Si3C [111]

Comparison with experimental mesurements

Several experimental measurements of thermal conductivity in IPNFs and NWs have been reported in the literature. In this section, we compare works from Ju [START_REF] Ju | Phonon Scattering in Silicon Films with Thickness of Order 100 Nm[END_REF], Liu [START_REF] Liu | Thermal Conductivity Measurements of Ultra-Thin Single Crystal Silicon Layers[END_REF], Li [START_REF] Li | Thermal Conductivity of Individual Silicon Nanowires[END_REF] and Hochbaum [START_REF] Hochbaum | Enhanced Thermoelectric Performance of Rough Silicon Nanowires[END_REF] with results from our Monte Carlo simulations.

In Figure 32 (a) and (b), we compare the calculated and measured conductivities in Si3C IPNF and NW, respectively. In all these experiments, the device length L is about 1 µm. In the case of IPNF, the simulated thermal conductivities fit the experimental data of [START_REF] Liu | Thermal Conductivity Measurements of Ultra-Thin Single Crystal Silicon Layers[END_REF] and slightly underestimate the conductivity from [START_REF] Ju | Phonon Scattering in Silicon Films with Thickness of Order 100 Nm[END_REF]. For nanowires, our results are close to [START_REF] Li | Thermal Conductivity of Individual Silicon Nanowires[END_REF] ones and the thermal conductivity relationship with the width are relatively well reproduced, except for the shortest width W = 22nm. Hochbaum's results [START_REF] Hochbaum | Enhanced Thermoelectric Performance of Rough Silicon Nanowires[END_REF] show a conductivity that is an order of magnitude lower. They explain this difference by the strong roughness on the NW boundaries due to the fabrication method. However, even our fully diffusive model fails to reproduce these results. [START_REF] Liu | Thermal Conductivity Measurements of Ultra-Thin Single Crystal Silicon Layers[END_REF], [START_REF] Ju | Phonon Scattering in Silicon Films with Thickness of Order 100 Nm[END_REF], [START_REF] Li | Thermal Conductivity of Individual Silicon Nanowires[END_REF] and [START_REF] Hochbaum | Enhanced Thermoelectric Performance of Rough Silicon Nanowires[END_REF]).

In the case of Li's measurements on NW [START_REF] Li | Thermal Conductivity of Individual Silicon Nanowires[END_REF], we further compare the temperature dependence of the thermal conductivity in Figure 33. Only three widths W are considered: 37nm, 56nm and 115nm, as the 22nm device is too small to be correctly simulated using the bulk materials parameters. Our Monte Carlo simulation reproduces correctly the evolution of the conductivity with the temperature. The differences may be due to other sources of thermal resistance in the experimental setup, such as contact resistance. 70 50 100 150 200 250 300 [START_REF] Li | Thermal Conductivity of Individual Silicon Nanowires[END_REF].

Temperature T [K] Thermal conductivity κ [W m -1 K -1 ] W = 115 nm W = 56 nm W = 37 nm

Thermal transport modeling at solid-solid interfaces

This chapter studies the thermal transport at interfaces between semi-conducting materials. One of the goals of this thesis is to characterize the properties of polytype interfaces. Another one is to implement an advanced interface model in a Full-band Monte Carlo simulator (MC) based on data from Molecular Dynamics (MD) simulations, which accurately model the position of atoms at the interfaces.

After presenting the studied interfaces, we determine the simulation parameters in MD to estimate the interface thermal conductance (ITC) at a Si/Ge interface.

The implementation of Diffusive Mismatch Model (DMM) to describe the phonon transmission at interfaces in the MC simulator is then tested, and the resulting ITC is compared to the semi-analytical (SA) DMM formula. Using both methods, we investigate the ITC of several interfaces, and their temperature dependence.

Finally, we study some spectral properties of the heat flux at the interfaces.

Simulated interfaces

Several interfaces between Si and Ge in two different phases are considered. As previously, the cubic phases are refered as '3C' and the hexagonal phases '2H'. A lattice orientation of [100] is used when nothing else is specified, '2Hx' and '2Hz' correspond to [10-10] and [0001], respectively.

First the interfaces between cubic Si & Ge in the same lattice orientation are investigated:

• Si3C/Ge3C (stands for Si3C within a [100] orientation in contact with a Ge3C within the same orientation)

• Si3C/Ge3C [110] • Si3C/Ge3C [111]
Then the polytype interfaces:

• Si3C/Si2Hx • Si3C/Si2Hz • Ge3C/Ge2Hx, refered as "aligned" polytype interfaces (see Figure 34 In the case of DMM, pseudo-interfaces between identical materials are also investigated:

• Si3C/Si3C • Si3C/Si3C [110] • Si3C/Si3C [111]
The mismatched polytype Ge interface is the kind of interfaces mainly observed in the Ge polytype nanowires fabricated by [START_REF] Vincent | Novel Heterostructured Ge Nanowires Based on Polytype Transformation[END_REF].

In Figure 35 (d), the influence of the length L of the two materials on each side of the interface is investigated. From 2 to 8 unit cells, the ITC is inversely propotionnal to the length. This dependence is due to the periodicity of the superlattice, and is related to a significant ballistic transport between the interfaces. The ITC is almost constant for L higher than 10 cells.

NEMD parameters

The NEMD method is described in §3.2.5. Our average temperature are higher in NEMD (T =1000K) than in the EMD (T =400K) in order to reduce the mean free path of phonons. This allows to sudy systems with smaller lengths between the thermostats and interfaces leading to lower computational times.

The default simulation parameters chosen for our NEMD simulations are:

• a relaxation time of 0.5 ns,

• a thermalization time with thermostats of 5 ns,

• a data collection (500,000 timesteps of 0.5fs) time of 1 ns, • a simulated system width W of 4 unit cells, • a simulated system length of 32 unit cells,

• an average temperature of 1000K,

• a temperature bias ∆T of 100K between thermostats.

The related temperature profile between thermostats is plotted for T =1000 K and ∆T = (a) 10K, (b) 50K and (c) 100K in Figure 36. For ∆T = 10K, we see that the temperature fluctuations are so strong with respect to the temperature gradient at the interface, that it leads to unphysical results. Temperature biases higher than 50 K are required in this structure. The ITC is extracted by using the temperature regression shown in red lines to compute the effective temperature bias at the interface. The ITC for several values of the simulation parameters in NEMD are presented in Figure 37.

In Figure 37 (a), we see the evolution of ITC with the temperature difference between the thermostats. The results reach a plateau above 100K.

In Figure 37 (b), we see that a simuation time longer than 0.4ns (800,000 timesteps) is required.

In Figure 37 (c), several box widths W are tested. Contrary to the case of EMD (cf. Figure 35 (c)), the thermal conductance is higher for small W . For the following NEMD simulations, a width of 4 unit cells is selected (N.B. 8 unit cells for EMD).

In Figure 37 (d), the length L of materials on each side is changed. This is the most important point as it directly impacts the computation time. As the ballistic transport is important in these nanometer scale devices, a strong dependence of length L can be observed. At the diffusive limit a plateau appears.

In [START_REF] Landry | Thermal Boundary Resistance Predictions from Molecular Dynamics Simulations and Theoretical Calculations[END_REF], the diffusive limit is reached in a structure with 400 unit cells between thermostats at 500K. As their system is not periodic in the transport direction, it would correspond in our case to a structure having 800 unit cells for each material. Due to the limitation of our computational ressources, the longest simulated devices have 512 unit cells. With 512 unit cells at 1000K, the diffusive limit is not reached in our simulations. This is consistent with the previously shown MC results of nanofilms (cf. §4).

In bulk NEMD simulations, the thermal conductivy is extrapolated by using the Mathiessen theoretical model presented in §2. In our NEMD simulation the ITCs are L-dependent. However, as far as we know, no L-depend model have been developped for the ITC. This is an open issue as the NEMD method is used in the next section to estimate the spectral transmission of phonons at interfaces.

Spectral distribution of the interfacial heat flux

This section studies the spectral distribution of the interfacial heat flux computed by using NEMD simulations.

The spectral NEMD heat flux is calculated from the forces and velocities of atoms around the interface. This method is detailed in [START_REF] Sääskilahti | Role of Anharmonic Phonon Scattering in the Spectrally Decomposed Thermal Conductance at Planar Interfaces[END_REF] in the case of an interface between simple materials (i.e. between Argon and an artificially heavy Argon). It is applied here to interfaces between Si and Ge.

The spectral interfacial heat flux Q(ω) is the sum of the heat fluxes between every pair of atoms i and j across the interface between the materials A and B. The cut-off distance of interatomic potentials limits the number of pairs of atoms included in the calculation.

First, an inter-atomic correlation function K ij is defined as:

K ij (t 1 -t 2 ) = 1 2 ⟨ ⃗ F ij (t 1 ) • (⃗ v i (t 2 ) + ⃗ v j (t 2 )) ⟩ (66)
where ⃗ F ij is the force applied from atom j on atom i, and ⃗ v their velocities. The averaging is performed over the simulation time.

The spectral heat flux can then be calculated from its Fourier transform Kij (ω) by using:

Q i→j (ω) = 2 A Re[ Kij (ω)] (67)
Finally, the spectral conductance can be recovered:

G(ω) = Q(ω) ∆T = 1 ∆T ∑ i∈A,j∈B Q i→j (ω) (68)
Figure 38 shows the resulting normalized spectral heat flux at a Si3C/Ge3C and a mismatched Ge3C/Ge2H interfaces.

These spectral distributions are significantly different. The polytype Ge interfaces have two strong peaks due to the acoustic and optical modes. The heat flux is mainly transmitted via the acoustic modes of phonons. The transmission at the Si3C/Ge3C interface is more uniform. Nevertheless, the second peak is higher and the contribution of optical modes are more important that the acoustic ones.

They both have a similar cutoff in angular frequency around 70 rad.THz. The angular frequency range is limited by the available phonon states in Ge. However, in the case of Si3C/Ge3C, a small spectral heat flux is observed above the cutoff, around 85 rad.THz, which is attributed to anharmonic interactions at the interface. In MD simulations, a detailed atomic model of the interfaces is implemented.

The Full Band DMM model (which is used in MC simulations) provides a simpler physical model of the interface. In Figure 39, we compare the interfacial heat flux of Si3C/Ge3C interfaces computed by using these two methods. The DMM result shows three peaks while only two peaks are present in MD. The DMM spectral flux is higher at lower frequencies. Besides, the DMM cutoff frequency slightly different than in MD as the phonon dispersion relation in MD is not exactly the same. As expected, there is no contribution above the cutoff within the DMM method. Figure 40 shows the the DMM spectral transmission coefficients in a Si3C/Ge3C interface, and a polytype Si3C/Si2Hx interface. As in the previous figure, the frequency cutoff is around 60 rad.THz for Si3C/Ge3C, while the polytype Si interface cutoff is around 95 rad.THz. Since the transmission coefficients are calculated by using the DMM impedance of materials, all the simulated polytype interfaces have their transmission coefficients oscillating around t = 1 2 . For Si3C/Ge3C, it is notable that the first and third peaks in the spectral heat flux correspond to high transmission values, while the middle peak is only related to a low transmission t = 0.2.

ITC by using MC simulations

To study heterojunctions and superlattices by using Full Band Monte Carlo simulation of phonons, the interface transmission is derived from DMM as explained in §3.1.8.2. This DMM spectral transmission is used to validate our implementation of the heterostructure transmission.

The first test is performed with a Si3C/Si3C interface i.e. a constant transmission of t A→B = t B→A = 1 2 , and a Si3C/Si2Hx interface. An heterojunction without any thermostat (i.e. without any phonon injection) is simulated with different initial temperatures on each side of the interface. As expected, after a transient regime, the temperature profile relaxes to the same equilibrium value, and the net heat flux at the interface becomes null.

Then thermostats are added to recover the same configuration as those used in the two previous chapters, i.e. with a hot thermostat and a cold thermostat along the X axis. The temperature profile at steady-state for a Si3C/Si3C heterojunction (HJ) is shown in Figure 41. The red and blue diamonds show the temperature of the thermostats, and the grey vertical line represents the position of the interface.

For small devices, i.e. L = 100nm, the transport is in the intermediate transport regime between ballistic and diffusive ones. As expected, we observe a temperature drop at the contacts with the thermostats, and also at the inner interface. The ITC is calculated as a first attempt by simply dividing the average heat flux in the device and the temperature difference between the two sides of the interface.

In Figure 42, the ITC is calculated in heterojunctions of different lengths L by the mentioned approach. The ITC is expected to be independent of L. This occurs in the ballistic regime but the ITC seems to be around twice the values by using the SA formula. Besides, when L increases, the ITC decreases below the expected value. To recover the expected result, the temperature difference to consider must be different. Only the population of incident phonons on the interface must be included in the computation of the temperature. Indeed, when a strong outof-equilibirum transport regime occurs, this last approach gives a significantly different temperature drop.

In order to calculate a correct ITC, we simulate very small HJ (L = 1nm) without any phonon-phonon scattering mechanism. The incident heat flux on each face is directly the heat flux injected from thermostats in which temperatures are known. The thermostat temperature can thus be used to calculate the temperature drop at the interface ∆T .

The results are given in Table 7 for several interfaces at 300K, and compared to results from the SA formula. They are mostly equivalent for all studied interfaces and orientations except for Si3C/Ge3C interfaces where differences up to 10% are present.

For some interfaces (marked wiht '*'), the specular boundaries on the Y and Z axes are substituted by fully-diffusive boundaries due to technical limitations. Besides, increasing ∆T does not modify significantly the computed ITC. As permutting the hot and cold thermostats modifies the ITC by about 1%, no significant rectification of interface is observed.

Materials

SA ITC G W m -2 K -1 MC ITC G W m -2 K -1 Si3C/
In order to estimate the ITC in more complex devices, a specific method to directly calculate the temperature of the incident heat flux at the interfaces should be implemented in the future.

The thermal conductivity κ of HJ is computed as in CPNF. The evolution of κ with the device length L is plotted in Figure 43, with crosses. We find that it follows closely the results of SA formula for CPNF thermal conductance (continuous line, see §2.4) by using L 2 (which is the actual length of the materials in contact) instead of L.

10 -9 10 -8 10 -7 10 -6 10 -5 10 -4

Temperature dependence of the interface thermal conductance

We compare in Figure 44 our EMD and NEMD ITC results for several temperatures. The EMD ITC are independant on temperatures. They are on average equal to G = 650 M W m -2 k -1 , slightly above the results from [START_REF] Chalopin | Thermal Interface Conductance in Si/Ge Superlattices by Equilibrium Molecular Dynamics[END_REF].

However, the NEMD ITC increases when the average temperature increases. This phenomenon was has been previously observed in [START_REF] Landry | Thermal Boundary Resistance Predictions from Molecular Dynamics Simulations and Theoretical Calculations[END_REF].

As NEMD results are performed in short devices (in order to reduce simulation times), the computed conductances are low. But as seen in Figure 37 (d) the diffusive limit is higher than our EMD results. After having studied Si/Ge in the [100] lattice orientation, we compare it to the [110] and [111] orientations. In Figure 45, the temperature dependence of these three orientations are calculated by using EMD. In these simulations, the sizes of the simulation box are 8 × 8 × (32 + 32), 8 × 6 × (22 + 22) and 6 × 6 × (18 + 18) unit cells.

There is almost a 50% increase of ITC between the [100] and [110] In Figure 46, we compare the conductances of the two Ge polytype interfaces to the Si/Ge [100] one, computed by using EMD. We find that the mismatched polytype ITC is on the same order of magnitude as Si/Ge interfaces, about 50% higher. When stacked in the polytype nanowires, this low ITC value should lead to low a thermal conductivity of the device that could be interesting for thermoelectric applications. The aligned polytype Ge interface is however about 10 times more conductive. 300 400 500 600 700 800 900 10 9 The temperature dependence of the reference Si3C/Ge3C ITC calculated by using the SA DMM formula is shown in Figure 47 for three lattice orientations. These results have been presented in [START_REF] Larroque | Phonon Transmission at Si/Ge and Polytypic Ge Interfaces Using Full-Band Mismatch Based Models[END_REF]. As the DMM is a model based on the impedence of the materials, the results are quite different from those of the atomistic MD methods. While the [110] ITC is again higher than the [100] and [111] ones, they all have the same temperature dependence. The ITCs increase when the temperature increases but reach a plateau above 600K. Since the DMM model for the transmission can be seen as a strong approximation, our future work will implement the transmission coefficients computed by using MD in our Full Band MC simulator.

10 10 Temperature T [K] Interface thermal conductance G [W m -2 K -1 ] Si3C/Ge3C [100] Ge3C/Ge2H "mismatched" Ge3C/Ge2H "aligned"

Conclusion

Throughout this thesis, several thermal properties of nanostructures are studied. Two complementary simulation methods are used: Molecular Dynamics, which considers the position of atoms and interatomic potentials, and Monte Carlo simulation for phonons. Previous results on bulk thermal conductivity in Si and Ge are reproduced as well as the thermal conductance of Si/Ge interfaces. An original MC simulator with a full-band model of the material's dispersion relation and scattering rates computed by DFT have been developed during this thesis. As these parameters are based on ab initio calculations, this code can be used to study a large class of materials or crystalline phases. Models for rough boundaries and interface between materials are also implemented.

Both methods give insight on the thermal conductivity of the hexagonal phases of Si and Ge. The computed thermal properties of both materials are in agreement with previous theories and measurements. The thermal conductivity is lower for hexagonal phases as compared to the cubic ones.

It has been shown that even if spectral semi-analytical models can estimate satisfactorily the cross plane thermal conductivity in quasi-ballistic transport in nanofilms, they are disappointing when the rough interfaces are dominant as in the case of nanofilms with in-plane configuration or in nanowires. Besides, the contribution of each phonon branch to the heat flux is complex and highly anisotropic. This is also depicted by the polar and azimuthal descriptions of heat fluxes in nanofilms and nanowires. The introduction of rough boundaries impacts mostly acoustic modes (transverse acoustic for Si3C) and tends to focus the heat flux in the main transport direction.

Finally, polytype interfaces can exhibit low thermal conductivities, in the same order of magnitude as Si/Ge interfaces. The Molecular Dynamics methods provide a spectral decomposition of the interfacial heat flux, which is compared to results from the Diffusive Mismatch Model. In future works, more realistic models of interfaces computed by using MD will be implemented in our Monte Carlo simulator. This will allow the complete simulation of polytype nanowires and provide deep insight of the heat transport in these nanodevices.

A Synthèse en français

A.1 Introduction

La conversion thermoélectrique est un sujet de recherche très actif. Ses principales applications sont liées à la récupération d'énergie pour les systèmes d'alimentation électrique autonomes ainsi qu'au refroidissement des microcontrôleurs. Actuellement, les systèmes thermoélectriques courants reposent sur des composés rares et toxiques comme le Tellurure de Bismuth, le Tellurure de Plomb, etc. Leur remplacement par le Silicium et le Germanium, plus abondants et largement utilisés par l'industrie microélectronique, serait intéressant. Cependant, leurs propriétés thermoélectriques intrinsèques sont médiocres, en particulier près de la température ambiante.

Afin d'améliorer la conversion thermoélectrique des appareils à base de Si et de Ge, les nanotechnologies offrent de nouvelles voies d'optimisation. En effet, la nanostructuration permet d'améliorer les propriétés intrinsèques des nanostructures de plusieurs ordres de grandeur.

Dans ce cadre, cette thèse étudie les propriétés thermiques des nanofils polyphasés en Si et Ge récemment fabriqués. Ils sont constitués d'un empilement quasi périodique de phases cubiques et hexagonales de Si ou Ge. La présence de nombreuses interfaces polyphasées devrait être bénéfique pour les propriétés thermoélectriques. Les caractérisations expérimentales de ces nanofils étant encore très difficiles, le développement en parallèle d'outils numériques spécifiques est d'une grande importance pour le développement de cette technologie.

Au cours de cette thèse, un simulateur de Monte Carlo original dédié au transport des phonons a été développé. Il s'agit d'un simulateur "full band" et 3D dans l'espace réel capable de fournir des informations détaillées sur le transport thermique dans des nanodispositifs complexes. Comme les paramètres "semiempiriques" requis sont calculés à l'aide de données ab initio, il est possible de modéliser une grande classe de matériaux et de phases, même s'ils n'ont pas été caractérisés expérimentalement. Par ailleurs, des études réalisées à l'aide de la Dynamique Moléculaire sont présentées afin de paramétrer la modélisation d'interface semi-transparente dans le code MC.

Les différents modèles théoriques utilisés sont présentés dans le §A.2. Les méthodes de simulation Monte Carlo et de Dynamique Moléculaire sont dans le §A.3. Dans le §A.4, nos méthodes sont validées en étudiant les propriétés de matériaux massifs et de couches fines. Dans le §A.5, l'effet de surfaces rugueuses sur des nanostructures est présenté. Enfin, dans le §A.6, les propriétés d'interfaces entre différentes phases de Si et Ge sont étudiées.

A.2 Modèles théoriques

Dans ce chapitre, nous détaillons les bases théoriques utilisées dans cette thèse. Les deux premières sections présentent les dispositifs étudiés et les différents types de modélisation de matériaux. Les sections suivantes portent sur la modélisation du transfert thermique et présentent des approches semi-analytiques pertinentes.

A.2.1 Dispositifs simulés

Dans cette thèse, différents types de nanostructures sont étudiés : nanofilms, nanofils (NW) et hétérojonctions (HJ). Les nanofilms sont considérés dans des configurations où le flux thermique est perpendiculaire (CPNF) ou transverse (IPNF) aux interfaces. Des faces spéculaires sont introduites pour simuler les dimensions supposées infinies (respectivement 2, 1, 0 et 2 dimensions pour CPNF, IPNF, NW et HJ). Ces dispositifs sont paramétrés par leur longer L (distance entre les thermostats, alignés selon l'axe de transport thermique X) et leur largeur W .

A.2.2 Modèles de matériaux

Pour résoudre l'équation du transport thermique dans un système à l'état solide, il faut connaître au préalable la relation de dispersion et les taux de diffusion.

Le mouvement des phonons dans un matériau est perturbé par plusieurs mécanismes de diffusion, qui peuvent être caractérisés par leur fréquences d'interactions. Habituellement, la principale source de résistance thermique est due aux diffusions phonon-phonon.

Il est possible d'avoir une représentation analytique de la relation de dispersion et des taux de diffusion en utilisant plusieurs approximations. Cependant dans cette thèse, une description complète ("full band") de la zone de Brillouin (BZ) est utilisée, c'est-à-dire que les relations de dispersion et les taux de diffusion sont calculés numériquement pour un grand nombre de vecteurs d'ondes dans le volume complet de la BZ. Deux méthodes ont été utilisées pour calculer ces paramètres: l'ABCM [START_REF] Larroque | Étude Théorique de L'anisotropie Du Transport Thermique Dans Les Nanostructures à Base de Silicium et de Germanium[END_REF]) et la DFT [START_REF] Togo | Distributions of Phonon Lifetimes in Brillouin Zones[END_REF].

La DFT est une méthode ab initio, qui est donc pertinente pour étudier des matériaux pas encore (ou peu) caractérisés expérimentalement. Les paramètres que nous avons utilisés ont été calculés par [START_REF] Chaput | Ab Initio Based Calculations of the Thermal Conductivity at the Micron Scale[END_REF].

Dans les deux cas, la BZ est discrétisée en N vecteurs d'ondes, avec respectivement N = 31 × 31 × 31 = 29791 et N = 31 × 31 × 19 = 18259 pour les phases cubiques et hexagonales. La pulsation, la vitesse de groupe et (avec la DFT) le taux de diffusion phonon-phonon sont calculés pour chaque vecteur d'onde ⃗ q et mode m.

En raison de la nature discrète de notre description de la BZ, la définition des états d'iso-énergie doit être relâchée, entrainant des fluctuations sur la pulsation des phonons. La conservation de la pulsation est limitée par le pas de discrétisation sur la pulsation ∆ω.

A.2.3 Modèles de faces rugueuses

Dans les nanostructures, il faut tenir compte de la contribution des faces rugueuses sur la résistance thermique. En effet, cette contribution peut devenir plus forte que celle des mécanismes de diffusion phonon-phonon.

Pour cela, le modèle de Soffer considère qu'une particule peut subir lors d'une collision soit une réflexion spéculaire, soit une réflexion diffusive [START_REF] Soffer | Statistical Model for the Size Effect in Electrical Conduction[END_REF]. Cela est déterminé par une probabilité de réflexion spéculaire dépendante de l'angle d'incidence du phonon, d'un paramètre de rugosité de surface ∆, et d'un paramètre de corrélation spatiale de la rugosité.

Dans le cas où la corrélation spatiale de la rugosité de surface est négligeable, la probabilité de réflexion spéculaire à une face rugueuse est : Le modèle SA balistique évalue la conductivité thermique dans des nanofilms (CPNF) infiniment courts, utilisant le formalisme de Landauer.

p specular = e -(2cos(θ)∆|⃗ q|) 2 (
κ balistique = L V s ∆T ⎛ ⎜ ⎜ ⎝ ∑ état j |⃗ vj •⃗ n|>0 ℏω j |⃗ v j • ⃗ n|f BE (ω j , T chaud ) + ∑ état j |⃗ vj •⃗ n|<0 ℏω j |⃗ v j • ⃗ n|f BE (ω j , T froid ) ⎞ ⎟ ⎟ ⎠ (A2)
Le modèle SA diffusif est dérivé de l'équation de transport de Boltzmann, et correspond à la conductivité d'un matériau massif.

k diffusif = V s ∑ état j ℏω j |⃗ v j • ⃗ n| 2 1 λ j ∂f BE ∂T (ω j , T eq ) (A3)
où λ j est le taux de diffusion phonon-phonon pour l'état j.

Un modèle SA Matthiessen est d'abord utilisé pour estimer la conductivité de dispositifs dans un régime transitif entre les limites balistique et diffusive.

1

κ Matthiessen = 1 κ balistique + 1 κ diffusif (A4)
Des modèles plus complexes basés sur une décomposition spectrale de la conductivité sont généralement plus fiables. Nous avons défini des modèles spectraux pour trois dispositifs : SA CPNF, SA IPNF et SA NW, utilisant le modèle de Soffer pour les faces rugueuses (cf. §A.2.3).

Pour ces trois modèles, le taux de diffusion des phonons dans un état j est modifié tels que:

λ CPNF,j = λ j + |⃗ v j • ⃗ n| L (A5) λ IPNF,j = λ j + |⃗ v j • ⃗ n| L + |⃗ v j • ⃗ n ⊥,1 | W (A6) λ NW,j = λ j + |⃗ v j • ⃗ n| L + |⃗ v j • ⃗ n ⊥,1 | W + |⃗ v j • ⃗ n ⊥,2 | W (A7) où L = L 2 et W = W ( 1 1-p spéculaire -1 2
) .

Leurs conductivités thermiques κ CPNF , κ IPNF , et κ NW sont alors calculées avec la formule diffusive (Eq. A3) en adaptant les taux de diffusions.

A.2.5 Modèles de conductance thermique d'interface

La conductance thermique d'interface (ITC) est liée à la chute de température au bord de l'interface entre deux matériaux, quand un flux thermique est imposé.

Le modèle utilisé dans cette thèse pour définir une interface semi-transparente est une adaptation du Diffusive Mismatch Model (DMM) pour notre modèle de matériaux "full band" [START_REF] Larroque | Phonon Transmission at Si/Ge and Polytypic Ge Interfaces Using Full-Band Mismatch Based Models[END_REF]. La DMM considère une interface de mauvaise qualité, de telle sorte que tous les phonons en collision avec l'interface subissent une diffusion. Les phonons peuvent ensuite être soit transmis, soit réfléchit, selon une probabilité de transmission t. La BTE décrit l'évolution temporelle de la distribution de phonons f j (⃗ r, ⃗ q, t) dans l'espace de phase, où j est l'état des phonons, ⃗ r la position dans l'espace réel et ⃗ q leur vecteur d'onde. (diffusion, collisions, etc.).

G DMM (ω, T ) = ∂f BE ∂T (ω, T )I A (ω)t A→B (ω) (A8) où G DMM (ω,
∂f j ∂t + ⃗ v • ⃗ ∇f = ∂f ∂t ⏐ ⏐ ⏐ ⏐ interactions , (A9) où ⃗ v j est la vitesse de groupe (⃗ v j = ⃗ ∇ ⃗ q ω j , ω j est la pulsation) et ∂f ∂t ⏐ ⏐ ⏐ interactions est le terme d'interactions
Nous résolvons la BTE pour le transport de phonons avec une méthode particulaire Monte Carlo [START_REF] Péraud | Efficient Simulation of Multidimensional Phonon Transport Using Energy-Based Variance-Reduced Monte Carlo Formulations[END_REF].

Dans cette approche stochastique, les trajectoires d'un grand nombre de particules semi-classiques sont choisies au hasard. Les trajectoires de chaque particule sont décrites comme une séquence de vols libres (trajectoires linéaires dans l'espace réel sans changement dans l'espace réciproque) et d'événements de diffusion instantanés (pas de changement dans l'espace réel mais avec une modification de l'état dans l'espace réciproque). Le mécanisme de diffusion qui met fin à un vol libre peut être soit une diffusion phonon-phonon, soit une collision avec une face du dispositif. Pour chaque particule, l'état initial, la durée de chaque vol libre, ainsi que le type et l'effet de chaque événement de diffusion sont choisis au hasard en fonction des taux de diffusion appropriés.

A.3.1.2 Techniques d'optimisation

Pour réduire le nombre de particules simulées et donc la puissance de calcul requise, les particules simulées ne sont pas des phonons mais des paquets de phonons. Ces paquets ont un nombre de phonons variable, tels que leur énergie totale E p est toujours constante. Cela permet de simplifier grandement le traitement des diffusions élastiques et garantit dans tous les cas la conservation d'énergie [START_REF] Klitsner | Phonon Radiative Heat Transfer and Surface Scattering[END_REF]. Le paramètre E p définit la précision de la simulation.

Puisque seule la déviation de la distribution de phonons par rapport à la distribution d'équilibre est pertinente pour étudier le transport thermique, seulement celle-ci est simulée. Une température de référence T 0 est définie, proche des températures dans les dispositifs simulés, et la distribution d'équilibre à cette température est soustraite. Cela permet de réduire de manière significative le nombre de particules à simuler, et de réduire la variance des résultats [START_REF] Péraud | Efficient Simulation of Multidimensional Phonon Transport Using Energy-Based Variance-Reduced Monte Carlo Formulations[END_REF].

A.3.1.3 Température effective

La température de chaque section du dispositif simulé est mise à jour périodiquement, puisque les taux de diffusion phonon-phonon dépendent de la température locale. Pour cela, la relation énergie-température est pré-calculée pour chaque matériau. L'énergie totale des particules dans une section est sommée, puis la relation énergie-température est inversée numériquement pour obtenir la température locale.

A.3.1.4 Algorithme

En premier lieu, la distribution initiale de particule est générée dans l'ensemble du dispositif. La simulation entre ensuite dans deux boucles : l'une pour chaque pas de temps, et la deuxième pour chaque particule. À chaque pas de temps, la température locale est mise à jour dans chaque section du dispositif, et de nouvelles particules sont injectées par les thermostats. Ensuite, chaque particule est déplacée pendant son temps de vol libre, qui peut être est interrompu par une interaction phonon-phonon ou une collision avec les bords du dispositif.

En considérant que les interactions sont des processus poissoniens (cf. [START_REF] Jacoboni | The Monte Carlo Method for Semiconductor Device Simulation[END_REF]), au début d'un vol libre, le temps avant la prochaine interraction t interaction est déterminé aléatoirement :

t interaction = - ln(n random ) λ j (T c ) (A10)
où n random est un nombre aléatoire uniforme dans l'intervalle ]0;1] et λ j (T c ) est le taux d'interaction, dans notre cas le taux de diffusion phonon-phonon d'une particule dans l'état j à la température T c .

Bien que la diffusion phonon-phonon soit un processus à trois phonons (les ordres supérieurs étant ignorés), nous la traitons comme un processus à deux phonons [START_REF] Lacroix | Monte Carlo Transient Phonon Transport in Silicon and Germanium at Nanoscales[END_REF]. Ceci simplifie grandement l'algorithme, et comme chaque particule a une énergie fixe, la conservation de l'énergie est exacte. Cette approximation peut être considérée dans le cas où un grand nombre de phonons ont un grand nombre de collisions, où le comportement moyen est alors correct. Quand une particule subit une diffusion phonon-phonon à la fin de son vol libre, son état précédent est perdu [START_REF] Mazumder | Monte Carlo Study of Phonon Transport in Solid Thin Films Including Dispersion and Polarization[END_REF] et remplacé par un nouvel état. Cet état est sélectionné aléatoirement selon la distribution d'équilibre, pondérée par le taux de diffusion de chaque état (cf. [START_REF] Lacroix | Monte Carlo Transient Phonon Transport in Silicon and Germanium at Nanoscales[END_REF]).

p j ∝ λ j (T c )n eq j (T c ) (A11)
Les thermostats injectent un flux de particule constant, dépendant de leur température.

Les méthodes EMD utilisent un formalisme de Green-Kubo pour déterminer les propriétés thermiques, alors que les méthodes NEMD sont similaires à des mesures expérimentales de dispositifs.

Ces méthodes ont été décrites dans des travaux précédents :

• EMD massif [START_REF] Sellan | Size Effects in Molecular Dynamics Thermal Conductivity Predictions[END_REF][START_REF] Schelling | Comparison of Atomic-Level Simulation Methods for Computing Thermal Conductivity[END_REF]Volz and Chen 2000) • EMD aux interfaces [START_REF] Chalopin | Thermal Interface Conductance in Si/Ge Superlattices by Equilibrium Molecular Dynamics[END_REF] • NEMD massif [START_REF] Schelling | Comparison of Atomic-Level Simulation Methods for Computing Thermal Conductivity[END_REF] • NEMD aux interfaces [START_REF] Landry | Thermal Boundary Resistance Predictions from Molecular Dynamics Simulations and Theoretical Calculations[END_REF]

A.4 Matériaux massifs et nanofilms

A.4.1 Transition entre régimes balistique et diffusif

Les différents modèles SA ont été comparé dans des dispositifs de différentes longueurs L. Pour les dispositifs très longs, tous les modèles (sauf le modèle balistique) convergent vers la valeur diffusive de conductivité thermique qui correspond à la conductivité du matériau massif. De même, à l'exception du modèle diffusif, ils convergent asymptotiquement vers la conductance du modèle balistique dans les dispositifs courts. Cela confirme que les modèles sont équivalents à ces limites.

En comparant aux résultats de simulation MC, on trouve que la transition entre ces limites n'est pas correctement exprimée par le modèle Matthiessen (jusqu'à 60% de différences observées). Cependant, le modèle SA CPNF est très proche des résultats MC sur toute la gamme de longueurs. Cela indique qu'un traitement spectral de la conductivité est nécessaire pour représenter de manière précise la transition entre ces régimes. 10 -9 10 -8 10 -7 10 -6 10 -5 10 -4 10 0

A.4.2 Phases hexagonales de Si et Ge

Dans la

10 1 10 2 Longueur du dispositif L [m] Conductivité thermique κ [W m -1 K -1 ] Si3C Si2Hx Si2Hz

A.5 Nanostructures rugueuses

Ce chapitre considère l'effet des faces rugueuses, introduites dans les nanofilms IPNF, et les nanofils NW.

A.5.1 Effets dus à la géométrie des dispositifs

La dégradation de la conductivité thermique dans des dispositifs longs (donc diffusifs) est liée au nombre de faces rugueuses présentes. Cependant, dans les appareils ultra-courts dans lesquels le transport thermique est balistique, l'évolution des conductivités thermiques sont les mêmes. Dans le cas de IPNF et NW de très grandes largeurs (W > 1µm), ou de paramètres de rugosité très faibles (∆ < 0.1nm), la conductivité tend vers la valeur limite des CPNF. Il est notable que contrairement aux CPNF, les modèles semi-analytiques (SA), correspondent pas précisément aux simulations MC pour les IPNF et les NW, bien que les tendances soient les mêmes.

A.5.3 Comparaison avec des mesures expérimentales

Dans le cas des mesures de Li sur NW [START_REF] Li | Thermal Conductivity of Individual Silicon Nanowires[END_REF], nous comparons la dépendance en température de la conductivité thermique sur la figure A5. Seulement trois largeurs W sont considérées : 37nm, 56nm et 115nm, car les dispositifs de 22nm sont trop fins pour être simulés correctement en utilisant les paramètres des matériaux massifs. Notre simulation Monte Carlo reproduit correctement l'évolution de la conductivité avec la température. Les différences peuvent être dues à d'autres sources de résistance thermique dans les mesures expérimentales, comme des résistances de contact. [START_REF] Li | Thermal Conductivity of Individual Silicon Nanowires[END_REF].

A.6 Transport thermique aux interfaces solide-solide

Ce chapitre étudie le transport thermique aux interfaces entre des matériaux semi-conducteurs. L'un des objectifs de cette thèse est de caractériser les propriétés des interfaces polyphasées. Un autre consiste à mettre en oeuvre un modèle d'interface avancé dans un simulateur de Monte Carlo "full band" basé sur des données issues de simulations de dynamique moléculaire (MD), qui modélisent avec précision la position des atomes aux interfaces.

A.6.1 Distribution spectrale du flux thermique aux interfaces

Le flux thermique spectral est calculé en NEMD à partir des forces et des vitesses des atomes autour de l'interface [START_REF] Sääskilahti | Role of Anharmonic Phonon Scattering in the Spectrally Decomposed Thermal Conductance at Planar Interfaces[END_REF]. Dans le cas d'interface devrait réduire significativement la conductivité du dispositif, ce qui le rendrait intéressant pour des applications thermoélectriques. Par contre, l'interface polyphasée "alignée" a une ITC environ 10 fois plus grande.

A.7 Conclusion

Dans cette thèse, plusieurs propriétés thermiques de nanostructures sont étudiées. Deux méthodes de simulation complémentaires sont utilisées : la Dynamique Moléculaire, qui prend en compte la position des atomes et des potentiels interatomiques, et des simulation Monte Carlo de phonons. Les résultats précédents sur la conductivité thermique en Si et Ge massifs sont reproduits ainsi que la conductivité thermique d'interfaces Si/Ge. Un simulateur Monte Carlo original, utilisant un modèle "full band" de la relation de dispersion et des temps de relaxation calculés par DFT, a été développé au cours de cette thèse. Comme ces paramètres sont basés sur des calculs ab initio, ce code peut être utilisé pour étudier une grande classe de matériaux ou de phases cristallines. Des modèles de faces rugueuses et d'interfaces entre matériaux sont également étudiés.

Les deux méthodes de simulation donnent un aperçu de la conductivité thermique des phases hexagonales de Si et Ge. Les propriétés thermiques calculées des deux matériaux sont en accord avec les théories et mesures précédentes. La conductivité thermique est plus faible pour les phases hexagonales que pour les phases cubiques correspondantes.

Nous avons démontré que même si les modèles semi-analytiques spectraux peuvent estimer de manière satisfaisante la conductivité thermique dans le cas de transport quasi-balistique dans des nanofilms, elles sont décevantes lorsque l'effet des faces rugueuses est dominant comme dans le cas des nanofils. De plus, la contribution de chaque mode de phonon au flux thermique est complexe et hautement anisotrope. Les distributions angulaires du flux thermique dans les nanofilms et les nanofils le montrent également. L'introduction de faces rugueuses a un impact principalement sur les modes acoustiques et tend à concentrer le flux thermique dans la direction principale du transport.

Pour finir, les interfaces polyphasées en Ge peuvent présenter de faibles conductivités thermiques, du même ordre de grandeur que les interfaces Si/Ge. Des méthodes de Dynamique Moléculaire fournissent une décomposition spectrale du flux thermique aux interfaces, qui est comparée aux résultats du Diffusive Mismatch Model. Dans des travaux futurs, des modèles plus réalistes d'interfaces paramétrés à l'aide de simulations de Dynamique Moléculaire seront implémentés dans notre simulateur Monte Carlo. Ceci permettra la simulation complète des nanofils polyphasés et fournira une connaissance approfondie du transport thermique dans ces nanodispositifs.

R ésum é :

Les g én érateurs thermo électriques convertissent directement l' énergie thermique en énergie électrique. Ils pourraient devenir de plus en plus utiles à des fins de r écup ération d' énergie et font l'objet de recherches actives. Cependant, les meilleurs mat ériaux thermo électriques sont rares et Abstract : Thermoelectric generators are able to directly convert heat into electrical energy. They could have a great potential in terms of energy harvesting, but unfortunately, the best thermoelectric materials are rare and pollutant. Silicon and Germanium would be attractive materials if their thermoelectric efficiency were improved. For this purpose, nanostructuring is a possible route, for instance via the introduction of rough boundaries or interfaces between materials. Recently, polytype nanowires (composed of a sequence of cubic and hexagonal phases of Si and Ge) have been fabricated, but the experimental characterization of such complex nanostructures with exotic materials is challenging. In this thesis, we study the details of thermal transport in nanostructures with numerical simulations. An original Monte Carlo method is developed, with a full band ab initio description of materials. It includes models for the rough boundaries and the solid-solid interfaces.

Molecular Dynamics simulations are also performed to characterize the properties of interfaces. We confirm that the hexagonal phases of Si and Ge have lower thermal conductivity than their cubic counterparts. The full band model shows a strong anisotropy in the heat flux. Usual semi-analytical models failed to reproduce the thermal conductivity of simulated nanostructures with rough boundaries. Besides, those boundaries tend to focus the heat flux in the main direction of the nanostructure. Finally, some polytype interfaces can have an interfacial conductance almost as low as Si/Ge interfaces, and thus could improve significantly the thermoelectric efficiency of polytype nanowires. The presented Monte Carlo method could easily be used with a wide range of materials, and it can model arbitrarily complex nanostructures. In the future, the results from Molecular Dynamics simulation will be used to parametrize a more realistic model of solidsolid interfaces.

Universit é Paris-Saclay

Espace Technologique / Immeuble Discovery Route de l'Orme aux Merisiers RD 128 / 91190 Saint-Aubin, France

  3 EMD bulk conductivity at 1000K in cubic and hexagonal phases of Si and Ge, with Tersoff potential. . . . . . . . . . . . . . . . . 4 SA bulk conductivities at 1000K and 300K in cubic and hexagonal Si. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 Heat flux contribution per mode in 1µm long CPNF, (a) in Si3C, (b) Si2Hx, and (c) in Si2Hz. . . . . . . . . . . . . . . . . . . . . . 6 Heat flux contribution per phonon mode in NW with L = 1µm and W = 100nm, (a) in Si3C, (b) Si2Hx, and (c) in Si2Hz. . . . . 7 Comparison of ITC at 300K from DMM SA and ballistic MC for several interfaces. . . . . . . . . . . . . . . . . . . . . . . . . . . .List of Figures1 Evolution of the maximum ZT over time. Reproduced from (Heremans et al. 2013). . . . . . . . . . . . . . . . . . . . . . . . 2 Germanium polytype nanowire (Vincent et al. 2014). . . . . . . . 3 Schema of the simulated devices: (a) cross-plane nanofilm, (b) inplane nanofilm, (c) nanowire and (d) heterojunction. Cold (hashed blue) and hot (hashed red) thermostats. Rough boundaries (green). Solid-solid interface (hashed purple). . . . . . . . . . . . . . . . . 4 Cartography of the angular frequency ω in the BZ. Wave vectorsq i,j are in [10 13 2π/m]. Schema of the BZ in (a) Si3C, (b) Si2H. Iso-energies in (c) the (100) plane of Si3C, (d),(e),(f) in the (10-10), (11-20) and (0001) planes of Si2H, for the first phonon mode. . . 5 Cartography of phonon-phonon scattering rates λ in the BZ. Wave vectors q i,j are in [10 10 Hz]. Iso-energies in (a) the (100) plane of Si3C, (b),(c),(d) in the (10-10), (11-20) and (0001) planes of Si2H, for the first phonon mode. . . . . . . . . . . . . . . . . . . 6 Schema of a possible phonon trajectory. . . . . . . . . . . . . . . 7 Phonon energy density-temperature relation in Si3C. . . . . . . . 8 Main Monte Carlo algorithm. . . . . . . . . . . . . . . . . . . . . 9 Estimation of thermal relaxation time τ from SA formula and MC simulations. . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 Time evolution of the cumulated phonon frequency change during phonon-phonon scattering events of all simulated particles in a nanowire. Parameters: L = 1µm, W = 100nm, ∆ = 0.5nm. . . . 11 Time evolution of the cumulated phonon frequency change during rough boundary scattering events of all simulted particles in a nanowire. Parameters: L = 1µm, W = 100nm, ∆ = 0.5nm. . . . 12 Temperature profiles in cross-plane nanofilms, simulated with varying reference temperature T 0 , and other parameters equal. . 13 Evolution of the thermal conductivity as a function of the upper bound of the integral on correlation time t corr . Each black line shows the result of an independent simulation. The red line shows the average, and the error bars are confidence interval at 95%. . 14 Inverse of thermal conductivity as a function of the inverse of the distance between thermostats, for Si at 1000K. Results from NEMD simulations (blue line). The red and green lines are the extrapolations at 1 L = 0, for the average and the maximum and minimum errors, respectively. . . . . . . . . . . . . . . . . . . . . 15 Diagram of the simulation box for simulating an interface in NEMD. 16 Temperature profile of the simulation box along the transport direction of a Si/Ge interface with NEMD. . . . . . . . . . . . . 17 Thermal conductivities in CPNF as function of film length, in Si3C. 18 Thermal conductivity in CPNF as function of device length L, in Si3C (blue), Si2Hx (green) and Si2Hz (red). . . . . . . . . . . . . 19 Spectral distribution of the phonon energy for each mode in a 1µm long CPNF. (a) Si3C (b) Si2H. Crosses for MC. Lines for SA CPNF. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 Cumulative spectral distribution of the heat flux in a 1µm long CPNF. (a) Si3C, (b) Si2Hx and (c) Si2Hz. Crosses for MC simulations, continuous lines for SA CPNF model. . . . . . . . . 21 Spherical coordinates, showing the polar angle θ and the azymuthal angle ϕ. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 Angular density of state and angular heat flux in Si3C, Si2Hx and Si2Hz, in a 1µm long CPNF. (a), (c), and (e) are the aDOS, (b), (d), and (f) are the angular heat fluxes, in Si3C, Si2Hx and Si2Hz, respectively . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 Angular heat flux as function of the polar θ and azymuthal angles ϕ in Si3C, Si2Hx, and Si2Hz. . . . . . . . . . . . . . . . . . . . . 24 (a) Angular density of state and (b) angular heat flux in an artificial isotropic material, in a 1µm long CPNF. . . . . . . . . . 25 Thermal conductivity as function of device length in Si3C for CPNF , IPNF and NW . . . . . . . . . . . . . . . . . . . . . 26 Thermal conductivity as function of device width in Si3C for CPNF , IPNF and NW . . . . . . . . . . . . . . . . . . . . . 27 Thermal conductivity as function of the surface roughness parameter in Si3C for CPNF , IPNF and NW . . . . . . . . . . . . 28 Cumulative spectral distributions of the heat flux in a 1µm long CPNF and NW for each phonon mode. . . . . . . . . . . . . . . . 29 Heat flux angular distribution in Si2Hx NW vs. CPNF, L = 1µm and W = 100nm. . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 Thermal conductivity as function of device length in Si3C for NW and NW FD . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 Thermal conductivity κ as function of device length L in Si3C for NW FD in [100], [110] and [111] lattice orientations. . . . . . . . 32 Thermal conductivity in Si3C in (a) IPNF and (b) NW, as function of device width with L = 1µm. Lines are MC simulations, crosses experimental results from

  . . . . . . . . . . . . 34 Diagram showing the position of atoms forming the (a) "aligned" and (b) "mismatched" Ge polytype interfaces. Blue and red dots are from the cubic and hexagonal phases respectively. . . . . . . 35 By using EMD, ITC at the Si3C/Ge3C interface vs. (a) number of independent runs, (b) simulated time, (c) width W , and (d) length L. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 Temperature profiles in the Si3C/Ge3C superlattices along the direction normal to the interfaces, at three imposed temperature differences ∆T : (a) 10K, (b) 50K and (c) 100K. White and grey backgrounds represent Si and Ge, respectively. Hashed areas represent the thermostats. Red bold lines are linear interpolations used to extrapolate the temperature at the interfaces. . . . . . . 37 By using NEMD, Si3C/Ge3C ITC vs. (a) temperature bias ∆T , (b) simulated time, (c) width W , and (d) length L. . . . . . . . . 38 Spectral distributions of the interfacial heat flux at Si3C/Ge3C and Ge3C/Ge2H mismatched interfaces from NEMD simulations. 39 Spectral distributions of the interfacial heat flux at Si3C/Ge3C interfaces from NEMD simulations and DMM. . . . . . . . . . . 40 Spectral transmission coefficients at Si3C/Ge3C and polytype Si3C/Si2Hx interfaces from DMM. . . . . . . . . . . . . . . . . . 41 Temperature profile in a Si3C/Si3C HJ, with a DMM interface at x = 50nm. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 Evolution of the ITC calculated by using irrelevant temperature drops with the device length L. . . . . . . . . . . . . . . . . . . . 43 Thermal conductivity in Si3C/Si3C HJ as function of device length. 44 Thermal conductances of Si/Ge strained interfaces as a function of temperature. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 Interface thermal conductances of Si/Ge strained interfaces for several lattice orientations at different temperatures using EMD. 46 Interface thermal conductance of polytype Ge interfaces and Si/Ge [100] interface as a function of temperature using EMD. . . . . . 47 Si/Ge ITC from SA ABCM. . . . . . . . . . . . . . . . . . . . . .

Figure 1 (

 1 Figure 1 (reproduced from Heremans et al. (2013)) shows the evolution of the figure of merit ZT of several materials between 1940 and 2010. Constant progress was made until 1970 when ZT of bulk materials reached a plateau. The

Figure 3 :

 3 Figure 3: Schema of the simulated devices: (a) cross-plane nanofilm, (b) inplane nanofilm, (c) nanowire and (d) heterojunction. Cold (hashed blue) and hot (hashed red) thermostats. Rough boundaries (green). Solid-solid interface (hashed purple).

Figure 4 :

 4 Figure 4: Cartography of the angular frequency ω in the BZ. Wave vectors q i,j are in [10 13 2π/m]. Schema of the BZ in (a) Si3C, (b) Si2H. Iso-energies in (c) the (100) plane of Si3C, (d),(e),(f) in the (10-10), (11-20) and (0001) planes of Si2H, for the first phonon mode.

Figure 5 :

 5 Figure 5: Cartography of phonon-phonon scattering rates λ in the BZ. Wave vectors q i,j are in [10 10 Hz]. Iso-energies in (a) the (100) plane of Si3C, (b),(c),(d) in the (10-10), (11-20) and (0001) planes of Si2H, for the first phonon mode.

]Figure 7 :

 7 Figure 7: Phonon energy density-temperature relation in Si3C.

Figure 8 :

 8 Figure 8: Main Monte Carlo algorithm.

Figure 9 :

 9 Figure 9: Estimation of thermal relaxation time τ from SA formula and MC simulations.

Figure 10 :

 10 Figure 10: Time evolution of the cumulated phonon frequency change during phonon-phonon scattering events of all simulated particles in a nanowire. Parameters: L = 1µm, W = 100nm, ∆ = 0.5nm.

Figure 11 :

 11 Figure 11: Time evolution of the cumulated phonon frequency change during rough boundary scattering events of all simulted particles in a nanowire. Parameters: L = 1µm, W = 100nm, ∆ = 0.5nm.

Figure 12 :

 12 Figure 12: Temperature profiles in cross-plane nanofilms, simulated with varying reference temperature T 0 , and other parameters equal.

Figure 13 :

 13 Figure 13: Evolution of the thermal conductivity as a function of the upper bound of the integral on correlation time t corr . Each black line shows the result of an independent simulation. The red line shows the average, and the error bars are confidence interval at 95%.

Figure 15 :

 15 Figure 15: Diagram of the simulation box for simulating an interface in NEMD.

Figure 16 :

 16 Figure 16: Temperature profile of the simulation box along the transport direction of a Si/Ge interface with NEMD.

Figure 17 :

 17 Figure 17: Thermal conductivities in CPNF as function of film length, in Si3C.

Figure 19 :Figure 20 :

 1920 Figure 19: Spectral distribution of the phonon energy for each mode in a 1µm long CPNF. (a) Si3C (b) Si2H. Crosses for MC. Lines for SA CPNF.

Figure 21 :

 21 Figure 21: Spherical coordinates, showing the polar angle θ and the azymuthal angle ϕ.

Figure 22 :

 22 Figure 22: Angular density of state and angular heat flux in Si3C, Si2Hx and Si2Hz, in a 1µm long CPNF. (a), (c), and (e) are the aDOS, (b), (d), and (f) are the angular heat fluxes, in Si3C, Si2Hx and Si2Hz, respectively

Figure 23 :

 23 Figure 23: Angular heat flux as function of the polar θ and azymuthal angles ϕ in Si3C, Si2Hx, and Si2Hz.

Figure 24 :

 24 Figure 24: (a) Angular density of state and (b) angular heat flux in an artificial isotropic material, in a 1µm long CPNF.

Figure 25 :

 25 Figure 25: Thermal conductivity as function of device length in Si3C for CPNF , IPNF and NW .

Figure 26 :

 26 Figure 26: Thermal conductivity as function of device width in Si3C for CPNF , IPNF and NW .

Figure 27 :

 27 Figure 27: Thermal conductivity as function of the surface roughness parameter in Si3C for CPNF , IPNF and NW .

Figure 28 :

 28 Figure 28: Cumulative spectral distributions of the heat flux in a 1µm long CPNF and NW for each phonon mode.

WFigure 29 :

 29 Figure 29: Heat flux angular distribution in Si2Hx NW vs. CPNF, L = 1µm and W = 100nm.

  10 -9 10 -8 10 -7 10 -6 10 -5 10 -4

Figure 30 :

 30 Figure 30: Thermal conductivity as function of device length in Si3C for NW and NW FD .

Figure 31 :

 31 Figure 31: Thermal conductivity κ as function of device length L in Si3C for NW FD in [100], [110] and [111] lattice orientations.

Figure 32 :

 32 Figure 32: Thermal conductivity in Si3C in (a) IPNF and (b) NW, as function of device width with L = 1µm. Lines are MC simulations, crosses experimental results from[START_REF] Liu | Thermal Conductivity Measurements of Ultra-Thin Single Crystal Silicon Layers[END_REF],[START_REF] Ju | Phonon Scattering in Silicon Films with Thickness of Order 100 Nm[END_REF],[START_REF] Li | Thermal Conductivity of Individual Silicon Nanowires[END_REF] and[START_REF] Hochbaum | Enhanced Thermoelectric Performance of Rough Silicon Nanowires[END_REF]).

Figure 33 :

 33 Figure33: Thermal conductivity in NW as function of temperature, compared with exp. results from[START_REF] Li | Thermal Conductivity of Individual Silicon Nanowires[END_REF].

  (a)) • Ge3C [00-1]/Ge2H [0-110], refered as "mismatched" polytype interfaces (see Figure 34 (b))

Figure 35 :

 35 Figure 35: By using EMD, ITC at the Si3C/Ge3C interface vs. (a) number of independent runs, (b) simulated time, (c) width W , and (d) length L.

Figure 36 :

 36 Figure 36: Temperature profiles in the Si3C/Ge3C superlattices along the direction normal to the interfaces, at three imposed temperature differences ∆T : (a) 10K, (b) 50K and (c) 100K. White and grey backgrounds represent Si and Ge, respectively. Hashed areas represent the thermostats. Red bold lines are linear interpolations used to extrapolate the temperature at the interfaces.

Figure 37 :

 37 Figure 37: By using NEMD, Si3C/Ge3C ITC vs. (a) temperature bias ∆T , (b) simulated time, (c) width W , and (d) length L.

Figure 38 :

 38 Figure 38: Spectral distributions of the interfacial heat flux at Si3C/Ge3C and Ge3C/Ge2H mismatched interfaces from NEMD simulations.

Figure 39 :

 39 Figure 39: Spectral distributions of the interfacial heat flux at Si3C/Ge3C interfaces from NEMD simulations and DMM.

Figure 40 :

 40 Figure 40: Spectral transmission coefficients at Si3C/Ge3C and polytype Si3C/Si2Hx interfaces from DMM.

Figure 41 :

 41 Figure 41: Temperature profile in a Si3C/Si3C HJ, with a DMM interface at x = 50nm.

Figure 42 :

 42 Figure 42: Evolution of the ITC calculated by using irrelevant temperature drops with the device length L.

Figure 43 :

 43 Figure 43: Thermal conductivity in Si3C/Si3C HJ as function of device length.

Figure 44 :

 44 Figure 44: Thermal conductances of Si/Ge strained interfaces as a function of temperature.

Figure 45 :

 45 Figure 45: Interface thermal conductances of Si/Ge strained interfaces for several lattice orientations at different temperatures using EMD.

Figure 46 :

 46 Figure 46: Interface thermal conductance of polytype Ge interfaces and Si/Ge [100] interface as a function of temperature using EMD.

Figure 47 :

 47 Figure 47: Si/Ge ITC from SA ABCM.

Figure A1 :

 A1 Figure A1: Conductivité thermique de CPNF en fonction de leur longueur L, en Si3C (bleu), Si2Hx (vert) et Si2Hz (rouge).

Figure A2 :

 A2 Figure A2: Intégrale de la distribution spectrale du flux thermique pour chaque mode dans un CPNF de 1µm de long, en Si3C, avec le modèle SA CPNF (ligne continue) et des simulations MC (croix).

Figure A3 :

 A3 Figure A3: Distributions angulaires (a) de la densité d'état, (b) du flux thermique, dans un NW de 1µm de long en Si3C.

A. 5 . 2

 52 Impact sur les distributions spectrales et angulaires du flux thermique La Figure A4 (a) compare la distribution spectrale du flux dans un CPNF (ligne continue) et un NW (pointillés) de même dimension (L = 1µm et W = 100nm), en Si3C, avec un paramètre de rugosité ∆ = 0.5nm. On observe que la réduction du flux par les faces rugueuses degrade plus fortement le transport des modes acoustiques. La proportion de flux transporté par les modes optiques est donc plus grande dans les NW que dans les CPNF, bien que toujours très faible. La Figure A4 (b) compare la distribution angulaire du flux selon l'angle azimutal ϕ dans un CPNF (ligne continue) et un NW (pointillés), en Si2Hx, avec les mêmes paramètres que précédemment. Comme vu dans le chapitre précédent, le flux en CPNF est transporté par les orientations cristallines ayant une grande densité d'état. Cependant, l'introduction de faces rugueuses réduit principalement le flux dans les orientations hors de la direction du transport (par exemple [21-30]). Il y a donc une concentration du flux thermique.

  Figure A4: (a) Intégrale de la distribution spectrale du flux thermique pour chaque mode, dans un CPNF (ligne continue) et un NW (pointillés) de 1µm de long, en Si3C, depuis des simulations MC. (b) Distribution angulaire du flux selon l'angle azymuthal ϕ dans un CPNF (ligne continue) et un NW (pointillés), en Si2Hx, depuis des simulations MC

Figure A5 :

 A5 FigureA5: Conductivité thermique de NW en Si3C en fonction de la température, comparés aux mesures expérimentales[START_REF] Li | Thermal Conductivity of Individual Silicon Nanowires[END_REF].

Figure A7 :

 A7 Figure A7: Conductance thermique d'interface Si3C/Ge3C et d'interfaces Ge polyphasées "alignée" et "non alignée", en fonction de la température, en EMD.

  polluants. Le Silicium et le Germanium seraient des mat ériaux int éressants si leur efficacit é thermo électrique était am élior ée. Pour ce faire, la nanostructuration est une voie possible, par exemple en introduisant des faces rugueuses ou de nouvelles interfaces semitransparentes. R écemment, des nanofils polyphas és (compos és d'une alternance de phases cubiques et hexagonales de Si et Ge) ont ét é fabriqu és, mais la caract érisation exp érimentale de nanostructures aussi complexes comprenant des mat ériaux exotiques peut être difficile. Dans cette th èse, nous étudions en d étail le transport thermique dans des nanostructures avec des simulations num ériques. Une m éthode Monte Carlo originale a ét é d évelopp ée, avec une description "full band" des mat ériaux. Elle inclut des mod èles pour les faces rugueuses et les interfaces entre mat ériaux. Des simulations de Dynamique Mol éculaire sont également effectu ées pour caract ériser les propri ét és des interfaces. Nous confirmons que les phases hexagonales de Si et Ge ont une conductivit é thermique inf érieure à celle des phases cubiques correspondantes. Le mod èle "full band" montre que le flux thermique est fortement anisotrope. Des mod èles semi-analytiques habituels n'ont pas pu reproduire la conductivit é thermique des nanostructures simul ées avec des faces rugueuses. De plus, ces faces ont tendance à concentrer le flux de chaleur dans la direction principale de la nanostructure. Enfin, certaines interfaces polyphas ées peuvent avoir une conductance thermique presque aussi faible que les interfaces Si/Ge, et pourrait ainsi am éliorer significativement l'efficacit é thermo électrique des nanofils polyphas és. La m éthode Monte Carlo pr ésent ée peut facilement être utilis ée pour étudier une large gamme de mat ériaux, et elle est capable de mod éliser des nanostructures arbitrairement complexes. A l'avenir, les simulations en Dynamique Mol éculaire seront utilis ées pour param étrer un mod èle plus r éaliste d'interfaces. Title : Theoretical study of thermoelectric properties in nanostructures Keywords : thermoelectricity, simulation, Monte Carlo, Molecular Dynamics, thermal transport

Table 1 :

 1 Si bulk thermal conductivity at 1000K from EMD and NEMD.

	Capinski Exp. Ruf NEMD Ter. NEMD SW NEMD SW Sellan NEMD SW Schelling EMD Ter. EMD SW EMD SW Sellan EMD SW Schelling	44 45.8 49.7 ± 1.2 52.7 ± 1.7 40 ± 8 65 ± 16 77.2 ± 21.3 60 ± 17 60 ± 12 62 ± 16
	Exp. Oz NEMD Ter. NEMD SW EMD Ter. EMD SW	κ W m -1 K -1 15 19.7 ± 0.5 28.5 ± 0.4 32.9 ± 8.0 82.8 ± 27.4

Table 2 :

 2 Ge bulk thermal conductivity at 1000K from EMD and NEMD.Our MD results are consistent with previous MD simulations. According to our error estimations, they overestimate the experimental values. Besides, Ge results within the EMD method using Stillinger Weber potential appear disappointing and only the Tersoff potential will be used in the following for Ge study.In Table3, the thermal conductivity at T = 1000K of cubic and hexagonal phases of Si and Ge are compared by using the EMD method with Tersoff potential.

	Si 77.2 ± 21.3 32.4 ± 7.2 Ge Hexagonal 72.0 ± 16.8 31.6 ± 11.2 Cubic

Table 3 :

 3 EMD bulk conductivity at 1000K in cubic and hexagonal phases of Si and Ge, with Tersoff potential.

Table 4

 4

	. The hexagonal phase is studied in two lattice orientations: [10-10] (Si2Hx) and [0001] (Si2Hz).

Table 4 :

 4 SA bulk conductivities at 1000K and 300K in cubic and hexagonal Si.

Table 5 :

 5

	(b)	1 (TA) 2 (TA) 3 (LA) 4 (TO) 5 (TO) 6 (LO) Mode Heat flux contribution 1 20.0% 2 23.6% 3 17.5% 4 17.8% 5 14.4% 6 1.2% 7 1.8% 8 2.8% 9 0.22% 10 0.22% 11 0.27% 12 0.08%	flux contribution 29.5% 41.5% 25.2% 4.2% -0.3% -0.1% (c) Mode Heat flux contribution 1 25.3% 2 28.4% 3 12.6% 4 19.1% 5 16.6% 6 -3.8% 7 2.6% 8 -0.32% 9 -0.02% 10 -0.01% 11 -0.33% 12 -0.28%

Heat flux contribution per mode in 1µm long CPNF, (a) in Si3C, (b) Si2Hx, and (c) in Si2Hz.

Table 6 :

 6 Heat flux contribution per phonon mode in NW with L = 1µm and

	(b)	Heat flux contribution 37.0% 39.9% 18.9% 4.8% -0.4% -0.2% Mode Heat flux contribution 1 (TA) 2 (TA) 3 (LA) 4 (TO) 5 (TO) 6 (LO) 1 18.4% 2 20.7% 3 20.0% 4 15.7% 5 14.5% 6 1.9% 7 2.7% 8 4.6% 9 0.32% 10 0.48% 11 0.40% 12 0.18% (c) Mode Heat flux contribution 1 23.5% 2 24.3% 3 14.3% 4 26.9% 5 17.6% 6 -7.8% 7 3.8% 8 -1.2% 9 0.01% 10 -0.09% 11 -0.53% 12 -0.68%

Table 7 :

 7 Comparison of ITC at 300K from DMM SA and ballistic MC for several interfaces.

	Si3C Si3C/Si3C [110] Si3C/Si3C [111] Si3C/Si2Hx Si3C/Si2Hz Si3C/Ge3C Si3C/Ge3C [110] Si3C/Ge3C [111]	536 527 526 512 485 245 239 249	531 528 524 * 515 480 266 * 241 * 272 *

4 Modèles de transport thermique et conductivité thermique

	A1)
	où ⃗ q est le vecteur d'onde incident, θ l'angle d'incidence, et ∆ le paramètre de rugosité.
	Ce modèle a permis à Kazan et al. (2010) de reproduire les mesures expérimen-tales de Hochbaum et al. (2008)
	A.2.

La conductivité thermique d'un matériau peut être calculée par une formule semi-analytique (SA) à partir de la relation de dispersion et des taux de diffusion "full band" (cf. §A.2.2). Des modèles similaires ont été développés par

Mingo et al. (2003)

et

[START_REF] Aksamija | Anisotropy and Boundary Scattering in the Lattice Thermal Conductivity of Silicon Nanomembranes[END_REF]

, l'originalité de cette thèse est d'utiliser des taux de dispersion "full band" calculés en DFT.

3 Simulation Monte Carlo et Dynamique Moléculaire

	A.3.1 Méthode Monte Carlo
	A.3.1.1 Principes
	La méthode de Monte Carlo est une approche stochastique pour résoudre l'équation de transport de Boltzmann (BTE). Dans cette thèse, nous présentons un nouveau simulateur pour le transport de phonons en MC, qui résout la BTE en 3D dans l'espace réel et l'espace réciproque (approche "full band"). La relation de dispersion des phonons et les taux de diffusion phonon-phonon, sont paramétrés dans toute la zone Brillouin en utilisant les calculs de DFT ab initio (cf. §A.2.2).
	Ce chapitre présente les méthodes de simulations numérique utilisées dans cette thèse pour étudier le transport thermique dans des nanostructures.
	D'abord, notre méthode Monte Carlo (MC) originale développée partiellement lors de cette thèse est détaillée, puis les méthodes de Dynamique Moléculaire (MD) sont rappelées.

T ) est la décomposition spectral de l'ITC à la température T , f BE est la distribution d'équilibre de Bose Einstein et I A est le flux de phonons à l'interface.

A.

Table A1 ,

 A1 la conductivité thermique des phases cubiques et hexagonales de Si et Ge à T = 1000K sont calculées par EMD. Dans les deux cas, les phases hexagonales ont des conductivités inférieures aux phases cubiques correspondantes.

	Si 77.2 ± 21.3 32.4 ± 7.2 Ge Hexagonale 72.0 ± 16.8 31.6 ± 11.2 Cubique

Table A1 :

 A1 Conductivité thermique à 1000K par EMD, pour des phases cubiques et hexagonales de Si et Ge.Ces matériaux sont ensuite étudiés avec des simulations MC et le modèle SA CPNF. La FigureA1montre l'évolution de la conductivité avec la longueur du dispositif dans des CPNF en Si cubique (Si3C) et Si hexagonal dans les orientations cristallines [10-10] (Si2Hx) et [0001] (Si2Hz). Les courbes ont les mêmes caractéristiques, et confirment que la conductivité des phases hexagonales est inférieure à celle de la phase cubique pour toutes les longueurs. À la limite diffusive, κ diffusif = 138, 100 et 74W m -1 K -1 respectivement pour Si3C, Si2Hx et Si2Hz. De plus, on observe une anisotropie entre les deux phases hexagonales, puisque la conductivité diffusive est supérieure de 26% en Si2Hx comparé à Si2Hz.

Acknowledgements

This work is supported by a public grant overseen by the French National Research Agency (ANR) as part of the "Investissements d'Avenir" program (Labex NanoSaclay, reference: ANR-10-LABX-0035).

Selection of MD simulation parameters

The accuracy of ITC computation by using MD methods relies on a good selection of the relevant simulation parameters. The following MD simulations uses the Stillinger-Weber potential. Indeed, the Tersoff potential used for studying bulk materials (see §4), could not easily be implement by using the spectral analysis that we have choosen (cf. §6.3, the method distinguishes the contribution of atoms at the interface that is not directly consistent with the Tersoff approach). All simulations only consider the steady state.

EMD parameters

The EMD method is reminded in §3.2.3. The main simulation parameters for EMD and their default values are:

• 10 independent runs, • a simulated time of 0.5 ns, • a data collection time (1,000,000 timesteps of 0.5 fs) of 0.5 ns, • a simulated box width W of 8 unit cells, • a simulated box length L of 32 unit cells of material on each side of the interface, • an equilibrium temperature T of 400 K.

In Figure 35, it is presented the ITC values for several deviations of the previous simulation parameters from their default value.

In Figure 35 (a), the ITC is calculated by averaging over an increasing number of independent runs. The values reach a plateau after 5 runs.

In Figure 35 (b), an increasing total simulation time is used for averaging each independent run. At least 400,000 timesteps are required in this case.

In Figure 35 (c), the influence of the width of the simulation box is investigated. A sudden drop of the ITC is observed in boxes smaller than 6 unit cells. In such devices smaller than the wavelength of some important phonons, the heat transfer cannot be simulated correctly.

Aux faces lisses, les particules subissent une réflexion spéculaire : la composante normale à la face du vecteur d'onde est inversée.

Lors d'une collision avec une face rugueuse, la particule a une probabilité de subir soit une réflexion spéculaire, soit une réflexion diffusive [START_REF] Soffer | Statistical Model for the Size Effect in Electrical Conduction[END_REF]. Dans le cas d'une réflexion diffusive, un nouvel état est sélectionné aléatoirement dans l'isoénergie concernée, ayant une vitesse positive selon la normale à la face. La probabilité de sélection de ces états est pondérée par la probabilité de réflexion spéculaire et la composante normale de la vitesse (conformément à la Loi de Lambert) :

Lors d'une collision avec une face semi-transparente, la particule a une probabilité de subir soit une transmission diffusive, soit une réflexion diffusive, selon une probabilité de transmission t A→B calculée par DMM. L'implémentation est similaire aux faces rugueuses, le nouvel état après diffusion étant sélectionné dans l'isoénergie concernée, et ayant une orientation de vitesse cohérente. De la même manière, la probabilité de sélection de ces états est pondérée par la probabilité de réflexion et la composante normale de la vitesse :

A.3.2 Dynamique Moléculaire

Les méthodes de Dynamique Moléculaire (MD) sont des simulations classiques des trajectoires d'atomes. Comme chaque position atomique est déterminée lors de l'initialisation, il est possible de modéliser avec précision les interfaces entre matériaux par ces techniques.

Dans une simulation MD, chaque atome a une masse, une position et une vitesse définies. La deuxième équation de Newton dépendante du temps est résolue numériquement pour obtenir les trajectoires atomiques. Les forces instantanées des interactions entre atomes sont calculées à partir de potentiels semi-empiriques. Enfin, toutes les propriétés thermodynamiques pertinentes peuvent être obtenues à partir d'une analyse statistique du système.

Ces techniques sont assez polyvalentes mais ont leurs limites. Tout d'abord, ces méthodes exigent beaucoup de puissance de calcul, donc seuls de petits systèmes de quelques milliers d'atomes sont gérables. Afin d'éviter de diffuser les phonons aux limites, des conditions périodiques sont généralement appliquées, mais la taille du système par rapport au libre parcours des phonons doit être prise en compte. Enfin, nous devons garder à l'esprit qu'il s'agit d'une simulation classique qui ne peut rendre compte d'effets quantiques. Les simulations MD sont donc limitées à des régimes de températures élevées supérieures à la température de Debye du matériau, qui peut être relativement élevée dans le cas des semiconducteurs standards.

Des méthodes à l'équilibre (EMD) et hors équilibre (NEMD) ont été utilisées dans les travaux suivants, pour étudier des matériaux massifs et des interfaces. • interface "alignée" Ge3C [100] // Ge2H [10-10],

• interface "non alignée" Ge3C [00-1] // [0-110], correspondant aux interfaces observées dans les nanofils polyphasée de Vincent et al. (2014).

Nous constatons que l'ITC de l'interface polyphasée "non alignée" est du même ordre de grandeur que l'ITC d'interfaces Si/Ge (supérieure d'environ 50%). Dans une nanostructure possédant un grand nombre de cette interface, la faible ITC