Abstract

Bamboo-like SiC nanowires (NWs) have specific geometric shapes, which have the potential to suppress thermal conductivity by phonon boundary scattering. In this work, phonon transport behaviors in the 3C-SiC, 4H-SiC, and 6H-SiC crystal lattices are studied by the Monte Carlo (MC) method, including impurity scattering, boundary scattering, and Umklapp scattering. Phonon relaxation times for Umklapp (U) scattering for the above three SiC polytypes are calculated from the respective phonon spectra, which have not been reported in the literature. Diffuse boundary scattering and thermal rectification with different aspect ratios are also studied at different temperatures. It is found that the thermal conductivities of the bamboo-like SiC polytypes can be lowered by two orders of magnitude compared with the bulk values by contributions from boundary scattering. Compared with bamboo-like 4H-SiC and 6H-SiC NWs, 3C-SiC has the largest U scattering relaxation rate and boundary scattering rate, which leads to its lowest thermal conductivities. The thermal conductivity in the positive direction is larger than that in the negative direction because of its lower boundary scattering relaxation rate.

1. Introduction

SiC is one of the most important third-generation semiconductors in the fields of electronic devices and high-temperature components. Particularly in harsh conditions, SiC devices are much better than Si and GaAs devices [1]. Though SiC has hundreds of stable polytypes, the most commonly used ones are 3C-SiC, 4H-SiC, and 6H-SiC [2]. In recent years, there have been many experimental and theoretical efforts invested in the development of SiC nanostructures [35]. Researchers [68] have proposed that thermal rectification effects are dependent on the size of the device and the asymmetric geometry. If phonons can be managed to transmit information as electrons and photons, thermal devices will also have more practical applications in phononic devices and the field of thermoelectric conversion.

To determine the boundary and interfacial effects of nanoscale structures, some groups have studied core-shell NWs [9], variable cross-section NWs [10, 11], and rough NWs [12, 13] in recent years. In thermal studies of SiC, Ni et al. [14] and Chantrenne and Termentzidis calculated thermal conductivities and frequency density of states of (DOS) 3C-SiC and 6H-SiC NWs with constant cross-sections by the nonequilibrium molecular dynamics (NEMD) method [15]. Bamboo-like SiC NWs were synthesized by a vapor-liquid-solid mechanism and chemical vapor deposition [1623]. Termentzidis et al. predicated the thermal conductivities of 3C-SiC and 2H-SiC NWs with constant and variable cross-sections via the NEMD method [24], in which the physical model was similar to that of bamboo-like NWs. Zianni appraised the impact of diameter modulation on the thermal conductivity reduction in core-shell silicon nanowires [25]. By a kinetic theory model, Zianni et al. [2628] investigated transport properties of phonon and electron in diameter-modulated nanowires and pointed out the amount of disorder which suppresses thermal conductance deeply according to phonon transmission coefficients.

Approaches to investigate nanoscale heat conduction generally fall into two main categories: analytic and numerical methods. Analytic methods incorporate the Boltzmann transport equation (BTE), Green’s theorem, and so on. The most widely used numerical methods are the NEMD and MC methods. In the MD method, macroscopic properties can be obtained by calculating the transmission function of every atom. However, if the system contains a large number of atoms, the computations will be immensely expensive. Furthermore, when the simulated temperature is lower than the Debye temperature, the BTE will be invalid for the exchange between classical statistics and quantum statistics, and the results should be revised by quantum corrections. In contrast to the MD method, the MC method treats atomic thermal vibrations as quantum particles and regards atomic thermal interactions as phonon-phonon scattering. Since the calculation of the transmission function is not required in MC simulations, the computational complexity is reduced to a moderate extent. The transmission parameters and distributions of phonons can be gained statistically without quantum corrections at low temperatures.

Based on the Debye model, Peterson first applied the MC method to solve thermal transport in a one-dimensional model by simplifying the phonon group velocities and polarization [29]. Most later MC methods were derived from Peterson’s work and extended its computational accuracy and geometric complexity [3037]. Using realistic phonon spectra, the accuracy of the MC method can be improved significantly [3843]. In particular, sophisticated nanostructures have attracted more interest for studying the effects of phonon ballistic transport [44, 45] and interfacial scattering [46]. Through MC simulations, researchers proposed that the surface roughness is the main factor leading to ultralow and amorphous-limit thermal conductivity [47, 48]. Bong and Wong studied the impact of anisotropy scattering using the MC method and suggested that forward scattering has a more obvious effect on the thermal conductivity than backward scattering [49]. Rickman et al. studied thermal transport across the grain boundary and suggested that one can engineer thermodynamic and transport properties in materials by inducing interfacial layering transitions via changes in temperature or pressure [50].

2. Monte Carlo Method of Phonon Transport

In the presented work, the entire MC simulation process is composed of four steps: initialization, drift, scattering, and statistics. For the initialization, the number of phonons with angular frequency at a given temperature can be calculated by the Bose-Einstein distribution:where denotes the Planck constant and represents the Boltzmann constant.where is the total number of phonons in volume . The phonon spectra are divided into 1000 equal-sized intervals. The first term on the right-hand side of (2) is the total number of transverse acoustic (TA) phonons, and the second term represents the number of longitudinal acoustic (LA) phonons. and denote the DOS of the LA phonons and TA phonons, respectively. is the th frequency interval. The phonon group velocity can be obtained by The temperatures of each cell can be counted by their energies, and this energy can be obtained bywhere is the material volumic energy and represents the phonon polarization.

After initialization, every phonon has a drifting time step . should be less than the average relaxation time for a phonon to avoid missing phonon scattering events. In this work, is set at 1.0 ps. Based on Matthiessen’s formula, the average relaxation time can be calculated by where is the phonon impurity scattering relaxation time [31]. denotes the normal scattering relaxation time, and is the Umklapp scattering relaxation time. is the boundary scattering relaxation time. Given the value of , the average scattering probability can be calculated from Whether a scattering event occurs should be judged by comparing the value of and a generated random number. In this work, is not included in (6). The conventional approach to solve phonon boundary scattering is composed of an indirect method and a direct method. Introduced in Ziman’s work [51] and developed by scholars [52, 53], the indirect method has a simplified , which assumes that the lateral surfaces of the simulated structure are flat. In the direct method, boundary scattering only depends on whether its position exceeds the boundary. In our study, the direct method is adopted. If a boundary scattering event occurs, the phonon will undergo diffuse reflection. The thermal flux and thermal gradient can be gathered by the synchronous ensemble method at each time in Figure 1. The vertical solid lines denote the time point . The interval time between two adjacent dashed lines is time step . The distance between two adjacent solid lines is , where is set to 102–103 in the simulations. To simplify the explanations in Figure 1, is set to 4. The transverse lines represent the time course of phonon 1 and phonon 2. Solid circles denote impurity scattering or U scattering, and hollow circles denote boundary scattering.

At each time point , all sample phonons should be relabeled, and their velocities and polarization should be inherited. The energy variation of each cell needs be counted to calculate the thermal flux. The energy-temperature list may be calculated from (4) in advance. Using the bisection method, the temperature of each cell can be determined by the energy-temperature list. Given the thermal flux and temperature gradient, thermal conductivities can be worked out by Fourier’s law.

3. Models for Phonon Transport in Bamboo-Like NWs

As shown in Figure 2, 3C-SiC has a zinc blende crystal structure with an ABC stacking sequence, while 4H-SiC and 6H-SiC have hexagonal crystal structures with ABCB and ABCACB stacking sequences. In Figure 3, hollow circles in the bamboo-like schematic drawing denote phonons distributed in the bamboo-like NWs, and the right-top inset is a scanning electron microscopy (SEM) image [21]. The thermal flux flows along the positive direction, indicated by the arrow symbol, from the hot part to the cold part. If we swap the positions of the hot part and cold part, the thermal flux will flow along the negative direction. The right-bottom inset is a bamboo joint schematic drawing confined by three parameters: , and . represents the bamboo joint’s top width, is its height, and denotes its bottom width. is defined as the aspect ratio to modulate the shape of the bamboo joint. In simulations, given , we can change the values of and to investigate the phenomena of phonon diffuse boundary scattering and thermal rectification effects.

In this work, three kinds of phonon scattering are considered: impurity scattering, U scattering, and boundary scattering. To consider impurity scattering, we use the same equation as [26] to calculate its relaxation rates. Since U scattering obeys the conservation of energy and momentum, it can be neglected in the simulations. The combined phonon relaxation time for impurity scattering and U scattering can be given by Matthiessen’s formula, as in (5), and the possibility of boundary scattering is only dependent on its locations in the simulated physical model.

3.1. Phonon U Scattering Rate for 3C-SiC, 4H-SiC, and 6H-SiC

From a theoretical point of view, thermal vibrations are often quantized as phonons, obeying the Boltzmann distribution. The phonon energy can be obtained from the spectra of 3C-SiC, 4H-SiC, and 6H-SiC, shown in Figure 4. Since optical phonons have much lower velocities than acoustic ones, we neglect the optical branches, following the previous literature [49, 54], in our simulations. Compared with 4H-SiC and 6H-SiC, 3C-SiC has a higher acoustic cut-off frequency, as shown in Figure 4. of 3C-SiC is along the direction, while 4H-SiC and 6H-SiC have the same expressions along the direction, where is the length of the axial vector [55].

When a phonon undergoes a U scattering event with two other phonons, their energies and wave vectors will be changed, as shown in the first Brillouin zone in Figure 5. For example, T1 is one of the degenerate acoustic transverse branches of 3C-SiC, and T2 is the translational curve of T1. is the acoustic longitudinal branch, and is the mirroring curve of at the symmetry axis of . P2 and P3 are intersection points of curve T2 and . Point P2 belongs to curve T2, while P3 belongs to T3. P1 and P2 represent two acoustic phonons with and , while P3 represents a longitudinal phonon with . They satisfy the U scattering requirement and , where denotes the reciprocal lattice vector. There are several types of U scattering events: , , and . The expressions of the phonon U scattering relaxation rate for the former two U scattering events are given aswhere is the Grüneisen parameter, for 3C-SiC [56], and for 6H-SiC and 4H-SiC [57]. is Planck’s constant, is the Boltzmann constant, and is the density of SiC. is the local temperature. is the phase velocity, and is the group velocity. is the optimum radius of the integral surface in the process of U scattering [58]. For the event , (8) can be employed to calculate :Figure 6 shows that the relaxation rate for is much larger than those for and , since the number of transverse acoustic phonons participating in U scattering is much larger than that of longitudinal acoustic phonons, according to the Bose-Einstein distribution. Figure 7 shows the U scattering combined relaxation rates for 3C-SiC, 4H-SiC, and 6H-SiC. Below 200 K, the U scattering relaxation rates decrease sharply with increasing temperature. Above 400 K, the U scattering relaxation rates for the three SiC polytypes converge to the orders of 10−11–10−12 gradually. 3C-SiC has the largest U scattering relaxation rate, while 4H-SiC has a closer value to that of 6H-SiC above 100 K. From this, it is can be speculated that 6H-SiC has a larger thermal conductivity than 3C-SiC or 4H-SiC, which is in accordance with experiments [59].

3.2. Phonon Boundary Scattering in Bamboo-Like NWs

While transporting through bamboo-like NWs, a phonon may encounter a large number of boundary scattering events. The border of the bamboo-like NWs should be formulized first. If a phonon’s coordinates surpass the border, it will be reflected diffusely by the boundary because of lattice defects and surface stresses. The direction of its velocity will be changed along the border randomly following equation:where and represent the components of the velocity along the -axis and -axis. and are generated random numbers. is a coordinate transformation matrix decided by the equation of the boundary border. A sample phonon’s flying trajectories are recorded as hollow circles in Figure 8 for a period of 100 . When the sample phonon knocks at the boundary, it will be reflected diffusely.

4. Results and Discussion

To verify our MC model and determine the value of the impurity scattering parameter Bi, the bulk thermal conductivities of the three SiC polytypes were calculated and are given in Figure 9. The circle symbols, triangle symbols, and rectangle symbols represent the experimental values [60]. The curves were obtained by our MC simulations of impurity scattering and U scattering. Below 100 K, the thermal conductivities fluctuate drastically at low temperatures and reach peak values at 80–100 K. Above 600 K, the thermal conductivities of the three polytypes converge to 0.5–2 W/cmK. In this case, the impurity scattering parameter Bi is fit to 5 × 10−48 from the experiments, which is fixed in our simulations. We employ Fourier’s law to calculate the thermal conductivities by the average cross-sectional area of the bamboo junction. As shown in Figure 9, the discrepancies between the experiments and the simulations are considerable. This gap may result from two sources: one is from the approximate value of Bi, which is uncertain in the experiments [60], while the other is perhaps from the Grüneisen parameters in (7) and (8), which are also approximate values in [56, 58]. Though there are considerable gaps between the simulations and experiments, the overall trends match well with experiments. To some extent, this MC model is helpful to investigate thermal transport in nanoscale structures qualitatively.

4.1. Density of States (DOS) of Phonon Frequency and Velocity in 3C-SiC

From a phonon’s DOS, distributions of the phonon frequency and velocity can be obtained directly, which is beneficial to analyze U scattering and boundary scattering in bamboo-like NWs. Figures 10 and 11 show the DOS of the phonon frequency and phonon velocity of 3C-SiC at 200 K and 600 K. It can be seen in Figure 10 that, with increasing temperature, the ratio of phonons with high frequencies increases. Figure 11 indicates the variation in the phonon group velocity. The average phonon group velocity equals 4276.3 m/s at 600 K and equals 4441.5 m/s at 200 K, which manifests that a higher temperature leads to a lower phonon group velocity. We also calculated the DOS of the phonon frequency and phonon velocity of 4H-SiC and 6H-SiC. The results show that the three SiC polytypes have similar trends of phonon frequency and group velocity.

4.2. Phonon Transport in Bamboo-Like 3C-SiC NWs

Boundary effects play an important role in nanoscale thermal transport. It is worth researching the thermal conductivities and thermal rectifications of bamboo-like NWs quantitatively. As Figure 12 shows, the thermal conductivities of bamboo-like NWs with different sizes are simulated along the positive direction. is fixed at 100 nm, and the thickness of the simulated NWs is set to 50 nm. The thermal conductivities of 3C-SiC with different aspect ratios are presented over temperatures from 100 K to 500 K. Since the number of phonons with high velocities decreases as the temperature rises, the boundary scattering relaxation rate also decreases with increasing temperatures. Hence, in contrast to bulk conditions, the thermal conductivities at high temperatures are larger than those at low temperatures in bamboo-like NWs. Additionally, the thermal conductivity becomes steady and irrelevant at aspect ratios above 2.3. In Figure 13, we investigate the dependence of the thermal conductivity on , with values of 100 nm, 150 nm, and 200 nm. When increasing from 150 nm to 200 nm, the gap between the thermal conductivities is smaller than that when going from 100 nm to 150 nm. To further study boundary scattering, we trace 5000 sample phonons and calculate their average boundary relaxation time, which is shown in the inset of Figure 13. The phonon boundary scattering relaxation time for  nm is very close to that of  nm and apparently larger than that of  nm, which affirms exactly the effect of boundary scattering.

Moreover, the thermal conductivities in the negative direction are investigated in Figures 14 and 15. As shown in Figure 14, the thermal conductivities increase with temperature and reach a steady state when the aspect ratio exceeds 2.4. The thermal conductivities in the negative direction are significantly lower than those in the positive direction. In Figure 15, the bamboo-like NW with  nm has the lowest thermal conductivities, as seen in Figure 13. The gap between NWs with  nm and  nm is larger than that in the positive direction. The inset figures in Figures 13 and 15 show that the boundary scattering relaxation time in the negative direction is shorter than that in the positive direction, which gives rise to a higher possibility for boundary scattering.

Thermal conductivity trends to maintain steady. The insets also indicate the discrepancies of boundary relaxation time between different aspect ratios decrease with the increasing value of aspect ratios. That is, the contribution of phonon boundary scattering to thermal resistances is restricted by the aspect ratio. Moreover, thermal conductivities of 4H-SiC and 6H-SiC in positive direction and negative direction are simulated, and the calculated results are used to thermal rectification in the next section.

4.3. Thermal Rectification of Bamboo-Like 3C-SiC, 4H-SiC, and 6H-SiC NWs

The thermal rectification effect is similar to the regime of the electronic diode, which exists in bamboo-like SiC NWs because of the different boundary scattering rates in the positive and negative directions. The efficiency of thermal rectification can be obtained from the following equation [61, 62]:where and represent the thermal conductivities in the positive and negative directions, respectively. Figure 16 shows the thermal conductivities and rectification efficiencies in the negative and positive directions at 300 K. With increasing aspect ratio, the thermal conductivities in the positive and negative models for 3C-SiC decrease first and then converge to 0.3 W/mK. increases from zero with increasing aspect ratio below 1.1 and then decreases to zero with a further increase in the aspect ratio. The trends are caused by the dependency of phonon boundary scattering on the aspect ratio. As approaches , the aspect ratio approaches zero, which would represent a constant cross-section nanowire. In this case, the rectification efficiency should also approach zero, and the resulting thermal conductivities in the forward and reverse directions can be compared to constant cross-section values. In addition, as the aspect ratio approaches infinity, the bamboo features also diminish, which also may result in near-zero rectification efficiencies.

In Figure 17, we also calculate , , and for 4H-SiC. Compared with 3C-SiC, 4H-SiC has higher and , which is consistent with the deduction from Figure 7. To study the contributions of boundary scattering, the phonon boundary relaxation time for 3C-SiC, 4H-SiC, and 6H-SiC was calculated and is given in Figure 18. The phonon boundary scattering relaxation time (inversely proportional to the boundary scattering relaxation rate) for 3C-SiC is shorter than that for 4H-SiC, resulting in a lower thermal conductivity. Namely, thermal resistances arising from boundary scattering in the bamboo-like NWs are more predominant for 3C-SiC than for 4H-SiC or 6H-SiC. The reason mainly relies on the differences between their phonon spectra, which leads to distinctive phonon group velocities and polarization.

In comparison to 3C-SiC and 4H-SiC, 6H-SiC has the highest and , as shown in Figure 19, due to the lowest U scattering rate and lowest boundary scattering rate. It can also be seen that the three studied SiC polytypes have different peak values of , and the highest is less than 0.5. Thus, by means of modulating the shape of the bamboo-like NWs, the thermal conductivities can be dropped by 200–500 times, compared with the bulk values. However, when the aspect ratio increases or decreases away from the value corresponding to the peak , the thermal rectification efficiency becomes ambiguous, and the trend of is shaped as a parabola going downwards.

In the previous works [9, 10, 24, 61], most of the authors have similar predictions that thermal conductivity decreases with the increasing aspect ratio. Different from their opinions, we suggest that thermal conductivities in positive direction and negative position do not increase with the aspect ratio at all times. When the aspect ratio increases above 2.4, thermal conductivity trends to maintain steady. The insets of Figures 13 and 15 indicate the discrepancy of boundary relaxation time between different aspect ratios decrease with the increasing value of aspect ratios. That is, the contributions of phonon boundary scattering to thermal resistances are restricted by the aspect ratio.

5. Conclusion

We have presented MC investigations of bamboo-like SiC NWs for 3C-SiC, 4H-SiC, and 6H-SiC on phonon transport. Particularly, phonon U scattering relaxation rates and boundary scattering relaxation time are calculated separately in this work, providing qualitative theoretical supports to future MC studies. It is observed the special shape of the bamboo-like boundary plays an essential role in reducing the thermal conductivities as previous reports [2528]. Particularly at low temperatures, the thermal conductivities are suppressed distinctly. Thermal conductivities in the positive direction are higher than those in the negative direction, which is because the boundary scattering relaxation rate in the positive direction is lower than that in the negative direction. Increasing the aspect ratio from 0.1 to 3.5, the rectification efficiency increases from zero at first and then decreases to zero with further increase in the aspect ratio, which demonstrates that the geometric shape of the bamboo-like junction has limited effect on tailoring the thermal rectification efficiency.

Conflicts of Interest

The author declares that they have no conflicts of interest.

Acknowledgments

Z. Wang would like to acknowledge the financial support from the National Natural Science Foundation of China (51106043 and 51373048) and the basic research project of Henan Provincial Department of Education Science and Technology (13A470176 and 2014GGJS-059).