Numerical study of thermal conductivity based on phosphorene anisotropy: Including [110] direction and related phosphorus nanotubes

Materials Today Communications 22 (2020) 100814

Contents lists available at ScienceDirect

Materials Today Communications journal homepage: www.elsevier.com/locate/mtcomm

Numerical study of thermal conductivity based on phosphorene anisotropy: Including [110] direction and related phosphorus nanotubes

T

Fan Zhua, Hang Yina,*, Ning Weib, Jing Wanc a

College of Water Conservancy and Civil Engineering, Shandong Agricultural University, Tai’an, 271018, China Jiangsu Key Laboratory of Advanced Food Manufacturing Equipment and Technology, Jiangnan University, Wuxi 214122, China c Shanghai Institute of Applied Mathematics and Mechanics, Shanghai Key Laboratory of Mechanics in Energy Engineering, Shanghai University, Shanghai 200072, China b

A R T I C LE I N FO

A B S T R A C T

Keywords: Thermal conductivity Phosphorene Anisotropy Phonon Nanotube

Thermal transport anisotropy is proved as an intrinsic property of phosphorene in armchair (AC) and zigzag (ZZ) direction. Thermal conductivities along the [110] direction (ST) and its perpendicular direction (PS) are systematically studied using molecular dynamics (MD) simulations. Thermal transport tendency from large to small is in the order of ZZ, PS, ST, and AC. The accuracy of MD in describing anisotropy is proved by acoustic phonon dispersion, and the discrepancy of specific thermal conductivity from the first-principle calculation is mainly due to the mismatch of low-frequency optical phonons. There is an enhancement of thermal conductivity with uniaxial tensile strain along ST and PS direction. Potential energy distribution through atom vibration is chosen to describe strain effect from different calculation methods. Nanotubes rolling from phosphorene have a gentle decrease of thermal conductivity with finite length along the same direction. These results shed light on direction-dependent thermal transport properties of phosphorene along with the applicability in MD simulation.

1. Introduction

according to first-principle calculation. Zhu et al. [31] calculated the thermal conductivities along AC and ZZ directions with density functional theory (DFT) by solving Boltzmann transport equation (BTE), and the values were κAC = 24.3 W/mK and κZZ = 83.5 W/mK for 10 μm phosphorene at 300 K. Zhang et al. [29] and Hong et al. [17] both studied monolayer phosphorene thermal conductivity with nonequilibrium molecular dynamics (NEMD), while their results also have significant differences. It is probably due to the selection of system length and simulation time for calculating heat flux [27]. Zhang et al. [29] also found that tensile strain can lead to a remarkable enhancement of thermal conductivity along AC direction, while the strain effect shows an opposite behavior with uniaxial tension in AC direction according to DFT calculation [25]. Recently, several studies [26,27,32] began to figure out the theory and method for eliminating such differences. They confirmed similar results separately with detailed theoretical discussion and strict simulation settings, but the differences in thermal conductivity values based on MD and DFT still cannot be neglected. The fitted thermal conductivities of phosphorene with infinite length in NEMD simulations are κAC = 20 W/mK and κZZ = 92 W/mK [27], while the DFT calculation claimed giant phonon anharmonicity can have a much lower result: κAC = 4.59 W/mK and κZZ = 15.33 W/ mK [32]. Fortunately, the study of thermal transport anisotropy is mainly focusing on the patterns that can be proved by the above

Phosphorene, owing to its intrinsic unique properties such as tunable direct bandgap, high hole mobility, negative Poisson’s ratio etc. [1–5], can be regarded as one of the remarkable members in two-dimensional (2D) materials’ family at present. Normally, it is representing for the monolayer or few-layered structure exfoliated from bulk black phosphorus [6–9]. Considering the puckered (or hinge-like) crystal structure of phosphorene, its physical properties usually vary with lattice direction showing distinct anisotropy [10–13]. The unit cell of phosphorene as simple orthorhombic lattice is shown in Fig. 1, and the [100] or [010] direction can be named as armchair (AC) or zigzag (ZZ) direction according to the configuration of crystal boundary. In particular, the AC and ZZ directions occupy the prior electrical and thermal transport direction separately, which suggests a potential application in designing thermoelectric device with large thermoelectric figures of merit (ZT) in AC direction [14–16]. As a result, it provides a hot topic of phosphorene in discussing thermal transport anisotropy [17–29]. With the limitation of experimental technique [19,21,30], thermal conductivity of monolayer phosphorene can only be examined through numerical methods up to date. For example, Fei et al. [14] first reported thermal conductivities of monolayer phosphorene κAC = 11.8 W/mK and κZZ = 30 W/mK for each direction with empirical expression ⁎

Corresponding author. E-mail address: [email protected] (H. Yin).

https://doi.org/10.1016/j.mtcomm.2019.100814 Received 17 September 2019; Received in revised form 26 November 2019; Accepted 26 November 2019 Available online 30 November 2019 2352-4928/ © 2019 Elsevier Ltd. All rights reserved.

Materials Today Communications 22 (2020) 100814

F. Zhu, et al.

Fig. 1. Top view of phosphorene with the arrows indicating armchair (AC), zigzag (ZZ) and stair (ST) directions parallel to [100], [010] and [110] directions respectively in simple orthorhombic lattice. PS is representing for the direction perpendicular to stair. Simple monoclinic lattice of 2D phosphorene based on AC and ST directions is also shown with insert map.

and each of the directions can be set parallel to heat flux in NEMD simulation. As for PNTs rolled along the three crystal orientations of phosphorene, we choose AC-PNT and ST-PNT (Fig. 2(d)) for further discussion that have much higher thermal stability than ZZ-PNT [36,42,43]. All the MD simulation processes adopted in this work are chosen to use the open-source code LAMMPS [44]. The phosphorus atom interaction through covalent bonds in phosphorene and PNTs is described with Stillinger-Weber (SW) many-body potential parameterized by Jiang [45]. Phonon dispersion curves generated by SW potential perform similar trend with DFT calculation [46], which suggests a reliable pattern of thermal transport anisotropy with MD simulation. It is also proved [27] that thermal conductivities calculated with Jiang’s parameterized SW potential are in a much reliable range than the SW potential parameterized by Xu et al. [47] For all simulation models, periodic boundary conditions are applied along the thermal flow direction in phosphorene and PNTs. A few layers of atoms at both ends are fixed during simulation, which confines temperature gradient only generated in simulation box with finite length. At the start of simulation, 500 ps relaxation is applied with a constant temperature of 300 K. The time step for integration is 0.5 fs. After that, heat source and cold sink are directly applied with NoséHoover thermostat at 310 K and 290 K respectively, while the other blocks are performed in the micro-canonical ensemble (NVE). The NEMD process is last for 4.5 ns in total, and the temperature of the system comes to a stable state after 1.5 ns. As a result, a constant heat flux J is obtained. According to Fourier’s law, thermal conductivity κ can be calculated:

discussion with a similar trend (anisotropy ratio between ZZ and AC direction basically larger than 2.2) [26]. Based on the pristine phosphorene, recent studies [22,33–35] also focused on the thermal transport properties of in-plane and out-of-plane heterostructure with graphene. As shown in Fig.1, the green atoms are parallel to the [110] direction of phosphorene in simple orthogonal lattice. According to the layout, the word “stair” is chosen to name the direction (ST direction), [36] and together with its perpendicular direction named as PS direction cutting from pristine phosphorene can form orthogonal superlattice for later discussion. Based on AC and ST direction, the unit cell can form single monoclinic lattice for 2D phosphorene. As far as we know, the mechanical and electrical properties along ST direction were studied before [36,37], but its thermal transport property remains unknown. Moreover, the production of monolayer phosphorene nanoribbon with single crystal orientation was recently reported with liquid dispersion method [38], which provide the possibility to form nanotube by rolling-up nanoribbon according to the related experimental and theoretical research [39–41]. Hao et al. [42] investigated the thermal conductivity of armchair black phosphorus nanotubes (AC-PNT, rolling from phosphorene along AC direction) through NEMD simulation, and the fitted value of infinite length is 44.8 W/mK. However, there is no direct comparison of thermal conductivity variation from phosphorene to PNT with the same simulation settings, and the thermal transport properties of other types of PNT (i.e. ST-PNT [36], rolling from phosphorene along ST direction, and PS direction is parallel to tube axis) should also be investigated for a further discussion. In this work, a series of MD simulations are used to study the thermal transport anisotropy of phosphorene and compared the results with above mentioned two kinds of PNTs. Together with lattice dynamic discussion, a clear understanding of thermal transport properties along different directions is based on phonon dispersion relation, and the differences between MD and DFT calculations can also be explained. Meanwhile, the influence of size and strain along ST and PS directions are further studied.

κ=−

J A (∂T ∂l)

(1)

where ∂T/∂l is the temperature gradient along the direction of heat flux, and A is the cross-sectional area of the corresponding model. As for phosphorene, Aphos=w t, where w is the width of phosphorene and t is monolayer thickness taken as balanced layer separation distance 0.524 nm. As for PNT, which can be treated as a hollow cylinder, the cross-sectional area is calculated with APNT=πD t, where D is the diameter of PNT. All the results in this work are generated with averaged values from five times MD simulation. The evaluation of phosphorene thermal anisotropy can be discussed in the view of lattice dynamics. Both simple orthorhombic and simple monoclinic lattice of 2D phosphorene (Fig. 1) are used for phonon spectra calculation. In the corresponding reciprocal space, AC, ZZ and PS direction can be described. According to fluctuation-dissipation theory, atom displacements at a specific temperature of 300 K during MD simulation can be recorded to form a dynamical matrix. Thus, phonon dispersion curves in this work are established with the dynamical matrix [48]. For comparison, we also calculate the phonon

2. Models and methods Fig. 2(a) gives the schematic of NEMD simulation models. Both phosphorene and PNTs are divided into 20 blocks along the heat flux direction, and the 3rd block and 18th block are set as the heat source and cold sink respectively. The first 0.5 nm at both ends labeled with green atoms are fixed during simulation. In order to study the thermal transport anisotropy, two types of phosphorene are conducted. One is the pristine phosphorene with its edge along AC and ZZ direction, and the other one is cutting from pristine phosphorene along ST and PS direction. The local layout for each direction is shown in Fig. 2(c). It can ensure that the in-plane boundaries of phosphorene are all orthotropic 2

Materials Today Communications 22 (2020) 100814

F. Zhu, et al.

Fig. 2. (a) Sketch maps of NEMD process in phosphorene and PNT. The red and blue blocks are representing for the heat source and cold sink separately. The boundary atoms labeled with green are fixed during the simulation. (b) The typical temperature distribution of NEMD in phosphorene and PNT. (c) Detailed layout of two types of phosphorene. In each NEMD process, heat flux is set along the four indicated directions. (d) Top and side view of AC-PNT and ST-PNT. Heat flux is along the tube axial direction.

density of states (PDOS) at frequency ω by taking the Fourier transform of velocity autocorrelation function (VACF)[49]:

D (ω) =

1 2π

(t ) v (0) 〉 iωt e dt ∫ 〈〈vv (0) v (0) 〉

(2)

3. Results and discussion At first, thermal transport anisotropy is discussed combine with size effect. Phosphorene size in this work is determined by the minimum unit cell along PS direction for periodic boundary condition [36]. The width of phosphorene remains 23 nm as a uniform value for all related simulations due to the minor effect of thermal conductivity [29], while the length is varying from 23 nm to 200 nm. In Fig. 3, thermal conductivities along different directions show a great difference at the same length, and for each direction, there is a rising trend with increasing length. The higher thermal conductivity for longer phosphorene is mainly due to long-wave phonon contribution [26]. When comparing with the simulation results with Müller-Plathe approach provided by Zhang et al. [29] at the same length, thermal conductivities in this work are a little higher than the previous work. For example, at the length of 30 nm, the thermal conductivities are κAC = 3.148, κST = 5.672, κPS = 6.995 and κZZ = 8.072 W/mK, respectively. In Zhang’s research, the results are κAC = 2.866 and κZZ = 6.857 W/mK. At the length of 100 nm, the corresponding values are κAC = 6.691, κST = 13.992, κPS = 17.902, and κZZ = 21.2 W/mK. It is also higher than the results of κAC = 5.759 and κZZ = 17.01 W/mK [29]. These differences are mainly originated from the actual length for thermal transport is not equal to system length. In the Müller-Plathe approach with periodic boundary, the effective length is half of the system length [27]. As for the fixed end in this work, the length between the heat source and cold sink is taken 70 percent of system length. With equal valid length, thermal conductivity matches well from both works. Normally, thermal conductivities with finite length are also used to

Fig. 3. Thermal conductivities versus length of phosphorene along different directions.

extrapolate the thermal conductivity (κ∞) of infinite system length. The effective length (L) of phosphorene should larger than effective phonon mean free path (λ). With inverse thermal conductivity versus inverse length fitting linearly, the fitted equation can be taken as the following form [50]:

1 1 ⎛λ = + 1⎞ κ κ∞ ⎝ L ⎠

(3)

However, the value of κ∞ is directly affected by the length variation range, and the appropriate phonon mean free path generated from different studies still remains some discrepancy [27,28,51]. The deviation of the fitting values from previous MD simulations are partially originated with the above reason. Moreover, thermal conductivities from MD simulations [17,27,29] are larger than some of the DFT results 3

Materials Today Communications 22 (2020) 100814

F. Zhu, et al.

Fig. 4. (a) Linear fitting relationships of inverse thermal conductivity versus inverse length along each direction. (b) Definition of direction angles for ST, PS and inplane PeP bond with AC direction. (c) Relationship of thernmal conductivities versus direction angles from AC direction.

[14,18,32]. In this work, linear fitting equations are established with length larger than 80 nm (Fig. 4(a)), and thermal conductivities κ∞ calculated from fitting equation intercept are κAC = 23.98, κST = 50.76, κPS = 88.5 and κZZ = 109.17 W/mK, respectively. The results can be used to describe the thermal transport anisotropy pattern of phosphorene. Apparently, the prior direction for thermal transport is still ZZ direction, and its anisotropy ratio with AC direction is 4.55. Both the thermal conductivity and anisotropy ratio match well with the recent results from multiple MD simulations [27]. The order of thermal conductivity from large to small is ZZ, PS, ST and AC, which has the same pattern of direction angle between AC direction (Fig. 4(b)&(c)). In other words, there is an increasing trend of thermal conductivity with direction angle. The results suggest that the inherent direction-dependent thermal transport property of phosphorene can have a possibility in adjusting heat flux much precisely. Known as thermal transport anisotropy of phosphorene is mainly attributed to the orientation-dependent phonon dispersion rather than phonon-phonon scattering [14,21], phonon dispersion curves are calculated based on MD process (PDC-MD, Fig. 5). The crystal lattice of monolayer phosphorene is usually taken as simple orthorhombic lattice from the definition of bulk black phosphorus, and the corresponding

phonon dispersion result is shown in Fig. 5(a). It is almost identical to the result from the software GULP [52] with the same SW potential [46]. When compared with DFT results (PDC-DFT) [14,18,28,31,32,53], acoustic branches also match well with each other. It is a common agreement [26] that phonons below the bandgap have dominant effects of thermal conductivity, especially for the acoustic branches can explain the differences between AC and ZZ direction. In the Brillouin zone of reciprocal space, Γ-X and Γ-Y paths are representing for ZZ and AC direction, respectively. For example, phonon group velocity of the longitudinal acoustic (LA) phonon branch at ZZ direction is 8.547 km/s, while the group velocity at AC direction is 4.198 km/s. These values are consistent with DFT calculation [14]. It can give a convincible proof for the accuracy of the MD approach in discussing thermal transport anisotropy. In general, all the optical phonon frequencies in PDC-MD are about 2 THz larger than PDC-DFT results. When compared PDOS curves calculated from Kong’s method [48] and VACF according to Eq. (2) in MD (PDOSKong and PDOSVACF in Fig. 5(a)), the bandgap and the phonon peaks in PDOSKong also shows a total blue shift from PDOSVACF. According to Kong’s method, the dynamical matrix for calculating phonon frequency is related to the second-order partial derivative of potential

Fig. 5. Phonon dispersion curves (PDC-MD) of phosphorene based on (a) simple orthorhombic lattice and (b) simple monoclinic lattice. Brillouin zone is marked through the insert map. The three acoustic branches are labeled with different colors. In (a), the orange line for TOz branch is taken from PDC-DFT results [18,31,32] for comparison. PDOS of pristine phosphorene calculated from Kong’s method [48] and VACF is also shown on the right. 4

Materials Today Communications 22 (2020) 100814

F. Zhu, et al.

originated from the phonon scattering enhancement of buckled phosphorene [29]. With the increase of tensile strain, thermal conductivities in both ST and PS directions are also increased. At the strain of 0.03, the increasing rates of thermal conductivity are 3.5 % (ST) and 1.2 % (PS), respectively. These values are larger than ZZ (0.5 %) and smaller than AC (7.5 %) from Zhang’s report [29]. According to the PDOS under tensile strain (Fig. 6(b)), although there is a little red shift at low-frequency peak, the phonon intensification under 4 THz probably indicates the phonon stiffing around B point at Γ-B path (Fig. 5(b)). In addition, the potential energy distribution is chosen to describe phonon anharmonicity with P atom vibration [32]. As shown in Fig. 7, one picked P atom is set to move along the uniaxial strain direction. The black square represents for the model of pristine phosphorene and the red circle represents for the phosphorene with 0.03 tensile strain. When the P atom moves close to the in-plane covalent P atom for each model (i.e. the atom moves left in AC and moves right in ST), there is an arresting potential energy increment in AC direction with tensile strain, while the minor increment in ST and PS direction can be detected with atom displacement under 0.03 nm. The enhancement of in-plane interaction can lead to a higher thermal conductivity that matches well with strain effect in the present study and Zhang’s result [29]. Furthermore, the asymmetry of potential energy distribution becomes higher in AC direction with tensile strain, while the potential energy remains symmetric in ZZ direction with or without strain. It can be used to describe the phonon anharmonicity that reduces the thermal conductivity in DFT calculation [32]. As a result, the decrease of thermal conductivity in AC direction with tensile strain reported by Ong et al. [25] can be explained. In ST and PS direction, phonon anharmonicity is reduced (take the green dash line for reference in Fig. 7(c)&(d)) under tensile strain, which could possibly lead to a positive effect of thermal conductivity in DFT calculation. In the model of PNT, thermal stability is directly affected by tube diameter [36,43], so diameter effect is taken as a prior factor for thermal conductivity discussion. In this process, the length remains 102 nm as uniform value, and the diameter varying from 2.0 nm to 7.5 nm for AC-PNT, along with 3.5 nm–8.7 nm for ST-PNT. The different range of the two kinds of PNT is based on the stability in NEMD. As shown in Fig. 8(a), the thermal conductivity increases with a larger diameter. When the diameter is beyond 3 nm and 5 nm for each model, the rising trend comes to a gentle stage. At the gentle rising stage, thermal conductivities of AC-PNT are larger than ST-PNT for about 4 W/mK. As a consequence, with the same diameter larger than 5 nm, its coupling effect with PNT length can be eliminated. Fig. 8(b) demonstrates the effect of PNT length varying from 30 nm to 200 nm, along with a similar diameter about 5.2 nm for both AC-PNT and ST-PNT. At the length of 150 nm, the thermal conductivity of ACPNT is 28.6 W/mK, while the value generated from the previous report [42] is about 23 W/mK with a diameter of 2.8 nm. The results agreed

energy over equilibrium positions. As a result, the PDC-MD accuracy is directly relied on the potential function, while the VACF with relationship to ensemble average of velocity shares similar variation range in PDOS with DFT [54]. It can be concluded that the thermal conductivity of phosphorene in using SW potential may have some disagreement with the pattern in PDOSVACF. Moreover, there is a soft transverse optical (TOz) branch (orange line in Fig. 5(a)) that can induce strong phonon anharmonicity for a much lower thermal conductivity in phosphorene[32]. The behavior of three optical branches under the bandgap in PDC-MD is all different from TOz branch in PDCDFT around the Γ point. Considering the present SW potential has one missing angle-angle crossing term that can underestimate the Poisson’s ratio of phosphorene [46], the discrepancy of actual thermal conductivity from MD is probably due to the ignoring of out-of-plane phonon anharmonicity. As for the simple monoclinic lattice shown in Fig. 1, Γ-X and Γ-B paths in the reciprocal space are representing for ZZ and PS direction, respectively (Fig. 5(b)). Phonon dispersion along Γ-X direction in Fig. 5(b) is identical to the same direction in Fig. 5(a), which can prove the validity of PDC-MD in simple monoclinic lattice. In the Γ-B direction, the phonon group velocity of LA branch around the center of Brillouin zone is about 17 km/s, which is almost the double of group velocity in Γ-X direction. It is mainly due to the PS direction has a small angle ∼4.7° with the in-plane PeP bond (the red bond in Fig. 4(b)), so that the phonons vibrate in the longitudinal direction is almost parallel to the in-plane PeP bond. When near the B point, the phonon frequency has a prominent decrease that finally leads to the thermal conductivity in PS direction lower than ZZ direction. The strain effect is normally regarded as one of the modulating methods for thermal conductivity [55,56] along with other physical properties [53,57,58], and the material’s residual strain is also a noticeable issue for precision fabrication of nano-devices. Phosphorene with puckered lattice structure can have the abnormal behavior of thermal conductivity increase under tensile strain, but the specific pattern is nearly opposite from MD and DFT calculation [23,25,29]. With uniaxial tensile strain along heat flux direction, Zhang et al. [29] suggest that thermal conductivity in AC direction has a sharp increase and the increasing tendency is much gentle in ZZ direction. Ong et al. [25] only find the increase in ZZ direction with the calculation based on DFT. In this work, we choose the monolayer phosphorene with the length of 52 nm to study strain effect along ST and PS directions. The strain rate is set as 0.001 ps−1. After the strain applied to the specific value, the above mentioned NEMD process is used to measure the thermal conductivity. Fig. 6(a) gives the variation of relative thermal conductivity with both compressive and tensile strain. When under compressive strain, the buckling behavior is easy to generate from 2D material with the small out-of-plane stiffness. The decrease of thermal conductivity is

Fig. 6. (a) Relative thermal conductivities along ST and PS direction with the effect of strain applied in the same direction of heat flux. The insert maps indicate the buckled structures under compressive strain −0.03 at ST and PS directions. (b) PDOS of phosphorene with tensile strain along ST and PS direction. 5

Materials Today Communications 22 (2020) 100814

F. Zhu, et al.

Fig. 7. Potential energy distribution with the displacement of P atom moves along each heat flux direction indicated by the red arrows. The green dash line indicates a harmonic distribution for reference.

respectively. It can be explained according to the PDOS in thermal transport direction (Fig. 9(b)). In the PDOS profile of phosphorene ZZ and AC-PNT, there exists some phonon softening at low-frequency range for AC-PNT. As for the comparison of phosphorene in PS direction and ST-PNT, the red shift of phonons in ST-PNT from related phosphorene is much obvious than AC-PNT.

well with each other when combining diameter effect. Known as heat flux is parallel to the tube axial direction, thermal conductivities of those two PNTs can be compared with the corresponding phosphorene in ZZ and PS direction. For example, at the length of 200 nm, κAC−PNT = 34.19 W/mK and κST-PNT = 24.96 W/mK. These numbers are a little lower than the phosphorene thermal conductivities: κZZ = 35.93 W/mK and κPS = 28.72 W/mK. According to phosphorene phonon dispersion, the flexural acoustic (FA) phonons have a little contribution in thermal transport, and PNT is theoretically bending from phosphorene. As a result, the LA phonons have the dominant effect of thermal conductivity in both phosphorene and PNT from the MD simulation. Hao et al. [42] proposed the theory of phonon restricted transport with different mean free path to explain a much larger differences from PNT to phosphorene. It is not convincible that both thermal conductivity and phonon mean free path are generated from different simulation results. When taken the reciprocal for linear fitting with lengths larger than 80 nm, the estimated thermal conductivities are κAC−PNT = 88.5 W/mK and κST-PNT = 54.64 W/mK (Fig. 9(a)). The descend ratios from corresponding phosphorene are 19 % and 38 %,

4. Conclusion In using the NEMD method, thermal conductivities of phosphorene and PNTs based on ST and PS direction are calculated and discussed along with AC and ZZ direction for comparison. The results show that thermal conductivity from large to small is in the order of ZZ, PS, ST and AC as the heat flux direction in monolayer phosphorene. According to phonon dispersion patterns, the accuracy in describing thermal transport anisotropic tendency with MD can be proved, while the deviation in describing out-of-plane phonon anharmonicity with MD approach makes thermal conductivity larger than DFT results. In addition, uniaxial tensile strain along both ST and PS direction have positive

Fig. 8. Thermal conductivities versus (a) diameter and (b) length with two kinds of PNTs. 6

Materials Today Communications 22 (2020) 100814

F. Zhu, et al.

Fig. 9. (a) Linear fitting relationships of inverse thermal conductivity versus inverse length for AC-PNT and ST-PNT. (b) PDOS of phosphorene and PNT for comparison in thermal transport direction. The diameters of PNTs are ∼5.2 nm.

effect of thermal conductivity in the same direction. At last, the thermal conductivity of AC-PNT is larger than ST-PNT, which depends on the axial direction layout when comparing with phosphorene. With finite length at nanoscale, there is only a gentle decrease of thermal conductivity from phosphorene to PNT. The present work provides a further understanding of direction-dependent thermal transport properties based on phosphorene, and helps out the evaluation of corresponding applicability in MD simulation.

[2] J. Qiao, X. Kong, Z.X. Hu, F. Yang, W. Ji, High-mobility transport anisotropy and linear dichroism in few-layer black phosphorus, Nat. Commun. 5 (2014) 4475. [3] V. Tran, R. Soklaski, Y. Liang, L. Yang, Layer-controlled band gap and anisotropic excitons in few-layer black phosphorus, Phys. Rev. B 89 (2014). [4] J.-W. Jiang, H.S. Park, Negative poisson’s ratio in single-layer black phosphorus, Nat. Commun. 5 (2014) 4727. [5] Y. Du, J. Maassen, W. Wu, Z. Luo, X. Xu, P.D. Ye, Auxetic Black Phosphorus: A 2D Material with Negative Poisson’s Ratio, Nano Lett. 16 (2016) 6701–6708. [6] L. Li, Y. Yu, G.J. Ye, Q. Ge, X. Ou, H. Wu, D. Feng, X.H. Chen, Y. Zhang, Black phosphorus field-effect transistors, Nat. Nanotechnol. 9 (2014) 372–377. [7] A. Castellanos-Gomez, L. Vicarelli, E. Prada, J.O. Island, K.L. Narasimha-Acharya, S.I. Blanter, D.J. Groenendijk, M. Buscema, G.A. Steele, J.V. Alvarez, H.W. Zandbergen, J.J. Palacios, H.S.J. van der Zant, Isolation and characterization of few-layer black phosphorus, 2d Mater. 1 (2014) 025001. [8] H.O. Churchill, P. Jarillo-Herrero, Two-dimensional crystals: phosphorus joins the family, Nat. Nanotechnol. 9 (2014) 330–331. [9] Z. Guo, H. Zhang, S. Lu, Z. Wang, S. Tang, J. Shao, Z. Sun, H. Xie, H. Wang, X.F. Yu, From black phosphorus to phosphorene: basic solvent exfoliation, evolution of Raman scattering, and applications to ultrafast photonics, Adv. Funct. Mater. 25 (2015) 6996–7002. [10] Q. Wei, X. Peng, Superior mechanical flexibility of phosphorene and few-layer black phosphorus, Appl. Phys. Lett. 104 (2014) 251915. [11] F. Xia, H. Wang, Y. Jia, Rediscovering black phosphorus as an anisotropic layered material for optoelectronics and electronics, Nat. Commun. 5 (2014). [12] H. Chen, P. Huang, D. Guo, G. Xie, Anisotropic mechanical properties of black phosphorus nanoribbons, J. Phys. Chem. C 120 (2016) 29491–29497. [13] L. Li, C. Feng, J. Yang, Tensile and compressive behaviors of prestrained singlelayer black phosphorus: a molecular dynamics study, Nanoscale 9 (2017) 3609–3619. [14] R. Fei, A. Faghaninia, R. Soklaski, J.A. Yan, C. Lo, L. Yang, Enhanced thermoelectric efficiency via orthogonal electrical and thermal conductances in phosphorene, Nano Lett. 14 (2014) 6393–6399. [15] J. Zhang, H.J. Liu, L. Cheng, J. Wei, J.H. Liang, D.D. Fan, J. Shi, X.F. Tang, Q.J. Zhang, Phosphorene nanoribbon as a promising candidate for thermoelectric applications, Sci. Rep. 4 (2014) 6452. [16] Q. Zeng, B. Sun, K. Du, W. Zhao, P. Yu, C. Zhu, J. Xia, Y. Chen, X. Cao, Q. Yan, Z. Shen, T. Yu, Y. Long, Y.K. Koh, Z. Liu, Highly anisotropic thermoelectric properties of black phosphorus crystals, 2d Mater. 6 (2019) 045009. [17] Y. Hong, J. Zhang, X. Huang, X.C. Zeng, Thermal conductivity of a two-dimensional phosphorene sheet: a comparative study with graphene, Nanoscale 7 (2015) 18716–18724. [18] A. Jain, A.J. McGaughey, Strongly anisotropic in-plane thermal transport in singlelayer black phosphorene, Sci. Rep. 5 (2015) 8501. [19] H. Jang, J.D. Wood, C.R. Ryder, M.C. Hersam, D.G. Cahill, Anisotropic thermal conductivity of exfoliated black phosphorus, Adv. Mater. 27 (2015) 8017–8022. [20] J.W. Jiang, Thermal conduction in single-layer black phosphorus: highly anisotropic? Nanotechnology 26 (2015) 055701. [21] S. Lee, F. Yang, J. Suh, S. Yang, Y. Lee, G. Li, H. Sung Choe, A. Suslu, Y. Chen, C. Ko, J. Park, K. Liu, J. Li, K. Hippalgaonkar, J.J. Urban, S. Tongay, J. Wu, Anisotropic inplane thermal conductivity of black phosphorus nanoribbons at temperatures higher than 100 K, Nat. Commun. 6 (2015) 8573. [22] T. Liang, P. Zhang, P. Yuan, S. Zhai, In-plane thermal transport in black phosphorene/graphene layered heterostructures: a molecular dynamics study, J. Chem. Soc. Faraday Trans. 20 (2018) 21151–21162.

Data availability The raw/processed data required to reproduce these findings cannot be shared at this time as the data also forms part of an ongoing study. CRediT authorship contribution statement Fan Zhu: Data curation, Formal analysis, Investigation, Methodology, Writing - original draft, Visualization. Hang Yin: Conceptualization, Formal analysis, Funding acquisition, Project administration, Resources, Writing - original draft. Ning Wei: Methodology, Writing - review & editing, Visualization. Jing Wan: Writing - review & editing, Visualization. Declaration of Competing Interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgements The authors acknowledge support from Shandong Provincial Natural Science Foundation, China (Grant No. ZR2017BA011) and the Scientific Research Foundation of Shandong Agricultural University (Grant No. 130/72130). The authors are also greatly thankful to the High Performance Computing Center of Shandong Agricultural University. References [1] A. Castellanos-Gomez, Black phosphorus: narrow gap, wide applications, J. Phys. Chem. Lett. 6 (2015) 4280–4291.

7

Materials Today Communications 22 (2020) 100814

F. Zhu, et al.

[40] K. Cai, L. Liu, J. Shi, Q.H. Qin, Self-assembly of a jammed black phosphorus nanoribbon on a fixed carbon nanotube, J. Phys. Chem. C 121 (2017) 10174–10181. [41] D. Pan, T.-C. Wang, C. Wang, W. Guo, Y. Yao, Self-assembled chiral phosphorus nanotubes from phosphorene: a molecular dynamics study, RSC Adv. 7 (2017) 24647–24651. [42] F. Hao, X. Liao, H. Xiao, X. Chen, Thermal conductivity of armchair black phosphorus nanotubes: a molecular dynamics study, Nanotechnology 27 (2016) 155703. [43] K. Cai, J. Wan, N. Wei, H. Cai, Q.H. Qin, Thermal stability of a free nanotube from single-layer black phosphorus, Nanotechnology 27 (2016) 235703. [44] S. Plimpton, Fast parallel algorithms for short-range molecular dynamics, J. Comput. Phys. 117 (1995) 1–19. [45] J.W. Jiang, Parametrization of stillinger-weber potential based on valence force field model: application to single-layer MoS2 and black phosphorus, Nanotechnology 26 (2015) 315706. [46] J.-W. Jiang, Y.-P. Zhou, Parameterization of Stillinger-Weber Potential for TwoDimensional Atomic Crystals, Handbook of Stillinger-Weber Potential Parameters for Two-Dimensional Atomic Crystals, IntechOpen, 2017. [47] W. Xu, L. Zhu, Y. Cai, G. Zhang, B. Li, Direction dependent thermal conductivity of monolayer phosphorene: Parameterization of Stillinger-Weber potential and molecular dynamics study, J. Appl. Phys. 117 (2015) 214308. [48] L.T. Kong, Phonon dispersion measured directly from molecular dynamics simulations, Comput. Phys. Commun. 182 (2011) 2201–2207. [49] J.M. Dickey, A. Paskin, Computer simulation of the lattice dynamics of solids, Phys. Rev. 188 (1969) 1407–1418. [50] P.K. Schelling, S.R. Phillpot, P. Keblinski, Comparison of atomic-level simulation methods for computing thermal conductivity, Phys. Rev. B 65 (2002). [51] G. Qin, Q.-B. Yan, Z. Qin, S.-Y. Yue, M. Hu, G. Su, Anisotropic intrinsic lattice thermal conductivity of phosphorene from first principles, J. Chem. Soc. Faraday Trans. 17 (2015) 4854–4858. [52] J.D. Gale, A.L. Rohl, The general utility lattice program (GULP), Mol. Simul. 29 (2003) 291–341. [53] B. Sa, Y.-L. Li, J. Qi, R. Ahuja, Z. Sun, Strain engineering for phosphorene: the potential application as a photocatalyst, J. Phys. Chem. C 118 (2014) 26560–26568. [54] R. Fei, L. Yang, Lattice vibrational modes and Raman scattering spectra of strained phosphorene, Appl. Phys. Lett. 105 (2014) 083120. [55] N. Wei, L. Xu, H.Q. Wang, J.C. Zheng, Strain engineering of thermal conductivity in graphene sheets and nanoribbons: a demonstration of magic flexibility, Nanotechnology 22 (2011) 105705. [56] M. Hu, X. Zhang, D. Poulikakos, Anomalous thermal response of silicene to uniaxial stretching, Phys. Rev. B 87 (2013). [57] A.S. Rodin, A. Carvalho, A.H. Castro Neto, Strain-induced gap modification in black phosphorus, Phys. Rev. Lett. 112 (2014) 176801. [58] J. Quereda, P. San-Jose, V. Parente, L. Vaquero-Garzon, A.J. Molina-Mendoza, N. Agrait, G. Rubio-Bollinger, F. Guinea, R. Roldan, A. Castellanos-Gomez, Strong modulation of optical properties in black phosphorus through strain-engineered rippling, Nano Lett. 16 (2016) 2931–2937.

[23] B. Liu, L. Bai, E.A. Korznikova, S.V. Dmitriev, A.W.-K. Law, K. Zhou, Thermal conductivity and tensile response of phosphorene nanosheets with vacancy defects, J. Phys. Chem. C 121 (2017) 13876–13887. [24] A. Maryam, G. Abbas, M. Rashid, A. Sattar, Directional mechanical and thermal properties of single-layer black phosphorus by classical molecular dynamics, Chinese Phys. B 27 (2018) 017401. [25] Z.-Y. Ong, Y. Cai, G. Zhang, Y.-W. Zhang, Strong thermal transport anisotropy and strain modulation in single-layer phosphorene, J. Phys. Chem. C 118 (2014) 25272–25277. [26] G. Qin, M. Hu, Thermal transport in phosphorene, Small 14 (2018) e1702465. [27] K. Xu, Z. Fan, J. Zhang, N. Wei, T. Ala-Nissila, Thermal transport properties of single-layer black phosphorus from extensive molecular dynamics simulations, Model. Simul. Mat. Sci. Eng. 26 (2018) 085001. [28] J. Zhang, H.J. Liu, L. Cheng, J. Wei, J.H. Liang, D.D. Fan, P.H. Jiang, J. Shi, Thermal conductivities of phosphorene allotropes from first-principles calculations: a comparative study, Sci. Rep. 7 (2017). [29] Y.-Y. Zhang, Q.-X. Pei, J.-W. Jiang, N. Wei, Y.-W. Zhang, Thermal conductivities of single- and multi-layer phosphorene: a molecular dynamics study, Nanoscale 8 (2016) 483–491. [30] Z. Luo, J. Maassen, Y. Deng, Y. Du, R.P. Garrelts, M.S. Lundstrom, P.D. Ye, X. Xu, Anisotropic in-plane thermal conductivity observed in few-layer black phosphorus, Nat. Commun. 6 (2015). [31] L. Zhu, G. Zhang, B. Li, Coexistence of size-dependent and size-independent thermal conductivities in phosphorene, Phys. Rev. B 90 (2014). [32] G. Qin, X. Zhang, S.-Y. Yue, Z. Qin, H. Wang, Y. Han, M. Hu, Resonant bonding driven giant phonon anharmonicity and low thermal conductivity of phosphorene, Phys. Rev. B 94 (2016). [33] Q.-X. Pei, X. Zhang, Z. Ding, Y.-Y. Zhang, Y.-W. Zhang, Thermal stability and thermal conductivity of phosphorene in phosphorene/graphene van der Waals heterostructures, J. Chem. Soc. Faraday Trans. 19 (2017) 17180–17186. [34] Y. Chen, Y. Zhang, K. Cai, J. Jiang, J.-C. Zheng, J. Zhao, N. Wei, Interfacial thermal conductance in graphene/black phosphorus heterogeneous structures, Carbon 117 (2017) 399–410. [35] X. Liu, J. Gao, G. Zhang, Y.-W. Zhang, Design of phosphorene/graphene heterojunctions for high and tunable interfacial thermal conductance, Nanoscale 10 (2018) 19854–19862. [36] Z. Zhao, H. Yin, K. Cai, W. Zhou, Mechanical stability of a nanotube from monolayer black phosphorus with the [110] direction as the tube’s circumference or generatrix, J. Chem. Soc. Faraday Trans. 20 (2018) 3465–3473. [37] X. Han, H.M. Stewart, S.A. Shevlin, C.R. Catlow, Z.X. Guo, Strain and orientation modulated bandgaps and effective masses of phosphorene nanoribbons, Nano Lett. 14 (2014) 4607–4614. [38] M.C. Watts, L. Picco, F.S. Russell-Pavier, P.L. Cullen, T.S. Miller, S.P. Bartus, O.D. Payton, N.T. Skipper, V. Tileli, C.A. Howard, Production of phosphorene nanoribbons, Nature 568 (2019) 216–220. [39] H.E. Lim, Y. Miyata, R. Kitaura, Y. Nishimura, Y. Nishimoto, S. Irle, J.H. Warner, H. Kataura, H. Shinohara, Growth of carbon nanotubes via twisted graphene nanoribbons, Nat. Commun. 4 (2013) 2548.

8