Abstract

Nonequilibrium molecular dynamics (NEMD) simulations are employed to gain an understanding of the effect of strain on the thermal conductivity of Si thin films. The analysis shows that the strain has an appreciable influence on the thermal conductivity of Si thin films. The thermal conductivity decreases as the tensile strain increases and increases as the compressive strain increases. The decrease of the phonon velocities and surface reconstructions generated under strain could explain well the effects of strain on the thermal conductivity of Si thin films.

1. Introduction

The Si thin films have shown important applications in many fields, such as microelectronic, optoelectronic, and microelectromechanical devices. Understanding and predicting thermal transport of Si thin films at the nanoscale are essential for improving the performance and reliability of these devices. Because the characterized size is down to the nanoscale, the heat transfer of Si thin films becomes unusual as compared to bulk materials. Researchers can design nanomaterials with certain thermal properties through probing the factors which affect thermal properties, so the thermal conductivity of Si thin films depending on the size of the system and defects have been researched by Boltzmann transport equation (BTE) based approaches and molecular dynamics (MD) simulations [15].

Nanoscale devices, particularly thin films, usually contain residual strain after fabrication [6], which may considerably affect the thermal transport properties of the material. Many researchers have attempted to simulate the effect of strain on the thermal properties of material using different methods. Bhowmick and Shenoy [7] derived the relationship of thermal conductivity as a function of temperature and strain using the Peierls-Boltzmann formulation; Zhang et al. [8] investigated how the thermal conductivity responded to the external tensile strain by nonequilibrium molecular dynamics (NEMD) simulations on functionalized bilayer graphene sheet; Wang and Shen [9] investigated the thermal properties of Ni/Al laminated structures using NEMD method; Xu and Li [10] presented a model that combined lattice dynamics and the phonon BTE to analyze strain effect on the cross-plane phonon thermal conductivity of silicon wire-germanium host nanocomposites. These studies have shown the significance of strain on the nanoscale thermal transport.

In this work, we carry out NEMD simulations to investigate the thermal properties of Si thin films. We find that strain has a strong influence on the thermal conductivity of thin films. Our work suggests that applying strain could serve as a means of tuning the thermal conductivity of nanoscale material.

2. Molecular Dynamic Simulation

The NEMD simulations of the Si thin films are performed using LAMMPS code. Figure 1 shows the simulation model in the present research. At the initialized condition, the positions of each atom are arranged according to the crystal lattice structure, and atomic velocities meet the Maxwell-Boltzmann distribution function. The thickness of the Si thin films chosen for this study is 10 Å, 15 Å, 20 Å, and 50 Å respectively, where Å is the lattice constant of Si. The periodic boundary conditions are applied in and direction, and the free boundary conditions are applied in direction. The cross-sectional area of is 4 Å   × 4 Å. The heating and cooling zones with the thickness of 3 Å create a temperature gradient in the system by controlling the energy given or taken from these regions. Adiabatic zones with the thickness of 2 Å are constructed to prevent the atoms escaping from the system, and the velocities of each atom in this region are 0.

Tersoff three-body potential model is employed to describe the interaction among different Si atoms in the simulation film that is reliable to describe the semiconductor atoms potential properties [11]. The Newton’s classical equations of motion for those atoms are solved with the Velocity-Verlet integration algorithm. The time step for the simulation is set as 1 fs. The simulations consist of two stages: the first stage is the constant-temperature simulation, in which the temperature is maintained at constant value with a coupling time of 2 × 106 MD steps; the second stage is the constant-energy one with a coupling time of 8 × 106 MD steps to ensure that the whole system reached an equilibrium state. The method proposed by Juncl and Jullien [12] is used to apply a specified heat flux by scaling the velocities of the atoms in the heating and cooling zones of the simulations. A specified amount of energy (set to 1% of ) is added to the heating zone and the same amount of energy is subtracted from the cooling zone at each time step.

The velocities of all the particles in the heating and cooling zones are scaled according to the following equation:where is the velocity of the center mass, which is set as zero at the beginning of our simulations. A scaling factor is calculated from the following equation:where is the kinetic energy of atoms in the hot or cold reservoirs.

Based on Fourier’s law of conduction, the thermal conductivity is given aswhere is the simulation step time and is the cross-sectional area. To obtain the temperature distribution, the system is divided into 20 planes along the direction, and the instantaneous temperature of each plane is calculated using the following equation:where is the mass of atoms, is the velocity of atoms, is the number of atoms in each plane, is Boltzmann constant, and .

For the analysis, 0~15% tensile and compressive strain is imposed on our NEMD simulation system by recalibrating the atomic coordinates. The rate is 0.005 nm per 105 steps. When each load applying process is finished, the thin films need to achieve equilibrium for a certain relaxation time.

3. Results and Discussion

All the simulations are implemented at 500 K. Figure 2 displays the calculated temperature distributions in the film with strain of −0.1, 0, and 0.1, respectively. The typical distributions have been obtained in all simulation runs. Both ends of the temperature profile show a nonlinear trend attributed to the strong scattering of heating and cooling zones. A linear temperature profile exists in the middle regions, and the gradient increases as the film strain increases.

As the film strain is changed in the range of −0.15~0.15, the thermal conductivity is collected in Figure 3. It is noted that the thermal conductivity decreases as the tensile strain increases and increases as the compressive strain increases. In addition, the strain has an even more significant effect on the thermal conductivity of the thin films with the thickness of 50 Å, which may suggest that the strain effect on the thermal conductivity is largely preserved in the nanoscale configuration over all the sizes. When the thickness of thin film is 10 Å, the thermal conductivity from its initial value of 4.65 W/m·K decreases by about 27% as the tensile strain increases to 0.1 and increases by 18% as the compressive strain increases to 0.1. Such trend is consistent with the earlier reported analyses of other nanostructures in literature [13, 14]. It is also observed from Figure 3 that the thermal conductivity of Si thin films increases with the increase in thickness. This result has been proven earlier by both experimental and computational works. The explanation of the size effect on the thermal conductivity is due to the discretisation of the wave vectors for small systems and the phonon scattering at the boundary surface [15, 16].

To explain the coupling between strain and thermal conductivity, we calculated the phonon density of states (PDOS) for different strain which can describe accurately the thermal transport property of the material. Figure 4 shows that the peaks of PDOS almost coincide for different strain at low frequencies, and they are more sensitive to the strain at high frequencies. Different from bulk material, high frequency phonons play a much more important role in thermal conduction of thin films. The peak moves towards the left when the tensile strain is applied, while it moves towards the right when the compressive strain is applied. The results imply that the tensile strain can shift the high frequency phonons into low frequency phonons, and this shift could slow down the phonon group velocities [17]. Under compressive strain, most phonons will be at a higher frequency which could enhance thermal transportation capacity [18].

In gas kinetic theory, the lattice thermal conductivity can be calculated using the following expression:where is the heat capacity per unit volume and and are the velocity of sound and mean free path of the phonons, respectively. It can be seen from (5) that the decline of phonon velocities due to the tensile strain will decrease the thermal conductivity of thin films.

On the other hand, the surface reconstructions of the simulation system contribute to the change in thermal conductivity. In Figures 5 and 6, the Si thin films with the thickness of 10 Å under different tensile and compressive strain exhibit different atomic configuration characteristics. Under small strain (Figures 5(a), 5(b), 6(a), and 6(b)), the surface disorder is significantly small, and there is no plastic deformation in the simulation system. As the applied strain is increased, dislocations could be generation in the Si thin films under large strain due to the fact that the atoms will recombine with each other to form new bonds. Figures 5(c), 5(d), 6(c), and 6(d) show that the disordered region expands to the whole conducting zones and the plastic deformation occurs gradually, which are stable agreeing with the dislocation simulations result that the partial dislocations are nucleated raising the contribution of extended full dislocations [19]. The simulation partial dislocations are even separated into two pieces in our simulation with the tensile strain of 0.15. The surface reconstructions can scatter phonons significantly which leads to a reduction of the thermal conductivity. For the compressive strain increase to 0.15, a convex distortion is formed near the center of the thin films, while this structure transformation does not reduce the thermal conductivity. It may be due to the fact that the surface effect is lower than the effect of phonon group velocity under compressive strain.

4. Conclusions

This study has performed an investigation into the thermal conductivities of Si thin films with varying strain. The simulation used the NEMD method to obtain a temperature gradient, and the thermal conductivities are then evaluated using the heat current and temperature gradient. The results have shown that the thermal conductivity of the Si thin films could be significantly decreased or increased by tensile or compressive strain. This is due to the mode-specific group velocities of phonons decreasing continuously from compressive strain to tensile strain. In addition, the surface disorder is significantly increased from compressive strain to tensile strain. Our demonstration of large changes in thermal conductivity driven by strain may point to a new direction of dynamic thermal management.

Competing Interests

The authors declare no competing interests.

Acknowledgments

The project was supported by the Fundamental Research Funds for the Central Universities of China (Grant no. DL12BB38).