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 [210], such as AlAs/GaAs, Si/Ge, Si/SiGe, and Si1−xGex/Si1−yGey, 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 [1114].

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 [1618] 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 107 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.

For a TA phonon, the same method can be used to solve the behavior of inelastic scattering as shown in Figure 2. Phonons with positions below q2 can pass through the Si/Ge interface without splitting. While a phonon in Si material is located at the region between q1 and q2, it will be reflected by the interface, unless it splits into a TA phonon p2 and a LA phonon p3 in Ge.

2.1.4. Renewal

At the end of one circulation, every sample phonon should be labeled again before the next circulation. The velocity and polarization of the sample phonon will be inherited. At this stage, we should count the energy variations of each cell to calculate thermal flux. In simulations, the energy list of temperatures is calculated from equation (6) in advance. Then, the temperature of each cell can be decided by its energy from this list by the bisection method. Given the thermal flux and temperature gradient, we can calculate thermal conductivities of the simulated superlattices from Fourier’s law.

3. Results and Discussions

Firstly, we investigate one period Ge/Si superlattice with cross section and a total length of 400 nm. This superlattice is discretized into eight equal layers with 50 nm thickness. The first layer serving as a hot part is set to 140 K, and the last layer serving as a cold part is set to 60 K. Furthermore, each layer is divided into five segments to calculate temperatures. The first 20 segments are set as Si material and the last 20 segments are set as Ge material. In this model, temperatures distribute linearly along the segment number. The results indicate the temperature difference and in the case of and . According to equation (8), the larger the , the more phonons are reflected by Ge/Si interface, causing the higher interfacial temperature jump.

If we swap the positions of the hot part and the cold part, heat flux will flow from Si to Ge material so-called Si/Ge. In Figure 3, while and , and . Since the phonon cut-off frequency of Ge is smaller than that of Si material in Figure 1, the longitudinal phonon from Ge can pass through the Ge/Si interface freely without inelastic splitting. Thus, and in Figure 3 are higher than those in Figure 4, and the Ge/Si superlattice has a lower interfacial thermal resistance than the Si/Ge superlattice.

Increasing the number of periods to five, the temperature distributions in the Ge/Si superlattice are shown in Figure 5. The total temperature jump between the hot part and the cold part equals to 61 K.

Reversing the direction of thermal flux as shown in Figure 6, we can get the result that the Si/Ge superlattice has a higher than the Si/Ge superlattice.

Further, to investigate the low temperature conditions, we set the temperature of the hot part to 100 K and set the cold part from 20 K to 50 K. As shown in Figure 7, with the increasing temperature gap , the average along thermal flux increases too. But, the variation trends of in the Si/Ge superlattice are different from those of Ge/Si. As to the Si/Ge superlattice in Figure 7, Si materials locate at the left sides of the first, the third, and the fifth interfaces, and Ge materials locate at the right sides. Therefore, phonons fly from Si to Ge which will undergo more interfacial scattering than flying from Ge to Si. on the first, the third, and the fifth interfaces in Si/Ge is higher than that on the second and the fourth interfaces. Conversely, of the Ge/Si superlattice has the opposite changes. Because of the stronger interfacial thermal resistances, the thermal fluctuations of the Si/Ge superlattice are much higher than those of the Ge/Si superlattice.

In Figure 8, thermal conductivities of Ge/Si superlattices with periods from one to five are calculated, which increase with the increasing temperatures from 30 K to 180 K and become flat at higher temperatures from 180 K to 200 K. It can be explained that the number of high-frequency phonons with low group velocities arises with temperatures, leading to the rise of interfacial transmissions by equation (8) and reducing the affections of interfacial thermal resistances. The Ge/Si superlattice with one period has the highest thermal conductivities in the simulations. Particularly, Ge/Si superlattices with three, four, and five periods have almost equal thermal conductivities at temperatures above 180 K. The results indicate that thermal conductivities can be suppressed effectively by means of increasing the number of interfaces at temperatures from 60 K to 180 K, but it will be invalid above 180 K. By tracing 105 sample phonons, we obtain the static results as shown in the inset of Figure 8, which depicts the convergence of thermal flux. It can be seen that the thermal flux of the Ge/Si superlattice with one period converges more rapidly and precisely than that of superlattices with two periods or three periods.

Thermal conductivities of Si/Ge superlattices are evaluated in Figure 9, which are lower than those of Ge/Si superlattices as a result of strong interfacial thermal resistances. Referring to the dispersions of Si and Ge, the interfacial phonon transmission of Ge/Si is larger than that of Si/Ge. Since the cut-off phonon frequencies of the LA branch and the TA branch in Si are much higher than those in Ge, phonons with high frequencies in Si have little probabilities of passing through the Si/Ge interfaces by inelastic scattering and most of them are reflected, which increases interfacial thermal resistances. For the stronger interfacial thermal resistance, the thermal conductivities of the Si/Ge superlattice are higher than those of the Ge/Si superlattice. Different from the Si/Ge superlattice, the curves of thermal conductivities increase slowly below 180 K and then converge at higher temperatures. Referring to the results, we may decrease thermal conductivities in the Si/Ge superlattice at the temperatures from 80 K to 180 K by increasing the number of interfaces. Tian et al. [27] predicted thermal conductivities of Si/Ge superlattices by means of the MD method and density functional theory (DFT). Compared with Tian’s work in Figure 10, our results are closer to the data calculated from the DFT method than those from the MD method. The trends of thermal conductivities in our work are similar with Tian’s work, and the gaps of thermal conductivity are mainly due to the different geometric sizes of the simulated models.

Furthermore, as shown in Figure 10, the efficiency of thermal rectifications is also calculated by where and represent thermal conductivities of Ge/Si superlattices and Si/Ge superlattices, respectively. Figure 10 shows the average values of converge to 0.38~0.61 at the temperature ranges from 80 K to 200 K.

4. Conclusion

In this work, a dedicated MC model based on phonon dispersion is proposed to investigate the thermal conductivity and thermal rectification dependencies of interfacial thermal resistances at low temperatures. From this model, the reason that why phonons transmit from Si to Ge material sustain more interfacial resistances than flying from Ge to Si can be explained. According to the results, thermal conductivities of Si/Ge and Ge/Si superlattices can be modulated by increasing or decreasing the number of their interfaces at temperatures below 180 K. But, this solution is no longer valid at higher temperatures. Besides, the thermal rectification between the Si/Ge and Ge/Si superlattices is found to vary from 0.38 to 0.61, which can be used to control phonon transmissions in composite materials made up of Si and Ge materials, such as superlattice films and superlattice nanowires at low temperatures.

Data Availability

All data included in this study are available upon request by contact with the corresponding author.

Conflicts of Interest

The authors declare that there is no conflict of interest regarding the publication of this paper.

Acknowledgments

Z. Wang would like to acknowledge the financial support by the subsidy scheme for the Fundamental Research Funds for the Henan Provincial Colleges and Universities in Henan University of Technology (2017QNJH26) and the basic research project of Henan Province Department of Education Science and Technology (19A460012). This study was partly supported by the National Basic Research Program of China (51775170) and the Young Backbone Teachers in Colleges and Universities in Henan (2014GGJS-059).