Journal of Nanomaterials

Volume 2019, Article ID 5862979, 9 pages

https://doi.org/10.1155/2019/5862979

## Thermal Transport in Silicon-Germanium Superlattices at Low Temperatures

College of Mechanical and Electrical Engineering, Henan University of Technology, Zhengzhou 450007, China

Correspondence should be addressed to Zan Wang; moc.361@zweurt

Received 7 May 2019; Revised 20 August 2019; Accepted 4 September 2019; Published 15 October 2019

Academic Editor: Balachandran Jeyadevan

Copyright © 2019 Zan Wang et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### Abstract

Interfacial thermal resistances between heterogeneous materials are still a challengeable subject since the mechanism to explain it quantitatively is not clear in spite of its importance. We propose a Monte Carlo (MC) model to study phonon interfacial elastic and inelastic scattering behaviors for superlattices composed of Si and Ge materials, which substantially reduces the amount of computations. In particular, below Debye temperatures, the molecular dynamics (MD) solution is not precise enough for semiconductors because of quantization errors. In this work, thermal conductivities and thermal rectifications of Si/Ge and Ge/Si superlattices with different periods are investigated separately at temperatures below 200 K.

#### 1. Introduction

In nanoscale heterostructures, both boundary effect and interfacial effect are the predominant features. Superlattice films are often synthesized by two different materials and form periodic structures. Owing to different lattice constants between two materials, their phonon dispersions are not the same and the phenomenon of lattice mismatch exists, which causes the special thermal properties and electronic characters [1]. Interfacial thermal resistances in superlattices [2–10], such as AlAs/GaAs, Si/Ge, Si/SiGe, and Si_{1−x}Ge_{x}/Si_{1−y}Ge_{y}, have attracted keen interests of many researchers. Rezgui et al. [2] presented an enhanced ballistic-diffusive equation to solve the heat transport across germanium-silicon interface and emphasized the importance of phonon-boundary interactions and surface roughness effect. Gordiz and Henry [3] suggested more than 15% of the interface conductance arising from less than 0.1% of the modes between 12 and 13 THz. Hahn et al. [4] found that thermal boundary resistances for different heat flow direction generated a rectification factor of the Si/Ge interface. Hijazi and Kazan [5] presented a predictive Boltzmann model for the cross-plane thermal conductivity in superlattices. Si-Ge superlattices doped by SiGe and Ge nanodots were also studied theoretically, and the results manifested the contributions of nanodots to limiting the thermal conductivities [11–14].

Fiedler et al. [15] presented a combined experimental and theoretical investigation of the thermoelectric properties of p-doped Si-Ge superlattices and revealed that doping increases their thermal conductivity about 50%. Recently, Singh et al. [19] investigated the effect of phonon dispersion on thermal transport across Si/Ge interfaces and studied the interactions between optical and acoustic phonons. Si-Ge superlattice nanowires were investigated by researchers [16–18] by MD method.

In this work, we apply the MC method to explore phonon transport in Si/Ge and Ge/Si superlattices and study the properties of interfacial phonon scatterings including elastic scattering and inelastic scattering. This study is aimed at developing an approximate MC model to predict the effective thermal conductivity of Si/Ge and Ge/Si superlattices.

The behaviors of phonon transport in nanoscale structure are often studied by analytical and numerical methods. Analytical methods mainly include Boltzmann transport equation (BTE), Green function, and lattice dynamics [20]. Numerical methods mainly include MC method and MD method. Because of the complexities of phonon scattering events, especially in special shaped nanoscale structures and heterostructures, it is difficult to get accurate analytical solutions by analytical methods without simplifying the analytical equation massively. Instead of analytical methods, basing on the principle of statistics, numerical methods are more suitable for simulating random scattering events in the process of phonon transport in complex conditions. In the method of MD, the simulated physical model is kept consistent with the actual crystal structure. Considering atom reciprocal potential such as Stillinger-Weber potential, all the kinetic properties and trajectories of each atom in the physical model should be calculated by Newton’s second law, which will consume a lot of computing time [21, 22]. Besides, below Debye temperature, the MD method is invalid for the exchanges between classical statistics and quantum statistics, which limits its applications at low temperatures [23]. Basing on phonon dispersion distribution, Bose-Einstein distribution, and Matthiessen formula, the MC method is a kind of probability model, in which phonon’s transport properties such as velocity and frequency are just decided by random selections.

#### 2. Methodology and MC Simulation

In this work, the MC method is chosen to study the mechanism of phonon transport in Si/Ge superlattice for two reasons. One is that temperatures of the simulated model are below Debye temperature, which is beyond the applications of the MD method. The other is that the simulated model with a size of contains nearly 10^{7} atoms, which will cost a lot of computational time.

##### 2.1. The Process of MC Simulation

The simulated superlattice is composed of several equivalent cells. Heat flux flows from the leftmost cell defined as the hot part to the rightmost cell defined as the cold part. In simulations, phonons drift from the hot part to the cold part by imposing temperature gradient. One full simulated circulation can be divided into four stages: initialization, drift, interfacial scattering, and renewal.

###### 2.1.1. Initialization

As a semiclassical equation, BTE has been successfully used to investigate the transmissions of an ensemble of particles that interact with each other via short-range forces. The behaviors of phonon transport in superlattices can be modeled using the BTE, which can be written as [24, 25] where is the distribution function of an ensemble of phonons, and is the phonon group velocity. The left side of equation (1) denotes the changes of the distribution function due to drift and causes the departure from the equilibrium. The right side of equation (1) restores the equilibrium distribution by phonon scatterings and can be expressed as [24, 25] where is the relaxation time of phonon impurity scattering, denotes the phonon normal scattering relaxation time, and is the phonon Umklapp scattering relaxation time. is the phonon boundary scattering relaxation time, and is the phonon interfacial scattering relaxation time.

At the beginning of initialization, the number of phonons, velocities, positions, frequencies, and polarizations should be created by the MC method. The number of phonons with angular frequency at the temperature of can be calculated from the Bose-Einstein distribution as [24, 25] where , denotes the Planck constant. represents the Boltzmann constant. For silicon or germanium materials, three acoustic branches including one longitudinal acoustic branch and two transverse acoustic branches in the phonon dispersions should be taken into account. The total phonon number in the th spectral interval with a given volume at temperature can be calculated from [24, 25] where is the number of transverse acoustic (TA) phonons and represents the number of longitudinal acoustic (LA) phonons. denotes the density of state of LA phonons and represents the density of state of TA phonons. is the number of spectral intervals and is the th spectral interval. In this work, is set to 1000. Phonon’s group velocity and position can be obtained as follows [24, 25]: where represents the average frequency of the th spectral interval, and denotes the phonon wave vector according to the th spectral interval (in one dimensional model, can be written as ). The temperature of each cell can be obtained by energies they have, and this energy can be calculated by [24, 25] where is the total energy of the material with a volume , represents phonon’s polarization, and the temperature of each cell can be got from the bisection method.

###### 2.1.2. Phonon Drift

After being initialized, phonons drift at the interval of time space . Based on the Matthiessen formula, the total lifetime of phonon can be calculated as

In this work, we set to 1.0 ps and stop the simulation while thermal flux converges to within 1%. Whether one phonon is involved in an event of interfacial scattering is decided by its coordinate position and transmission coefficient.

###### 2.1.3. Interfacial Phonon Scattering and Transmission Coefficient

The regime of interfacial phonon scattering is very complex and still not be fully understood. There are two major problems that should be solved before figuring out the physical picture. One is the computation of phonon interface transmission, and the other is to dispose interfacial phonon scattering. At present, the main analytical methods are AMM and diffuse mismatch model (DMM), and the subsequent models are mostly based on these two classical models.

For the classical models such as AMM and DMM, the contributions of interfacial quality are neglected to simplify the complexities. Basing on Chen’s formula [26], we introduce the roughness factor into equation (8) to make the model closer to the reality. represents the quality of the interface and is caused by the lattice mismatch, lattice imperfection, the synthetic process, etc. where is the phonon transmission coefficient from material 1 to material 2, and and represent thermal capacities of material 1 and material 2, respectively. and indicate phonon group velocities. From MC simulations, , , , and can be calculated from a phonon’s dynamic information.

Whether a phonon can pass through the interface or not is determined by two steps in this work. The first step is judged by equation (8), which is a simplified model used by previous studies. If equation (8) is not be satisfied, the phonon will be reflected directly. If the first step is satisfied, we should use the second step to further judge whether the phonon satisfies elastic or inelastic scattering conditions. If not, the phonon will be also reflected. That is, a phonon passing through the interface should not only satisfy equation (8) but also meet the conditions of the second step. From equation (8), we can judge preliminarily whether one phonon is able to pass through or not, but the processes of interface phonon scattering are still difficult to dispose. Phonon transport between two different materials is assumed as a movement between their phonon dispersion curves. That is, one phonon with a given frequency on the dispersion curves of material 1 will pass through the interface elastically, if the corresponding point with the same frequency can be found on the dispersion curves of material 2. Otherwise, it will be reflected elastically or be scattered inelastically by the interface. To aid comparisons and explanation for interfacial scattering, we put the dispersions of Si and Ge in one figure, which are calculated separately. Figure 1 shows the mechanisms of transporting from Si to Ge material for a LA phonon labeled as p1 on the dispersion of Si. The TA1 curve of Ge is obtained by shifting TA curve vertically to the extreme location, and its lowest phonon frequency equals to the cut-off frequency of LA. We get the position of q1 on LA curve of Si, which has the same frequency with the highest frequency on the TA1 curve. That is, TA1 and q1 are fixed in this model. If the LA phonon p1 is judged to transmit the interface successfully by equation (8), it must be scattered into two TA phonons located at p2 and p3 separately since there does not exist one point having the same frequency with p1 on the curves of TA or LA of Ge. The positions of p2 and p3 are determined by shifting the TA curve of Ge and random selections by . Subtracting the frequency to , the value of can be obtained. The position of p can be selected randomly by , where is a random number ranging from 0 to 1. Then, drawing a horizontal line through p, p2 is obtained from the LA curve of Ge and q3 is also obtained from the -axis. Shifting the TA curve to TA2 with a dash curve, the lowest phonon frequency of TA2 equals to . Drawing a horizontal line through p1, p3 can be obtained from the shifted TA2 curve. p3 has the same wave vector with p3 on the TA curve of Ge. In this case, LA phonon p1 in Si is split into the LA phonon p2 and the TA phonon p3 in Ge, . It should be noted, the random number plays an important role in solving interfacial phonon scattering, which is generated uniformly by the procedure. With the increasing number of sample incident phonons and random number , the simulated system converges to an equilibrium state.