# Unexpectedly Large Tunability of Lattice Thermal Conductivity of Monolayer Silicene via Mechanical Strain

## Abstract

Strain engineering is one of the most promising and effective routes toward continuously tuning the electronic and optic properties of materials, while thermal properties are generally believed to be insensitive to mechanical strain. In this paper, the strain-dependent thermal conductivity of monolayer silicene under uniform bi-axial tension is computed by solving the phonon Boltzmann transport equation with force constants extracted from first-principles calculations. Unlike the commonly believed understanding that thermal conductivity only slightly decreases with increased tensile strain for bulk materials, it is found that the thermal conductivity of silicene first increases dramatically with strain and then slightly decreases when the applied strain increases further. At a tensile strain of 4%, the highest thermal conductivity is found to be about 7.5 times that of unstrained one. Such an unusual strain dependence is mainly attributed to the dramatic enhancement in the acoustic phonon lifetime. Such enhancement plausibly originates from the flattening of the buckling of the silicene structure upon stretching, which is unique for silicene as compared with other common two-dimensional materials. Our findings offer perspectives of modulating the thermal properties of low-dimensional structures for applications such as thermoelectrics, thermal circuits, and nanoelectronics.

###### pacs:

## I Introduction

Two-dimensional (2D) materials have been extensively studied in the past decade because of their novel physical and chemical properties Novoselov et al. (2004); Seol et al. (2010); Yan et al. (2014) and potential applications.Renteria et al. (2014); Tao et al. (2015) For example, it has been found that graphene has extremely high thermal conductivity,Balandin et al. (2008) which has great potential in the applications including electronic cooling and composite materials. Silicene is the silicon counterpart of graphene and another typical 2D material with a honey-comb lattice structure. Compared to graphene, silicene is more compatible with silicon-based semiconductor devices and therefore has greater potential in nanoelectronic applications. Silicene has also been found to have opened a tunable band gap when a transverse electric field is applied.Drummond et al. (2012); Ni et al. (2012); Kaloni et al. (2013) Monolayer silicene has been successfully fabricated on substrates such as Ag(110),Aufray et al. (2010) Ir(111),Meng et al. (2013) and Ag(111)Lalmi et al. (2010) surfaces. Recently, Tao et al. have demonstrated silicene transistors operating at room temperature.Tao et al. (2015) Although the performance is still moderate and the lifetime of this transistor is only a few minutes, it has attracted significant research interest in silicene based devices.Nguyen (2014); Lian and Ni (2015); Schwierz et al. (2015)

On the other hand, the intrinsic physical properties of silicene, such as lattice thermal conductivity, have been an active area of research. Although the thermal conductivity of silicene has not been measured in experiments due to the difficulty of synthesizing free standing silicene, several numerical simulations have predicted the thermal conductivity of silicene and the results at 300K range from 5 to 69 W/mK.Li and Zhang (2012); Ng et al. (2013); Pei et al. (2013); Hu et al. (2013); Xie et al. (2014); Gu and Yang (2015) Most of the numerical simulations are based on classical molecular dynamics and the discrepancy of results mainly arises from the different interatomic interaction potentials used. Notably, first-principles-based lattice dynamics predicted that the thermal conductivity of silicene is in the range of 20-30 W/mK, Xie et al. (2014); Gu and Yang (2015) which should be more reliable. In our previous first-principles calculations,Xie et al. (2014) the thermal conductivity of 9.4 W/mK at 300 K was not refined due to the small cutoff used for the anharmonic force constant calculation and not imposing the acoustic sum rule.Li et al. (2012) Our new calculation with larger cutoff and denser -mesh gives a value of around 25 W/mK, which is consistent with other literature.Gu and Yang (2015) Despite recent efforts to describe the properties of unstrained silicene, in real applications, nanoscale devices usually contain residual strain after fabrication. Bhowmick and Shenoy (2006) It is thus important to investigate possible strain effects on the property of silicene. It was found that a mechanical tensile strain less than 5% could tune the electronic structure of siliceneYan et al. (2015) and larger tensile strain (7.5%) could induce a semimetal-metal transition.Liu et al. (2012) On the other hand, using first-principles it has been demonstrated that the silicene structure remains buckled even when 12.5% tensile strain is applied.Liu et al. (2012); Lian and Ni (2013)

In comparison to the structural and electronic properties, the strain effect on the lattice thermal conductivity of silicene is less investigated. Pei et al.Pei et al. (2013) and Hu et al. Hu et al. (2013) investigated the effect of uniaxial strain on the thermal conductivity based on the classical non-equilibrium molecular dynamics method. Pei et al. studied tensile strain up to 12% and concluded that the thermal conductivity first increases slightly (around 10% increment) and then decreases with an increased amount of tensile strain. Hu et al. found that the thermal conductivity of silicene sheet and silicene nanoribbon experiences monotonic increase by a factor of 2 with tensile strain up to 18%. The modified embedded-atom method (MEAM)Baskes (1992) and original Tersoff potentialTersoff (1989) were used in their simulations, respectively. However, both potentials are developed for bulk silicon, so directly applying those potentials to the new 2D silicene structure is questionable. For example, the Tersoff potential cannot even reproduce the buckled structure of silicene and the MEAM potential seems to overestimate the buckling distance. It is well known that the interatomic potential directly determines the quality of classical molecular dynamics simulation. Therefore, in order to precisely predict the strain effect on the lattice thermal conductivity of silicene and identify the underlying mechanism, it is necessary to calculate the lattice thermal conductivity of silicene under different strains using a more accurate method.

In this paper, the strain dependent thermal conductivity of monolayer silicene is calculated based on single mode relaxation time approximation (RTA) and iterative solution of the Boltzmann Transport Equation (BTE), where the harmonic and anharmonic force constants are determined using first-principles calculations. The contributions of different modes under different strains are analyzed. The governing mechanisms are analyzed and compared with other materials.

## Ii Methods and Simulation Details

From the solution of the BTE, the lattice thermal conductivity is obtained as Li et al. (2014)

(1) |

where and denote the , , or direction. , , , and are Boltzmann constant, temperature, the volume of the unit cell, and the number of points in the first Brillouin zone respectively. The sum goes over phonon mode that consists of both wave vector and phonon branch . is the equilibrium Bose-Einstein distribution function. is the reduced Planck constant. is the phonon frequency and is the phonon group velocity in direction. The last term is expressed in Ref. 29 as

(2) |

where is the phonon RTA lifetime. is a correction term that eliminates the inaccuracy of RTA by iteratively solving BTE. When is fixed to be equal to zero, the RTA result for thermal conductivity is obtained. Equation (1) can be rearranged with the expression for volumetric phonon specific heatTurney et al. (2009) and the RTA result for thermal conductivity becomes

(3) |

When the isotope scattering, boundary scattering and impurity scattering are ignored, the intrinsic three-phonon RTA lifetime is computed as the inversion of the intrinsic scattering rateLi et al. (2014)

(4) |

where and denote the second and third phonon mode that scatter with phonon mode . and are the intrinsic three-phonon scattering rates for absorption processes and emission processes respectively. Both iterative method and RTA are used to predict thermal conductivity in our calculation. The method is discussed in detail in other literature.Li et al. (2014); Omini and Sparavigna (1995); Turney et al. (2009)

First-principles calculations were carried out using the VASP package.Kresse and Furthmüller (1996) In all our calculations, we used the projector augmented-waves methodBlöchl (1994) and the Perdew-Burke-Ernzerhof (PBE) exchange and correlation.Perdew et al. (1996) A large energy cutoff of 400 eV was chosen. A vacuum spacing of 15 Å was used to prevent interactions between layers. The electronic stopping criteria was eV. The hexagonal symmetry was enforced during the geometry optimization. A hexagonal primitive unit cell was first generated and optimized with 30301 -mesh for electronic integration, and then a 551 supercell was built and re-optimized with a 661 -mesh until the modulus of the force acting on each atom was less than eV/Å. For unstrained silicene, the external pressure in plane was -0.02 kB after supercell optimization. The supercell was then used to compute the force constants required for the phonon dispersion calculation. We used the Phonopy packageTogo et al. (2008) to compute and diagonalize the dynamical matrix and obtain the phonon dispersion curve. The anharmonic force constants were extracted using the code from ShengBTE packageLi et al. (2014) called thirdorder.py. For this calculation, up to the fourth nearest neighbors were considered. Thirdorder.py also applies the sum rules to the anharmonic force constants. Finally, we use the ShengBTE package to compute the thermal conductivity with the harmonic and anharmonic force constants. A 1011011 -mesh, that samples the first Brillouin zone, and a thickness of 4.2 Å, which is twice the van der Waals radius of silicon, were considered. Isotope scattering with the natural isotopic distribution of silicon was considered when solving the BTE, as implemented in ShengBTE. We also tested all the strained cases without isotope scattering and the results were quite similar. For the strained structures we proceeded in the same way. They were however optimized with a fixed strained lattice constant. Strained structures were generated by stretching the optimized unit cell by a certain percentage , where is the optimized lattice constant of unstrained silicene and is that of strained silicene.

## Iii Results and discussion

### iii.1 Structure

The optimized lattice constant of unstrained silicene is 3.87 Å. Figure 1 shows the Si-Si bond length and buckling distance of silicene as a function of strain. The buckling distance first decreases significantly with increasing strain from 0% to 4%, and then decreases slowly from 4% to 6%, and finally stays almost unchanged from 6% to 10%. The small fluctuation from 6% to 10% strain can be attributed to numerical uncertainty. On the other hand, the Si-Si bond length keeps increasing when strain becomes larger. This result is similar to previous first-principles calculations.Liu et al. (2012) It is known that bonding is weaker in silicene than that in graphene because of the longer bond distance. The bonding will be dehybridized into -like bonding,Cahangirov et al. (2009); Şahin et al. (2009) so silicene cannot have a complete planar structure as graphene, even with large strain. We also perceive that the ratio of buckling distance to bond length (nearest-neighbor distance) keeps decreasing with larger strain, which means that the structure becomes more planar under larger strain. The structures of unstrained silicene and silicene under 10% strain are shown in the inset of Fig. 1. It can be clearly seen that strain will reduce the buckling distance and result in longer atomic bond length.

### iii.2 Dispersion

Figure 2 shows the phonon dispersion curve in high-symmetry directions. The phonon dispersion curve plays a crucial role in computing the correct thermal conductivity.Paulatto et al. (2013) Our result of phonon dispersion curve for unstrained silicene is similar to other first-principles calculations. Drummond et al. (2012); Scalise et al. (2012); Şahin et al. (2009) Since silicene has a small buckling, its structure does not have reflectional symmetryLindsay et al. (2010) across the plane. As a result, the vibrational pattern of flexural acoustic mode is not purely out-of-plane, as demonstrated in our previous work.Zhang et al. (2014) Hereafter we denote this flexural acoustic phonon branch as FA branch to avoid confusion with the purely out-of-plane acoustic (ZA) branch in graphene. Another direct consequence of the buckling is that the FA branch is not quadratic near the zone center (). Instead, it has a large linear componentGu and Yang (2015) and therefore a well-defined group velocity (Fig. 2). In addition, in the range of 3-6 THz, the longitudinal acoustic (LA) branch and flexural optical (FO) branch have an avoided crossing.Delaire et al. (2011) This is again because the LA and FO modes have the same symmetry.Dove (1993) This is different from graphene where the out-of-plane optical (ZO) and LA branch cross at about 25 THz (834 ).Yan et al. (2008)

Comparing with the dispersion of unstrained silicene, the dispersion curves of optical phonon modes overall shift downward when the applied tensile strain increases. This shift is mainly due to the reduction of the material stiffness under tensile strain, which is similar to other bulk and low dimensional materials.Li et al. (2010); Picu et al. (2003) It is interesting to note that the frequency gap at the avoided crossing is reduced when the tensile strain increases. We believe this can be attributed to the fact that the structure of strained silicene is becoming more planar when a larger strain is applied, as we discussed earlier.

In order to quantify the acoustic phonon group velocity near the zone center, we calculated the group velocities at point for three acoustic modes. The group velocity of FA mode increases monotonically from 1075.7 m/s for unstrained case to 3287.0 m/s for 10% strain. This is related to the increased stiffness in out-of-plane direction when we apply tensile strain in the in-plane directions. The group velocity of LA mode decreases monotonically from 9548.3 m/s for unstrained case to 7184.6 m/s for 10% strain, which is due to the increase in the Si-Si bond length and the weakened Si-Si atomic interaction in the in-plane directions. For transverse acoustic (TA) mode, the group velocity first increases from 5629.7 m/s for unstrained case to 5804.4 m/s for 4% strain, and then slightly decreases to 5235.7 m/s for 10% strain. It should be noted that the strain dependence of zone center phonon group velocity of silicene is quite different from that of bulk silicon, in which the group velocity only changes slightly under 3% strain.Parrish et al. (2014) The strain dependence is quite similar to graphene, where the slope of ZA branch increases with strain and LA branch decreases with strain.Fugallo et al. (2014)

### iii.3 Thermal conductivity

The thermal conductivity values of silicene at 300K under different tensile strains are computed with both RTA and iterative method, as shown in Fig. 3. The results of unstrained silicene are 25.5 and 27.9 W/mK for RTA and iterative method, respectively. This is similar to the result of Gu et al., Gu and Yang (2015) in which the thermal conductivity of unstrained silicene is predicted to be around 25 W/mK. We note that RTA and iterative method give a similar trend in the strain dependence of thermal conductivity of silicene. Both methods predict that thermal conductivity first increases significantly and then decreases slightly. The highest thermal conductivity value (195.3 W/mK from RTA) appears at 4% strain and is about 7.5 times that of unstrained case. Such a significant increase is quite anomalous. Usually the thermal conductivity of a bulk material, such as silicon,Parrish et al. (2014) diamond,Broido et al. (2012) and argon,Chernatynskiy and Phillpot (2013) would only slightly decrease under tensile strain. For graphene, a similar method predicted that the thermal conductivity only slightly increases with 4% strain.Fugallo et al. (2014); Kuang et al. (2015) In some other empirical molecular dynamics based calculations, it was even found that the thermal conductivity of graphene decreases with tensile strain.Li et al. (2010) With all the strains we considered, the maximum thermal conductivity occurs at 4% strain, which is the turning point where buckling distance of silicene does not decrease significantly with strain any more. This implies that the thermal conductivity of silicene have a strong correlation with the buckling distance. It should be noted that at 7.5% strain silicene is predicted to become metallic,Liu et al. (2012) where electrons may also contribute to the total thermal conductivity. In our calculation, we considered the contribution of phonons to the thermal conductivity only.

In the subsequent discussions, we choose the RTA result to explain this anomalous behavior of silicene because phonon lifetime is well-defined in RTA scheme. No extrapolation of thermal conductivity with respect to -mesh size is performed here because a denser -mesh would produce slightly higher thermal conductivity and similar strain-dependence (we tested up to 2012011). It should be pointed out that in Ref. 21 it was claimed that the thermal conductivity of silicene would diverge with the sample size. Similar conclusions have been drawn in some other literature for graphene. However, there has been strong debate on the possible divergence of thermal conductivity of graphene. For example, Fugallo et al.Fugallo et al. (2014) argue that the thermal conductivity of graphene will converge when the simulated sampling size goes up to 1 mm. In their work, exact phonon BTE is solved and first-principles calculations are employed to extract harmonic and anharmonic interatomic force constants (IFCs). Barbarino et al. Barbarino et al. (2015) also reach the same conclusion with approach-to-equilibrium molecular dynamics simulations for graphene sample of 0.1 mm in size. With a finite -mesh we sampled, we actually exclude those extremely long wavelength acoustic phonon modes, which are believed to be responsible for the possible divergence of the thermal conductivity.Bonini et al. (2012) For real applications a finite sample has to be used, and the wrinkles and defects are generally unavoidable, so the sample cannot have diverged thermal conductivity. Since we are focusing on the strain effect of silicene, we choose a finite -mesh to avoid this issue.

To understand the anomalous strain dependence of thermal conductivity, we first decompose the thermal conductivity contributions into different phonon modes, and the results are also plotted in Fig. 3. Different phonon modes are sorted by their frequencies. The lowest branch is taken as the FA branch while the highest three branches are taken as optical branches. The other two branches are then TA or LA branch. This is a simple and commonly used method to sort different phonon modes in the entire Brillouin zone. We have carefully checked the symmetry of eigenvectors and find this algorithm is reliable for almost all the phonons except at a few high symmetry points. It is evident that the acoustic branches give the dominant contribution over the full strain range, while contribution from optical phonons is in the range of 4% - 22%. The percentage of optical phonon contribution is similar to silicon nanowires.Tian et al. (2011) Figure 3 also shows that LA and TA modes contribute more than half of the total thermal conductivity and govern the trend of strain dependence of thermal conductivity. Especially when the strain is smaller than 4%, TA and LA modes contribute more than 65% of the total thermal conductivity. This is quite different from graphene, for which it is believed that the pure out-of-plane ZA mode has a major contribution to the total thermal conductivity. Seol et al. (2010); Lindsay et al. (2010, 2011) The thermal conductivity contributed by FA mode first increases with strain up to 6% and then slightly decreases. At zero strain, the contribution from FA mode is only about 5% while at 6% strain the contribution from FA mode increases up to around 45%.

### iii.4 Heat capacity, group velocity, and lifetime

From Eq. (3) we know that thermal conductivity is related to volumetric phonon specific heat (heat capacity) , group velocity , and phonon lifetime . In order to find out the dominant factor for the anomalous strain dependence of thermal conductivity, we substitute each of the three terms for unstrained silicene with the value of strained silicene. For example, when heat capacity is replaced, the thermal conductivity is calculated as

(5) |

where the superscript denotes the applied strain in strained silicene and denotes unstrained silicene. The results for three cases are plotted in Fig. 4. We see that the calculated thermal conductivity changes significantly when lifetime is replaced with the value of strained silicene. At 4% strain, the highest value is about 7 times that of unstained case. This shows that the unusual strain dependence of thermal conductivity is mainly due to the change in phonon lifetime. From 0% to 6% strain the thermal conductivity with lifetime replaced changes dramatically, indicating that lifetime is the dominant factor in this range of strain. From 6% to 10% strain the thermal conductivity with lifetime replaced changes about 17.8%, while the thermal conductivity with group velocity or heat capacity replaced changes -12.7% and -5.7% respectively. In this range of strain, these changes are on the same order of magnitude. The three competing factors balance in this range, so the change in the thermal conductivity is small.

### iii.5 Lifetime

To further quantify the lifetime variation under different strains, we plotted the frequency dependent phonon lifetime in Fig. 5. Since acoustic phonons are the dominant heat carriers in the thermal transport in silicene, we only show the lifetimes for acoustic phonons. In addition, those negligible phonon modes whose aggregate contribution to thermal conductivity is less than 0.1% are excluded to reduce the amount of data points. From Fig. 5 it can be seen that, except for a few phonon modes, the acoustic phonon lifetimes of strained silicene are consistently significantly larger than that of unstrained case. The top panel in Fig. 5 indicates that the overall FA phonon lifetimes keep increasing for strain from 0% to 6% but decrease a little for strain from 6% to 10%. The major heat carriers whose aggregate contribution to overall thermal conductivity is more than 50% are those low-frequency acoustic phonons with frequencies under 2.8 THz. The bottom panel shows that the lifetime of TA/LA phonons with frequencies lower than 2.8 THz would overall increase when silicene is strained from 0% to 4%, and then decrease afterwards. The transition from increased to slightly decreased lifetime occurs in the range of 4-6% strain for all the important acoustic phonon modes, which is consistent with the strain-dependent thermal conductivity. It should be noted that such a significant change in phonon lifetime with tensile strain is quite unusual. In Ref.47, the strain-dependent phonon lifetime for solid argon and silicon was calculated using a similar approach. The phonon lifetimes of both materials show quite small strain dependence (changes are within 300%). Here for silicene, the phonon lifetime of the majority of low frequency FA phonons increases by two to three orders of magnitude with strain and the variation of TA/LA lifetime is also one to two orders of magnitude.

### iii.6 Weighted phase space and scattering channel

To understand such a dramatic change of acoustic phonon lifetime, we also calculated the phase space defined in Ref. 56 and the “weighted phase space” defined in Ref. 57 (results not shown). The weighted phase space is an expression of frequency-containing factors that quantifies the phonon scattering probability for a particular dispersion curve. The weighted phase space was used to successfully explain the ultralow thermal conductivity of filled skutterudite YbFeSb.Li and Mingo (2015) However, unlike YbFeSb whose low thermal conductivity is mainly attributed to the allowed phase space for scattering, the variation of weighted phase space of silicene under different strains cannot fully explain the large variation of phonon lifetime. We also calculated thermal conductivity of using the harmonic force constants of 4% strain and anharmonic force constants with 0% strain, the result is only about twice the unstrained thermal conductivity. Together with the calculation result of weighted phase space, we find that the variation of neither harmonic nor anharmonic force constants alone can fully account for the variation of phonon lifetime.

We further investigate the different scattering channels to quantify the importance of different phonon modes in the scattering processes. Figure 6 shows the scattering rates of acoustic phonons along to M direction, and only major scattering channels that has a large contribution to overall scattering rate are included. Scattering rates for emission processes are multiplied by 1/2 to avoid counting the same process twice. Note that along to M direction the FA, TA, and LA branches can be easily separated, so the scattering channels for the three branches are plotted separately in Fig. 6.

Figure 6 (a,b,c) show the scattering rates of different scattering channels of FA phonons for unstrained (0%), 4%, and 10% strained silicene respectively. Note that the legends are the same for these three subfigures while the scales for the -axis are not the same. It can be seen that the total scattering rates of FA phonons decrease orders of magnitude from 0% to 4% but decrease only a little from 4% to 10%. In unstrained silicene, dominant scattering channels are the scattering among the FA modes (i.e., FA+FA FA and FAFA+FA processes). However for strained silicene, either 4% or 10%, the dominant scattering channels become FA+FATA/LA, FA+FAO and FA+TA/LAO, where O indicates optical phonons. Our data also show that scattering rates of FA+TA/LAO processes also decrease from 0% to 4% strain. We attribute the decrease of scattering rates for processes involving odd number of FA phonons to the change of buckling distance. For graphene, because of the reflectional symmetry of the structure, the phonon scattering process involving odd number of ZA phonons is not allowed, leading to a very small scattering rate of ZA mode Lindsay et al. (2010). Our observation of FA mode is in line with their discussion on graphene: when the strain is larger, the silicene structure becomes more planar, so the scattering processes involving odd number of FA modes (especially the scattering process involving 3 FA modes) decreases.

As we noted before, the dramatic change of thermal conductivity from 0% to 4% is mainly due to the in-plane modes (TA and LA). Therefore, we need to further check the scattering channel of TA and LA modes. Figure 6 (d,e,f) plot the scattering rates of TA phonons for 0%, 4%, and 10% strained silicene respectively. Similar to the trend of FA phonons, the total scattering rates reduces significantly from 0% to 4% strain. From 4% to 10% strain, scattering rates still decrease slightly in the high frequency region but increase in low frequency. For unstrained silicene, dominant scattering channels are TA+FATA/LA and TAFA+FA. For silicene under 4% tensile strain, dominant scattering channels are TA+OO, TA+TA/LAO and TAFA+FA. For silicene under 10% strain, dominant scattering channels become TA+TA/LAO and TA+FATA/LA. The variation of scattering rates for different scattering channels for TA phonon modes are the most complicated. The scattering rates of the common dominant process TAFA+FA become smaller with larger strain. For TA+TA/LAO process, the scattering rates increase. For TA+OO, the scattering rates first become larger from 0% to 4%, and then become smaller from 4% to 10%. For TA+FATA/LA, the trend is opposite to previous one: scattering rates first become smaller then become larger. Overall, we find that with larger strain the scattering with FA phonon becomes weaker while the scattering with optical phonon becomes stronger. These competing mechanisms result in the change of the total scattering rates.

Figure 6 (g,h,i) at the bottom are the scattering rates of LA phonons for 0%, 4%, and 10% strained silicene respectively. The total scattering rates of LA phonons decrease monotonically from 0% to 4% and then to 10% strain, but we also note that the total scattering rate of LA phonon does not change as much as FA and TA modes. For different strains LA+TA/LAO is the dominant scattering channel for all the three cases, presumably due to the relatively higher frequency of LA mode. At the low frequency regime of the unstrained case the LAFA+FA is dominate, but at larger strain this channel is becoming less important.

From the analysis of the scattering channel, we have a few observations. First, the scattering process among FA mode is significantly reduced with larger tensile strain, which is due to the reduced buckling distance and more planar structure. This could explain the significant enhancement of phonon lifetime for FA modes as seen in Fig. 5. We should note that this is not the major reason for the enhanced thermal conductivity from 0% to 4% strain, because FA modes have relatively small contributions to thermal conductivity in this range. Second, the scattering rates of TA modes decrease significantly from 0% to 4% strain, mainly due to the reduced scattering with FA modes. This is the major factor that the thermal conductivity of silicene increases in this range. Third, the LA phonon scattering rates do not change significantly under different strains and thus are not responsible much to the large tunability of thermal conductivity. Lastly, we should note that the above analysis is based on the to M region of the Brillouin zone. To distinguish TA and LA phonon is not always possible for any point, so it is not safe to conclude that the enhanced thermal conductivity is mainly due to TA modes. We rather believe that the enhanced thermal conductivity is mainly due to the enhanced lifetime of both LA and TA modes because they scatter less with FA mode when strain is applied.

## Iv Summary

In conclusion, we performed first-principles calculations to predict the lattice thermal conductivity of silicene under strain. Phonon BTE is solved both in the RTA scheme and iteratively in our prediction. Both methods yield a similar trend in the change of thermal conductivity with respect to tensile strain. It is shown that within 10% tensile strain the thermal conductivity of silicene first increases dramatically and then decreases slightly. The maximum thermal conductivity was found when 4% tensile strain was applied, and the value was about 7.5 times that of the unstrained case. Such a dramatic change is quite unusual for solid materials, and could be used as a thermal switch together with thermal diodes to build thermal circuits. This trend is mainly due to the strain-dependent phonon lifetime, which is related to the variations of both harmonic and anharmonic force constants under strain. FA phonon lifetimes increase significantly under tensile strain because the structure becomes more planar. This leads to a large increase of their contribution to overall thermal conductivity, but is not the major reason for the significant change of overall thermal conductivity within 4% strain. The significant enhancement of thermal conductivity from 0% to 4% strain is mainly due to the reduced scattering of TA and LA phonons with FA phonons. Our result suggests that other 2D materials with intrinsic buckling may have similar strain dependence of thermal conductivity, which is left for further investigation.

###### Acknowledgements.

This work was supported by the National Natural Science Foundation of China (Grant No. 51306111) and Shanghai Municipal Natural Science Foundation (Grant No. 13ZR1456000). M.H. acknowledges the support by the Deutsche Forschungsgemeinschaft (DFG) (project number: HU 2269/2-1). Simulations were performed with computing resources granted by HPC () from Shanghai Jiao Tong University. The authors thank Dr. Wu Li for answering questions regarding ShengBTE.### References

- K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, Y. Zhang, S. V. Dubonos, I. V. Grigorieva, and A. A. Firsov, Science 306, 666 (2004).
- J. H. Seol, I. Jo, A. L. Moore, L. Lindsay, Z. H. Aitken, M. T. Pettes, X. Li, Z. Yao, R. Huang, D. Broido, N. Mingo, R. S. Ruoff, and L. Shi, Science 328, 213 (2010).
- R. Yan, J. R. Simpson, S. Bertolazzi, J. Brivio, M. Watson, X. Wu, A. Kis, T. Luo, A. R. Hight Walker, and H. G. Xing, ACS Nano 8, 986 (2014).
- J. D. Renteria, D. L. Nika, and A. A. Balandin, Applied Sciences 4, 525 (2014).
- L. Tao, E. Cinquanta, D. Chiappe, C. Grazianetti, M. Fanciulli, M. Dubey, A. Molle, and D. Akinwande, Nature Nanotechnology 10, 227 (2015).
- A. A. Balandin, S. Ghosh, W. Bao, I. Calizo, D. Teweldebrhan, F. Miao, and C. N. Lau, Nano Letters 8, 902 (2008).
- N. D. Drummond, V. Zólyomi, and V. I. Fal’ko, Physical Review B 85, 075423 (2012).
- Z. Ni, Q. Liu, K. Tang, J. Zheng, J. Zhou, R. Qin, Z. Gao, D. Yu, and J. Lu, Nano Letters 12, 113 (2012).
- T. P. Kaloni, Y. C. Cheng, and U. Schwingenschlögl, Journal of Applied Physics 113, 104305 (2013).
- B. Aufray, A. Kara, S. Vizzini, H. Oughaddou, C. Léandri, B. Ealet, and G. L. Lay, Applied Physics Letters 96, 183102 (2010).
- L. Meng, Y. Wang, L. Zhang, S. Du, R. Wu, L. Li, Y. Zhang, G. Li, H. Zhou, W. A. Hofer, and H.-J. Gao, Nano Letters 13, 685 (2013).
- B. Lalmi, H. Oughaddou, H. Enriquez, A. Kara, S. Vizzini, B. Ealet, and B. Aufray, Applied Physics Letters 97, 223109 (2010).
- N.-T. Nguyen, Micro and Nanosystems 6, 205 (2014).
- C. Lian and J. Ni, Phys. Chem. Chem. Phys. 17, 13366 (2015).
- F. Schwierz, J. Pezoldt, and R. Granzner, Nanoscale 7, 8261 (2015).
- H.-p. Li and R.-q. Zhang, EPL (Europhysics Letters) 99, 36001 (2012).
- T. Y. Ng, J. Yeo, and Z. Liu, International Journal of Mechanics and Materials in Design 9, 105 (2013).
- Q.-X. Pei, Y.-W. Zhang, Z.-D. Sha, and V. B. Shenoy, Journal of Applied Physics 114, 033526 (2013).
- M. Hu, X. Zhang, and D. Poulikakos, Physical Review B 87, 195417 (2013).
- H. Xie, M. Hu, and H. Bao, Applied Physics Letters 104, 131906 (2014).
- X. Gu and R. Yang, Journal of Applied Physics 117, 025102 (2015).
- W. Li, L. Lindsay, D. A. Broido, D. A. Stewart, and N. Mingo, Physical Review B 86, 174307 (2012).
- S. Bhowmick and V. B. Shenoy, The Journal of Chemical Physics 125, 164513 (2006).
- J.-A. Yan, S.-P. Gao, R. Stein, and G. Coard, Physical Review B 91, 245403 (2015).
- G. Liu, M. S. Wu, C. Y. Ouyang, and B. Xu, EPL (Europhysics Letters) 99, 17010 (2012).
- C. Lian and J. Ni, AIP Advances 3, 052102 (2013).
- M. I. Baskes, Physical Review B 46, 2727 (1992).
- J. Tersoff, Physical Review B 39, 5566 (1989).
- W. Li, J. Carrete, N. A. Katcho, and N. Mingo, Computer Physics Communications 185, 1747 (2014).
- J. Turney, E. Landry, A. McGaughey, and C. Amon, Physical Review B 79, 064301 (2009).
- M. Omini and A. Sparavigna, Physica B 212, 101 (1995).
- G. Kresse and J. Furthmüller, Computational Materials Science 6, 15 (1996).
- P. E. Blöchl, Physical Review B 50, 17953 (1994).
- J. P. Perdew, K. Burke, and M. Ernzerhof, Physical Review Letters 77, 3865 (1996).
- A. Togo, F. Oba, and I. Tanaka, Physical Review B 78, 134106 (2008).
- S. Cahangirov, M. Topsakal, E. Aktürk, H. Åahin, and S. Ciraci, Physical Review Letters 102, 236804 (2009).
- H. Şahin, S. Cahangirov, M. Topsakal, E. Bekaroglu, E. Akturk, R. T. Senger, and S. Ciraci, Physical Review B 80, 155453 (2009).
- L. Paulatto, F. Mauri, and M. Lazzeri, Physical Review B 87, 214303 (2013).
- E. Scalise, M. Houssa, G. Pourtois, B. v. d. Broek, V. Afanasâev, and A. Stesmans, Nano Research 6, 19 (2012).
- L. Lindsay, D. A. Broido, and N. Mingo, Physical Review B 82, 235428 (2010).
- X. Zhang, H. Xie, M. Hu, H. Bao, S. Yue, G. Qin, and G. Su, Physical Review B 89, 054310 (2014).
- O. Delaire, J. Ma, K. Marty, A. F. May, M. A. McGuire, M.-H. Du, D. J. Singh, A. Podlesnyak, G. Ehlers, M. D. Lumsden, and B. C. Sales, Nature Materials 10, 614 (2011).
- M. T. Dove, Introduction to lattice dynamics (1993).
- J.-A. Yan, W. Y. Ruan, and M. Y. Chou, Physical Review B 77, 125401 (2008).
- X. Li, K. Maute, M. L. Dunn, and R. Yang, Physical Review B 81, 245318 (2010).
- R. C. Picu, T. Borca-Tasciuc, and M. C. Pavel, Journal of Applied Physics 93, 3535 (2003).
- K. D. Parrish, A. Jain, J. M. Larkin, W. A. Saidi, and A. J. H. McGaughey, Physical Review B 90, 235201 (2014).
- G. Fugallo, A. Cepellotti, L. Paulatto, M. Lazzeri, N. Marzari, and F. Mauri, Nano Letters 14, 6109 (2014).
- D. A. Broido, L. Lindsay, and A. Ward, Physical Review B 86, 115203 (2012).
- A. Chernatynskiy and S. R. Phillpot, Journal of Applied Physics 114, 064902 (2013).
- Y. Kuang, L. Lindsay, and B. Huang, Nano Letters 15, 6121 (2015).
- G. Barbarino, C. Melis, and L. Colombo, Physical Review B 91, 035416 (2015).
- N. Bonini, J. Garg, and N. Marzari, Nano Letters 12, 2673 (2012).
- Z. Tian, K. Esfarjani, J. Shiomi, A. S. Henry, and G. Chen, Applied Physics Letters 99, 053122 (2011).
- L. Lindsay, D. A. Broido, and N. Mingo, Physical Review B 83, 235428 (2011).
- L. Lindsay and D. A. Broido, Journal of Physics: Condensed Matter 20, 165209 (2008).
- W. Li and N. Mingo, Physical Review B 91, 144304 (2015).