Study of phonon transport across several Si/Ge interfaces using full-band phonon Monte Carlo simulation
11 Study of phonon transport across several Si/Ge interfaces using full-band phonon Monte Carlo simulation
N. D. Le , B. Davier , P. Dollfus , J. Saint-Martin Université Paris-Saclay, CNRS, Centre de Nanosciences et de Nanotechnologies, 91120, Palaiseau, France Department of Mechanical Engineering, The University of Tokyo, Tokyo 113-856, Japan
Abstract
A Full Band Monte Carlo simulator has been developed to consider phonon transmission across interfaces that are perpendicular to the heat flux. This solver of the Boltzmann transport equation which does not require any assumption on the shape the phonon distribution can naturally consider all phonon transport regimes from the diffusive to the fully ballistic regime. Hence, this simulator is used to study single and double Si/Ge heterostructures from the micrometer scale down to the nanometer scale i.e. in all phonon transport regime from ballistic to fully diffusive. A methodology to estimate the thermal conductivities and the thermal interfaces is presented. I Introduction
The management of the heat transfer at the nanoscale is crucial for the design and optimization of devices especially in thermoelectric [1] and nanoelectronics [2]. However, common macroscopic models such as heat Fourier’s formalisms [3] which assumes local equilibrium is not yet relevant in nanodevices that are shorter than the mean free path of phonons [4]. Thus, in complex nanostructures, the phonon distributions could significantly differ from the Bose Einstein statistics, especially if thermal interfaces are involved. Several analytical models are able to capture the out of equilibrium transport regime of heat in homogenous systems [5]–[8]. Besides, the modeling of interface thermal resistance also called the Kapitza’s resistance [9] is mainly based on the Acoustic Mismatch Model (AMM) [10] and the Diffusive Mismatch Model (DMM) [11] and their extensions [12], [13][14], [15]. Only recently an analytical modeling of heat transfer in heterostructure based on a set of a few parameters relevant at the nanometer scale was presented [16]. Nevertheless, the most versatile approaches to investigate the thermal transfers at the nanoscale are the use of advance numerical simulations. Among them, Molecular Dynamics (MD) simulations that calculate the classical trajectory of atoms have been widely applied [17]–[19]. These atomistic and time dependent simulations naturally consider the anharmonicity of phonons via semi-empirical parameters that can now be derived via ab-initio methods [20], [21]. However, the application of MD remains limited to systems containing thousands of atoms, i.e. with a typical size in the order tens of nanometers. The Non Equilibrium Green’s Functions framework [22]–[24] is another approach that is able to consider mechanical waves at the nanoscale and thus their possible reflection [23], [25]. This method that is not time dependent can nevertheless investigate systems larger than those possible in MD but the implementation of the anharmonicity requires a strong cost in term of computational burden [26]. The stochastic particle Monte Carlo (MC) method [27]–[29] is an efficient and accurate method to solve the time dependent Boltzmann’s Transport Equation (BTE) much beyond the linear approximation. This semi-classical framework is relevant from the nano to the microscale even in complex geometries in which thermal interface are involved [30]–[32]. The phonon dispersion could use a simple parabolic approximation or considered a full band description [33]. Moreover, accurate descriptions of scattering mechanisms and in particular phonon/phonon scatterings can be implemented at the particle level considering scattering rates that could be semi-empirical or derived from ab-initio method [34] [35]. In the present work, phonons transport in heterostructures have been investigated in all phonon transport regimes by using our Full-band MC simulation and the transition between the different phonon transport regimes is investigated. The implementation of the interface modelling is detailed in section II. Results computed by using MC simulations in both simple and double Si/Ge heterostructures are presented in section III. II Modeling of heterostructures
II.1.
Definitions based on semi-spherical temperature concept
This work, in accordance with our previous works [16], [33], is based on the concept of pseudo hemispherical (or directional) temperatures. The definition of the effective/ballistic thermal conductivities 𝜅 𝑒𝑓𝑓 / 𝜅 𝑏𝑎𝑙𝑙 , the total/ballistic thermal conductances and the interface thermal conductance G int are related to these temperatures: Hemispherical temperatures T + /T - are defined as the temperature of phonon sub-populations with a positive/negative (oriented toward the cold/hot thermostat, respectively [16]) giving the same energy density E but with an equilibrium distribution. Thus, it satisfied the following relationship: 𝐸(𝑇 + /𝑇 − ) = ∑ ℏ𝜔 𝑠 𝑓 𝐵𝐸 (𝜔 𝑠 , 𝑇 + /𝑇 − ) 𝑠𝑡𝑎𝑡𝑒 𝑠, 𝑣 𝑥𝑠 >0/𝑣 𝑥𝑠 <0 (1) where f BE is the Bose-Einstein distribution function and 𝜔 𝑠 the phonon pulsation of state s . It should be mentioned that the common approach is the use of a local pseudo temperature considering all the phonon population [33]. In a 1D homogeneous system of length L in contact with a hot thermostat at a hot temperature T H and a cold one at a cold temperature T C , the thermal flux Q can be expressed (in all transport regimes) by using these different expressions: 𝑄 = 𝐺 𝑡𝑜𝑡𝑎𝑙 . Δ𝑇 𝑐𝑜𝑛𝑡𝑎𝑐𝑡 = 𝜅 𝑒𝑓𝑓 Δ𝑇 𝑐𝑜𝑛𝑡𝑎𝑐𝑡 𝐿 = 𝜅 𝑏𝑎𝑙𝑙𝑖𝑠𝑡𝑖𝑐 Δ𝑇 𝑙𝑜𝑐𝑎𝑙 𝐿 = 𝐺 𝑏𝑎𝑙𝑙𝑖𝑠𝑡𝑖𝑐 Δ𝑇 𝑙𝑜𝑐𝑎𝑙 (2) where the two temperature differences are defined as follows: Δ𝑇 𝑙𝑜𝑐𝑎𝑙 (𝑥) = 𝑇 + (𝑥) − 𝑇 − (𝑥) ∆𝑇 𝑐𝑜𝑛𝑡𝑎𝑐𝑡𝑠 = 𝑇 𝐻 − 𝑇 𝐶 . For a homogeneous structure in which a temperature bias ∆𝑇 𝑐𝑜𝑛𝑡𝑎𝑐𝑡𝑠 is applied and working at (pseudo) temperature 𝑇̅ , its thermal conductivity 𝜅 𝑒𝑓𝑓𝑒𝑐𝑡𝑖𝑣𝑒 is given by: 𝜅 𝑒𝑓𝑓𝑒𝑐𝑡𝑖𝑣𝑒 = Q∆𝑇 𝑐𝑜𝑛𝑡𝑎𝑐𝑡𝑠 = Ω(2𝜋) ∑ ℏ𝜔 𝑠 |𝑣 𝑠,𝑥 |𝜆 𝑚𝑓𝑝,𝑠 . 𝜕𝑓 𝐵𝐸 𝜕𝑇 (𝜔 𝑠 , 𝑇̅) 𝑠𝑡𝑎𝑡𝑒 𝑠 (3) where Ω is the volume of one phonon state in the first Brillouin zone, 𝜆 𝑚𝑓𝑝,𝑠 is the mean free path of mode s and the phonon group velocity of the state s along the x -direction which is the heat transport direction. The associated conductance 𝐺 is given by 𝐺 = 𝜅 𝑒𝑓𝑓𝑒𝑐𝑡𝑖𝑣𝑒 𝐿 . In a ballistic system 𝜆 𝑚𝑓𝑝,𝑠 =L/2 and thus the ballistic conductivity 𝜅 𝑏𝑎𝑙 is given by: 𝜅 𝑏𝑎𝑙 = 𝐿 Ω(2𝜋) ∑ ℏ𝜔 𝑠 |𝑣 𝑠,𝑥 | 𝜕𝑓 𝐵𝐸 𝜕𝑇 (𝜔 𝑠 , 𝑇̅) 𝑠𝑡𝑎𝑡𝑒 𝑠 (4) Eq. (3) and (4) are used to compute the semi-analytical results presented in this paper. In heterostructures, interfaces have their own thermal conductance. The interface thermal conductance is here defined as follows: 𝐺 𝑖𝑛𝑡 = 𝑄Δ𝑇 𝐼 = 𝑄𝑇 + (𝑥 − 𝜀) − 𝑇 − (𝑥 + 𝜀) (5) where the local temperature at interface Δ𝑇 𝐼 = 𝑇 + (𝑥 − 𝜀) − 𝑇 − (𝑥 + 𝜀) is the difference between the hemispherical temperatures on each side of the interface. II.2.
Full-band Monte Carlo simulations 1.
Full-band description of phonons
In this study, two materials have been investigated: cubic Si and cubic Ge. Performing Full-Band Monte Carlo simulations requires the knowledge of the phonon dispersion, phonon group velocity and phonon scattering rate. The phonon dispersion in both cubic Si and Ge have been calculated by using the Adiabatic Bond Charge Model (ABCM)[36] [37]. The first Brillouin zone is discretized into a set of 29791 (31×31×31) wave vectors. The phonon band structure in cubic Ge is shown in Fig. 1a. The details about the fitting method and phonon dispersion in cubic Si are given in Ref. [38].
Fig. 1 (a).Phonon dispersion in cubic Ge fitted from ABCM method. (b) Cartography of phonon lifetime in cubic Ge. The colorbar indicates the density of states.
The phonon-phonon scattering rates in cubic Si have been computed by using the Density Functional Theory (DFT) as previously done in[35], [39]. For Ge, the same spectral distribution as in Si is used (as Si and Ge have the same crystallographic structure) but a proportionality factor on their amplitude is applied to fit the Ge experimental thermal conductivity at ambient temperature. The spectral phonon lifetimes in Ge used in this work, which are the inverse of the phonon scattering rate, are shown in Fig. 1b. Monte Carlo simulation for phonon
The particle Monte Carlo (MC) method is a stochastic method to solve the time dependent Boltzmann Transport Equation (BTE) without any assumption on the shape of the phonon distribution functions. The full-band phonon dispersion and scattering rates in the materials are used as inputs for the Monte Carlo simulation. At each simulation step, the values of the pseudo and hemispherical temperatures at each position in the structure are updated according to the local phonon distributions. The details of the Monte Carlo simulation for phonon employed in this work is presented in Ref.[33]. Several kinds of boundaries are considered in the MC simulation: semi-infinite and rough boundaries. Semi-infinite boundaries are defined by using periodic boundary conditions. In that purpose a couple of opposite external faces (i.e. that are located face to face) must be defined a priori. At these boundaries, when a particle exits from one face a particle with exactly the same properties is created on the opposite (and associated) face. In the case of a rough boundary, when a phonon collides to the boundary, two kinds of elastic reflection can occur: specular or diffusive reflection. The probability of a specular reflection for an incident phonon of wave vector 𝒒 and incident angle 𝜃 is provided by the Soffer’s model [40]: 𝑝 𝑠𝑝𝑒𝑐𝑢𝑙𝑎𝑟 (|𝒒 |, cos 𝜃 ) = exp(−(2 cos 𝜃 Δ|𝒒 |) ) (6) where Δ is an empirical parameter characterizing the roughness of the interface. In the case of specular reflection, the normal component to the boundary 𝑣 ⊥ of the reflected phonon wave vector is changed: 𝑞 𝑗 ′ ,⊥ = −𝑞 𝑗,⊥ (that is not obvious in the case of a full-band description if the reflection surface is not oriented along a high symmetry plane of the material) while its component parallel to the boundary is conserved: 𝑞 𝑗 ′ ,∕∕ = 𝑞 𝑗,∕∕ . It should be mentioned that as a specular reflection has no effect on the heat flux parallel to the interface, it can also be used to model semi-infinite boundaries. In the case of a diffusive reflection (with a probability of 𝑠𝑝𝑒𝑐𝑢𝑙𝑎𝑟 .), the reflected phonon wave vector is randomly selected on an iso-energy surface because the phonon energy is conserved and satisfies the orientation condition: 𝑣 𝑗 ′ ,⊥ 𝑣 𝑗,⊥ < 0 (reflection). To avoid non-physical accumulation of phonons near the rough boundary, the probability of the final state among the iso-energy states is non-uniform. Indeed, according to the Lambert’s law, the probability of the state j’ with a reflected angle 𝜃 𝑟 is proportional to the normal component of the group velocity to the boundary [33]: 𝑝 𝑑𝑖𝑓𝑓𝑢𝑠𝑖𝑣𝑒 (𝐪 𝑗′ , 𝜃 𝑟 ) ∝ |𝐯 𝑗 ′ (𝐪 𝑗 ′ )| cos 𝜃 𝑟 (7) Heterojunction modeling in Monte Carlo simulation
Fig. 2. Schema of phonon diffusion mechanisms at an interface
As schematized in Fig. 2, the interface between materials A and B is modeled by using the Full band version of the Diffusive Mismatch Model (DMM) [12]. Even if this approach does not capture all the complex physic of the heat transfer occurring at the interface [41] and in particular inelastic phenomena are neglected, this approach is easy to implement and provides reasonable prediction of interface thermal conductance [42]. In DMM it is assumed that all transmitted or reflected phonons have undergone a diffusive scattering at the interface (no specular reflection occurs) and the interface scattering is assumed to be elastic. For an incident phonon of wave vector 𝐪 𝑗 and angular frequency 𝜔 colliding this interface from the A-side, the spectral transmission probability 𝑇 𝐴→𝐵 (𝜔) is given by [42]: 𝑇 𝐴→𝐵 (𝜔) = 𝐼
𝐵,𝐧 (𝜔)𝐼
𝐴,𝐧 (𝜔) + 𝐼
𝐵,𝐧 (𝜔) (8) where 𝐼 𝐴 𝐵⁄ ,𝐧 (𝜔) is expressed by: 𝐼 𝐴,𝐧 (𝜔) = 12 . Ω(2𝜋) ∑ 𝛿(𝜔 𝑠 − 𝜔)|𝐯 𝑠 ∙ 𝐧| 𝑠𝑡𝑎𝑡𝑒𝑠 𝑠 𝑖𝑛 𝐵𝑍 𝐴 (9) Where 𝐧 is a unit vector perpendicular to the interface. The summation is performed over all the phonon states s in the materials A in its entire first Brillouin zone 𝐵𝑍 𝐴 . Hence, when a phonon collides the interface, the reflection (phonon remains in the same material) or transmission (the phonon is transferred to the other side of the interface) process are selected randomly according to the transmission probability presented in Eq.(9). Then, the wave vector and the mode of the final phonon are randomly selected (independently from the incident state) among the isoenergy states j in the final material with a the non-uniform probability 𝑃 𝑗 to satisfy the Lambert’s cosine law i.e.: 𝑃 𝑗 (𝐪 𝑗′ , 𝜃 𝑟 ) ∝ |𝐯 𝑗 ′ (𝐪 𝑗 ′ )| cos 𝜃 𝑟 . Investigated structures
In this work, homogenous and heterogenous nanostructures of length L along the X axis are investigated. In all structures, the thermal flux is along the X axis. The cells located at the two extrema have their external face placed in contact with a hot thermostat of temperature 𝑇 𝐻 and a cold thermostat of temperature 𝑇 𝐶 , respectively. All nanostructures have square cross-sections and a cubic mesh is employed. First, all structures are uniformly meshed along the X-axis into 20 equally sized cells. Then, the meshing of the heterojunctions is refined about the interfaces and the thermostats. The length of these refined cells located 5 nm around the interfaces are shorter or equal than 1 nm. Fig. 3. Homogeneous nanostrucrures: (a) Cross-plane nanofilm (CPNF). (b) In-plane nanofilm (IPNF). (c) Nanowires (NW). Red/blue faces are hot/cold thermostats, respectively. Transparent external faces: periodic boundaries. Green faces: rough boundaries.
Fig. 4. Heterogeneous nanostructures: (a) Simple heterojunction. (b) Double heterojunction. Red/blue faces are hot/cold thermostats. Transparent external faces: periodic boundaries. Yellow faces: DMM interfaces.
Three configurations of homogeneous nanostructures illustrated in Fig. 3 are studied: • Cross-plane nanofilms (CPNF) in which periodic boundary conditions are applied along both Y and Z directions to describe semi-infinite boundaries. • In-plane nanofilms (IPNF) in which periodic boundary conditions are applied for XZ planes to describe semi-infinite boundaries and diffusive boundary conditions are applied to XY planes to represent external faces (rough boundaries). • Nanowires (NW) in which the diffusive boundary conditions are applied to all XY and XZ planes. In addition, two different kinds of heterogeneous nanostructures in the cross-plane configuration schematized in Fig. 4 are studied: • Simple heterojunctions in which a DMM heterojunction (in yellow) is located in the middle of the heterojunction ( 𝑥 = 𝐿 2⁄ ) joining two homogeneous bars. • Double heterojunctions in which three equally long homogeneous bars are joined to form two DMM interfaces at positions 𝑥 = 𝐿 3 ⁄ and 𝑥 = 2𝐿 3⁄ . III
Results
III.1.
Ge homogeneous nanostructures
As shown previously for Si in Ref [33], Fig. 5 shows the variation of the thermal conductivity of Ge in CPNF, IPNF and NW as a function of the nanostructure length. When the nanostructure length L becomes smaller than a few tens of nanometers which is the order of magnitude of the effective phonon mean free path in these structures, the phonon transport regime tends to be more and more ballistic. The evolution of thermal conductivities of the three nanostructures have the same behaviors: they tend asymptotically to the ballistic thermal conductivity of Ge (dashed line), which is a linear function of L . In long nanostructures i.e. much longer than the phonon mean free path, the diffusive transport regime occurs, and the thermal conductivities tend to become L independent. In CPNF, the thermal conductivity of bulk Ge of 60 W.m -1 .K -1 at 300 K [43] tends to be recovered. In nanostructures, the rougher interfaces there are in the system, the lower the thermal conductivities. Thus, at a given system length, NWs exhibit a lower conductivity than the IPNFs which have a lower conductivity that the CPNFs. Fig. 5. Thermal conductivity of Ge homogeneous nanostructrues of the three homogeneous nanostructures as a function of the nanostructure length L. a) b)
Fig. 6. The ratio 𝜅 𝑒𝑓𝑓𝑒𝑐𝑡𝑖𝑣𝑒 𝜅 𝑏𝑎𝑙𝑙𝑖𝑠𝑡𝑖𝑐 ⁄ as a function of the length L. (a) Si nanostructures (b) Ge nanostructures. To quantify the transition between different phonon transport regimes, the evolution of the Knudsen number 𝐾 𝐷 is widely used. 𝐾 𝐷 , characterizing the degree of ballisticity of the transport, is defined as the ratio between the effective conductivity 𝜅 𝑒𝑓𝑓𝑒𝑐𝑡𝑖𝑣𝑒 (cf. definition of Eq. 3) and the ballistic one 𝜅 𝑏𝑎𝑙 (cf Eq. 4). It depends on the nanostructures as follows: 𝐾 𝐷 = 𝜅 𝑒𝑓𝑓𝑒𝑐𝑡𝑖𝑣𝑒 𝜅 𝑏𝑎𝑙 ⁄ = {𝜅 𝐶𝑃𝑁𝐹 𝜅 𝑏𝑎𝑙 ⁄ 𝑓𝑜𝑟 𝐶𝑃𝑁𝐹𝜅 𝐼𝑃𝑁𝐹 𝜅 𝑏𝑎𝑙 ⁄ 𝑓𝑜𝑟 𝐼𝑃𝑁𝐹 𝜅 𝑁𝑊 𝜅 𝑏𝑎𝑙 ⁄ 𝑓𝑜𝑟 𝑁𝑊 𝐾 𝐷 computed by using our MC code is presented in Fig. 6 for different nanostructures of length L varying from 1 nm to 1 μm made of Si (a) and Ge (b). As expected , the value of 𝐾 𝐷 decreases monotonically from 1 in ultra-short systems where a fully ballistic transport regime occurs down to 0 in ultra-long systems where a fully diffusive transport regime occurs. The intermediate transport regime is defined here as the regime where K D values are comprised between 0.8 (L<10 nm) and 0.1 ( L >1µm). III.2.
Simple Si/Ge heterostructure
Fig. 7 shows the profiles of the pseudo temperature T and the two hemispherical temperatures T + and T - in a 100 nm-long Ge/Si heterojunction. The Ge bar is in contact with a hot thermostat at temperature T H = 302 K and the Si bar is in contact with a cold thermostat at temperature T C = 298 K. Along the two homogenous bars, the temperatures vary continuously whereas there is an abrupt drop of the temperatures at the interface. More surprisingly, in our MC simulations all the considered temperatures do not evolve linearly along the bars and there is a significant curvature of the temperature profile near the interface in the Ge side which has a lower thermal conductivity than Si. Thus, a coarse mesh around the interface could lead to a significant overestimation of the temperature drop at the interface and an underestimation of the related interface thermal conductance. For this reason, the refinement of the mesh in the vicinity of the interfaces is mandatory. Fig. 7. Temperature profiles in a 100 nm – long Ge/Si heterojunction
In a 20 nm – long Si/Ge heterojunction, the heat flux computed by using our Monte Carlo simulations are plotted in Fig. 8 as for different temperature difference (T H -T C ) between thermostats while keeping the working temperature of the heterojunction at 300 K. The Si/Ge heterojunction is considered for both configurations: (i) T H :Ge/ T C :Si in which the Ge bar is in contact with the hot thermostat while the Si bar is in contact with the cold thermostat; and (ii) T H : Si/ T C :Ge in which the Si bar is in contact with the hot thermostat while the Ge bar is in contact with the cold thermostat. Fig. 8 a clearly shows a linear relationship between the heat flux density and the temperature difference (T H -T C ) for both thermal polarizations. Using linear regressions, the following relationships are obtained: i. T H :Ge/T C :Si : 𝑄 = (191.30 MW. m −2 . K −1 )(𝑇 𝐻 − 𝑇 𝐶 ) + 14.40 MW. m −2 ii. T H :Si/T C :Ge : 𝑄 = (191.35 MW. m −2 . K −1 )(𝑇 𝐻 − 𝑇 𝐶 ) − 16.24 MW. m −2 A total thermal conductance of 191 MW.m -2 .K -1 , independent of the temperature difference configuration is obtained in the 20 nm–long heterojunction. Due to the numerical nature of the Monte Carlo method, the heat flux density at zero temperature bias does not vanish, but it takes a nonzero value which is around -/+16 MW.m -2 in this 20 nm long device. We assumed that this value corresponds to the maximum numerical resolution 𝛿𝑄 of our simulation. Besides, in an MC algorithm the interval of 95% confidence of the flux called 𝐼𝐶 𝑄 can be reduced as small as expected if the number of samples is large enough. Finally, the error bars for heat flux density are defined here as 𝐸𝑟𝑟𝑜𝑟𝑏𝑎𝑟 𝑄 = 𝑚𝑖𝑛{𝛿𝑄, 𝐼𝐶 𝑄} . Similarly, the error bars for the interface thermal conductance are defined by summing: 𝛿𝐺 𝑖𝑛𝑡 𝐺 𝑖𝑛𝑡 = 𝛿𝑄𝑄 + 𝛿(∆𝑇 𝐼 )∆𝑇 𝐼
0 where the estimation value of 𝛿(∆𝑇 𝐼 ) is 𝛿𝑇 = 𝛿𝑄 𝐺 𝑡𝑜𝑡 ⁄ . (a) (b) Fig. 8. (a) Heat flux density and (b) interface thermal conductance in a 20 nm – long Si/Ge heterojunction for different temperature differences (T H -T C ). Fig. 8 b presents the interface thermal conductance as a function of the temperature difference (T H -T C ) . For temperature difference greater than or equal to 2K, the interface thermal conductance converges to a value of about 240 MW.m -2 . K -1 that is consistent with the experimental value of Ref. [26], [44]. Although the error bars are large at low temperature bias, due to the numerical fluctuation of the Monte Carlo method, the value of 𝐺 𝑖𝑛𝑡 is still within the error bars, indicating that 𝐺 𝑖𝑛𝑡 is independent on the temperature difference (𝑇 𝐻 − 𝑇 𝐶 ) . . In Ref. [45] a rectification of 5% was estimated by using MD at the Si/Ge interfaces. This effect could be due to the asymmetric phonon distributions on each side of the interface, phenomenon captured by our approach, and/or the inelastic scatterings at the interfaces which are not included here. As there is an overlap of the error bar extension in the two thermal configurations, we can conclude here, but a small thermal rectification less than 5% can be observed at the Si/Ge interface in our MC simulations. A set of Si/Ge heterojunction of length L varying from 1 nm to 1 μm with thermostat temperatures of T H = 302 K and T C = 298 K are considered. The heat flux density of the heterojunctions as a function of the heterostructure length L is plotted in Fig. 9 a . In ultra-short heterostructure ( L<
10 nm), the heat flux density Q approaches the value of 1 GW.m -2 in the ballistic regime. In long heterostructures ( L >200 nm) in which the heat transport is diffusive, the heat flux density tends to be proportional to consistently with the classical Fourier’s law. The related thermal conductances G int are plotted in Fig. 9 b . As the values of G int remains within the error bars, the interface thermal conductance has values of 240 MW.m -2 K -1 , for both temperature difference configuration and for all length L . Once again, no significant thermal rectification is observed with respect to our error bars of around 5%. It is noted that the width of the error bars is larger for long heterostructures than for short heterostructures. Another interesting parameters of an interface is its Kapitza length l defined as l = 𝜅 𝑑𝑖𝑓𝑓𝑢𝑠𝑖𝑣𝑒 /G int where 𝜅 𝑑𝑖𝑓𝑓𝑢𝑠𝑖𝑣𝑒 is the thermal conductivity. It corresponds to the length of a homogenous bar made of material exhibiting the same thermal conductance as one interface. Here, the Kapitza lengths are 540 nm and 250 nm in Si and Ge, respectively. 1 (a) (b) Fig. 9. (a) Heat flux density and (b) interface thermal conductance of Si/Ge heterojunction as a function of the heterostructure length L (T H =302 K, T C =298 K). According to our Monte Carlo algorithm, the total thermal flux density passing through the double heterojunction is 239.75±21.37 MW.m -2 . The calculated values for the two interface thermal conductances are G =271.34±31.78 MW.m -2 .K -1 (at x=100 nm) and G =234.95±27.52 MW.m -2 .K -1 (at x=200 nm). With respect to the value extracted in simple heterostructure (of 240 MW.m -2 .K -1 ), the values extracted in double heterostructure remain within the error bar of the method. Fig. 10. Temperature profile in Ge/Si/Ge double heterojunction.
In addition, the effective thermal conductivities of the Ge and Si bars are 17.25±1.71 W.m -1 .K -1 , 46.96±5.85 W.m -1 .K -1 and 17.29±1.71 W.m -1 .K -1 , for the left Ge bar (0≤x≤100 nm), the middle Si bar (100 nm≤x≤200 nm), and the right Ge bar (200 nm≤x≤300 nm), respectively. These results match the value obtained in 100 nm-long CPNF and shown in Ref [33] and Fig 5, respectively. These results illustrate the consistence of our simulation methodology even if complex geometries are simulated. It is important to note that these consistencies between our MC results in terms of interface conductances and thermal conductivities are obtained thanks to the use of the refinement of the 2 mesh near the interface and the used of the hemispherical temperature are mandatory. Otherwise, the computed interface thermal conductance appears artificially size dependent even using DMM transmission coefficients. Finally, in Fig. 11, the total conductance of the simple and double heterostructures computed by using our MC simulations, are plotted as a function of the device length L . In ultra-small devices where the ballistic transport occurs, the total conductance in heterostructures is governed by the Si/Ge interface conductance. In long device, much longer than the Kapitza length of the interface, the effect of the interfaces become negligible. The main interest of the presented MC approach is to capture the heat transfer in the intermediate regime even in complex geometry such as double heterostructures of few tens of nanometers long. Fig. 11. Total conductance of simple and double Si/Ge heterojunctions cf. Fig 4 of various lengths L. IV Conclusion
In this paper, our Full Band Monte Carlo simulation considering phonon transmission across interfaces that are perpendicular to the heat flux have been presented. The phonon transmissions are computed by using a Full band version of the DMM formalism. These simulations have been used to study single and double heterostructures made of Si/Ge interfaces. A clear methodology to estimate the numerical resolution of the method has been developed. By considering the hemispherical temperature, i.e. by distinguishing the phonons by the direction of their displacement and by carefully refining the mesh around the interfaces in order to capture the curvature of the temperature profile, we have calculated thermal interface conductance that are not size dependent or heterostructure dependent. 3 This original Monte simulator providing a deep insight of the phonon transport in particular in terms of spectral properties and transient responses can be used to study complex geometries with many interfaces of different types: rough/specular, reflecting, or semi-transparent interfaces. V References [1] W. Liu, Q. Jie, H. S. Kim, and Z. Ren, “Current progress and future challenges in thermoelectric power generation: From materials to devices,”
Acta Materialia , vol. 87, pp. 357–376, Apr. 2015, doi: 10.1016/j.actamat.2014.12.042. [Online]. Available: https://linkinghub.elsevier.com/retrieve/pii/S1359645414009689. [Accessed: 25-Nov-2020] [2] “THE INTERNATIONAL ROADMAP FOR DEVICES AND SYSTEMS:MORE MOORE Report 2017,” IRDS, 2017 [Online]. Available: https://irds.ieee.org/images/files/pdf/2017/2017IRDS_MM.pdf [3] J. B. J. baron Fourier,
Théorie analytique de la chaleur . F. Didot, 1822. [4] D. G. Cahill et al. , “Nanoscale thermal transport. II. 2003–2012,”
Applied Physics Reviews , vol. 1, no. 1, p. 011305, Mar. 2014, doi: 10.1063/1.4832615. [Online]. Available: http://aip.scitation.org/doi/10.1063/1.4832615. [Accessed: 26-Sep-2019] [5] J. Kaiser, T. Feng, J. Maassen, X. Wang, X. Ruan, and M. Lundstrom, “Thermal transport at the nanoscale: A Fourier’s law vs. phonon Boltzmann equation study,”
Journal of Applied Physics , vol. 121, no. 4, p. 044302, Jan. 2017, doi: 10.1063/1.4974872. [Online]. Available: http://aip.scitation.org/doi/10.1063/1.4974872. [Accessed: 17-Jul-2020] [6] J.-P. M. Péraud and N. G. Hadjiconstantinou, “Extending the range of validity of Fourier’s law into the kinetic transport regime via asymptotic solution of the phonon Boltzmann transport equation,”
Phys. Rev. B , vol. 93, no. 4, p. 045424, Jan. 2016, doi: 10.1103/PhysRevB.93.045424. [Online]. Available: https://link.aps.org/doi/10.1103/PhysRevB.93.045424. [Accessed: 07-Sep-2020] [7] A. T. Ramu and J. E. Bowers, “A Generalized Enhanced Fourier Law,”
Journal of Heat Transfer , vol. 139, no. 3, p. 034501, Mar. 2017, doi: 10.1115/1.4034796. [Online]. Available: https://asmedigitalcollection.asme.org/heattransfer/article/doi/10.1115/1.4034796/449448/A-Generalized-Enhanced-Fourier-Law. [Accessed: 17-Jul-2020] [8] J. Ordonez-Miranda, R. Yang, and J. J. Alvarado-Gil, “A constitutive equation for nano-to-macro-scale heat conduction based on the Boltzmann transport equation,”
Journal of Applied Physics , vol. 109, no. 8, p. 084319, Apr. 2011, doi: 10.1063/1.3573512. [Online]. Available: http://aip.scitation.org/doi/10.1063/1.3573512. [Accessed: 02-Sep-2020] [9] P. Kapitza, “The study of heat transfer in helium II,”
J. Phys.(USSR) , vol. 4, no. 1–6, pp. 181–210, 1941. [10] W. A. Little, “The transport of heat between dissimilar solids at low temperatures,”
Can. J. Phys , vol. 37, no. 49, p. 334, 1959. [11] E. T. Swartz and R. O. Pohl, “Thermal boundary resistance,”
Reviews of modern physics , vol. 61, no. 3, p. 605, 1989 [Online]. Available: http://journals.aps.org/rmp/abstract/10.1103/RevModPhys.61.605. [Accessed: 25-Feb-2016] [12] P. Reddy, K. Castelino, and A. Majumdar, “Diffuse mismatch model of thermal boundary conductance using exact phonon dispersion,”
Applied Physics Letters , vol. 87, no. 21, p. 211908, 2005, doi: 10.1063/1.2133890. [Online]. Available: 4 http://scitation.aip.org/content/aip/journal/apl/87/21/10.1063/1.2133890. [Accessed: 26-Feb-2016] [13] J. Larroque, P. Dollfus, and J. Saint-Martin, “Phonon transmission at Si/Ge and polytypic Ge interfaces using full-band mismatch based models,”
Journal of Applied Physics , vol. 123, no. 2, p. 025702, Jan. 2018, doi: 10.1063/1.5007034. [Online]. Available: http://aip.scitation.org/doi/10.1063/1.5007034. [Accessed: 26-Sep-2019] [14] S. Merabia and K. Termentzidis, “Thermal conductance at the interface between crystals using equilibrium and nonequilibrium molecular dynamics,”
Phys. Rev. B , vol. 86, no. 9, p. 094303, Sep. 2012, doi: 10.1103/PhysRevB.86.094303. [Online]. Available: https://link.aps.org/doi/10.1103/PhysRevB.86.094303. [Accessed: 11-Mar-2020] [15] G. Chen, “Thermal conductivity and ballistic-phonon transport in the cross-plane direction of superlattices,”
Phys. Rev. B , vol. 57, no. 23, pp. 14958–14973, Jun. 1998, doi: 10.1103/PhysRevB.57.14958. [Online]. Available: https://link.aps.org/doi/10.1103/PhysRevB.57.14958. [Accessed: 01-Oct-2019] [16] B. Davier et al. , “Revisiting thermal conductivity conductance at the nanoscale Homogenous system,” arXiv , 2020. [17] Y. Chalopin, K. Esfarjani, A. Henry, S. Volz, and G. Chen, “Thermal interface conductance in Si/Ge superlattices by equilibrium molecular dynamics,”
Physical Review B , vol. 85, no. 19, p. 195302, May 2012, doi: 10.1103/PhysRevB.85.195302. [Online]. Available: http://link.aps.org/doi/10.1103/PhysRevB.85.195302. [Accessed: 10-Dec-2014] [18] D. P. Sellan, E. S. Landry, J. E. Turney, A. J. H. McGaughey, and C. H. Amon, “Size effects in molecular dynamics thermal conductivity predictions,”
Physical Review B , vol. 81, no. 21, p. 214305, Jun. 2010, doi: 10.1103/PhysRevB.81.214305. [Online]. Available: http://link.aps.org/doi/10.1103/PhysRevB.81.214305. [Accessed: 09-Dec-2015] [19] E. S. Landry and A. J. H. McGaughey, “Thermal boundary resistance predictions from molecular dynamics simulations and theoretical calculations,”
Physical Review B , vol. 80, no. 16, pp. 165304–1, Oct. 2009, doi: 10.1103/PhysRevB.80.165304. [Online]. Available: https://link.aps.org/doi/10.1103/PhysRevB.80.165304. [Accessed: 11-Sep-2017] [20] B. Qiu and X. Ruan, “Molecular dynamics simulations of lattice thermal conductivity of bismuth telluride using two-body interatomic potentials,”
PHYSICAL REVIEW B , p. 6. [21] T.-Q. Duong, C. Massobrio, G. Ori, M. Boero, and E. Martin, “Thermal resistance of an interfacial molecular layer by first-principles molecular dynamics,”
The Journal , p. 6, 2020. [22] N. Mingo, L. Yang, D. Li, and A. Majumdar, “Predicting the Thermal Conductivity of Si and Ge Nanowires,”
Nano Letters , vol. 3, no. 12, pp. 1713–1716, Dec. 2003, doi: 10.1021/nl034721i. [Online]. Available: http://pubs.acs.org/doi/abs/10.1021/nl034721i. [Accessed: 01-Jul-2015] [23] F. Mazzamuto et al. , “Enhanced thermoelectric properties in graphene nanoribbons by resonant tunneling of electrons,”
Phys. Rev. B , vol. 83, no. 23, p. 235426, Jun. 2011, doi: 10.1103/PhysRevB.83.235426. [Online]. Available: http://link.aps.org/doi/10.1103/PhysRevB.83.235426 [24] A. Alkurdi, S. Pailhès, and S. Merabia, “Critical angle for interfacial phonon scattering: Results from ab initio lattice dynamics calculations,”
Applied Physics Letters , vol. 111, no. 9, p. 093101, Aug. 2017, doi: 10.1063/1.4997912. [Online]. Available: http://aip.scitation.org/doi/10.1063/1.4997912. [Accessed: 11-Sep-2017] [25] J. Carrete et al. , “Phonon transport across crystal-phase interfaces and twin boundaries in semiconducting nanowires,”
Nanoscale , vol. 11, no. 34, pp. 16007–16016, 2019, doi: 5 10.1039/C9NR05274G. [Online]. Available: http://xlink.rsc.org/?DOI=C9NR05274G. [Accessed: 09-Dec-2020] [26] Y. Lee et al. , “Quantum treatment of phonon scattering for modeling of three-dimensional atomistic transport,”
PHYSICAL REVIEW B , p. 6, 2017. [27] S. Mazumder and A. Majumdar, “Monte Carlo Study of Phonon Transport in Solid Thin Films Including Dispersion and Polarization,”
Journal of Heat Transfer , vol. 123, no. 4, p. 749, 2001, doi: 10.1115/1.1377018. [Online]. Available: http://HeatTransfer.asmedigitalcollection.asme.org/article.aspx?articleid=1445087. [Accessed: 13-Jul-2016] [28] D. Lacroix, K. Joulain, and D. Lemonnier, “Monte Carlo transient phonon transport in silicon and germanium at nanoscales,”
Physical Review B , vol. 72, no. 6, Aug. 2005, doi: 10.1103/PhysRevB.72.064305. [Online]. Available: http://link.aps.org/doi/10.1103/PhysRevB.72.064305. [Accessed: 05-Jul-2016] [29] J.-P. M. Péraud, C. D. Landon, and N. G. Hadjiconstantinou, “Monte carlo methods for solving the boltzmann transport equation,”
Annual Review of Heat Transfer
Front. Energy Res. , vol. 6, p. 28, May 2018, doi: 10.3389/fenrg.2018.00028. [Online]. Available: http://journal.frontiersin.org/article/10.3389/fenrg.2018.00028/full. [Accessed: 21-Oct-2019] [31] M.-S. Jeng, R. Yang, D. Song, and G. Chen, “Modeling the Thermal Conductivity and Phonon Transport in Nanoparticle Composites Using Monte Carlo Simulation,”
J. Heat Transfer , vol. 130, no. 4, p. 042410, 2008, doi: 10.1115/1.2818765. [Online]. Available: http://HeatTransfer.asmedigitalcollection.asme.org/article.aspx?articleid=1449099. [Accessed: 18-Dec-2020] [32] T. Hori, J. Shiomi, and C. Dames, “Effective phonon mean free path in polycrystalline nanostructures,”
Appl. Phys. Lett. , vol. 106, no. 17, p. 171901, Apr. 2015, doi: 10.1063/1.4918703. [Online]. Available: http://aip.scitation.org/doi/10.1063/1.4918703. [Accessed: 18-Dec-2020] [33] B. Davier et al. , “Heat transfer in rough nanofilms and nanowires using full band ab initio
Monte Carlo simulation,”
Journal of Physics: Condensed Matter , vol. 30, no. 49, p. 495902, Dec. 2018, doi: 10.1088/1361-648X/aaea4f. [34] L. Yang and A. J. Minnich, “Thermal transport in nanocrystalline Si and SiGe by ab initio based Monte Carlo simulation,”
Scientific Reports
Ab initio based calculations of the thermal conductivity at the micron scale,”
Appl. Phys. Lett. , vol. 112, no. 3, p. 033104, Jan. 2018, doi: 10.1063/1.5010959. [36] W. Weber, “Adiabatic bond charge model for the phonons in diamond, Si, Ge, and α-Sn,”
Physical Review B , vol. 15, no. 10, pp. 4789–4803, 1977, doi: 10.1103/PhysRevB.15.4789. [37] A. Valentin, J. Sée, S. Galdin-Retailleau, and P. Dollfus, “Study of phonon modes in silicon nanocrystals using the adiabatic bond charge model,”
J. Phys.: Condens. Matter , vol. 20, no. 14, p. 145213, Apr. 2008, doi: 10.1088/0953-8984/20/14/145213. 6 [38] J. Larroque, P. Dollfus, and J. Saint-Martin, “Full-Band modelling of phonons in polytype Ge and Si,”
J. Phys.: Conf. Ser. , vol. 906, p. 012007, Oct. 2017, doi: 10.1088/1742-6596/906/1/012007. [39] A. Togo, L. Chaput, and I. Tanaka, “Distributions of phonon lifetimes in Brillouin zones,”
Phys. Rev. B , vol. 91, no. 9, p. 094306, Mar. 2015, doi: 10.1103/PhysRevB.91.094306. [40] S. B. Soffer, “Statistical Model for the Size Effect in Electrical Conduction,”
Journal of Applied Physics , vol. 38, no. 4, pp. 1710–1715, Mar. 1967, doi: 10.1063/1.1709746. [41] S. Sadasivam et al. , “Thermal transport across metal silicide-silicon interfaces: First-principles calculations and Green’s function transport simulations,”
Phys. Rev. B , vol. 95, no. 8, p. 085310, Feb. 2017, doi: 10.1103/PhysRevB.95.085310. [Online]. Available: https://link.aps.org/doi/10.1103/PhysRevB.95.085310. [Accessed: 15-Feb-2021] [42] J. Larroque, P. Dollfus, and J. Saint-Martin, “Phonon transmission at Si/Ge and polytypic Ge interfaces using full-band mismatch based models,”
Journal of Applied Physics , vol. 123, no. 2, pp. 025702–025702, 2018, doi: 10.1063/1.5007034. [43] C. J. Glassbrenner and G. A. Slack, “Thermal Conductivity of Silicon and Germanium from 3°K to the Melting Point,”
Phys. Rev. , vol. 134, no. 4A, pp. A1058–A1069, May 1964, doi: 10.1103/PhysRev.134.A1058. [44] M. Goto et al. , “Ultra-low thermal conductivity of high-interface density Si/Ge amorphous multilayers,”
Appl. Phys. Express , vol. 11, no. 4, p. 045202, Apr. 2018, doi: 10.7567/APEX.11.045202. [Online]. Available: https://iopscience.iop.org/article/10.7567/APEX.11.045202. [Accessed: 11-Dec-2020] [45] R. Rurali, X. Cartoixà, and L. Colombo, “Heat transport across a SiGe nanowire axial junction: Interface thermal resistance and thermal rectification,”