Approaching the alloy limit of thermal conductivity in single-crystalline Si-based thermoelectric nanocomposites: A molecular dynamics investigation

Single-crystalline Si-based nanocomposites have become promising candidates for thermoelectric applications due to their prominent merits. Reducing the thermal conductivity κ without deteriorating the electrical properties is the key to improve their performance. Through non-equilibrium molecular dynamics simulations, we show that κ of single-crystalline Si-based nanocomposites can be reduced to the alloy limit by embedding various nanoinclusions of similar lattice constants but different lattice orientations or space symmetries with respect to the matrix. The surprisingly low κ is mainly due to the large acoustic phonon density of states mismatch caused by the destruction of lattice periodicity at the interfaces between the nanoinclusions and matrix, which leads to the substantial reduction of phonon group velocity and relaxation time, as well as the enhancement of phonon localization. The resulting κ is also temperature-insensitive due to the dominance of boundary scattering. The increase in thermal resistance induced by lattice structure mismatch mainly comes from the nanoinclusions and the channels between them and is caused by the enhanced boundary scattering at the interfaces parallel to the heat flux. Approaching the alloy limit of κ with potentially improved electrical properties by fillers will remarkably improve ZT of single-crystalline Si-based nanocomposites and extend their application.

I n recent years, bulk semiconductor nanocomposites, i.e., bulk semiconductors with nanoinclusions, have become promising candidates for large-scale thermoelectric applications due to their prominent performance, low cost and structural stability [1][2][3][4][5] . The efficiency of thermoelectric materials is often measured by a dimensionless figure of merit ZT 5 S 2 sT/(k L 1 k e ), where S is the Seebeck coefficient, s is the electrical conductivity, T is the temperature, and k L and k e are the lattice and electronic contributions to the thermal conductivity k, respectively. With the nano-size inclusions of a characteristic dimension comparable to the phonon mean free path (MFP) or wavelength, nanocomposites can strongly suppress the phonon transport with less deterioration in electrical properties through classical and quantum confinement effects, therefore leading to a much enhanced ZT. Remarkable ZT values or improvements over the bulk materials have been reported for AgPb m SbTe 21m systems 1 (ZT 5 2.2 at 400 K), In 0.53 Ga 0.47 As/ErAs nanocomposites 2 (200% improvement), p-type BiSbTe nanocomposites 3 (ZT 5 1.4 at 373 K) and PbTe-SrTe systems 5 (ZT 5 ,2.2 at 915 K), which are comparable with those of low-dimensional systems [6][7][8][9][10][11] . More importantly, the nano-bulk combination provides an opportunity to engineer some materials that are unsuitable for thermoelectric applications in their bulk form, such as singlecrystalline Si (ZT < 0.01 at 300 K due to its high k L 5 148 w/mK 12 ), as low-cost high-performance thermoelectrics. For single-crystalline Si-based thermoelectric materials, a crystalline matrix is preferred to achieve high charge mobility. However, k of a single crystal is generally much higher than that of its alloy, which is widely believed to be the lowest attainable k in crystalline solids (''alloy limit'') 2 . Nanopores have been recently introduced to reduce k of Si, leading to significantly enhanced ZT. Tang et al. 13 reported a ZT 5 0.4 at 300 K for nanopatterned holy Si films by decreasing k to ,2 W/mK. Yu et al. 14 and Hopkins et al. 15 also found that k of Si nanomeshes could be as low as 1.9 W/mK and 6.8 W/mK respectively without compromise on electrical properties at room temperature. Very recently, Yang et al. 16,17 predicted a ZT as high as 0.76 at 300 K for n-type three-dimensional Si phononic crystals, in which k can be reduced by a factor up to 500 with little deterioration in the electrical properties compared with those of bulk Si.
With a similar structure, Si-based nanocomposites have rarely been experimentally investigated for thermoelectric applications. The open question is whether k of Si-based nanocomposites can be further reduced so as to approach the alloy limit.
The key to answer this lies in a fundamental understanding of thermal transport in nanocomposites. A few simulation works have been conducted to investigate k of nanocomposites using Boltzmann transport equation (BTE) 18 , Monte Carlo (MC) 19 and molecular dynamics (MD) simulations 20 . These investigations showed that k depends on both the volumetric fraction and size of the nanoinclusions and proposed to reduce k using a high interface density. Some theoretical models have also been developed for the quantitative prediction of the effects of nanoinclusions on k. Kim et al. 21 proposed an approximate analytical solution to estimate the phonon scattering cross section of polydispersed spherical nanoparticles. Mingo et al. 22 calculated k of the ''nanoparticle-in-alloy'' nanocomposites and predicted a potential 5-fold increase in ZT at room temperature. All these investigations demonstrate the great potential of decreasing k of Si-based nanocomposites, which is generally attributed to the enhanced interfacial scattering. However, further efforts are still needed to explore the microscopic origin of the strong interfacial scattering and the strategy to approach the alloy limit in Si-based nanocomposites. It is of importance to illuminate how the various dissimilarities between the nanoinclusions and matrix as well as the design parameters influence k of Si-based nanocomposites and the phonon transport properties such as group velocity and relaxation time. These phonon transport details are highly desired for understanding the k reduction mechanism and are helpful for the design of nanocomposites.
In this work, we investigate k of single-crystalline Si-based nanocomposites using nonequilibrium molecular dynamics (NEMD) simulations and conduct a systematic phonon transport analysis in both frequency domain and real space domain. Various nanoinclusions of similar lattice constants but different mass densities, lattice orientations or space symmetries from those of the matrix are considered. k of samples with various nanoinclusion sizes and pitch spacings are calculated. The influence of different nanoinclusions on k of the nanocomposites is revealed. We find that lattice orientation or space symmetry mismatches can even reduce k of nanocomposites to the alloy limit by depressing the group velocity and relaxation time, as well as enhancing phonon localization. To further the understanding of k reduction in these nanocomposites, the local real-space thermal transport and resistance distribution in the Ge/Si nanocomposites are illuminated and the correlation between the enhanced scattering and the surface vibration states is discussed.

Results
The representative configurations of Si-based nanocomposites are shown in Fig. 1(a). The cylindrical nanoinclusions with a square cross section d P 3 d P are embedded in the matrix and separated by the spacing d S . The nanocomposites are constructed by embedding nanoinclusions into bulk Si of [100] orientation (parallel to the Z axis) at the designed positions. Eight types of nanoinclusions are considered, including Ge100 (Ge nanocrystal with orientation [100] parallel to the Z axis), Ge110, Ge111, Si110, Si111, Si 0.5 Ge 0.5 alloy (Fd3 _ m group), NiSi 2 nanocrystal (Fm3 _ m group) and pores. Ge, SiGe and NiSi 2 are good candidates because their lattice constants are very close to that of Si and the electron transport is expected to be little affected.
We calculated k using NEMD with a typical simulation domain shown in Fig. 1(b) and the modified embedded atom method (MEAM) potentials 23 , which have been successfully applied to describe a wide range of materials, in particular semiconductors such as Si, Ge and their alloys [24][25][26][27] . It is noted that for Si and Ge the MEAM potentials can well reproduce the acoustic phonon transport properties but overestimate the frequencies of optical phonons, which can be addressed by 2NN MEAM potentials 28,29 . Considering that optical phonons contribute little (,5% for bulk Si) to the thermal conductivity 30,31 , the MEAM potentials can still lead to reasonable accuracy for these materials. We first verified the accuracy of these potentials by calculating the bulk k of various samples at 300 K. The NEMD simulations give 145 6 20 W/mK for Si, 71 6 11 W/mK for Ge, 4.9 6 0.6 W/mK for Si 0.5 Ge 0.5 and 7.3 6 0.7 W/mK for NiSi 2 , respectively. These values are in good agreement with the experimental results (148 W/mK for Si 12 , 59.9 W/mK for Ge 12 , 4.6 W/mK for Si 0.5 Ge 0.5 32 and ,10 W/mK for NiSi 2 33 ), indicating the suitability of these potentials for the targeted investigations.
Figures 2(a) and (b) show k of all samples with infinite length in the direction perpendicular to nanoinclusions at 300 K, as a function of d S and d P , respectively, in comparison with that of bulk Si 0.5 Ge 0.5 . For nanoporous structures, to evaluate their applicability as thermoelectric materials, k eff based on the effective cross section should be used due to the existence of pores, which is related to the porosity w and commonly used k A based on the total cross section though A similar treatment has also been employed by Hopkins et al. 15 . It is shown that k of all nanocomposites are significantly reduced with respect to that of bulk Si and the reduction depends on both the configuration and the type of nanoinclusions. Figure 2(a) shows that k of the samples increase almost linearly with the increasing d S , indicating that important phonon modes are in the ballistic transport regime within this d S range. This is also supported by the investigations of spectral contributions to k 30,31,35,36 , which found that phonon modes with MFPs from tens of nanometers to a few micrometers contribute significantly to k of bulk Si while those with MFPs , 10 nm have negligible contributions. Compared with other samples, k of Ge100/ Si and SiGe/Si grow much faster with the increasing d S . Figure 2(b) shows that k of Ge100/Si and SiGe/Si decrease significantly with the increase of d P while those of other samples almost remain constant. The comparison between Figs. 2(a) and 2(b) indicates that k is more sensitive to d S . For the sample Ge100/Si, k decreases by 33.5 W/mK with an increase of the volumetric fraction of Ge by 18.8% when d S decreases; however, by varying d P k only decreases by 16.4 W/mK although the volumetric fraction of Ge increases by 31.3%. The later analysis will show that this is due to the strong scattering at the interfaces parallel to the heat flux. Because the intrinsic MFPs of important phonon modes in the matrix are much longer than those of nanoinclusions, the conductance of channel regions will be more strongly affected by d S compared with the influence of d p on that of nanoinclusions and thus leads to more significant variation of k in the nanocomposites. For fixed configurations, nanoporous Si exhibits the lowest k for most cases; however, k of some nanocomposites can be as low as that of nanoporous Si. Considering that some nanoinclusions can even improve electrical properties, Si-based nanocomposites might achieve a better thermoelectric performance than nanoporous Si. Furthermore, appropriate patterning can even make k of Si-based nanocomposites approach that of bulk Si 0.5 Ge 0.5 alloy, which probably results in an optimal ZT. Also, the type of nanoinclusions significantly influences k. The thermal conductivity of Ge100/Si is larger than that of SiGe/Si, which is reasonable when considering the lower k of SiGe. Unexpectedly, most other samples show significantly lower k compared with SiGe/Si even though the inclusions have a higher bulk k. The effect of nanoinclusions on k is generally attributed to their size and density differences between the nanoinlcusions and matrix 22 . However, it seems that our results cannot be simply explained by these factors. For instance, for the same configuration with d S of 1.1 nm and d P of 1.1 nm, k (6.9 W/mK) of Ge110/ Si is less than one quarter of that (30.9 W/mK) of Ge100/Si although k of Ge is isotropic. Actually, the embedding of nanoinclusions also introduces lattice constant mismatch and interatomic potential mismatch except for density difference. In MD simulations, we can individually investigate the effect of each factor although these factors usually couple together. It is found that the lattice constant mismatch is typically , 1% for the relaxed samples. To evaluate the effect of lattice constant mismatch, we modified the potential of Ge to eliminate this mismatch while keeping other properties the same. The simulations show that k of Ge100/Si with the modified potential is very close to that of the original one (k 5 30.9 W/mK), indicating a negligible influence of the small lattice constant mismatch. This agrees with the previous observation that k changes little when the strain caused by lattice constant variation is less than 1% in bulk Si 37 . We then investigated the effects of density difference and potential mismatch on k by modifying the atomic mass and potential of the nanoinclusions in the Ge100/Si. In Sample A, the potential for Ge is replaced by that of Si while in Sample B the atomic mass of Ge is changed to that of Si. It is found that k of the former (31.2 W/mK) changes little while the latter exhibits a k value close to that of bulk Si. We further changed the atomic mass of nanoinclusions and found that k decreases with enhanced density difference, agreeing with previous investigations 22 . Based on the above analysis, one can conclude that density difference is the main reason that causes the reduced k in Ge100/Si nanocomposites while interatomic potential and small lattice constant mismatches considered here have little effect.
Although density difference accounts for the reduction of k in Ge100/Si, the mechanism that leads to the further reduction of k in Ge110/Si and Ge111/Si is still unclear. Considering the isotropic thermal properties in Ge, there should be little difference for k between them. Therefore, the difference in k is probably due to the lattice orientation mismatch. We note that the samples can be divided into two groups based on the lattice structures of nanoinclusions: lattice-structure-matched nanocomposites, including Ge100/ Si and SiGe/Si, and lattice-structure-mismatched ones, i.e., the lattice orientation and/or group symmetry of the nanoinclusions are different from those of the matrix. The significant difference between k of these two groups indicates that these mismatches may be crucial for the further reduction of k.
Due to the mismatches in lattice orientation and/or group symmetry, the vibrational properties of interfacial atoms may be affected. Figure 3(a) shows the PDOS of interfacial Si and Ge atoms in Ge100/ Si, SiGe/Si, Ge110/Si, Si110/Si and NiSi 2 /Si, with respect to those of bulk Si and Ge. For all the samples, there are significant PDOS mismatches between the interfacial atoms, accounting for the interfacial scattering and thus the reduced k. For example, the interfacial Si (Ge) atoms in Ge100/Si and SiGe/Si vibrate at higher frequencies than those in the bulk, which has also been reported in the study of other solid-solid interfaces 29,38,39 . These higher frequency phonons only exist close to the interface and cannot participate in the energy transport across interface without phonon scatterings, because they are significantly mismatched with phonon modes far from the interface. The new observations mainly lie in the lower frequency regime. For Ge100/Si and SiGe/Si, the acoustic PDOS of interfacial Si (Ge) atoms below 5 THz almost overlap with those of the bulk although their optical PDOS are significantly suppressed. However, in Ge110/Si, Si110/Si and NiSi 2 /Si, the acoustic PDOS of the interfacial Si and Ge atoms significantly shift to the lower frequency and the optical peaks of PDOS are further reduced with respect to those in Ge100/Si and SiGe/Si. Such a shift is caused by the destruction of lattice periodicity near the interfaces. For the interfaces in Ge100/Si and SiGe/Si, only element disorder arises and lattice periodicity is preserved, which mainly affects optical phonons but has much less influence on acoustic phonons. However, for the interfaces in Ge110/Si, Si110/ Si and NiSi 2 /Si, there is also lattice disorder due to the mismatches in lattice orientation or space symmetry except for the element disorder, which therefore remarkably disturbs acoustic phonons and has stronger influence on optical phonons. Because acoustic phonons dominate the thermal transport in Si and Ge, the acoustic PDOS mismatch can strongly suppress the acoustic phonon transport and thus reduces k to the alloy limit. The PDOS mismatch in Si110/Si is as strong as that in Ge110/Si, indicating that the lattice disorder almost dominates the PDOS mismatch. If there is a large difference between the lattice arrangements of interfacial atoms, atomic reconstruction usually happens 38 and modifies the lattice periodicity of interfacial atoms. As shown in Fig. 3(b), there is a significant atomic reconstruction for the interfacial atoms in Ge110/Si while another two samples almost remain their initial lattice positions. In Ge100/Si and NiSi 2 /Si, the lattice arrangements of nanoinclusions are very similar to that of the matrix, which makes the interfacial atoms only deviate slightly from their initial positions. However, the difference in interfaces between Ge110/Si and NiSi 2 /Si indicates that acoustic PDOS mismatch may also arise even when there is no significant atomic reconstruction. This is because the vibrational properties of an atom are related to its interactions with surrounding atoms and atomic reconstruction is just one possible way to change the coupling.
To clarify the effects of lattice orientation mismatch on the phonon transport in Si-based nanocomposites, we calculated the phonon group velocities V g and relaxation times t of Ge100/Si, Ge110/Si, and Si110/Si, in comparison with those of bulk Si, as shown in Figs. 4(a) and (b). All the samples used in the phonon spectrum analysis have the same d S of 1.1 nm and d P of 1.1 nm. The group velocities were obtained by numerically differentiating the phonon dispersions extracted from EMD simulations 40 . The relaxation times were obtained by fitting the autocorrelation function of the total energy of each phonon mode with an exponential function 34,41 . Obviously, group velocities of Ge100/Si, especially for optical modes, are greatly reduced with respect to those in bulk Si, which can be attributed to the lower group velocities in nanoinclusions and the zone folding effects, as reported in nanoporous structures 34,42 . When the frequency approaches 0, group velocities in Ge100/Si are comparable to those in bulk Si but they rapidly decrease with the increasing frequency and remain relatively low values at higher frequencies.
The reduction of relaxation times in these nanocomposites is also significant over most frequency range. More importantly, the group velocities and relaxation times of most phonons in Ge110/Si and Si110/Si are further reduced with respect to those of Ge100/Si. This can be ascribed to the enhanced boundary scattering caused by the acoustic PDOS mismatch and enhanced optical PDOS mismatch. One can find much shorter relaxation times in Ge110/Si from 1 to 5 THz, agreeing well with the frequency range where acoustic  PDOS mismatch happens. Due to the same configuration and the dominance of boundary scattering, Ge110/Si and Si110/Si exhibit similar group velocities and relaxation times.
It has been reported that interfaces/surfaces may result in phonon localizations in Si nanotubes 43 and Ge-Si core-shell nanowires 44 . To understand the localization effect, the participation ratios of phonon modes in these samples were calculated, as shown in Fig. 4(c). The participation ratios of most modes in the nanocomposites are significantly reduced with respect to those in bulk Si, indicating that the introduction of nanoinclusions makes phonons more localized. Compared with Ge100/Si, the participation ratios of many modes in Ge110/Si are significantly lower, indicating that lattice-orientation-mismatch leads to the enhancement of localization effect. The localization generally arises at surfaces/interfaces and can be enhanced by disorder 43 . As mentioned above, latticeorientation-mismatch induces lattice disorder at the interfaces, which therefore enhances the localization. Similar behavior has been observed in nanoporous Si 45 .
Therefore, one can conclude that the further reduction of k in Ge110/Si and Si110/Si compared to Ge100/Si is mainly caused by the lattice orientation mismatch, which results in the significant reduction of group velocities and relaxation times, as well as the enhancement of phonon localizations. Similar phenomena are also found for Ge111/Si and Si111/Si. Compared with Si110/Si and Si111/Si, Ge110/Si and Ge111/Si show slightly lower k, which is likely due to the additional scatterings caused by element mismatch. Obviously, boundary scattering dominates the thermal transport in these nanocomposites, making their k values similar for the same configurations. The dominance of the boundary scattering also explains why k of latticestructure-mismatched nanocomposites are almost independent on d P .
Similarly, the acoustic PDOS mismatch caused by different space symmetries also greatly suppresses the phonon transport in nanocomposites. As shown in Fig. 2, NiSi 2 /Si shows a substantially lower k than that of SiGe/Si whose nanoinclusions have the same space symmetry as the matrix, although k of NiSi 2 is even higher than that of SiGe. The acoustic PDOS mismatch also significantly depresses group velocities and relaxation times in NiSi 2 /Si, as shown in Figs. 5(a) and (b). Over most frequency range, especially between 1 and 5 THz, NiSi 2 /Si nanocomposite exhibits significantly lower group velocities and shorter relaxation times with respect to those in SiGe/Si. The enhancement of phonon localization in NiSi 2 /Si is obvious, as marked by the significantly lower participation ratios of most modes compared with those in SiGe/Si. Thus, space-symmetrymismatch accounts for these variations and the further reduction of k in NiSi 2 /Si.
The above analysis shows that the lattice orientation and space symmetry mismatches can lead to significant acoustic PDOS mismatch near interfaces and thus further reduce k of nanocomposites to the alloy limit by depressing group velocities and relaxation times, as well as enhancing phonon localizations. For thermoelectric applications, metal silicides may be preferable fillers because of their better electrical properties. Several investigations suggest that nanoinclusions may enhance the power factor with respect to that of the matrix, e.g., the power factor can be improved in InGaAs embedded with ErAs that can act as electron donor 2,22 . Such a low k combined with even improved electrical properties may make Si-based nanocomposites a promising candidate for thermoelectric applications. In addition, the lattice structure mismatches can significantly suppress low frequency phonons that are generally difficult to block, providing a new strategy to suppress acoustic phonon transport.
Next, we investigated the effects of lattice structure mismatch on the temperature dependence of k. Figure 6 shows k of Ge100/Si and Ge110/Si with d S of 2.2 nm and d P of 1.1 nm between 300 K and 1100 K. Significantly different temperature dependences are observed. k of Ge110/Si is insensitive to temperature while that of Ge100/Si notably decreases as the temperature rises. It is found that the MFPs of most modes in Ge110/Si are less than or comparable with the interatomic distance and those phonons with MFPs significantly larger than the interatomic distance only contribute , 6% to k. The dominant modes with short MFPs mainly limited by the strong boundary scattering are hardly affected by temperature. Therefore, the application temperature range for these nanocomposites will be significantly broadened. Similar temperature dependence of k was also observed in nanoporous SiGe alloy 46 .
The spectrum analysis provides an overall frequency-domain explanation of the effect of lattice orientation and space symmetry mismatches on the thermal transport. Understanding the temperature distribution and thermal resistance related to the configurations will provide more details about local thermal transport in the real space and help with the design and optimization of nanocomposites. We therefore calculated the temperature distributions and profiles for the Ge100/Si and Ge110/Si with d S 5 5.4 nm and d P 5  Figure 7(a) shows the temperature distribution of Ge100/Si for the heat flux along the Z direction, indicating an obvious nonuniformity. Figure 7(b) further shows the temperature variations along the normalized coordinates X* 5 0.125, 0.250 and 0.500 in Ge100/Si and Ge110/Si. The temperature distributions and profiles show that a temperature jump arises at the interfacial regions. For the line away from the nanoinclusion (X* 5 0.125), temperature decreases monotonically. But for the lines passing through the interface (X* 5 0.250) and nanoinclusion center (X* 5 0.500), abnormal temperature gradients are observed. The temperature slope for the line X* 5 0.500 is negative in the Border regime, showing a clear hot region in front of the nanoinclusions and a cold region behind them, which are caused by the phonon energy accumulation and depletion, respectively. This phenomenon probably results from phonon ballistic transport through channels, as explained by Yang et al. 18 . Similar behavior has also been observed in nanoporous structures 34 . The overall temperature profiles along the heat flow direction are shown in Fig. 7(c). It is clear that the total temperature drop for Ge110/Si is significantly larger than that in Ge100/Si, corresponding to a lower k in the former. The temperature slope of Ge110/Si in the intermediate region is more than two times larger than that of Ge100/Si, indicating increased thermal resistance due to the lattice orientation mismatch.
Because the heat flux distribution is also nonuniform in nanocomposites, thermal resistances of the thermal boundary phase (TBP), nanoinclusions between TBP and channel regions were then calculated to quantitatively evaluate the influence of lattice orientation mismatch, as shown in Fig. 7 the matrix side and half UC on the nanoinclusion side because the PDOS of atoms within half UC near the interfaces are significantly modified. From the values of temperature difference DT, heat flow J and cross section area A c , the thermal resistance R can be determined by R 5 DTA c /J. The thermal resistance of the Border region is a summation of those in the left and right Borders. The TBP thermal resistance R TBP includes the contributions from the two TBPs. Clearly, lattice orientation mismatch strongly increases the thermal resistance of the nanoinclusion and channel regions while only slightly enhances R TBP . In the nanoinclusion and channel regions, the enhanced boundary scattering due to the PDOS mismatch significantly decreases the MFP and thus results in increased thermal resistances. Similarly, it has been reported that the strong boundary scattering determines the thermal resistance of channel regions that dominate the thermal transport in nanoporous Si 34 . In both samples, R TBP is important while R TBP of Ge110/Si is only 23% larger than that of Ge100/Si. The variation of R TBP in both samples can be rationalized by the diffuse mismatch model (DMM) 47 , which provides relatively good agreements with experiments at room or higher temperatures 38,48 . When the PDOS near the interface are used, DMM predicts a 20% larger R TBP for Ge110/Si in comparison with that in Ge100/Si, agreeing well with the NEMD simulation results. Similarly, the thermal resistances of different regions in NiSi 2 /Si are enhanced due to the space symmetry mismatch.

Discussion
In summary, we have calculated k of single-crystalline Si-based nanocomposites embedded with various nanoinclusions using NEMD simulations and conducted a phonon transport analysis in both frequency domain and real space domain. Specifically, we designed nanoinclusions with a lattice orientation or space symmetry different from that of the matrix. It was found that k of nanocomposites with lattice-structure-mismatched nanoinclusions are significantly lower than those of the lattice-structure-matched ones and can even approach the alloy limit. A comparison between the PDOS of interfacial atoms and bulk shows a significant acoustic PDOS mismatch and enhanced optical PDOS suppression due to the destruction of lattice periodicity in nanocomposites, which therefore introduce strong scatterings for acoustic and optical phonons respectively. Further phonon spectrum analysis shows that the nanocomposites with lattice-structure-mismatched nanoinclusions exhibit notably lower group velocities, shorter relaxation times and stronger phonon localizations. Specifically, the low frequency phonons are greatly suppressed due to the destruction of lattice periodicity, providing a new strategy to block long-MFP phonons.
Additionally, the lattice structure mismatch leads to a temperature-insensitive k in Si-based nanocomposites from 300 K to 1100 K, which can be ascribed to the very short MFPs of dominant phonon modes caused by strong boundary scattering. The local thermal transport details were further illuminated by exploring the temperature distributions and profiles in the nanocomposites, showing that with the introduction of lattice structure mismatch, the thermal resistance of interfacial regions only slightly increases while those of nanoinclusions and channels between them are remarkably improved. This is due to the enhanced boundary scattering at the interfaces parallel to the heat flux.
Approaching the alloy limit of k combined with even improved electrical properties will significantly improve the ZT of singlecrystalline Si-based nanocomposites, which thus may become a promising bulk thermoelectric candidate. Especially, the insensitive temperature dependence of k in these nanocomposites will therefore extend their applications to a wide temperature range.

Method
NEMD simulation. The NEMD simulations were performed in LAMMPS 49 using the direct method 50 , with the periodic boundary condition applied to the X, Y directions and the fixed boundary condition to the Z direction. The MEAM potential parameters for Si and Ge are from Refs. 26 and 27, respectively. The cross potential parameters for Si-Ge and Ni-Si are from Refs. 24 and 25, respectively. The nanoinclusions are embedded in the intermediate region of the simulation domain, which is separated with the heat source/sink by a buffer layer to eliminate the nonlinear effect. Initially, the system was equilibrated at a given temperature within the NPT (constant number of particles N, pressure P and temperature T) ensemble with a Nose-Hoover thermostat for 300 ps and then in the NVE (constant number of particles N, volume V and energy E) ensemble for 100 ps. The two ends of the simulation domain, with a thickness of one unit cell (UC), were in contact with Langevin thermostats to form the heat source and sink. The temperature of the heat source and sink was set as T H 5 T 0 1 D/2 and T C 5 T 0 2 D/2, respectively, where T 0 is the mean temperature of the system and D is the temperature difference. After the system reached a steady state, a linear temperature profile was achieved in the intermediate region for 3-6 ns. The time-averaged temperature profile in the intermediate region was linearly fitted to obtain the temperature gradient, which was used to calculate the thermal conductivity according to Fourier's law where J is the heat flow and A c is the cross section area. To eliminate the finite-size effects, for each configuration, a series of simulations with variable domain lengths (L Z ) ranging from 65.2 nm to 273.7 nm were performed. The thermal conductivity in bulk limit k was obtained by fitting to 1/k A -1/L Z relation under the gray approximation and extrapolating to L Z R ' 50 , where a good linear correlation was found.
Participation ratio. Phonon localization can be quantitatively characterized by the participation ratio p l defined for each eigenmode l as 51 where N is the total number of atoms and e ia,l is the ath eigenvector component of eigenmode l for the ith atom. The participation ratio measures the fraction of atoms participating in a given mode and effectively indicates the localized modes with O(1/ N) and delocalized modes with O(1).
Diffuse mismatch model. According to the diffuse mismatch model, the phonon boundary resistance R b between two materials 1 and 2 is given by where v m is the cutoff frequency of the softer material, f p~½ exp ( hv=k B T){1 {1 is the Bose-Einstein (equilibrium) distribution, D p is the PDOS of interfacial atoms and v g is the phonon group velocity. t 1 R 2 is the transmission coefficient from the material 1 to 2 and can be calculated based on the phonon mode l by t 1?2 (v)~P l v g,2,l D p,2,l (v) P l v g,1,l D p,1,l (v)z P l v g,2,l D p,2,l (v) :