Molecular Dynamics Simulations on Influence of Defect on Thermal Conductivity of Silicon Nanowires

A non-equilibrium molecular dynamics simulation method is conducted to study the thermal conductivity (TC) of silicon nanowires (SiNWs) with different types of defects. The impacts of defect position, porosity, temperature, and length on the TC of SiNWs are analyzed. The numerical results indicate that SiNWs with surface defects have higher TC than SiNWs with inner defects, the TC of SiNWs gradually decreases with the increase of porosity and temperature, and the impact of temperature on the TC of SiNWs with defects is weaker than the impact on the TC of SiNWs with no defects. The TC of SiNWs increases as their length increases. SiNWs with no defects have the highest corresponding frequency of low-frequency peaks of phonon density of states; however, when SiNWs have inner defects, the lowest frequency is observed. Under the same porosity, the average phonon participation of SiNWs with surface defects is higher than that of SiNWs with inner defects.


INTRODUCTION
A one-dimensional material system is the smallest size structure that can be used for efficient electronic and optical transmission (Xia et al., 2003;Hsu et al., 2015) and has unique heat transfer properties . It is considered to be an ideal choice for semiconductor materials. Silicon nanowire (SiNW) is a typical one-dimensional material, and there are many mature methods to obtain this material, such as laser ablation, thermal evaporation, chemical vapor deposition, template method, solution method, and hydrothermal method (Yang et al., 2014). It has not only stable semiconductor characteristics of macroscopic silicon material but also field emission (Frederick et al., 1999), thermal conductivity (TC) (Zhang and Zhang, 2013), visible photoluminescence (Feng et al., 2000), and quantum confinement effects (Kobayashi and Hiramoto, 2008) that are different from bulk silicon materials. Therefore, SiNWs can be widely applied in microelectronic device, nanophotonic devices, optoelectronic devices, and new energy (Boukai et al., 2008;Laik et al., 2008;Chen et al., 2018;Anufriev and Nomura, 2019;Zhang et al., 2019).
The TC is an important thermoelectric performance parameter of the material, and the research on the TC of SiNWs plays an important role in its application, especially in the field of microelectronics. Therefore, many investigations have been conducted to explore the heat transfer performance of SiNWs. The experimental measurement of the TC of a single SiNW at low temperature indicates that it is much lower than the bulk value  and shows linear temperature dependence alternative to the conventional T3 dependence (Bourgeois et al., 2007;Chen et al., 2008). By adjusting the structures, the TC of SiNWs can be tuned. For instance, the fishbone SiNWs allow more flexibility in controlling the heat conduction as compared with pristine SiNWs (Maire et al., 2018). Molecular dynamics simulation is a feasible and reasonable way to learn the mechanism of TC of SiNWs and the impact of some influencing factors that cannot be considered by the experiment. The numerical results obtained by Volz and Chen (1999) indicate that the TC of SiNWs is about two orders of magnitude lower than that of silicon bulk materials, and this trend has also been confirmed in experiments. The atomistic computation conducted by Markussen et al. (2008) shows that ultrathin pristine SiNWs with diameters ranging from 1 to 5 nm at room temperature have highly anisotropic heat conductance. Chen et al. (2011) found that the TC of SiNWs decreases with its surface-to-volume ratio (SVR) and is linearly related to SVR for all cross-sectional geometry when the crosssectional area is > 20 nm 2 . Mingo et al. (2003) compared three methods for predicting TC of SiNWs vs. temperature and found the method employing full phonon dispersion relations of the material is better than the traditional Callaway and Holland models. Through simulations, Donadio and Galli (2010) found that the presence of a crystalline/amorphous interface has a significant influence on the TC of thin SiNWs. Based on the non-equilibrium molecular dynamics method using Stillinger-Weber potential model and Nose-Hoover thermostat, Wang et al. (2009) predicted the TC of SiNWs. From the results, we can see that the reciprocal of the TC has a linear relationship with the length, and the TC dependence on temperature is consistent with the experimental results.
In addition, regulating the TC of SiNWs is of great significance to its application. Besides the abovementioned factors, such as temperature and size, introducing defects into SiNWs is also an important method. At room temperature, by doping isotopes, the TC of SiNWs can reduce exponentially . The TC of SiNWs with roughened surfaces is lower compared with that of the pristine ones, and its value is smaller under a smaller roughness wavelength and larger amplitude (Liu and Chen, 2010). For SiNWs, the TC can decrease close to two orders of magnitude due to the comprehensive influence of surface ripples and core defects (Zago et al., 2012). By raising the percentage of random vacancy defects, linearly decreased TC of SiNWs is observed in the results obtained from SW potential, and non-linearly decreased TC is observed in the results based on Tersoff interatomic potential. It indicates that vacancy defects can promote the efficiency of thermoelectric energy conversion of SiNWs (Shahraki and Zeinali, 2015). A simulation conducted by Zhan et al. (2014) indicated that the stacking faults can greatly reduce the TC of the SiNW, and the extrinsic stacking faults can induce a larger reduction than the twin boundary and intrinsic stacking faults. Yang et al. (2014) compared and analyzed the TC of silicon materials under different porosities, and the TC decreased with the increase of the porosity. Romano et al. (2012) computed the effective TC of porous SiNWs using the phonon Boltzmann transport equation and found the correlation between the phononic TC, the porosity, and surface-to-volume ratio.
In summary, the defects have a significant influence on the regulation of TC of SiNWs. However, research on the regulation mechanism of defects on the TC of SiNWs, especially the impact of spatial defects, is still limited. To this end, based on a nonequilibrium molecular dynamics simulation method, this paper conducts a numerical study that explores the defect location, porosity, and other factors impacting TC of SiNWs along with the impact mechanism.

Construction of Defective SiNWs
The SiNWs can be obtained by arranging the Si atoms regularly and repeatedly. As shown in Figure 1A, all SiNWs in the current study contain six Si spores and have a size of 3.258 nm in the x and y directions. To investigate the impact of porosity and temperature, SiNWs in the z direction have 200 Si spores and the length of 108.6 nm. In addition, SiNWs with lengths of 48.87, 70.59, and 152.04 nm are also constructed. Holes are dug into the inside or surface of the complete SiNWs to construct two types of defective SiNWs.
The SiNWs with surface defects and inner defects are shown in Figures 1B-E, respectively. In these figures, the Φ and d s represent the diameter of the sphere defect and the distance between two adjacent sphere defects. In the current study, d s is 5.43 nm, and the Φ and d s are all the same for SiNWs with inner and surface defects. For the SiNWs with inner defects, the sphere center coordinates of the sphere defect can be expressed as (1.629, 1.629, n × 5.43), where n is an integer from 1 to 20. Correspondingly, for the SiNWs with surface defects as shown in Figures 1C-E, center coordinates of the sphere defect in up and low surfaces, and the sphere center coordinates are (0, 1.629, n × 5.43) and (3.258, 1.629, n × 5.43), where n is an integer from 1 to 20 as well.
Porosity, ε, can be changed by adjusting the diameter of the sphere defect and is defined as where V defects is the volume of defects and V total is the volume of complete SiNWs. In this study, Φ of 4, 6, 8, 10, 12, 13, 14, and 15 Å are applied, and the corresponding porosities are 0.58, 1. 96, 4.64, 9.07, 15.68, 19.94, 24.91, and 30.64%, respectively.

Simulation Method
The non-equilibrium molecular dynamics simulation method of velocity-exchanging method (Müller-Plathe, 1997) is applied to determine the TC of various SiNWs in this paper. The SiNWs are divided into 54 layers with the single-layer thickness of about 2 nm in the z direction, and the number of atoms in each layer is about 1,000 to ensure the accuracy and authenticity of statistics.
Based on the velocity-exchanging method, the hot and cold areas are established, and the highest atom velocity in the cold area is exchanged with the lowest atom velocity in the hot area. During initialization, periodic boundary conditions are selected in the direction of heat flow (the z direction), and free boundary conditions are selected in the other two directions. The interaction between Si atoms is described by the Tersoff potential function (Tersoff, 1989;Volz and Chen, 1999). A simulation time step is set as 0.5 fs, and the initial system temperature is set as 500 K. During the simulation, the system will first run 2 million steps under the NVT system to reach steady-state equilibrium, and then the system will run 400 million steps under the NVE and perform velocity exchanging. When the system is stable, the average temperature, T, of each layer can be obtained by following the formula (2) according to the Boltzmann energy equalization theorem.
where N, m, and v are the atom number in the statistical area, the atomic mass, and the atomic velocity, respectively. k B is the Boltzmann constant. After determining the temperature of each layer, the TC of SiNWs can be calculated by where v h is the velocity in the hot area, v c is the velocity in cold area, t is the running time, S is the system cross-sectional area, and z is the coordinate position of the particle.
In the current study, the TC of SiNWs under each porosity and temperature is simulated several times by adjusting the random number in the speed initialization. Note that the thermal conductivities presented in the figures are the average values.

RESULTS AND DISCUSSIONS
The TC of Defective SiNWs Figure 2A shows the TC of two defective SiNWs with various porosities obtained by formula (3). As the porosity increases, the TC of both defective SiNWs with surface and inner defects decreases. When the porosity is up to 30%, the TC of defective SiNWs approaches 1.14 W/m•K and the impact of porosity tends to be flat. Under the same porosity and temperature of 500 K, the TC of SiNWs with surface defects is always higher than that of SiNWs with inner defects. In addition, by comparing these TC at 300 and 500 K in the Figure 2A, under the same porosity, TC at 300 K is always higher than that at 500 K.
On the macroscale, the TC is a physical parameter of the material which does not change with the size of the material. However, for the low-dimensional materials, the TC of the material is related to the length of the material due to the size effect. For this reason, the TC of three types of SiNWs with different lengths is analyzed. As shown in Figure 2B, the TC of SiNWs with no defects increases significantly as the length increases, while that of SiNWs with surface and inner defects increases slightly with the increase in length, and the increasing trend is gentle.
As mentioned above, the TC of SiNWs varies in temperature. The TC of SiNWs with no defects and the TC of SiNWs with surface and inner defects, both having a porosity of 15.68% when the temperature is 300-800 K, are studied. As indicated in Figure 2C, the TC of SiNWs with no defects decreases significantly with increasing temperature and is much higher than that of SiNWs with defects, while the TC of SiNWs with surface and inner defects decreases slightly with increases in temperature. It also can be seen that, when the temperature stays the same, the TC of SiNWs with surface defects is always higher than that of SiNWs with inner defects.

Phonon Density of States of Defective SiNWs
The phonon density of states reflects the frequency distribution of phonons in a material. When phonons are concentrated in the low-frequency region, it indicates that the thermal resistance of the material is relatively high. Therefore, the difference in the TC of the material can be effectively explained by analyzing the phonon density of states, which can be calculated as (Sokhan et al., 2000;Pang et al., 2014).
where n is the modulus of the lattice vibration in the interval from ω to (ω + ω). Since the total modulus is equal to the total degrees of freedom, and if a simple crystal has N atoms where ω m is the maximum frequency. The n can be calculated as where dq is the vertical distance between the two equal frequency planes and ds is the area element. By incorporating Equations (5) and (6) into Equation (4), the phonon density of states can be obtained by the following equation (7) g Frontiers in Energy Research | www.frontiersin.org Figure 3 shows the total of local phonon densities of states of SiNWs with no defects and SiNWs with surface defects and inner defects both with porosity of 15.68%. It can be seen that the corresponding frequency of peaks of phonon density of states for SiNWs with no defects is the highest and that for the SiNWs with inner defects is the lowest. Therefore, SiNWs with no defects have the maximum TC. Furthermore, as shown in Figure 3, atoms at the typical characteristic positions are selected to reflect the typical local phonon density of states of the SiNW. For the SiNWs with no defects, the peak of phonon density of states of selected atomic is close to that of the total density of states, and the frequency of the peak is about 15.33 THz. The peak frequency of phonon density of states of SiNWs with surface defects is relatively lower than that of SiNWs with no defects. In addition, a few peaks arise in the low-frequency region. Therefore, it is explained that the TC of SiNWs with surface defects is lower than that of SiNWs with no defects, which is mainly caused by surface scattering of phonons. However, obvious peaks at the low-frequency region are observed for SiNWs with inner defects, which deviate obviously by comparing with SiNWs with no defects and SiNWs with surface defects. However, low frequency leads to a decrease in its TC. It gives the reason why the TC of SiNWs with inner defects is always lower than that of SiNWs with surface defects under the same porosity.

Phonon Participation
Phonon participation, P (Xia et al., 2003), during transportation process can be applied to characterize the localization of the vibration mode, of which the size reflects the participation degree of atoms in the phonon mode during heat transfer. For example, during the heat transfer process, if only one atom in the system with the number of N participates in the energy transfer process, the value of P will be 1/N. If all the atoms in the system participate in the energy transfer process, the value of P will be 1. Therefore, phonon participation of the material can effectively and intuitively reveal the difference in the TCs of SiNWs with different defects. The phonon participation can be calculated as follows: where N and u are the total number and velocity of atom. Figure 4 shows the phonon participation of SiNWs with no defects and two defective SiNWs with porosity of 15.68%. It can be seen in Figure 4 that the average phonon participation of SiNWs with no defects during the heat transfer process is about 0.4828, which is higher than that of SiNWs with surface defects (the average value is 0.3886) and SiNWs with inner defects (the average value is 0.3448). These results verify that the TC of SiNWs with no defects is higher than that of SiNWs with defects, and the TC of SiNWs with surface defects is higher than that of SiNWs with inner defects under the same porosity.

CONCLUSIONS
In this paper, a non-equilibrium molecular dynamics simulation method of velocity-exchanging method is used to investigate the TC of defective SiNWs with different types of defects. The influences of defect position, porosity, length, and temperature on the TC of SiNWs are discussed, and the regulation mechanism of defect position on the TC of SiNWs is revealed. The main conclusions are as follows.
(1) Under the same porosity, SiNWs with surface defects have higher TC than SiNWs with inner defects. The TC of SiNWs gradually decreases with the increase of porosity. If the porosity is > 30%, the TC of defective SiNWs tends to be flattened by the porosity.
(2) At the same temperature, the TC of SiNWs with no defects is higher than that of SiNWs with defects, and the TC of SiNWs with surface defects is higher at the same porosity. With the increase of temperature, the TC of SiNWs decreases, and the impact of temperature on the TC of SiNWs with no defects becomes more obvious. (3) Affected by the scale effect, the TC of SiNWs increases as their length increases. (4) SiNWs with no defects have the highest corresponding frequency of low-frequency peaks of phonon density of states; however, when SiNWs have inner defects, the lowest frequency is observed. (5) The phonon participation of SiNWs with no defects is much higher than that of SiNWs with defects. When porosity stays the same, the average phonon participation of SiNWs with surface defects is higher than that of SiNWs with inner defects.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
HL and XW conducted the simulations and prepared the article. XW and QR contributed to analysis of