Thermal conductivity of strained silicon: molecular dynamics insight and kinetic theory approach

TL;DR

This study combines molecular dynamics and kinetic theory to analyze how tensile and compression strains influence silicon's thermal conductivity, providing insights into phonon behavior under various conditions.

Contribution

It introduces a comprehensive approach using multiple interatomic potentials and kinetic theory to understand strain effects on silicon's thermal conductivity.

Findings

Thermal conductivity varies significantly with strain and temperature.

Different interatomic potentials yield consistent trends in thermal conductivity.

Kinetic theory calculations align well with molecular dynamics results.

Abstract

In this work, we investigated tensile and compression forces effect on the thermal conductivity of silicon. We used equilibrium molecular dynamics approach for the evaluation of thermal conductivity considering different interatomic potentials. More specifically, we tested Stillinger-Weber, Tersoff, Environment-Dependent Interatomic Potential and Modified Embedded Atom Method potentials for the description of silicon atom motion under different strain and temperature conditions. Additionally, we extracted phonon density of states and dispersion curves from molecular dynamics simulations. These data were used for direct calculations of thermal conductivity considering the kinetic theory approach. Comparison of molecular dynamics and kinetic theory simulations results as a function of strain and temperature allowed us to investigate the different factors affecting the thermal conductivity of strained silicon.

Tables

Table 1. Table 1: Stillinger-Weber potential parameters for silicon  Stillinger and Weber ( 1985 )

$ϵ$ (eV)	2.1683	$\sigma$ (Å)	2.0951
$a$	1.80	$\lambda$	21.0
$\gamma$	1.20	$\cos {\theta }_0$	-0.333333333333
$A$	7.049556277	$B$	0.6022245584
$p$	4.0

Table 2. Table 2: EDIP potential parameters for silicon  Justo  et al. ( 1998 )

$A$ (eV)	7.982173	$B$ (Å)	1.5075463
$a$ (Å)	3.121382	$c$ (Å)	2.5609104
$\alpha$	3.1083847	$\beta$	0.0070975
$\eta$	0.2523244	$\gamma$ (Å)	1.1247945
$\lambda$ (eV)	1.453310	$\mu$	0.6966326
$u_1$	1.1247945	$\sigma$ (Å)	0.5774108
$u_2$	1.4533108	$Q_0$	0.2523244
$u_3$	0.6966326	$\rho$	1.2085196
$u_4$	1.2085196

Table 3. Table 3: MEAM potential parameters for silicon  Lee ( 2007 )

$A$	0.58	$B$ (dyne/cm2)	0.99$\cdot {10}^{12}$
$r_e$ (Å)	2.35	$E_c$ (eV)	4.63
${\beta }^{\left(0\right)}$	3.55	$t^{\left(0\right)}$	1.0
${\beta }^{\left(1\right)}$	2.57	$t^{\left(1\right)}$	1.8
${\beta }^{\left(2\right)}$	0.0	$t^{\left(2\right)}$	5.25
${\beta }^{\left(3\right)}$	7.5	$t^{\left(3\right)}$	-2.61
$C_{max}$	2.8	$C_{min}$	1.41

Table 4. Table 4: Tersoff potential parameters for silicon  Tersoff ( 1989 , 1990 )

$m$	3.0	$n$	0.78734
$R$ (Å)	2.85	$D$ (Å)	0.15
$c$	1.0039$\cdot {10}^5$	$d$	16.217
$A$ (eV)	1830.8	$B$ (eV)	471.18
${\lambda }_1$ (Å-1)	2.4799	$\beta$	1.1$\cdot {10}^{-6}$
${\lambda }_2$ (Å-1)	1.7322	$\cos {\theta }_0$	-0.59825
${\lambda }_3$ ( Å-1)	0.0

Table 5. Table 5: Parameters of fitting for phonon lifetime

	$A\left(s^{-1-\chi }{K}^{-\xi }\right)$	$\xi$ (units)	$\chi$ (units)	$B\left(K\right)$
compressed	$1.52321\cdot {10}^{-14}$	1.60416	1.08586	0
unstrained	$1.73548\cdot {10}^{-16}$	1.76149	1.03746	0
tensed	$2.87935\cdot {10}^{-18}$	1.91305	0.974069	0

Equations

$${\kappa }_{ij}=\frac{1}{V{k}_B{T}^2}{\int }_0^{\infty }dt⟨J_i(t)J_j(0)⟩$$

$$\kappa =\frac{1}{3V{k}_B{T}^2}{\int }_0^{t_c}dt{⟨\text{J}(t)\text{J}(0)⟩}_{t_s}$$

$${\tilde{\mathbf{u}}}_{k\alpha }(\mathbf{q})=\frac{1}{\sqrt{N_{UC}}}\sum _l{\mathbf{u}}_{lk\alpha }{e}^{-i\mathbf{q}{\mathbf{r}}_l}$$

$${\tilde{\mathbf{G}}}_{k\alpha ,k^{\prime }\beta }(\mathbf{q})=⟨{\tilde{\mathbf{u}}}_{k\alpha }(\mathbf{q}){{\tilde{\mathbf{u}}}^{\ast }}_{k^{\prime }\beta }(\mathbf{q})⟩$$

$${\mathbf{D}}_{k\alpha ,k^{\prime }\beta }(\mathbf{q})=\frac{k_{B}T}{M}[{\tilde{\mathbf{G}}}^{-1}(\mathbf{q}){]}_{k\alpha ,k^{\prime }\beta }$$

$$|{\mathbf{D}}_{k\alpha ,k^{\prime }\beta }(\mathbf{q})-{\delta }_{\alpha \beta }{\delta }_{k{k}^{\prime }}{\omega }^2(\mathbf{q})|=0$$

$${\tau }_q=\frac{1.41}{\sqrt{\frac{⟨A_q\ddot{A_q}⟩}{⟨A_q^2⟩}}\left(\frac{⟨{\ddot{A}}_q^2⟩⟨A_q^2⟩}{⟨A_q\ddot{A_q}{⟩}^2}-1\right)}$$

$${\kappa }_z(q,p)=C(q,p)v^2(q,p)\tau (q,p){\cos }^2[{\theta }_z(q)]$$

$$v=\frac{d\omega (q,p)}{dq}$$

$$C(q,p)=k_B{x}^2\frac{e^x}{V(e^x-1{)}^2}$$

$$x=\frac{\hslash \omega (q,p)}{k_{B}T}$$

$${\kappa }_z=\sum _q\sum _p{\kappa }_z(q,p)$$

$$\begin{gathered} \kappa =6\frac{k_B}{3\frac{a^3}{4}}\underset{0}{\overset{\infty }{\int }}d\omega (\frac{\hslash \omega }{k_{B}T}{)}^2\frac{e^{\frac{\hslash \omega }{k_{B}T}}}{(e^{\frac{\hslash \omega }{k_{B}T}}-1{)}^2}D(\omega )\times \\ \times \sum _p\xi (\omega ,p)v^2(\omega ,p)\tau (\omega ,p) \end{gathered}$$

$$\begin{array}{r}\xi (\omega ,p)\propto \frac{q^2(\omega ,p)}{v(\omega ,p)}\\ \sum _p\xi (\omega ,p)=1\end{array}$$

$${\tau }^{-1}(\omega ,p)={\tau }_U^{-1}(\omega ,p)+{\tau }_{BC}^{-1}(\omega ,p)+{\tau }_d^{-1}(\omega ,p)$$

$$\begin{gathered} {\tau }^{-1}(\omega ,p)={\tau }_U^{-1}(\omega ,p)= \\ =\{\begin{array}{l}{A}_T{\omega }^{{\chi }_T}{T}^{{\xi }_T}\exp (\frac{B_T}{T})\text{TA branches}\\ A_L{\omega }^{{\chi }_L}{T}^{{\xi }_L}\exp (\frac{B_L}{T})\text{LA branches}\end{array} \end{gathered}$$

$$\tau (\omega ,p)=A^{-1}{\omega }^{-\chi }{T}^{-\xi }\exp (-B/T)$$

$$U_{tot}=\sum _{i<j}{\varphi }_2(r_{ij})+\sum _{i<j<k}{\varphi }_3(r_{ij},r_{ik},{\theta }_{ijk}),$$

$${\varphi }_2=A_{ij}{ϵ}_{ij}\left[B_{ij}{\left(\frac{{\sigma }_{ij}}{r_{ij}}\right)}^{p_{ij}}\right]e^{\frac{{\sigma }_{ij}}{r_{ij}-a_{ij}{\sigma }_{ij}}},$$

$${\varphi }_3={\lambda }_{ijk}{ϵ}_{ijk}{\left[cos{\theta }_{ijk}-cos{\theta }_{0ijk}\right]}^2\times$$

$$\times e^{\frac{{\gamma }_{ij}{\sigma }_{ij}}{r_{ij}-a_{ij}{\sigma }_{ij}}}{e}^{\frac{{\gamma }_{ik}{\sigma }_{ik}}{r_{ik}-a_{ik}{\sigma }_{ik}}},$$

$$U_{tot}=\sum _{i<j}{\varphi }_2(r_{ij},Z_i)+\sum _{i<j<k}{\varphi }_3(r_{ij},r_{ik},Z_i),$$

$${\varphi }_2=A\left[{\frac{B}{r_{ij}}}^{\rho }-e^{-\beta Z_i^2}\right]e^{\frac{\sigma }{r_{ij}-a}},$$

$${\varphi }_3=e^{\frac{\gamma }{r_{ij}-a}}{e}^{\frac{\gamma }{r_{ik}-a}}h(\cos {\theta }_{ijk},Z_i),$$

$$Z_i=\sum _{m\ne i}f(r_{im}),f(r)=\{\begin{array}{ll}1& r<c\\ e^{\frac{\alpha }{1-x^3}}& c<r<a\\ 0& r>a\end{array},$$

$$h(l,Z)=\lambda \left[1-e^{Q(Z)(l+\tau (Z){)}^2}+\eta Q(Z)(l+\tau (Z){)}^2\right],$$

$$Q=Q_0{e}^{-\mu Z},\tau =u_1+u_2\left(u_3{e}^{-u_{4}Z}-e^{-2u_{4}Z}\right)$$

$$U_{tot}=\sum _i\left[F_i({\overline{\rho }}_i)+\frac{1}{2}\sum _{j\ne i}{S}_{ij}{\phi }_{ij}(R_{ij})\right],$$

$$F({\overline{\rho }}_i)=A{E}_c({\overline{\rho }}_i/{\overline{\rho }}_i^0)\ln ({\overline{\rho }}_i/{\overline{\rho }}_i^0),$$

$${\overline{\rho }}_i={\overline{\rho }}_i^{(0)}G(\Gamma ),$$

$$G(\Gamma )=\frac{2}{1+e^{-\Gamma }},\Gamma =\sum _{k=1}^3{t}_i^{(k)}{\left(\frac{{\overline{\rho }}_i^{(k)}}{{\overline{\rho }}_i^{(0)}}\right)}^2,$$

Full text

††preprint: JAP19-AR-02647

I Introduction

Tuning of thermal properties of various materials is one of the key demand in material research Volz et al. (2016); Cahill et al. (2003, 2014); Termentzidis (2017). First and foremost, such necessity arises because of continuous miniaturization of core components of various micro-devices. As a result, issues related to the improvement of heat management are more and more crucial. In this context any possibilities of increasing or lowering the thermal conductivity in semiconductor material are very important and can encounter a broad attention in several application fields.

Strain is an effective method for tuning the thermal property of various materials due to its flexibility. Strain can affect the thermal properties of a material by shifting its phonon frequencies Zi, Zhang, and Xie (1992); Sui and Herman (1993), thereby changing its heat capacity, phonon group velocities, and phonon lifetimes, all of which contribute to the lattice thermal conductivity. Many attempts have been done to study the thermal conductivity of solids under strain using experiments Andersson and Backstrom (1988); Goncharov et al. (2012); Murphy et al. (2014) and theoretical simulations Bhowmick and Shenoy (2006); Xu and Li (2009); Parrish et al. (2014). Ross et al. Ross et al. (1984) concluded that the thermal conductivity of semiconductors increases with compressive strain. Picu et al. Picu, Borca-Tasciuc, and Pavel (2003) studied strain effects on a model Lennard-Jones crystal using molecular dynamics simulations. They found that the thermal conductivity increases under compression and decreases under tension. It has been observed that the thermal conductivity decreases continuously when the strain changed from compressive to tensile for bulk silicon, silicon nanowires, and silicon thin films Li et al. (2010).

Besides, in addition to the control and the modification of material properties that can be expected, elastic stresses often arise as a result of technological processing of a crystalline solid, like: porous network formation Marty et al. (2006); Isaiev et al. (2015), nanostructuration Neimash et al. (2016); Fan et al. (2017), amorphization Newby et al. (2013); Orekhov and Stegailov (2018), nanoinclusion Korotchenkov et al. (2014); Kuryliuk et al. (2015). Additionally, other methods of heat fluxes control involving phononics membranes Midtvedt, Isacsson, and Croy (2014), dislocations Ni et al. (2014); Al-Ghalith, Ni, and Dumitrică (2016); Termentzidis et al. (2018); Wang et al. (2017), functionalization of surfaces by different coatings and shells Yang, Chen, and Dresselhaus (2005); Isaiev et al. (2014) lead to generation of strongly heterogeneous fields of elastic stresses.

For these reasons a particular emphasis on the understanding of the role of strain on thermal transport properties at the nanoscale is an urgent need. As a first step in elucidating this issue, one needs more physical insight regarding phonon transport in strained materials and silicon is a good starting point to link theory and experimentation. Silicon is a material of choice in nowadays micro- and nanoelectronics. Particularly, strained silicon is one of the best candidate for implementation in metal oxide semiconductors devices Munguía et al. (2007) because of the possibility of band gap tuning with straining. Furthermore, recent studies revealed mechanical-field control of electron spin qubits of donors impurities in silicon Mansir et al. (2018), which open possibilities of strained silicon application in spin-based quantum technologies Usman et al. (2018). The change of vibration states of the silicon is also an appropriate way to tune its thermal properties.

Molecular dynamic (MD) is an important tool for predicting the thermal properties of bulk semiconductors and nanostructures, including silicon. The accuracy of the results delivered by MD simulations depends critically on the reliability of interatomic potentials. Several semi-empirical potentials have been developed for Si. The most popular Si potentials were proposed by Stillinger and Weber (SW) Stillinger and Weber (1985) and Tersoff Tersoff (1988a). Other Si potential formats include the Environment-Dependent Interatomic Potential (EDIP) Justo et al. (1998), and the Modified Embedded Atom Method (MEAM) Baskes, Nelson, and Wright (1989) potentials. These original potentials were modified by many authors by slightly changing the analytical functions and improving the parametrization Tersoff (1988b, 1989, 1990); Lee and Baskes (2000); Lee et al. (2001). Several studies have been conducted to compare different empirical interatomic potentials for heat transport modeling in silicon Abs da Cruz et al. (2011); Howell (2012); El-Genk, Talaat, and Cowen (2018); however, no study has systematically compared potentials regarding strain effect on the thermal conductivity of bulk silicon. Therefore, the first goal of this work is to test the Tersoff III Tersoff (1988b, 1989, 1990), SW Stillinger and Weber (1985), EDIP Justo et al. (1998), and Second Nearest Neighbour (2NN) MEAM Lee and Baskes (2000); Lee et al. (2001) potentials for their ability to predict the the thermal conductivity of silicon under hydrostatic compressive and tensile strains. Interatomic parameters of each considered potentials are provided in the Appendix A.

In this study, we used equilibrium molecular dynamics (EMD) approach for simulation of thermal conductivity of strained silicon. Since thermal conductivity evaluated with molecular dynamics is sensitive to the interatomic interactions Abs da Cruz et al. (2011); Howell (2012); El-Genk, Talaat, and Cowen (2018), we chose the optimum potential for simulation of thermal conductivity of strained Si silicon. With the chosen potential, we calculated dispersion curves, density of states and phonon lifetimes for different temperatures and strains. However, the significant drawback of MD simulations is the classical description, which induce no energy quantization. Therefore, we applied additionally analytic approach based on the kinetic theory (KT) of gases (Chantrenne et al., 2005) for evaluation of the thermal conductivity with the input data calculated by molecular dynamics. Finally, we investigated correlations between the temperature dependence of thermal conductivity revealed by KT simulations with the data obtained by EMD.

II Simulation Methods

II.0.1 Formation of strained silicon and elastic properties calculation

All MD simulations were carried out with Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS). The initial system was set as a single crystal silicon with a diamond lattice structure. The initial lattice parameter was set to be equal to 5.43 Å. The $x$, $y$ and $z$ dimensions of the silicon slab as well as the simulation domain was set to be equal to 10 lattice cells. It has been shown that with these cell sizes, the system-size effects can be eliminated Sellan et al. (2010); Wang and Ruan (2017). Additionally, we tested the impact of the simulations domain size on the strain-stress curves, and found it insignificant (see Appendix B). The periodic boundary conditions were used in all directions to simulate a bulk crystal. The system was equilibrated under temperature $T$ in the Nosé-Hoover thermostat with the use of isotherm-isobaric (NPT) ensemble during 100 ps with time step equal to 1 fs to achieve equilibrium lattice parameter for each temperature and potential for the case of external pressure absence ($P$=0 Pa). The initial systems after equilibrium were used for thermal conductivity calculations of pristine Si structures.

Then, the pristine systems were slowly hydrostatically deformed in compression or tension. The rate of strain was set to be equal to $10^{-3}$ ps*-1*. After each 0.1 ps of deformation procedure the information regarding strain and induced stress in the system was collected. Such deformation was performed during 200 ps, and the final deformation was equal to $\pm 20\%$ to achieve stress-strain curve visualization. Additionally, for the strains equal to $\pm 1\%,\pm 2\%,\pm 3\%,\pm 4\%$, and $\pm 5\%$ the re-start files with atoms positions and velocities were collected. These files were used for further thermal conductivity evaluation.

Fig. 1 depicts the typical strain-stress curves of a silicon sample in compression/tensile tests for the four considered potentials at 300 K. As one can see from the figure, bulk modulus ($B$) can be estimated from the slope of a curve divided on 3 in the elastic deformation range (blue window in Fig 1). For all potentials $B$ values are the same and are equal to $B=100\pm 8$ GPa. The latter value matches well experimental data Hopcroft, Nix, and Kenny (2010); Mat . High frequency oscillations on strain-stress curves at large tensile deformations are related to material fracture. The abrupt changes of strain-stress curves obtained with Tersoff III and 2NN MEAM potentials at high compressive strain are caused by phase transformation in silicon under hydrostatic pressure Mizushima, Yip, and Kaxiras (1994); Yang, Li, and Q. Jiang” (2004).

II.0.2 Thermal conductivity evaluation

For the thermal conductivity evaluation, we used the Equilibrium Molecular Dynamics (EMD) approach, which is based on Green-Kubo formalism. In this framework, thermal conductivity can be evaluated from auto-correlation of heat flux as follows McGaughey and Kaviany (2004):

[TABLE]

where $V$ is the volume of the system, $k_{B}$ is Boltzmann’s constant, and $\left<J_{i}(t)J_{j}(0)\right>$ is the heat current autocorrelation function. In general, the conductivity is averaged across all three directions to obtain the scalar conductivity.

Since EMD is based on numerical simulations, the integration time as well as the averaging time in Eq. 1 were limited to finite values. Additionally, the thermal conductivity was averaged in different directions. Thus, Eq. 1 was used in the following form:

[TABLE]

where J is the heat flux vector, $t_{c}$ is a finite correlation time for which the integration is carried out, and $t_{s}$ is the sampling time over which the autocorrelation function is accumulated for averaging.

The atom structure input files of EMD calculation of pristine and strained samples are the ones described above. For thermal conductivity calculation, the systems were equilibrated again according to Nosé-Hoover thermostat during 1 ns, but with the use of canonical (NVT) ensemble. After equilibrium, we performed NVE integration during 25 ns. During this integration, the heat fluxes in different directions were collected and recorded. Thermal conductivities were extracted with the Eq. 2. The correlation time was chosen to be equal 10 ps, and the sampling time was 50 ps. It should be noted that we tested also sampling time equal to 100 ps to check stability of calculations. In this case, duration of NVE integration was equal to 100 ns. Particularly, thermal conductivity for Tersoff potential under 300 K with longer sampling time (100 ps) was found to be equal $\kappa$ = $229\pm 8$ W/(m K), and for shorter sampling time (50 ps) $\kappa$ = $222\pm 8$ W/(m K)). Thus, both values are in good agreement, and therefore for further calculation we used only the shorter sampling time equal to 50 ps. Since, thermal conductivity calculated based on Green-Kubo formalism is sensitive to initial conditions, the simulations were carried for 10 different seeds. The final value of thermal conductivity for each potential, strain and temperature was averaged among these seeds.

II.0.3 Calculation of phonon dispersion and density of states

To get more details about strain effect on the silicon thermal conductivity, the phonon density of states (DOS) and the phonon dispersions of silicon were calculated. Those calculations were performed using the FixPhonon module of LAMMPS. With this tool, dispersion curves can be obtained from the Green’s function formalism Kong et al. (2009); Kong (2011). In the latter, we consider the lattice vibrations at finite temperature $T$; the $k^{\rm th}$ basis atoms in the $l^{\rm th}$ unit cell are displaced from their equilibrium positions $\textbf{r}_{lk}$ by an amount of $\textbf{u}_{lk}$ (atom displacement). In the reciprocal space, the displacements are obtained from the Fourier transformation of the real space ones  Kong (2011) according to :

[TABLE]

where $N_{UC}$ is the total number of unit cells in the crystal, $\mathbf{r}_{l}$ is the equilibrium positions of $l^{\rm th}$ unit cell, $\alpha$ is the component of the atomic displacement of the $k^{\rm th}$ atom in Cartesian coordinates, and $\mathbf{q}$ is the wave vector. The Green’s function in reciprocal space given by :

[TABLE]

where superscript $\ast$ denotes the complex conjugate and $\langle\dots\rangle$ denotes the ensemble average. With the Green’s function method, the dynamical matrix can be expressed as:

[TABLE]

Here $M$ is the mass of Si atom. By solving the eigenvalues of dynamical matrix $\mathbf{D}$, we can get the frequencies of all phonon modes:

[TABLE]

where $\omega$ is the phonon frequency. Hence, the relations between $\mathbf{q}$ and $\omega$ are obtained. The phonon DOS curve is computed from the phonon dispersion by dividing the frequency range into many small segments and counting the number of states in each segment.

In this work, the dynamic matrix calculations were performed in the NVE ensemble with a Langevin thermostat. The systems were equilibrated for $0.5\cdot 10^{6}$ time steps, and then system is run in the NVE ensemble for $6\cdot 10^{6}$ time steps to record the atomic displacements and velocities. Using the eigenvalues of the dynamical matrices, the phonon DOS and the phonon dispersions were then calculated using the auxiliary post-processing code $phana$ Kong (2011).

II.0.4 Calculation of phonon relaxation times

The vibrational lifetime was estimated from the recently developed method based on the Monte Carlo-based moments approximation Dickel and Daw (2010a, b); Gao et al. (2014). Beginning with the harmonic force constant matrix, the normal modes are found, indexed by wave vector $q$ and branch $b$. The calculation involves ensemble averaging of products of normal mode amplitudes $A_{b}$ and accelerations $\ddot{A_{b}}$, which are obtained by projecting the atomic displacements and forces onto the normal modes. The mode lifetimes are obtained from Gao, Wang, and Daw (2015)

[TABLE]

Calculations of phonon relaxation times were carried out with $Jazz$ python wrapper for LAMMPS Gao, Wang, and Daw (2015); Jaz , implemented to calculate the lifetimes of vibrational normal modes.

III Elementary kinetic theory approach for phonons propagation

III.0.1 Thermal conductivity evaluation

In this section, we will describe elements of kinetic theory approach for phonons propagation to calculate the thermal conductivity of strained Si. For the calculation, we use the analytical model developed previously by P. Chantrenne et al Chantrenne et al. (2005). The model considers phonons to be particles that follow Bose-Einstein statistics. The thermal conductivity $\kappa$ in the direction $z$ associated with the phonons $(q,p)$ can be written as

[TABLE]

where $q$ is the wave vector, $p$ its polarization, and $v$ the group velocity determined from the dispersion curves:

[TABLE]

with $v$ the angular velocity, $\tau(q,p)$ is the phonon relaxation time due to the phonon-scattering phenomena, $\theta_{z}(q)$ is the angle between the wave vector $q$ and the direction $z$, and $C(q,p)$ is the specific heat per unit volume.

[TABLE]

[TABLE]

The total thermal conductivity is the sum of the individual contributions due to all the wave vectors $q$ and polarizations $p$:

[TABLE]

In our case, the system is isotropic and the size of it can be considered infinite. Thus, the sum in Eq. 12 becomes an integral. It is expressed in terms of angular frequency:

[TABLE]

where 6 is the number of polarizations, $\frac{a^{3}}{4}$ is the volume of a primitive cell, $D(\omega)$ is the phonon density of states ($\int\limits_{0}^{\infty}D(\omega)d\omega=1$), $\tau(\omega,p)$ is the phonon relaxation time, $\xi(\omega,p)$ are coefficients, that express the contributions of different polarizations to DOS.

[TABLE]

III.0.2 Relaxation times

In order to theoretically assess thermal conductivity, the last needed parameter is the phonon relaxation time. In general, several scattering mechanism are at play. According to Mathiessen rule, the resulting relaxation time (RT) is

[TABLE]

where $\tau^{-1}_{U}(\omega,p)$ is RT due to Umklapp scattering, $\tau^{-1}_{BC}(\omega,p)$ is RT due to the presence of system boundaries and $\tau^{-1}_{d}(\omega,p)$ is RT due to defect or impurity scattering. In our case, the system is considered as infinite (periodic boundary conditions set in EMD), so $\tau^{-1}_{BC}(\omega,p)=0$ , and the lattice is supposed to be perfect, so it has no defects. Thus,

[TABLE]

Besides, it should be mentioned that the implemented approach for calculations of relaxation times only gives a global relaxation time. Decomposition of different contribution is a challenging issue and is not the purpose of this work. Therefore, we fitted numerical data evaluated with Eq.7 by the following equation for all polarizations:

[TABLE]

The fitting was performed to overcome issues arised because of the significant scattering of data evaluated by MD (see Fig. 6). The fitting parameters are presented in Appendix C (Table 5).

IV Results and Discussion

In this section, the results of calculations of the thermal conductivity $\kappa$ of Si under hydrostatic strain $\varepsilon_{h}$ using SW, EDIP, 2NN MEAM and Tersoff III potentials at $T=300$K will be discussed. The strain-dependencies of the thermal conductivity predicted with the use of these potentials are plotted in Fig. 2. Firstly, it should be noted that thermal conductivity evaluated with MD simulations overestimates the experimental results for natural Glassbrenner and Slack (1964) and highly isotopically enriched Inyushkin et al. (2018) Si ($\kappa$=130-160 W/(m$\cdot$K) for 300 K). As discussed previously, MD is based on classical laws of motion and quantum phenomena are not taken into account. This behavior in well known, however in further analysis we will only focus on the dependence of the thermal conductivity versus applied strain as a main parameter for the chosen potentials reliability.

Fig. 2 shows that calculated thermal conductivity of unstrained SW silicon is about 393 W/m$\cdot$K. Under compression strain, the thermal conductivity of SW silicon rapidly decreases, which contradicts previously published MD simulation results for bulk silicon Li et al. (2010) and bulk argon Picu, Borca-Tasciuc, and Pavel (2003). For a tensile strain, the SW potential leads to non-monotonic dependence of thermal conductivity, which is not in agreement with other data already revealed by MD method Li et al. (2010); Picu, Borca-Tasciuc, and Pavel (2003). Therefore, the use of such potential is clearly inappropriate for the simulation of the thermal conductivity of silicon under hydrostatic strain.

Next, the value of $\kappa$ resulting from the EDIP calculations for unstrained silicon was found to be about 174 W/(m$\cdot$K). This potential gives also non-physical trends $\kappa(\varepsilon_{h})$: the thermal conductivity increases continuously from compression strain to tensile hydrostatic strain, which is in a sharp contrast to the expected dependence.

Despite, 2NN MEAM potential (see Fig. 2) gives a correct strain dependence of $\kappa$ with the increase of tensile strain, it gives an insignificant increases with the increase of compression strain. Therefore, the 2NN MEAM potentials may not be suitable for studying the strain effect on the thermal conductivity of silicon.

Finally, the calculated value of $\kappa$ for unstrained silicon with Tersoff III potential is about 229 W/(m$\cdot$K). This is in a good agreement with other MD results Lee et al. (2007). It can be noticed that the thermal conductivity of Tersoff silicon decreases continuously when the applied strain is changed from compressive to tensile. For a 4$\%$ tensile hydrostatic strain, there is a 37$\%$ reduction (from 229 W/(m$\cdot$K) to 144 W/(m$\cdot$K)) in thermal conductivity. On the other side, at a 4$\%$ compressive hydrostatic strain, the thermal conductivity increases 22$\%$ to 280 W/(m$\cdot$K). This result qualitatively corresponds to the experimental data obtained by Andersson and Backstrom Andersson and Backstrom (1988) for the thermal conductivity of silicon under hydrostatic pressure. According to their experiments, the thermal conductivity $\kappa$($P$) slowly increases at room temperature with hydrostatic pressure $P$ as $\kappa$/$\kappa_{0}$=1+0.004$P$, where $\kappa_{0}$ is the thermal conductivity of unstrained silicon and $P$ being measured in GPa. From Fig. 1, the strain of -4$\%$ corresponds to the hydrostatic pressure about of 15 GPa. Therefore, the expected increasing of the thermal conductivity according to this equation would be equal to 6$\%$, which is smaller than one calculated with the Tersoff potential. More importantly we notice that the general trend of $\kappa(\varepsilon_{h})$ variations obtained with Tersoff III potential, plotted in Fig. 2, is similar to the previously reported strain effects on the thermal conductivity of bulk silicon Andersson and Backstrom (1988); Li et al. (2010); Xu and Li (2009) and Si nanostructures Yang et al. (2014); Shahraki and Zeinali (2015).

The results shown in Fig. 2 indicate that only Tersoff III potential is reliable for thermal conductivity calculations of Si under hydrostatic strain, while SW, EDIP and 2NN MEAM empirical potentials fails in comparison with experimental Andersson and Backstrom (1988) and DFT results Parrish et al. (2014), especially for compressive strain. Such mismatching arises because these potentials were parametrized to reproduce physical properties of some stable structures of silicon. Thus, these potentials could not describe correctly the phase transformation considering the latter specific parametrization. Silicon undergoes a phase change from the diamond lattice structure to the $\beta$-Sn lattice structure at $T$=300 K when compressive stresses becomes greater than 12 GPa Yang, Li, and Q. Jiang” (2004). This transition, using MD and Tersoff potential, is observed for larger pressure (60 GPa) Mizushima, Yip, and Kaxiras (1994) than the one experimentally measured.

The first possible explanation for the differences between the potentials in a predicting the $\kappa$($\varepsilon_{h}$) dependencies of silicon under hydrostatic compression is the wide range of cutoff radius $r_{c}$ of each potential. The corresponding value $r_{c}$ is 3.77Å for SW potential, 3.12Å for EDIP potential, 4.3Å for 2NN MEAM potential, and 3.0Å for Tersoff III potential. Since these empirical potentials are optimized for a particular local environment of atomic configurations, they are generally not designed for conditions where the number of interacting neighbors can change abruptly. In fact, when the atomic structure is strongly compressed, the the number of neighboring atoms taken into account in the energy calculation may drastically change, leading to sharp energy variation. The fact that the Tersoff III potential is applicable over a relatively wide range of compressive strain is due to its short cutoff range.

Additionally, different functional forms and the number of fitted parameters also lead to various thermal conductivity behavior when modifying strain. For example, the SW potential is a linear combination of two- and three-body interaction terms that stabilize the diamond cubic structure of crystalline silicon. It has only eight independent parameters and it is fitted to reproduce few experimental properties of both crystallized and liquid silicon. The Tersoff III functional form includes many-body interactions thanks to a bond order term. This potential has 13 adjustable parameters, which are determined by heat of formation, cohesive energy, bulk modulus and the relative stability of various Si polytypes obtained from calculations and experiments. The EDIP potential has a functional form similar to that of Tersoff but slightly more complicated. Thirteen parameters are determined by fitting to a fairly small ab-initio database. MEAM is an angle dependent functional form that evolved out of the simpler embedded atom method. Because parameter values of Tersoff III potential were determined by fitting the elastic constants Tersoff (1988b), the potential is more likely to be applicable for the analysis of deformation behavior. The latter suggests that the Tersoff III potential would be a reliable choice for describing Si thermal transport properties under hydrostatic strain. So further consideration of this work will be carried out with this empirical potential.

Figure 3 shows the calculated thermal conductivity of bulk Si crystal with the use of Tersoff III potential under hydrostatic compressive and tensile strains at the temperature range 300-1100 K comparing to the results of the unstrained case. The results indicate that, for hydrostatic compression, the strain field increases the thermal conductivity and enables to affect the conductivity in the whole range of temperature from 300 to about 900 K, while the tensile hydrostatic strain decreases the thermal conductivity of silicon. For the temperatures over $\sim$900 K, the strain effects on the thermal conductivity becomes weak enough to be negligible.

Eventually, we present further analysis of the thermal conductivity variation with hydrostatic strain and temperature based on the Eq. 13. From this equation, we can consider the changes in phonon density of states, phonon group velocity, and phonon relaxation time to understand the effect of hydrostatic strain on the thermal conductivity of Si. Moreover Eq. 13 takes into account the energy distribution on the phonon mode as a function of the temperature level.

First, Fig. 4 shows the strain effect on the phonon DOS $D$($\omega$) for bulk and strained silicon at 3$\%$ compressive and tensile strains. The calculations with the Tersoff III potential for unstrained silicon satisfactorily reproduce both the structure of the experimental dependence $D$($\omega$) and the frequency positions of phonon modes. It is shown that the optical phonon modes show a red shift when the tensile strain is applied, while a blue shift is observed for the compressive strain. At a 3$\%$ tensile hydrostatic strain, the frequency of the optical phonon mode changes from 15.3 to 13.4 THz ($\Delta\omega$=-1.9 THz), while at a 3$\%$ compressive hydrostatic strain, the frequency shift is $\Delta\omega$=2.2 THz. The last result correlates with the experimental Raman spectra obtained by Weinstein and Piermarini Weinstein and Piermarini (1975) for silicon under the hydrostatic pressure. According to their experiments, the frequency $\omega_{LO}$ of the LO phonon mode in Si shifts with pressure $P$ as $\omega_{LO}$=15.57+0.16$P$-0.002$P^{2}$ [THz]. The spectral positions of the LO phonon bands obtained by using this relation for pristine and compressed at 3$\%$ silicon are shown by the straight vertical lines in Fig. 4. The expected frequency shift from the experimental data is about 1.1 THz and is smaller than one obtained from MD simulations with Tersoff III potential.

Consequently, optical phonons will have a lower or higher energy in comparison to the unstrained case when strain is applied. There are similar peak shifts for longitudinal acoustic (LA) phonons, indicating that the hydrostatic strain has a significant effect on them as well. In particular, the frequency shift for LA phonons peak in DOS is $\Delta\omega$=-0.9 THz (from 10.6 THz to 9.7 THz) and $\Delta\omega$=1.0 THz for 3$\%$ tensile and compressive hydrostatic strains, respectively. On the other hand, variations in transverse acoustic (TA) phonons are not significant compared to the LA and optical phonons. It is known that optical phonons contribute little to heat transfer Esfarjani, Chen, and Stokes (2011); Larkin and McGaughey (2013), so strain-induced transformations of phonon DOS may affect the thermal conductivity of silicon mostly due to frequency shift of the LA phonon modes.

In addition, our calculations show that temperature has little effect on the phonon DOS (not presented in Fig. 4). For example, the frequency shift for LA phonons peak in DOS is only $\Delta\omega$=-0.3 THz, when temperature changes from 300 to 1100 K.

The changes in phonon group velocity can be obtained from the dispersion curves, which are plotted in Fig. 5 for the typical hydrostatic strains ( -3$\%$, 0$\%$, and 3$\%$). The dispersion curves for unstrained Si are compared to experimental results, adopted from work[ Flensburg and Stewart (1999)]. Tersoff III overestimate TA modes by more than 40$\%$ for all large wave vectors, which leads to a logical increase in thermal conductivity compared to experimental one. As shown in Fig. 5, as the strain condition changes from compressive to tensile, transverse acoustic phonon branches change little, while the frequency of optical and longitudinal acoustic phonon branches decrease and whole range of frequency becomes narrowed. This agrees with the change of DOS presented in Fig 4. Decreasing of the maximum phonon frequency as the strain moves from compressive to tensile leads to decreasing the phonon group velocities.

The evaluated phonon group velocities are plotted in Fig. 5 by dotted line curves. As one can see, strain does not affect the TA phonon group velocity, while the LA one decrease as the applied strain is changes from compressive to tensile. These observations are consistent with our EMD thermal conductivity calculations presented above (Fig. 3). When the Si structure is under compression, LA phonon frequencies are larger as well as group velocities. It results in a larger phonon heat flux and thus in a bigger thermal conductivity as compared to pristine case. Oppositely, traction of Si structure reduces maximum frequencies and related group velocities and thus heat flux and thermal conductivity. In addition, our calculations show that the temperature has a negligible effect on the dispersion curves of silicon and, consequently, on the phonon velocities. This is confirmed, in Fig. 5 the dashed lines show the calculated dispersion curves for silicon at $T$=1100 K. For example, our calculations show (not presented in Fig. 5), that increase in temperature from 300 to 1100 K leads to changes in group velocity only $\Delta\upsilon\leq$2-4$\%$ for transverse acoustic phonons, and $\Delta\upsilon\leq$1$\%$ for longitudinal acoustic phonons in a wide range of wave vectors.

In Fig. 6, we present the frequency dependent phonon lifetimes in pristine and strained silicon at 300 K and 1100 K. As one can see, lifetimes $\tau$ varies insignificantly as the system moves from compressive to tensile hydrostatic strain. For example, the mean phonon lifetimes for LA phonons with frequency $\omega\simeq$10 THz are about $\tau$=7 ps, $\tau$=8 ps, and $\tau$=6 ps in pristine, compressed at -3$\%$, and tensed at 3$\%$ silicon, respectively. The strain induced changes in phonon lifetimes for TA modes are also insignificant, so reduction of phonon mean free path as the system moves from compressive to tensile hydrostatic strain is mostly associated with the decreasing of the group velocities described above. Additionally, Fig. 6 shows a significant (more than three times) reduction of the phonon lifetimes over almost whole frequency range when temperature increases from 300 K to 1000 K. This reduction results from increasing phonon scattering rate due to phonon-phonon interaction. Therefore, the trend shown on Fig 3 for the temperature dependence of the thermal conductivity is mainly associated with changes in phonon lifetimes. Finally, this analysis shows that the strain and temperature dependence of the thermal conductivity of silicon arises from the change of the phonon DOS, phonon group velocity, which monotonically decreases when strain goes from compressive to tensile, and the temperature dependence of the phonon relaxation time.

The comparison of MD data calculated with Tersoff III potential with the results of analytical simulation can be made from the Fig. 7. As one can see, there is qualitative correlations between results of both approaches. For analytical calculations the DOS and dispersion curves was directly taken from the MD simulations (Fig. 4 and Fig. 5, respectively), which were performed also with Tersoff III potential. The phonon lifetimes were firstly calculated with Eq. 7 for Tersoff III potential, and then, because of data scattering, they were fitted with Eq. 17. The calculated lifetimes with the fitting curves are presented on the Fig. 6. Fitting curves were used for analytical calculations of thermal conductivity with KT analytical approach.

The quantitative correlations is observed only in the temperature range from 500 to 900 K. Some mismatch in the low temperature range (below Debye temperature) can be explained by classical nature of MD simulations; i.e., the mean number of phonons does not depend on their vibration frequency in this case. Using quantum correction (see Appendix D) mostly give general trends. However, thermal conductivity could not be fully corrected with the latter approach Turney, McGaughey, and Amon (2009). Thus, for low temperature range, simulations approaches based on the KT simulations are more reliable. Some mismatches occurring at higher temperatures correspond to the increasing role of higher-order inharmonic phenomena in the strained silicon, because we get close to the melting temperature.

V Conclusions

In this study, we considered thermal transport in crystalline silicon under different magnitude of hydrostatic strains. Particularly, we considered Tersoff III, Stillinger-Weber (SW), Environment-Dependent Interatomic Potential (EDIP), and Second Nearest Neighbour (2NN) Modified Embedded Atom Method (MEAM) potentials for thermal conductivity evaluation of strained silicon with the use of equilibrium molecular dynamics approach. Firstly, the systems were created by expanding and compression of the crystalline silicon, during this step strain-stress curves were extracted for all potentials. All potentials give the similar value of the bulk modulus (slope of the strain-stress curve). However, only Tersoff potential gives dependence of the thermal conductivity on the strain, which corresponds well with experimental data. Thus, this potential was used further for calculation of dispersion curves, density of states and relaxation time for obtaining physical insight regarding the thermal conductivity behavior of strained silicon.

In addition, these data were used in analytic model based on the kinetic theory. In such way, the dependence of thermal conductivity on the temperature were calculated from analytical relations and correlated with the ones evaluated with equilibrium molecular dynamics approach. Both dependencies revealed good qualitative agreement. Nevertheless, some quantitative mismatches were observed for low and high temperatures. Mismatch at low temperatures may arise because quantum corrections (see Appendix D for more details on the latter correction). At the same time the discrepancies for high temperatures may arise as we get close to the melting temperature.

This results bring new insights about strained silicon thermal properties. In future work, the modeling of nanostructures, which exhibits such stained area could be of great interest to tailor phonon transport properties.

Acknowledgements.

This work has been partially funded by the CNRS Energy unit (Cellule Energie) through the project ImHESurNaASA. We want to acknowledge the partial financial support of the scientific pole EMPP of the University of Lorraine. The publication contains the results obtained in the frames of the research work ”Features of photothermal and photoacoustic processes in low-dimensional silicon-based semiconductor systems” (Ministry of Education and Science of Ukraine, state registration number 0118U000242). We wish to acknowledge Dr. Vladimir Lysenko (INL, INSA Lyon) for fruitful discussions.

Appendix A Potential models

A.1 Stillinger-Weber potential

The SW potential is written in the following form Stillinger and Weber (1985):

[TABLE]

[TABLE]

[TABLE]

The corresponding parameters for the SW potential are provided in table 1.

A.2 Environment-Dependent Interatomic Potential

The EDIP potential is written as Justo et al. (1998)

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

The corresponding parameters for the EDIP potential are provided in table 2.

A.3 Modified Embedded Atom Method Potential

In the 2NN MEAM formalism Lee and Baskes (2000); Lee et al. (2001), the total energy of a system is given by

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

The corresponding parameters for the 2NN MEAM potential are provided in table 3.

A.4 Tersoff Potential

The Tersoff potential describes the atomic interaction as follows Tersoff (1989, 1990):

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

The corresponding parameters for the Tersoff potential are provided in table 4.

Appendix B Impact of a cell number in simulation domain on stress-strain curves

The several sizes of simulation domain ($l$) were checked to find appropriate ones for further simulations. As an example, Fig. 8 presented the results of simulations of stress-strain curves for Tersoff III potentials Tersoff (1988a, 1989, 1990). The situation when the simulations domain has number of cells equal to 5, 10, 15, and 20 are demonstrated. In the Fig. 8, $a$ is the lattice parameter.

Appendix C Parameters of fitting of phonon lifetime

Phonon lifetimes evaluated with MD (Eq. 7) were fitted with the following equation:

[TABLE]

Parameters of the fitting are presented in the following table.

Appendix D Impact of quantum correction

Molecular dynamics operates with classical laws of motion to evaluate the systems with a large number of degrees of freedom. Thus, quantum phenomena are naturally excluded from the consideration in this case. One of the commonly used approach for taken them into account is based on the MD temperature rescaling to the temperature which includes quantum phenomena Soleimani et al. (2018). In this case, thermal conductivity may be corrected as follows Lee et al. (1991):

[TABLE]

where $\kappa_{qc}$ is the quantum corrected thermal conductivity, $\kappa_{MD}$ is the thermal conductivity evaluated with MD, $T$ is the temperature, and $T_{MD}$ is the classical temperature, which can be calculated as follows:

[TABLE]

where $T_{D}$ is the Debye temperature.

The Fig. 9 (left scale) demonstrates temperature dependence of the derivative $\frac{dT_{MD}}{dT}$. This derivative defines the ratio of quantum corrected thermal conductivity to MD one (Eq. 46). And as one can see from the figure, the quantum corrected thermal conductivity is less than MD one approximately to 10 % at the 300 K. Close to Debye temperature ($T_{D}$ = 658 K for silicon) the difference is less than 3 %, and only after 1000 K the value of the derivative is close to one. Thus, the impact of the quantum phenomena is significant for the temperatures less than $T_{D}$.

On the right scale of Fig. 9 we presented dependence of MD temperature on the quantum corrected one. We presented the line which corresponds to the case $T(T)$ with the dotted line.

The dependence of the quantum corrected in such way thermal conductivity is presented on the Fig. 10. In the correction we neglected change of Debye temperature in strained silicon. As one can see, that we achieve an excellent agreement of MD data with analytical approach. Nevertheless, the analytic model overestimate MD data for compressed silicon and underestimate for tensile one in the low temperature range. Such mismatch can be partially overcome with the correction of Debye temperature, which arise mainly with slightly increasing of sound velocity (see Fig. 1) in compressed silicon and decreasing in tensile one.

Finally, we should note that the use of the quantum correction in the presented above form is only illustrative, and thermal conductivity could not be fully corrected, see for example results of quantum and classical lattice dynamics simulations presented by Turney et al  Turney, McGaughey, and Amon (2009). Therefore, we decided to leave the part with the quantum correction only in the Appendix. Nevertheless, the presented results can be important for understanding of mismatch of kinetic theory analytic approach and MD data without quantum correction.