Strain effects on the interfacial thermal conductance of graphene/h-BN heterostructure

Fng Liu ,YouKun Gong ,Rui Zou ,Huiming Ning ,Ning Hu ,Yolu Liu ,Lingk Wu,Fuho Mo,Shoyun Fu,Chng Yn

a College of Aerospace Engineering,Chongqing University,Chongqing,400044,China

b School of Mechanical Engineering,Hebei University of Technology,Tianjin,300401,China

c State Key Laboratory of Reliability and Intelligence Electrical Equipment,Hebei University of Technology,Tianjin,300130,China

d National Engineering Research Center for Technological Innovation Method and Tool,Hebei University of Technology,Tianjin,300401,China

e College of Mechanical and Vehicle Engineering,Hunan University,Changsha,410082,China

f School of Chemistry,Physics and Mechanical Engineering,Queensland University of Technology (QUT),Brisbane,QLD,4001,Australia

Keywords:Graphene h-BN Heterostructure Strain engineering Interfacial thermal conductance NEMD

ABSTRACT Previous experimental and computational results have confirmed that the thermal conductivity of a twodimensional (2D) material can be considerably affected by strain.Numerous attention has been paid to explore the relevant mechanisms.However,the strain effects on the interfacial thermal conductance (ITC) of 2D heterostructure have attracted little attention.Herein,the non-equilibrium molecular dynamics (NEMD) simulations were conducted to the graphene/hexagonal boron nitride (GR/h-BN) heterostructure to investigate the strain effects on the ITC.Three types of strains were considered,i.e.,tensile strain,compressive strain,and shear strain.The results indicate that the strain can adjust the ITC for the GR/h-BN heterostructure effectively,and the strain loading direction also influences the ITC.Generally,the tensile strain reduces the ITC of the heterostructure,in addition to the BN-C system at small tensile strain;both the compressive strain and shear strain increase the ITC,especially at a small strain.For the NB-C system,it is more sensitive to the strain loading direction and the yx shear strain of 0.06 is the most effective way to strengthen the ITC.Our results also show that the out-of-plane deformation weakens the in-plane vibration of atoms,leading to a reduction of the interfacial thermal energy transport.

1.Introduction

Since graphene (GR) was fabricated by Geim and Novoselov [1],tremendous attention has been attracted due to its excellent material characteristics [2–7].In recent years,many single-layered 2D materials have been found,e.g.,molybdenum disulfide (MoS2),black phosphorus(BP),tungsten diselenide (WSe2) as well as the h-BN [8–12].Among these 2D materials,insulator h-BN possesses promising material characteristics,for instance,outstanding mechanical properties and large thermal conductivity,which result in receiving extensive attention[13–16].Similar to graphene,2D h-BN possesses a one-atomic layer structure [17].Meanwhile,h-BN is an appropriate candidate for producing heterostructure with GR,not only because it has a similar honeycomb structure with GR,but also because it has small lattice mismatches with GR [18].After forming a 2D in-plane GR/h-BN heterostructure,it obtains adjustable electronic and magnetic characteristics[19–21],which can be applied to manufacture 2D devices,e.g.,bandpass filter,field-effect transistor (FET),quantum tunnelling transistor,thermoelectric device,LED,solar cell,etc.[22–26].

Because of the formation of interfaces after combining the h-BN and GR,its interfacial thermal energy transport capacity becomes critical in determining the thermal characteristics of the GR/h-BN heterostructure.Relevant studies have confirmed that its interfacial thermal characteristics may impact by many factors,for instance,the system size,defect,temperature,doping,interface topography,and loading condition[27–31].For 2D materials,strain and deformation are common,especially,when they are assembled into nanodevices.Typically,two factors are thought to attribute to the strain and deformation,one from constraints imposed by device assembly,and the other from the application environments,e.g.,temperature or moisture changes [32].Recent investigations also indicate that strain engineering can effectively adjust the characteristics of materials,such as the optical and thermal properties[33–35].However,it should be noted that,for different 2D materials,the dependence of strain on the thermal characteristics is diverse.For some 2D materials,such as single-layered GR and MoS2,their thermal conductivities decrease with increasing tensile strain [36,37].In contrast,other 2D materials’ thermal conductivities increase with increasing tensile strain,such as Penta-SiN2[35].For these 2D materials,the diversity originates from the divergent strain dependences of the phonon-phonon scattering rates,heat capacity,and group velocity[32].

Our previous work indicated that doping in the GR/h-BN heterostructure can strengthen the ITC,especially the N atoms doping [29].However,some other factors,such as the strain effect and loading direction,were not considered in our previous work,which may adjust the ITC.Relevant studies regarding the influences of strain on this heterostructure systems’ ITC have been carried out.Such as,Hong et al.[27]applied tensile strain to the GR/h-BN system and obtained that its interfacial thermal resistance(ITR)increases with increasing strain.This tendency has also been found in graphene/silicene monolayer heterostructures[28].Similarly,some critical factors have not been considered in these previous studies,such as shear strain and loading direction.Inspired by the above works,in this paper,how strain engineering influences the ITC for the GR/h-BN heterostructure are systematically investigated.Based on this reason,three different types of strain are considered,i.e.,tensile strain,compressive strain,and shear strain as well as the loading direction.In order to analyze the influences of strain on the ITC of GR/h-BN heterostructure,the non-equilibrium molecular dynamics (NEMD) simulations are applied.We believe this work can provide a fundamental understanding of the strain effects on the ITC of the GR/h-BN heterostructure and provide inspirations for the applications as well as the control of 2D in-plane GR/h-BN heterostructure nanodevices,especially the thermoelectric devices[38,39].

2.Computational method and model

The computational results were obtained by using the NEMD simulations in the large-scale molecular dynamics simulations (LAMMPS)[40] software package.For the C,N,and B atoms in our models,the interatomic interactions were computed by the Tersoff potential[41,42].In the calculation,the periodic boundary conditions (PBC) were considered except thez-axis (out-of-plane direction and using free boundary condition),and thex-axis as shown in Fig.1(a)was chosen as the heat flux direction.Meanwhile,two ends of the model were constrained after applying the corresponding strain.The conjugate gradient(CG)algorithm was implemented to minimize the energy of the systems.To eliminate the residual stress of the systems,the isothermal-isobaric(NPT) ensemble was conducted with a time step of 1 fs,the temperature and pressure for the ensemble were set as 300 K and 0 atm,respectively,and the relaxation last 0.5 ns,the above relaxation processes are similar to our previous works [43–46].After the NPT ensemble,a strain was conducted to the system under the canonical(NVT)ensemble with the same temperature and then relaxed the system for another 0.5 ns.As shown in Fig.1(a),to conduct the NEMD simulations a heat source and a heat sink were introduced,respectively.Using the Langevin thermostat[47],the heat source and heat sink were maintained at 330 K and 270 K,respectively.After a relaxation of～10 ns under the microcanonical (NVE) ensemble,stable temperature distribution in thex-axis(i.e.,heat flux direction) was produced,and the heat flux and temperature distribution of the last 4 ns were collected.To obtain accurate results,10–15 independent simulations were conducted to obtain the average results.

Fig.1.Detailed diagram of the simulation system.(a) NEMD simulation model;two connection types at the interface: (b) BN-C system and (c) NB-C system.

For the in-plane GR/h-BN heterostructure,two interface types are commonly generated during the fabrication of heterostructure,one is the GR connected with h-BN by a zigzag edge,and the other is the GR connected with h-BN through an armchair edge.As the zigzag linking edge is more easily formed than the armchair linking edge in the heterostructure system during the fabrication[48,49],and the zigzag linking edge is usually applied to study the interfacial thermal characteristics of the GR/h-BN heterostructure [29,30,50].Therefore,we only study the interfacial thermal characteristics for the GR/h-BN system with the zigzag linking edge at the interface.The computational models can be seen in Fig.1(b)and 1(c),corresponding to two connection types.At the interface,the C atoms connect with N atoms,as demonstrated in Fig.1(b),this model is named as BN-C system.Similarly,the model in Fig.1(c),where the C atoms connect with B atoms,is named as NB-C system.

In this work,three different types of strains were considered for bothx-axis andy-axis to systematically investigate the strain effects on the ITC of the GR/h-BN heterostructure,i.e.,tensile strain,compressive strain,and shear strain,as presented in Fig.2.For all strains,the maximum strain was set as 0.12,which is less than the fracture strains of the GR and h-BN[51–54].

3.Results and discussion

The ITC is critical to the thermal properties of the heterostructure,which can be obtained by dividing the heat fluxJby the interfacial temperature difference ΔT,the detailed equation can be found in our previous reports[29,55].

It is worth noting that all the strains are uniaxial,the biaxial strains are out of consideration here.The tensile/compressive strains are defined by the following equation,

whereLordenotes the model's original size in the deformation direction,and ΔLrepresents the displacement in the deformation direction.

For small deformation,the shear strain can be calculated as whereLprepresents the model's original length which is perpendicular to the deformation direction(for instance,they-direction of the system for thexyshear deformation),and ΔLsdenotes displacement in the shear orientation(for instance,along thex-axis in thexyshear deformation).

To better understand the impact mechanisms of ITC under different strains,the analyses of phonon density of states (PDOS) or phonon vibration power spectrums (VPS) [55–58] for the system were also performed.The results of VPS or PDOS can be obtained by applying a Fourier transform to the VAF,i.e.,the atomic velocity autocorrelation function[55,56],

whereTmeans the total simulation time,ω represents the frequency,v0andvtare the atomic velocity at the time of 0 andt,respectively,and〈v0·vt〉 denotes the VAF.

Furthermore,in order to obtain the accurate quantification of the overlap for the VPS or PDOS,we calculated the value of overlapSfor the interfacial N,B,and C atoms,and theScan be obtained as follows [30,59],

whereDGRrepresents the PDOS of C atoms,Dh-BNmeans the sum PDOS of B and N atoms,and the relevantDcan be calculated by Eq.(3).

3.1.Uniaxial tensile strain

Fig.2.Schematic of different strain loadings.Tension: (a) and (c);compression: (b) and (d);shear: (e) and (f).

Fig.3.Variation of ITC under tensile strain.(a) Applied x-axial tensile strain;(b) applied y-axial tensile strain.

Fig.3 shows the variations of ITC under different tensile strains for the BN-C and NB-C systems,and these variations are consistent with the reference [60].We can know that the ITC of the BN-C system is much larger than that of the NB-C system regardless of whether the strain is utilized or not.The proposed explanation is that the strength of the C–N bond is larger than that of the C–B bond [30].To better understand the potential mechanisms,we compared the PDOS of the interfacial N,B,and C atoms.Compared with the NB-C system,we can see that the C atoms have more PDOS overlaps with the N and B atoms in the BN-C system,as shown in Fig.4,especially at low frequency (i.e.,within the scope of 0–5 THz).We also analyzed the overlap value of PDOS for the strain-free system and obtained that theSof the BN-C system and NB-C system are 0.0187 and 0.0181,respectively.The largerSmeans more thermal energy can be transported across the interface,resulting in larger ITC for the BN-C system.It should be noted that for the BN-C system with thex-axial tensile strain,the ITC increases at small strains(i.e.,less than 0.02),and a similar phenomenon can be found in Ref.[61].The possible reason is that after applying small tensile strain,the distancel1between atoms 1–4 increases;however,the angel θ of atoms 2-1-3(or 5-4-6)and distancel2between atoms 5–6(or 2–3)decrease,as shown in Fig.5,which slightly increase the thermal energy transport between atoms 5–6 and 2–3.For deeper insight into the increase of ITC for BN-C under small tensile strain,the overlap of PDOS for the BN-C system was analyzed and found that the overlap of PDOS increases at 5–7 THz and 40–43 THz along with the increases ofS,as shown in Fig.6(a).Beyond the strain of 0.02,the ITC decreases as thex-axial tensile strain increases.From Fig.6(b),we can see that at largex-axial tensile strains,the overlap of PDOS degenerates at high frequencies.Meanwhile,from Figs.4(a) and Fig.6(b),we can see that the G-bands of the C,B,and N atoms become softened when the strain increases,which means a redshift for the G-band and a reduction of phonon group velocities[27,62],leading to degradation for the thermal energy transport across the interface.

3.2.Uniaxial compressive strain

To study the effects of compressive strain on the ITC for the heterostructure system,thex-andy-axial uniaxial compressive strains were applied to the systems,respectively,as shown in Fig.2(b) and (d).The results show that both the ITCs of the BN-C and NB-C systems increase first,as shown in Fig.7,then the ITCs decrease with the compressive strain after reaching the maximum.We also found that under they-axial compressive strain the ITCs of the BN-C and NB-C systems usually larger than those of the strain-free BN-C and NB-C systems.However,at largexaxial compressive strain,the ITCs of the BN-C(at the strain of 0.12)and NB-C systems (at the strain of 0.10) become smaller than those of the corresponding strain-free systems.The possible reason is ascribed to the different deformation modes for the heterostructure systems,as presented in the inset of Fig.7(a)-7(b).

Because the in-plane GR/h-BN heterostructure has a single-layered 2D structure,after applying a small compressive strain to the heterostructure,a part of bonds between atoms are compressed,e.g.,for C atoms the bonds compressed from～1.425 Å to～1.375 Å,which increases the interactions between atoms and enhances the interfacial thermal energy transport capacity,resulting in the improvement of ITC.However,when the compressive strain increases,the compression of bonds does not continue.Without out-of-plane constraints,the heterostructure system produces an out-of-plane deformation,which decreases the in-plane vibration of atoms,leading to the reduction of ITC.At the same time,as shown in the inset of Fig.7(a)-7(b),the out-of-plane deformation increases with the compressive strain.We also observed that the out-of-plane deformation for the heterostructure system is dependant on the compressive direction.

3.3.Shear strain

Applying a shear strain to the freestanding GR/h-BN monolayer is challenging in experiments,to our best knowledge,the relevant studies have not been reported.In MD simulations a shear strain is easy to apply on the 2-D model.Therefore,the influences of shear strain on the ITC for the heterostructure system were also explored.In the current work,both thexyandyxshear strains were considered for the systems,as shown in Fig.2(e)and(f).It should be noted that we could not obtain stable results for the BN-C system with theyxshear strain,the possible reason is that the dramatic deformation in out-of-plane.Therefore,under theyxshear strain,we only discussed the effects of the shear strain on the ITC for the NB-C systems.As shown in Fig.8(a),the ITCs of the BN-C and NB-C systems increase and then decrease with thexyshear strain.The variations of ITC for the NB-C system under theyxshear strain are similar to the above discussion,as shown in Fig.8(b),i.e.,the ITC increases and then decreases with the strain.As under large shear strain,the heterostructure systems produce out-of-plane deformation,i.e.,wrinkle,as shown in Fig.9,which weakens the in-plane vibration of atoms,resulting in the reduction of ITC.

Fig.4.PDOS of interfacial atoms for heterostructure.(a) BN-C system;(b) NB-C system.

Fig.5.Atomics location variation under tensile strain.(a) Strain-free;(b) applied tensile strain.

Fig.6.Representative PDOS of the BN-C system with different x-axial tensile strains.(a) 0.02;(b) 0.06.The higher frequency peaks denote the G-band.

Fig.7.Variations of ITC under compressive strain.(a)Applied x-axial strain;(b)applied y-axial strain.The insets demonstrate the representative schematic of model deformation under different compressive strains,i.e.,0.02 and 0.10.

Fig.8.Shear strain effects on the ITC of heterostructure systems.(a) Applied xy shear strain.

For the BN-C heterostructure,compared with the corresponding strain-free system,the compressive strain andxyshear strain can enhance the ITC even at a large strain of 0.12,as shown in Figs.7 and 8(a).Meanwhile,we can see that the ITCs increase with the strain and then decrease after the critical value,as shown in Figs.7 and 8,the reason is that before the critical strain the system is mainly dominated by bond compression and small deformation.However,after the critical strain,the system is mainly dominated by large out-of-plane deformation,which reduces the interfacial thermal energy transport.For the NB-C system,by contrast,only they-axial compressive strain can enhance the ITC within the range of all strain values.Affect by the loading direction,the deformation of the NB-C system is different,as shown in Fig.9.Under this condition,as shown in Fig.8(a),thexyshear strain only enhances the ITC at a small strain(i.e.,＜0.06).However,we can see from Fig.8(b)that theyxshear strain can enhance the ITC at a larger strain (i.e.,＜0.10).It is worth noting that the error bars in Figs.7 and 8 are higher than Fig.3,the main reason is that the initial deformation location is random for the heterostructure systems under compressive strain and shear strain,which leads to the fluctuation of results.

We also compared the maximum ITC for the heterostructure system under different strain types and corresponding strain values,as shown in Table 1.For the BN-C system,both the compressive strain and shear strain have similar enhancement effects on the ITC,and the loading direction has little influence on the enhancement.For the NB-C system,by contrast,both strain type and loading direction have large effects on theenhancement of ITC.It is also found that theyxshear strain can effectively improve the NB-C system's ITC.

Table 1 Enhancement of ITC for the heterostructure systems under strain.

Fig.9.Out-of-plane deformation of NB-C heterostructure systems.(a) Applied xy shear strain;(b) applied yx shear strain.

4.Conclusion

To systematically study the effects of strain on the ITC for the GR/h-BN heterostructures(including two interface types,i.e.,BN-C system and NB-C system),three different types of strain were conducted to the systems.Through the NEMD simulations,we found that the ITC can be controlled by strain engineering.Generally,the tensile strains degrade the ITC of the GR/h-BN heterostructure.It should be noted that a small tensile strain (＜0.02) can slightly improve the ITC of the BN-C system.From the analyses,we can see that both the compressive strain and shear strain can enhance the ITC of the BN-C system.However,for the NB-C system,only they-axial compressive strain can enhance the ITC in the strain range of 0.0–0.12;thex-axial compressive strain,as well as the shear strain,just improve the ITC at a relatively small strain range.Compared with the BN-C system,the NB-C system is more sensitive to the strain loading direction,and theyxshear strain of 0.06 is the most effective way to enhance its ITC.The results also show that the out-ofplane deformation weakens the in-plane vibration of atoms,leading to a reduction of ITC.Our investigation may provide a fundamental understanding of the strain effects on the ITC of the GR/h-BN heterostructure and serves as theoretical guidance for the applications and control of thermoelectric devices.

Declaration of competing interest

The authors declare that they have no conflict of interest.

Acknowledgments

This research was funded by the National Natural Science Foundation of China (11902056,11632004,11902053,and U1864208),the National Key Research and Development Program of China(2018YFC1105800),the National Science and Technology Major Project(2017-VII-0011-0106),the Key Program for International Science and Technology Cooperation Projects of the Ministry of Science and Technology of China(2016YFE0125900),the Key Project of Natural Science Foundation of CQ CSTC(cstc2017jcyjBX0063),Science and Technology Planning Project of Tianjin(20ZYJDJC00030),Key Program of Research and Development of Hebei Province (202030507040009),the Fund for Innovative Research Groups of Natural Science Foundation of Hebei Province (A2020202002) and the Key Project of Natural Science Foundation of Tianjin (S20ZDF077),we also want to thanks the China Postdoctoral Science Foundation funded project (2019M653334 and 2020M680842).