Molecular dynamics simulations of strain-dependent thermal conductivity of single-layer black phosphorus

Wu Junwei　Tao Yi　Chen Chen　Chen Yuewen　Chen Yunfei

(1School of Mechanical Engineering, Southeast University, Nanjing 211189, China)(2Jiangsu Key Laboratory for Design and Manufacture of Micro-Nano Biomedical Instruments, Southeast University, Nanjing 211189, China)

Recent years have witnessed some advanced two-dimensional materials such as graphene, molybdenum disulfide (MoS2), hexagonal boron nitride (h-BN)and black phosphorus.Particularly, black phosphorus, one type of these new advanced materials, has exhibited various excellent characteristics.

SLBP has a puckered structure with each atom connected by three adjacent atoms through covalent bonds[1].Similar to graphene, it has two distinct directions, namely, the armchair direction and zigzag direction.SLBP has a direct band gap of 1.5 eV and higher carrier mobility[2].Those characteristics can be applied in nanoelectronic devices, such as field-effect transistors and photo-transistors[3-6].In addition to exploring the application in nanoelectronics[7], there is also much theoretical and experimental studies focusing on the possible device applied in thermoelectricity.Therein, it has been found that SLBP has a high figure of merit (ZT value)at room temperature[8].According to the first-principles calculations, it is greater than 1.0 at 300 K[8], and the extra tensile strain of 8% will make the value increase to 2.12 in the armchair direction[9].However, the ZT value is closely related to thermal conductivity and the relationship between them can be expressed asZ=S2α/κ, whereZ,S,α,κare the ZT value, the Seebeck coefficient, the electrical conductivity and the phononic thermal conductivity, respectively.Thus, in order to enhance the efficiency of thermoelectricity via increasing the ZT value, decreasing the phononic thermal conductivity can be seen as a direct and efficient approach.It is widely acknowledged that the presence of the tensile strain decreases the thermal conductivity of graphene, because the strain can enhance the phonon scattering and the phonon modes will shift down to the lower frequency[10].Consequently, more phonons will be active and the Umklapp scattering effect is prominent.According to the above information, exploring the specific relationship between tensile strain and thermal conductivity on SLBP nanoribbons extensively is important.

With the development of various experimental methods, the thermal conductivity of black phosphorus can be measured precisely.Although experimental observations are intuitive and effective, the obtained results are greatly influenced by the experimental operation skills.MD, a computational method, allows the microscopic dynamics of a molecular system to be investigated at atomic resolution and offers an effective tool to observe the transmission of the heart flux.In this paper, we conduct MD simulations on SLBP nanoribbons with different lengths to explore the thermal conductivity of bulk and the influence of the tensile strain on the thermal conductivity.

1　Model and Method

The structure of SLBP is depicted in Fig.1.The zigzag and armchair directions are placed along thexandydirections, respectively.In the unit cell with four atoms, the length of the structure along the zigzag direction isa1and that along the armchair direction isa2.The lattice constants of the unit cell area1=0.331 nm anda2=0.436 nm in the zigzag and armchair directions, which are in good agreement with other theoretical calculations[7].Moreover, the size of SLBP can be defined by its lengthLand widthW.In this paper, for a zigzag-edged black phosphorus nanoribbon, the number of unit cells along the zigzag edgeMhas been utilized to describe its lengthL.Meanwhile, the number of unit cells along the armchair edgeNclearly describes its widthW.Thus, the geometry of zigzag-edged black phosphorus nanoribbon can be defined by (M,N).The interlayer spacing 0.536 nm[11]of black phosphorus is taken as the thicknessh.In our study, all MD simulations results are performed by the LAMMPS package[12].The Stillinger-Weber (SW)potential parameterized by Jiang[13]is used to describe the covalent interaction between phosphorus atoms.Compared with the first-principles method, the SW potential provides a relative reasonable and reliable thermal conductivity value in both zigzag and armchair directions.Consequently, this potential is applied to explore the effects of the size and the strain on the thermal conductivity in this paper.

Fig.1　Configuration of single-layer black phosphorus

Thermal conductivity is obtained by non-equilibrium molecular dynamics (NEMD)using the LAMMPS package[14].For black phosphorus, thermal transport is anisotropic[15], and we, hence, calculate the thermal conductivity in both zigzag and armchair directions.Along the zigzag direction, five SLBP nanoribbons with different lengths are chosen, namely (20,50), (40,50), (60,50), (80,50), (100,50).Likewise, along the armchair direction, seven SLBP nanoribbons with different lengths are also chosen, namely (50,10), (50,20), (50,30), (50,40), (50,50), (50,100), (50,150), in order to explore the effects of the size on the thermal conductivity.To begin with, the structure is equilibrated at 300 K at a NVT ensemble for 106steps with the time step of 1 fs.Then we change the NVT ensemble to a NVE ensemble.A Berendsen thermostat is used to maintain the hot and cold source at 320 and 280 K with the time of 1 ns, respectively.The heat fluxJwill be given due to the exchange of kinetic energies.According to Fourier’s law of heat conduction, the thermal conductivity can be calculated as

κ=-J(A∂T/∂x)-1

(1)

whereAis the cross-sectional area of SLBP given asMha1(armchair direction)andNha2(zigzag direction).∂T/∂xis the thermal gradient obtained from the temperature curve by linearly fitting into the linear region.

When exploring the relationship between the strain and the thermal conductivity, the box size can be changed slowly to simulate the strain on SLBP nanoribbons[14].The strain is applied when the structure has been equilibrated at 300 K under a NVT ensemble for 1 ns.After the sizes of the box become stable, the temperatures of the hot and cold sources are maintained at 320 and 280 K, respectively, by a Berendsen thermostat with the same time.

2　Results and Discussion

Similar to all the nanomaterials, black phosphorus has a size-dependent thermal conductivity[16].When the sample size reaches nanometer scale, phonon group velocities and phonon MFP will be reduced due to phonon confinement, thus making thermal transport weak[17].Fig.2 shows the thermal conductivity of nanoribbons with different lengths along the armchair and zigzag directions at 300 K.The reciprocal of the thermal conductivity fitted by the inverse thermal conductivity depends on length.The reason for this phenomenon is the fact that the length of SLBP is not significantly longer than phonon MFP, and thus the phonon with MFP longer than the length of SLBP makes no contribution to thermal conductivity.The thermal conductivity of infinite SLBP nanoribbon can be calculated by using the following linear extrapolation[18]:

(a)

(b)

(2)

whereΛ,κandκare the MFP of infinite SLBP nanoribbon, thermal conductivity of film and bulk, respectively.As a result, the thermal conductivity of SLBP along the zigzag direction can be calculated as 26.55 W/mK.Moreover, the MFP of the infinite SLBP nanoribbon is also calculated as 47.9 nm based on Eq.(2).Using the same method, the estimated thermal conductivity of 5.76 W/mK in the armchair direction is obtained.Moreover, the MFP of the infinite SLBP nanoribbon can be calculated as 22.72 nm along the armchair direction, which is smaller than the previous result in the zigzag direction.

It demonstrates that the mechanical strain has a significant effect on the ZT value, and the extra tensile strain of 8% making the ZT value increase to 2.12 in the armchair direction[9].Furthermore, the thermal conductivity has a close relationship with the ZT value.For this reason, it is necessary to investigate the relationship between the mechanical strain and the thermal conductivity.

As a matter of fact, the phenomenon of strain reduction in thermal conductivity has been paid much attention to for flexible structures in a nanometer scale.As for familiar materials, both graphene and carbon nanotube possess this characteristic[19].We choose a (60, 50)SLBP nanoribbon with the size of 20 nm×20 nm to explore the relationship between the thermal conductivity and the strain.In order to maintain the stability of the structure in the simulation process, the strain is exerted along the armchair and zigzag directions, respectively, causing the nanoribbon length in two directions to vary from 0 to 9% and 0 to 3%, respectively.

In Fig.3, along the armchair direction, the effect of the tensile strain on the thermal conductivity must obviously not be overlooked.When the SLBP nanoribbon is under tensile strain, the relative thermal conductivity is greater than 1, proving that the tensile strain along the armchair direction has a positive influence on the thermal conductivity.However, the effect of the strain on the thermal conductivity along the zigzag direction is small.With the increasing strain along the zigzag direction, the change of the relative thermal conductivity is small.Furthermore, when the same strain is applied along two directions, the variation of the relative thermal conductivity along the armchair direction is always larger than that along the zigzag direction.This indicates that the relative thermal conductivity is more sensitive to the structural change along the armchair direction than that along the zigzag direction.

Fig.3　Calculated relative thermal concavity of SLBP along the armchair and zigzag directions as a function of strain at the temperature of 300 K

To further testify this strain-dependent thermal conductivity, the DOS for SLBP nanoribbons with different strains at 300 K should be calculated.In Fig.4(a), the larger strain value can lead to a higher peak value of DOS at the high frequency region along the armchair direction.On the contrary, larger strain value leads to a lower peak value of DOS at the low frequency region.Compared with the change of peak value of DOS along the armchair direction, the value increases at the low frequency region while it decreases at the high frequency region with the increased strain along the zigzag direction, as shown in Fig.4(b).According to the cut-off frequency of DOS, it should be noted that the black phosphorus Debye temperature is 649 K, which greatly exceeds the simulation temperature of 300 K.Hence, high frequency phonons are inactive and the thermal transport is dominated by low frequency phonons.According to the above mentioned, tensile strain can lead to increasing numbers of low frequency phonons and advance thermal transport along the zigzag direction.However, it reduces the thermal transport along the armchair direction.From the diagram of DOS, a reasonable explanation cannot be provided for the elevated relative thermal conductivity along the armchair direction.Undoubtedly, the thermal conductivity is dependent on three factors, namely the specific heat, the MFP of the phonon and the phonon group velocity.Moreover, both the specific heat and the phonon group velocity can be affected by DOS.In order to investigate the factors of increasing thermal conductivity thoroughly, phonon MFP has been fitted using Eq.(2).Fig.5 clearly depicts that the increasing strain can cause MFP increase along the armchair direction.When the strain reaches 9%, the MFP may be about 1.5 times as much as that without strain.The reason for this phenomenon is the puckered structure of black phosphorus.The strain along the armchair direction makes the length of vertical structure reduce owing to the positive Poisson’s ratio between armchair direction and out-of-plane direction[20].Then the structure of SLBP nanoribbon becomes flat, causing reduced phonon scattering and ascendant MFP.Thus, the increasing relative thermal conductivity along the armchair direction can be regarded as the effect of MFP, which plays a decisive role.

However, the phonon MFP decreases with the increasing strain along the zigzag direction.In Fig.5, MFP decreases about 3 nm when the applied strains increase from 0 to 3%.This indicates that the tensile structure along the zigzag direction can lead to more phonon scattering.According to the discovery of the negative Poisson’s ratio between zigzag direction and out-of-plane direction[20], strain applied along the zigzag direction can cause the length of vertical structure to increase.It can be speculated that strain applied along the zigzag direction aggravates phonon scattering and reduces MFP to about 3 nm.On the whole, both DOS and MFP have effect on thermal conductivity.Whereas applying the strain to the SLBP can slightly influence the values of DOS and MFP, the thermal conductivity, therefore, will not change obviously.

(a)

(b)

Fig.5　Comparison of phonon MFP with the strain of 0, 3%, 6% and 9% along the armchair direction and strain of 0,1%, 2% and 3% along the zigzag direction

3　Conclusion

In summary, MD simulation and the SW potential are applied to investigate influences of different lengths and tensile strain on the thermal conductivity of SLBP nanoribbons.Both the thermal conductivity and MFPs of SLBP are estimated to be 26.55 W/mK and 47.9 nm along the zigzag direction while 5.76 W/mK and 22.7 nm along the armchair direction.The simulated thermal conductivity for SLBP nanoribbon along the armchair direction increases a lot with the tensile strain compared with a slight increase in the zigzag direction.When the strain is applied along the armchair direction, the structure of the SLBP nanoribbon becomes more flat due to the positive Poisson’s ratio, resulting in less phonon scattering and a longer MFP.Furthermore, the MFP does play a pivotal role in the contribution to thermal conductivity.Conversely, the negative Poisson’s ratio causes more structures to be puckered and more phonon scattering along the zigzag direction.In this situation, the strain influences both the DOS and MFP values slightly and the thermal conductivity, therefore, will not change obviously.Our findings explore the relationship between strain and thermal conductivity reasonably and can provide a reliable method to estimate the MFP of black phosphorus.Moreover, the effect of tensile strain on phonon MFPs along two different directions has been investigated in detail.All these conclusions can inspire researchers to conduct a more in-depth study of the relationship between mechanical strain and thermal conductivity in the future.

[1]Wei Q, Peng X.Superior mechanical flexibility of phosphorene and few-layer black phosphorus[J].AppliedPhysicsLetters, 2014,104(25):251915.DOI:10.1063/1.4885215.

[2]Tran V, Soklaski R, Liang Y, et al.Layer-controlled band gap and anisotropic excitons in few-layer black phosphorus[J].PhysicalReviewB, 2014,89(23):817-824.DOI:10.1103/physrevb.89.235319.

[3]Li L, Yu Y, Ye G J, et al.Black phosphorus field-effect transistors.[J].NatureNanotechnology, 2014,9(5):372-377.DOI:10.1038/nnano.2014.35.

[4]Liu H, Neal A T, Zhu Z, et al.Phosphorene: An unexplored 2D semiconductor with a high hole mobility[J].ACSNano, 2014,8(4): 4033-4041.DOI:10.1038/nnano.2014.35.

[5]Qiao J, Kong X, Hu Z X, et al.High-mobility transport anisotropy and linear dichroism in few-layer black phosphorus[J].NatureCommunications, 2011,5:4475-4475.DOI:10.1038/ncomms5475.

[6]Hong T, Chamlagain B, Lin W, et al.Polarized photocurrent response in black phosphorus field-effect transistors[J].Nanoscale, 2014,6(15): 8978-8983.DOI:10.1039/c4nr02164a.

[7]Xiao J, Long M, Zhang X, et al.First-principles prediction of the charge mobility in black phosphorus semiconductor nanoribbons[J].JournalofPhysicalChemistryLetters, 2015,6(20): 4141-4147.DOI:10.1021/acs.jpclett.5b01644.

[8]Fei R, Faghaninia A, Soklaski R, et al.Enhanced thermoelectric efficiency via orthogonal electrical and thermal conductances in phosphorene[J].NanoLetters, 2014,14(11): 6393-6399.DOI:10.1021/nl502865s.

[9]Lv H Y, Lu W J, Shao D F, et al.Enhanced thermoelectric performance of phosphorene by strain-induced band convergence[J].PhysicalReviewB, 2014,90(8):085433.DOI:10.1103/physrevb.90.085433.

[10]Ma F, Zheng H B, Sun Y J, et al.Strain effect on lattice vibration, heat capacity, and thermal conductivity of graphene[J].AppliedPhysicsLetters, 2012,101(11): 111904.DOI:10.1063/1.4752010.

[11]Qin G Z, Yan Q B, Qin Z Z, et al.Hinge-like structure induced unusual properties of black phosphorus and new strategies to improve the thermoelectric performance[J].ScientificReports, 2014,4: 6946.DOI:10.1038/srep06946.

[12]Plimpton S.Fast parallel algorithms for short-range molecular dynamics[J].JournalofComputationalPhysics, 1995,117(1): 1-19.DOI:10.1006/jcph.1995.1039.

[13]Jiang J W.Parametrization of Stillinger-Weber potential based on valence force field model: Application to single-layer MoS2and black phosphorus[J].Nanotechnology, 2015,26(31): 315-706.DOI:10.1088/0957-4484/26/31/315706.

[14]Müller-Plathe F.A simple nonequilibrium molecular dynamics method for calculating the thermal conductivity[J].TheJournalofChemicalPhysics, 1997,106(14): 6082-6085.DOI:10.1063/1.473271.

[15]Jain A, McGaughey A J H.Strongly anisotropic in-plane thermal transport in single-layer black phosphorene[J].ScientificReports, 2015,5: 8501.DOI:10.1038/srep08501.

[16]Zhang Y Y, Pei Q X, Jiang J W, et al.Thermal conductivities of single-and multi-layer phosphorene: A molecular dynamics study[J].Nanoscale, 2016,8(1): 483-491.DOI:10.1039/c5nr05451f.

[17]Xu X, Pereira L F C, Wang Y, et al.Length-dependent thermal conductivity in suspended single-layer graphene[J].NatureCommunications, 2014,5: 3689.DOI:10.1038/ncomms4689.

[18]Yang F, Dames C.Mean free path spectra as a tool to understand thermal conductivity in bulk and nanostructures[J].PhysicalReviewB, 2013,87(3): 035437.DOI:10.1103/physrevb.87.035437.

[19]Korznikova E A, Baimova J A, Dmitriev S V.Effect of strain on gap discrete breathers at the edge of armchair graphene nanoribbons[J].EPL(EurophysicsLetters), 2013,102(6): 60004.DOI:10.1209/0295-5075/102/60004.

[20]Jiang J W, Park H S.Negative poisson’s ratio in single-layer black phosphorus[J].NatureCommunications, 2014,5: 4727.DOI:10.1038/ncomms5727.