Tuning the thermal conductivity of silicon carbide by twin boundary: a molecular dynamics study

Silicon carbide (SiC) is a semiconductor with excellent mechanical and physical properties. We study the thermal transport in SiC by using non-equilibrium molecular dynamics simulations. The work is focused on the effects of twin boundaries and temperature on the thermal conductivity of 3C-SiC. We find that compared to perfect SiC, twinned SiC has a markedly reduced thermal conductivity when the twin boundary spacing is less than 100 nm. The Si–Si twin boundary is more effective to phonon scattering than the C–C twin boundary. We also find that the phonon scattering effect of twin boundary decreases with increasing temperature. Our findings provide insights into the thermal management of SiC-based electronic devices and thermoelectric applications.


Introduction
Silicon carbide (SiC) is an excellent candidate material for high-power high-frequency applications in nano/microelectro-mechanical systems (N/MEMS) due to its wide electronic band gap, highly saturated electron drift velocity, excellent thermal and chemical stability [1,2], high creep and corrosion resistances [3][4][5][6][7]. However, the intrinsically high thermal conductivity of crystalline SiC leads to the reduction of figure of merit and unexpected thermoelectric performance, and hence limits its wider potential to be an efficient thermoelectric energy conversion material. To improve thermoelectric properties, preexisting defects or grain boundaries are often introduced into bulk material as phonon scattering barriers for tuning of thermal transport [8][9][10]. Among the practical thermal barriers, twin boundary (TB) is one of the coherent grain boundaries, which could enhance the figure of merit without a degradation in electrical conductivity by forming quantum wells [11,12].
The TBs in SiC have been synthesized in real applications. For instance, thin twin lamella were observed in twinned single crystals of 3C-SiC grown by the Baikov technique [13]. Stacking faults or TBs are usually produced parallel to the close-packed {1 1 1} plane [14]. It was proposed that the natural heterojunctions and superlattices with multilayer twins in SiC can facilitate phonon scattering [15]. Recently, many researches have been performed to study the phonon scattering at TB. Some atomistic simulations found that thermal resistance of TB is lower than other high-energy grain boundaries [10,[16][17][18]. The weak phonon scattering effect of TB have also been reported in twinned diamond [19]. On the other hand, recent experiments found marked reduction of thermal conductivity in Ge twinned NWs [20]. The efficient phonon blocking at TBs was also calculated in SiC nano-twinned superlattices NWs by Xiong et al [21] using MD simulations, where they ascribed the efficient phonon blocking to the combined effects of phonon scattering by TBs and surfaces. However, the isolated phonon scattering mechanism of TB in bulk SiC is still lacking.
In this paper, we perform non-equilibrium MD (NEMD) simulations to study the effects of two different types of TBs (i.e. Si-Si TB and C-C TB) on heat transport in bulk twinned SiC. We investigate the influences of temperature and sample length on the thermal conductivities of these twinned samples. The phonon scattering mechanism of TBs is also analyzed for understanding the reduced thermal conductivities of twinned 3C-SiC. Our results have implications for phonon engineering in other polytypes of SiC, i.e. the stacking fault band structures in 4H-SiC [19,22], for potential thermal management and thermoelectric applications.

Simulation methodology
As shown in figure 1, the model system of bulk 3C-SiC used in this work is a cubic supercell with a {1 1 1} TB in the middle. Figure 1(a) displays the atomic structure of the simulated supercell with [1 1 0], [1 1 2], and [1 1 1] orientations aligned with the x, y, and z directions, respectively. Figure 1(b) shows that there are six regions in the supercell: two fixed walls at outmost ends, two thermal transfer regions (being cold and hot respectively) in the middle, and two thermal bath regions. Such end-fixed model is physically equivalent to the double-length periodic model with periodic boundaries along all directions, as shown in figure 1(c). This periodic model can be considered as a simplified infinite system containing periodic TB, which is often called periodic supper-lattice. To reduce the computational cost, we use the 'end-fixed' model throughout this work, i.e. two outmost ends are fixed in the z direction and periodic boundary conditions applied in both x and y directions; the cross-section in each model is fixed at 4 × 6 unit cells along x and y directions, respectively. Figure 2 shows that the TB in the middle of the supercell has two types: one is the Si-Si bonded mirror interface (referred to as Si-Si TB) and the other is C-C bonded mirror interface (referred to as C-C TB).
MD simulations are conducted with LAMMPS [23]. The thermal conductivities of bulk SiC are calculated from NEMD simulation results. It is a 'direct approach' as experimental measurement and can facilitate the converge of our model systems, considering the fact that the heat current autocorrelation function (HCACF) is difficult to converge for systems consisting of structural defects such as TBs, grain boundaries, and disordered layers. Although the direct method may underestimate the thermal conductivity of our models due to their limited size [24,25], it is considered as an effective atomistic method to study the thermal transport properties of twinned systems.
To study the temperature effect on thermal transport of SiC, the simulation systems are first relaxed at different temper atures ranging from 100 K to 1500 K. The lattice parameter at each temperature is obtained by calculating the dimension of the equilibrated supercell, which is fully relaxed with an isothermal-and-isobaric (NPT) ensemble. Such NPT equilibration simulation is performed for 1 ns with the time step of 1 fs. Periodic boundary conditions are applied to the whole supercell during equilibration.
After equilibration, NEMD simulation is conducted in a constant volume and constant energy (NVE) ensemble. As experiments, we fix two outmost walls and prescribe a high temperature (i.e. T 0 + 50 K) in the hot bath region and a low temperature (i.e. T 0 − 50 K) in the cold bath region, respectively. After this temperature difference is imposed, NEMD is run for 0.5 ns with a time step 0.5 fs, so as to equilibrate heat transport from the hot to cold region. The thermal transfer region is divided into 160 slabs in z direction. The spatial average of temperature in each slab is calculated based on the kinetic energy of atoms therein. Then NEMD is run for another 1 ns to obtain the time average of temperature gradient along the z direction. For all the systems studied, the simulation time of 1.5 ns is long enough to reach the steady heat flux [26], such that the energy addition needed to keep high temperature in the hot region is approximately equal to the energy subtracted from the cold bath at each time step.
From the temperature gradient and resulting heat flux, the thermal conductivity, k, is calculated based on the Fourier Law: where J is the heat flux density normal to the cross-section, ∆T is the temperature difference between cold bath and hot bath regions, l is the heat transfer length; G is the temperature gradient, which can be derived by using the least-square method along the heat transfer direction. During NEMD simulations, periodic boundary conditions are applied in x and y directions while both ends in z direction are fixed. The atomic interaction in 3C-SiC is described by using an empirical Tersoff potential [27]. This many-body potential was originally developed for covalent crystalline solids and later adapted for different polymorphs of SiC. The total potential energy Φ in an N-atom system is given by where r ij is the interatomic distance between atom i and j; are the repulsive and attractive interaction functions, respectively; * b ij represents the many-body effect from atom k around atom i and j. We take the Tersoff function parameters for SiC from Porter et al [28]. These parameters have been validated by studying the predicted thermal conductivity and thermal expansion coefficients of 3C-SiC [29].

Figures 3(a) and (b)
show the typical temperature profile in the model system of Si-Si TB and C-C TB, respectively. Both profiles exhibit a sharp temperature gap at the TB. This gap at the Si-Si TB is larger than that at the C-C TB. For comparison, the temperature profile in a pristine SiC, as shown in figure 3(c), is continuous along the length direction. We normalize these temperature profiles by dividing their respective average temperature in the heat-transfer region and plot them in figure 3(d). Clearly, the normalized temperature gradient of pristine SiC is larger than that of twinned samples.
To study the effect of TB on thermal transport, we calculated temperature gaps across TB (∆T twin ), partial hot region (∆T hot ) and cold regions (∆T cold ). In this work, ∆T twin was defined by the difference between the interface temperatures, which were extrapolated from the temperature profiles in the hot and cold regions to the interface [30]. In all the cases studied, temperature profiles of separated thermal transfer regions are similar. ∆T hot is equal to ∆T cold and the length of hot region l hot is equal to the length of cold region l cold . The effective thermal conductivity k e can be expressed as where J e is the effective heat flux across the whole sample, = + l l l hot cold is the entire heat transfer length and the effective temperature gradient for the whole sample is twin . According to the widely used thermal circuit model, the total thermal resistance of the twinned sample can be defined as: The total resistance of pristine sample is  where J p is the effective heat flux across the pristine sample. In these equations, the exact measurement of the isolated interfacial temperature gap could be difficult due to the nonlinear temperature jump at TB. To consider this nonlinearity, we ascribe the interfacial (or Kapitza) resistance of TB as the difference between R e and R p , which can be derived as The above equations suggest that there exists a threshold thermal transfer length beyond which the temperature gap due to TB can be negligible. To find the threshold length, we calculate the effective heat flux (J e ), the effective thermal conductivity (k e ) and Kapitza resistance (R TB ) for samples with different lengths. Figure 4(a) presents the calculated heat flux as a function of sample length. Effective heat fluxes of all samples decrease with increasing thermal transfer length. It is also observed that heat fluxes of twinned samples are significantly less than those of pristine SiC, which confirms a strong phonon blocking mechanism at TB [21]. Hence, we conclude that besides sample length, heat flux is another important factor that controls the effective thermal conductivity. Under the same thermal potential between hot and cold bath regions, effective heat fluxes of samples with Si-Si TB are slightly less than those with C-C TB, indicating that phonon blocking at the Si-Si TB is slightly larger than that at the C-C TB. Figure 4(b) shows the effective thermal conductivity k e of samples with Si-Si TB and C-C TB and of a pristine sample as a function of sample length at 300 K. In all cases, effective thermal conductivities increase with sample length. For pristine samples, thermal conductivities saturate at length beyond 100 nm. This size dependence of thermal conductivity mainly results from two sources. One is boundary scattering between thermal bath and thermal transfer regions [31]. This scatting effect is pronounced for the sample length less than the bulk phonon mean free path (MFP) [24,25] and becomes less pronounced for longer cases. Another is phonon-phonon scattering during phonon transportation, which is proportional of the sample length. While for twinned systems, there is an additional TB phonon scattering, which accounts for the reduction of thermal conductivity k e from k p . According to figure 4(b), the effective thermal conductivities of both twinned samples show a similar trend as that of the pristine system. Both k e of twinned samples are markedly lower than the corresponding k p of pristine sample. Notably, the k e of Si-Si twinned samples are slightly less than that of C-C twinned ones. This is consistent with the typical temperature profile of the Si-Si twinned sample which has larger temperature gap at TB than C-C twinned samples. For example, in samples with the same length (l = 34.9 nm), the effective thermal conductivity of Si-Si twinned sample is 35.7 W mK −1 , which is ~8.2% less than C-C twinned sample and ~35.9% less than pristine sample. We note that, as the sample length increases to beyond 100 nm, the difference of k e between twinned and pristine samples decreases rapidly. This suggests a critical twinning-sensitive sample length at 100 nm below which thermal transport is sensitive to the presence of TB.
We further investigated size dependence of the interfacial thermal resistance of both twinned samples at different sample lengths, as shown in figure 4(c). The interfacial resistances presented here were calculated by equation (6) from the predicted k e of twinned samples and k p of pristine samples. It is known that total thermal resistances, including bath-transfer boundary resistance [31], conduction resistance and interfacial resistance, calculated by NEMD method are dependent on system size. Therefore, the size dependence of isolated Kapitza resistance needs to be calculated by keeping other resistances consistently. Here, we keep the domain size of corresponding twinned and pristine samples at the same length for each calculation of interfacial resistance. It is noted that the interfacial resistance of TB decrease with sample length increasing and can be negligible with sample length beyond 180 nm, indicating that phonon scattering at TB becomes insignificant as sample length increases to the critical sample length. Since the critical length of the endfixed model can be taken as the threshold spacing between TBs in the periodic model, we suggest a maximum twinningsensitive interval spacing (~180 nm) below which phonon engineering is effective by introducing TBs.
Next, we study temperature effects on the thermal conductivity of twinned SiC at a constant sample length of 21 nm. Figure 5 summarizes the effective thermal conductivities of both Si-Si twinned and C-C twinned samples as a function of temperature. For comparison, thermal conductivities of pristine SiC at different temperatures are also presented. All types of samples show a similar trend of deceasing thermal conductivities with increasing temperatures from 100 K to 1500 K, which is consistent with previous results in other solid materials [30,[32][33][34][35]. The reduction of thermal conductivities occurs primarily in the first stage from 100 K to 700 K. In this stage, large differences of thermal conductivities can be observed between pristine and twinned samples. However, as the temperature rises over 700 K, thermal conductivities of all samples begin to saturate and show little difference.
The temperature dependence of thermal conductivity can be explained by phonon kinetic theory, which describes the thermal conductivity as where C v is the specific heat at constant volume, v g is the phonon group velocity, Λ denotes the phonon mean free path and can be expressed as τ Λ = v g , and τ is the phonon relaxation time. For twinned samples, the total phonon relaxation time can be deduced from the Matthiessen rule: where τ u represents the contribution from phonon-phonon anharmonic (or umklapp) scattering, τ twin from phonon scattering at TBs, and τ b from bath-transfer boundaries. It is known that the phonon-phonon anharmonic scattering of 3C-SiC can be negligible at low temperatures [36]. We have demonstrated in our previous results that the boundary effect (τ b ) depends on sample length. Thus, at a short sample length or low temperatures, the phonon blocking at TB becomes dominant over the phonon-phonon scattering. As a result, the significant thermal resistance from TBs will lead to a large difference between twinned and pristine samples. As temperature increases, the vibrational energies of atoms in the system increase, leading to the enhanced phonon umklapp scattering along the sample. The enhancement of conduction resistance also occurs with increasing sample length. Thus, at a long sample or high temperatures, the conduction resistance will dominant over the interfacial resistance, leading to small difference between pristine and twinned samples.
We also calculated the partial thermal conductivity of the cold region and hot region of Si-Si twinned and C-C twinned samples at different temperatures, as shown in figure 5(b) and (c) respectively. For both twinned samples, the partial thermal conductivities of cold region (k cold ) and hot region (k hot ) are approximately equivalent to each other, indicating that phonon transport remains coherent after crossing TBs. It is also noted that the difference between partial and effective thermal conductivities of each twinned sample is only significant at low temperatures, confirming that phonon blocking of TB is effective at low temperatures.
To gain further insight into the effect of TB on the thermal conductivity, we compared phonon spectra of localized Si-Si twinned and C-C twinned samples (only including interface and its nearby 15 layer atoms) with that of pristine sample at the same sample size and orientation. The phonon spectra are computed by evaluating eigenvalues of system dynamic matrix (directly extracted from MD simulations) based on fluctuation-dissipation theory [37]. This method can capture the anharmonicity of phonons at finite temperature naturally [36]. The typical phonon density of states (PDOS) for each type of sample is shown in figure 6(a). It is clear that all samples show a typical two-peak PDOS curve, where the first low-frequency peak is primarily determined by Si atoms and the second by C atoms. Notably, the primary frequency band of pristine samples is broader than that of twinned samples. This phenomenon can also be observed from the comparison of phonon dispersion curves between the three samples, as shown in figure 6(b), where optical branches of pristine sample are higher than that of twinned samples. It seems that pristine sample could transfer wider phonon bands than twinned samples. This phonon spectral difference could explain why the thermal conductivity of pristine SiC is much higher than that of twinned samples. Figure 6(c) shows the phonon spectra of localized Si-Si twinned samples at different temperatures. All phonon spectra have two peaks at about 13 THz and 23 THz. With increasing temperature, the first peak shows 'blue shift' and the second peak shows 'red shift', leading to a narrower frequency band at elevated temperatures. As expected, the narrower band of PDOS results in the lower effective thermal conductivity at higher temperature. The similar temperature dependence of PDOS can also be found in C-C twinned sample. Figure 6(d) shows PDOS of Si-Si twinned samples as a function of sample length at 300 K. It is seen that longer samples have a broader frequency band with the first peak fixed at about 13 THz and the second peak showing 'blue shift'. The wider phonon band suggests an increased thermal conductance for longer sample, which is consistent with length dependence of thermal conductivity. On the other hand, we note that both PDOS peaks decrease with increasing sample length, which suggests a decreased thermal conductance with increasing sample length. According to figure 4, longer samples have larger thermal conductivities, implying that phonon blocking of TB are more significant than phonon-phonon scattering in cases at 300 K. This confirms that the interfacial resistance dominant over the conduction resistance at low temperatures.

Conclusions
We have studied the thermal transport properties of twinned 3C-SiC using non-equilibrium molecular dynamics simulations with the Tersoff potential. The thermal conductivities of twinned 3C-SiC samples are found to be much lower than that of pristine 3C-SiC. Particularly, the effective thermal conductivity of Si-Si twinned sample is less than that of the C-C twinned sample. Furthermore, the phonon blocking mechanism of twin boundary on the thermal conductivity depends sensitively on sample length. A threshold sample length is found to be around 100-180 nm, below which twin boundary can reduce the thermal conductivity of bulk SiC significantly. However, the phonon blocking effect of twin boundary becomes weak as temperature increases over 700 K. Nevertheless, these findings provide new insights into the phonon blocking mechanism for tuning thermal conductivity of 3C-SiC by twin boundaries.

Acknowledgments
The Central South University (Grant No. Kfkt2013-10) are appreciated. The authors gratefully acknowledge Prof Ting Zhu for his enlightening comments and revisions.