Thermal conductivities of phosphorene allotropes from first-principles calculations: a comparative study

Phosphorene has attracted tremendous interest recently due to its intriguing electronic properties. However, the thermal transport properties of phosphorene, especially for its allotropes, are still not well-understood. In this work, we calculate the thermal conductivities of five phosphorene allotropes (α-, β-, γ-, δ- and ζ-phase) by using phonon Boltzmann transport theory combined with first-principles calculations. It is found that the α-phosphorene exhibits considerable anisotropic thermal transport, while it is less obvious in the other four phosphorene allotropes. The highest thermal conductivity is found in the β-phosphorene, followed by the δ-, γ- and ζ-phase. The much lower thermal conductivity of the ζ-phase can be attributed to its relatively complex atomic configuration. It is expected that the rich thermal transport properties of phosphorene allotropes can have potential applications in the thermoelectrics and thermal management.

Single-layer black phosphorus, the so-called phosphorene, has emerged as a viable candidate in the field of two-dimensional atomic-layer materials due to its high carrier mobility and a thickness-dependent direct band gap. Besides phosphorene, several other stable layered phosphorus allotropes, with similar buckled or puckered honeycomb structure have been theoretically predicted by Zhu et al. 1,2 . It is suggested that layered phosphorus with distinct configuration possess different electronic properties, which can be further modulated by layer thickness, in-plane strain, and heterostructure assembling without energy penalty. Such tunable electronic properties of phosphorene allotropes are very beneficial for the optoelectronic and nanoelectronic applications. Although high performance field-effect transistors (FETs) and photovoltaic devices based on layered black phosphorus have been reported recently [3][4][5][6][7] , the anisotropic thermal transport properties of phosphorene may degrade the device reliability and performance since the low thermal conductivity along the armchair direction can lead to localized Joule heating in the confined system. On the other hand, the possibility to use phosphorene (α-phase) as thermoelectric material has been theoretically suggested [8][9][10] . For example, Liao et al. 10 showed that the power factor of phosphorene can reach as high as 70 μWcm −1 K −2 at appropriate carrier concentration. However, they predicted that the thermoelectric performance of α-phosphorene is poor at room temperature due to relatively higher lattice thermal conductivity, implying that suppressing thermal transport is an efficient way to enhance its thermoelectric efficiency. Although all-scale hierarchical architecturing and nanostructuring are efficient approaches to reduce the thermal conductivity of thermoelectric materials, the corresponding synthetic techniques are usually complicated and costly. As an alternative, high thermoelectric performance could be sought in materials with intrinsically low thermal conductivity 11 . It is thus quite necessary to investigate the thermal transport properties of various phosphorene allotropes, which could find potential applications in nanoelectronics and thermoelectric materials.
In this work, using phonon Boltzmann transport theory combined with first-principles calculations, we provide a comparative study on the thermal transport properties of five phosphorene allotropes. We demonstrate that thermal transport in the α-phosphorene exhibits strong orientation dependence, which is less obvious in the β-, γ-, δand ζ-phase. Moreover, we find that previously less studied ζ-phosphorene possesses much lower thermal conductivity, which is consistent with its relatively complex atomic configuration. Our theoretical work not only offers physical insight into the thermal transport in phosphorene allotropes, but also suggests that low thermal conductivity can be achieved in crystalline system constructed by light phosphorus atoms.

Computational Methods
Our theoretical calculations are performed by using a first-principles plane-wave pseudopotential formulation [12][13][14] as implemented in the Vienna ab-initio (VASP) 15 code. The exchange-correlation functional is in the form of Perdew-Burke-Ernzerhof (PBE) 16 with the generalized gradient approximation (GGA). For the structural   optimization, the energy convergence threshold is set to 1 × 10 −7 eV and the residual force acting on each atom is less than 10 −5 eV/Å. The cutoff energy for the plane-wave basis is set to be 500 eV, and uniform Monkhorst-Pack 17 k-mesh is applied to sample the Brillouin zone. To eliminate interactions between the phosphorene layer and its periodic images, we use a vacuum distance larger than 14 Å for the supercell geometry. The lattice thermal conductivity of phosphorene allotropes can be calculated by solving phonon Boltzmann transport equation (BTE) with an iterative self-consistent algorithm, as implemented in the so-called ShengBTE code [18][19][20] . In this approach, the thermal conductivity along the α direction can be calculated by: , Here → C q j , is the specific heat of the phonon mode with the wave vector → q and polarization j, , , is the corresponding phonon group velocity, and τ→ q j , is the self-consistent phonon relaxation time. Nq is the number of sampled q points in the Brillouin zone, and V is the volume of the unit cell (for low-dimensional system such as our studied phosphorene, however, the definition of a "volume" depends on the vacuum distance adopted, and we will come back to this point later). During the calculations of thermal conductivity, the only inputs are the harmonic (second-order) and anharmonic (third-order) interatomic force constants (IFCs) matrix, which can be extracted from first-principles calculations by using the finite displacement approach. To calculate the secondand third-order IFCs, a 6 × 6, 5 × 5, 5 × 4, 3 × 3, and 3 × 3 supercell with a uniform k-mesh of 3 × 3 × 1 is respectively used for the α-, β-, γ-, δand ζ-phosphorene. The phonon dispersion relation is obtained by using Phonopy package 21 with the harmonic IFCs as input. When dealing with the anharmonic ones, a cutoff distance of about 5.5 Å is employed, which has been confirmed by Jain et al. 22 to obtain converged thermal conductivity for both    Table 2. The Γ point group velocity of longitudinal acoustic phonons along the x-and y-direction for the five phosphorene allotropes.
It should be mentioned that the phonon BTE method has already been used to predict the lattice thermal conductivity of many bulk and low-dimensional structures such as transition-metal dichalcogenides, and the results are in good agreement with the experimental measurements [23][24][25] .

Results and Discussion
The crystal structures of phosphorene allotropes are schematically depicted in Fig. 1. Following the well-known notation for phosphorene allotropes 2 , we denote them as α-phosphorene ( Fig. 1(a)), β-phosphorene ( Fig. 1(b)), γ-phosphorene ( Fig. 1(c)), δ-phosphorene ( Fig. 1(d)), and ζ-phosphorene ( Fig. 1(e)), respectively. One can see that all the phases share the structural motif of threefold-coordinated P atoms and are characterized by a non-planar honeycomb atomic-layer. There are four atoms in the rectangular unit cell of αand γ-phosphorene, and eight atoms in δand ζ-phosphorene. For the isotropic β-phase, however, there are only two P atoms in the hexagonal unit cell. The optimum structural parameters for the five phosphorene allotropes are summarized in Table 1. It should be noted that our calculated results for the α-, β-, γand δ-phase are very close to previous first-principles calculations 2,26 . For the previously less studied ζ-phase, we see it has larger lattice constant and more atoms in the unit cell compared with those of other phases, suggesting that P atoms in the ζ-phosphorene are arranged in a more complicated way. Such behavior would set a crucial precondition for the ζ-phase to exhibit a lower thermal conductivity, as will be discussed in the following. To investigate the stability of the phosphorene phases, we have computed the energy difference (∆E) with respect to the most stable α-phosphorene by using the formula: where E is the total energy per atom of the phosphorene allotropes. As can be found from Table 1, the β-phase is almost as stable as the α-phase, while the γand δ-phosphorene have higher energies than that of the α-phase by 95 and 91 meV/atom, respectively. Although the energy difference between the ζand α-phosphorene is 135 meV/atom, we still expect that the ξ-phosphorene is a stable phase by considering the fact that ΔE is lower than the energy difference between the black and white phosphorus (160 meV/atom), which was also used as an reference to estimate the stability of nine new phosphorene polymorphs 27 . To further confirm the thermal  stability of the ζ-phosphorene, we have performed ab-initio molecular dynamics (MD). The system is modeled by a 4 × 4 × 1 supercell containing 128 atoms, and is simulated in a microcanonical ensemble for 1000 steps with a time step of 1.0 fs. Figure 2 shows the structural snap shots of the ζ-phosphorene at 300, 500, and 700 K. We see that the structure of ζ-phase has small fluctuations even at 700 K (additional MD with a longer simulation time gives similar results), which indicates that the ζ-phosphorene considered in our work is rather stable. Figure 3 plots the phonon dispersion relations of the five phosphorene allotropes. For the well-known α-, β-, γand δ-phase, our results are found to be consistent with recent works by using the finite displacement method 1, 2, 28, 29 . In addition, the calculated grüneisen parameters for the αand β-phosphorene (not shown here) agree well with those reported by Sun et al. 29 . All these results further confirm the reliability of our calculations. For the less studied ζ-phosphorene, we see there is no imaginary frequency in the phonon spectrum which suggests that it is also kinetically stable. As indicated in the figure, the relatively larger phonon gaps of α-(0.92 THz) and β-phosphorene (3.41 THz) suggests that their rates for three-phonon scattering are lower by considering the energy conservation. It is thus reasonable to expect that the αand β-phosphorene may have a relatively higher thermal conductivity, as will be discussed later. Moreover, we see that the phonon dispersion of α-phosphorene is highly asymmetric along the Γ-X and Γ-Y directions, which is caused by its puckered hinge-like structure 22,30 . In contrast, such behavior is less obvious for the other four phosphorene allotropes. In Table 2, we list the Γ point group velocities of longitudinal acoustic (LA) phonons for all the five phosphorene allotropes. We find that the group velocity of α-phosphorene along the x-direction is about two times of that along the y-direction, which means that the elastic constant will show significant orientation dependence. For the other four phosphorene allotropes, however, the differences between x-and y-direction are relatively small (it is identical for the β-phosphorene). Figure 4 shows the calculated lattice thermal conductivity (κ p ) of five phosphorene allotropes as a function of temperature in the range from 200 K to 800 K. As the value of κ p for two-dimensional system is closely related to the vacuum distance, we adopt the interlayer separation of the bulk counterparts from ref. 2, which is 5.30, 4.20, 4.21 and 5.47 Å for the α-, β-, γand δ-phase, respectively. For the less studied ζ-phase, we apply the same method as used in the ref. 2 to optimize the corresponding bulk structure, and get a interlayer separation of 4.89 Å. It can be seen that for all the five phosphorene allotropes, the κ p decrease with increasing temperature and roughly follow a T −1 dependence, indicating that the Umklapp process is the dominant phonon scattering mechanism in the temperature range considered. For the α-phosphorene, we see it indeed exhibits obvious anisotropic thermal transport. For example, the room temperature thermal conductivity along the x-direction is about 3.5 times higher than that along the y-direction. Such behavior is consistent with a recent work carried out by using first-principles calculations 22,28 and molecular dynamics simulations 31 , which confirm the reliability of our approach. We further find that the thermal conductivities of the other four phosphorene allotropes show less anisotropy, and their values are in the order of β-phase > δ-phase > γ-phase > ζ-phase. It should be noted that our calculated thermal conductivity of β-phosphorene is much close to those obtained by Zheng et al. 32 Table 3. Percentage contribution of different phonon branches to the room temperature thermal conductivity for the five phosphorene allotropes. but a bit different from that reported by Jain and McGaughey 22 . The discrepancy may be caused by the method to calculate the IFCs: Jain et al. used the density functional perturbation theory while we adopted the finite displacement technique. Meanwhile, the size of k-mesh and supercell used to calculate IFCs may also lead to such difference. Among the five investigated phosphorene allotropes, it is very hopeful to use β-phase to solve the thermal management issues in the phosphorene based FETs. In addition, a similar lattice structure and thermal conductivity to MoS 2 monolayer will make β-phosphorene a promising material to form superlattice structure with highly tunable electronic properties by controlling the layer thickness and stacking order. On the other hand, the thermal conductivity of ζ-phosphorene can be as small as 2.7 W/mK at 800 K, which is rare for crystalline system constructed by light atoms. Combined with its higher power factor of 0.0358 W/mK 2 along y-direction, the highest ZT value is predicted to be 2.1 for the p-type system, suggesting that ζ-phosphorene could be a promising high-performance thermoelectric material. It should be mentioned that the calculated thermal conductivity of α-phosphorene is much lower than that of graphene (2000~5000 W/mK). The reason is that the no-planar structure of α-phosphorene breaks the out-of-plane symmetry and promotes the phonon-phonon scattering of the out-of-plane acoustic (ZA) mode 22,30 . As can be seen from Table 3, the contribution of ZA branch to the thermal conductivity in the α-phosphorene is much smaller than that of graphene (76%) 34 . In addition, we see the contribution of ZA mode is smaller than that of transverse acoustic (TA) and longitudinal acoustic (LA) mode in the βand γ-phosphorene. Moreover, except for the β-phase, all the other four allotropes have relatively higher optical contribution, which may arise from the overlap of low frequency optical phonons with the longitudinal acoustic phonons (see Fig. 3). Such behavior will flatten the LA dispersion and induce suppression of thermal transport from LA phonons 35 . In this regard, neglecting the optical-phonon contributions will result in an underestimated thermal conductivity 9 . Surprisingly, the contribution from the optical branch is higher than 30% in the ζ-phosphorene, which can be attributed to its more complex atomic configuration. Similar behavior has also been found in the complex compounds such as SnSe and CoSb 3 skutterudite 36,37 .
To have a better understanding of the calculated thermal conductivity of the five phosphorene allotropes, we plot in Fig. 5 the normalized accumulative thermal conductivity at room temperature with respect to cutoff phonon mean free path (MFP). For all the phosphorene allotropes, one can see that heat is mainly carried by phonons with MFP in a broad range (in the orders of magnitude from 10 to 10 3 nm for the α-, β-, and ζ-phase, and from 10 to 10 4 nm for the γand δ-phase). To figure out which kind of phonons can have a significant effect, we give in Table 4 the MFP values corresponding to 50% κ P accumulation. It can be seen that such characteristic MFP for the βand δ-phosphorene have an order of magnitude of about 100 nm, which means that the lattice thermal conductivity could be effectively decreased if the sample size is smaller than this critical value. For the ζ-phosphorene, however, the major contribution comes from phonons with MFP of about 10~20 nm, which is similar to that found in silicene where phonons with MFP of 5~20 nm contribute more than 80% of the thermal conductivity 38 . It is interesting to note that the characteristic MFP of the γ-phase exhibits strong anisotropy, which offers extra flexibility to modulate its thermal conductivity along the y-direction.
To understand why the five phosphorene allotropes have quite different thermal conductivity, we first examine their group velocities of different phonon branches. As shown in Fig. 6, we see that the group velocities of optical phonons are comparable with that of acoustic phonons in the α-, γ-, δand ζ-phosphorene, which is consistent with the fact that optical phonons contribute considerably to the total thermal conductivity of these four phases.  Table 4. The MFP value (in unit of nm) corresponding to 50% κ P accumulation for the five phosphorene allotropes. For the acoustic phonons, we find that the group velocity of LA and ZA modes in the α-phase and LA and TA modes in the β-phase are relatively higher than those in the other three phases, which may result in a higher thermal conductivity of these two phases. We next focus on the phonon relaxation time (τ) of these phosphorene allotropes. As can be seen from Fig. 7, the β-phase has smaller relaxation time of optical phonon compared with those of the other phases. Moreover, the lowest relaxation time for acoustic phonons can be found in the ζ-phosphorene, especially for the TA and LA modes. As a result, the ζ-phase exhibits the lowest thermal conductivity among the five investigated phosphorene allotropes. To further estimate the phonon-phonon scattering, we plot in Fig. 8 the dimensionless total scattering phase space (the so-called P3 parameter). It is clear to find that the γand ζ-phosphorene possess larger scattering phase space. This means that there are more phase space allowing phonon-phonon scattering in the γand ζ-phosphorene, which can lead to a reduction of phonon relaxation time and thus the lattice thermal conductivity.

Conclusion
We demonstrate by phonon Boltzmann theory combined with first-principles calculations that thermal transport in the α-phosphorene exhibits strong orientation dependence, which is less obvious for the β-, γ-, δand ζ-phase. Moreover, we find that the previously less studied ζ-phosphorene exhibits considerably small thermal conductivity in a broad temperature range, which is very desirable for thermoelectric application. Detailed analysis of the characteristic MFP of these phosphorene allotropes suggests that the thermal conductivity of βand δ-phosphorene could be effectively modulated by controlling the sample size. In contrast to the general assumption, we find that the optical modes contribute significantly to the total thermal conductivity (except for the β-phosphorene) and thus cannot be ignored. On the experimental side, considering the fact that α-phosphorene can be synthesized using mechanical cleavage 3, 4, 39 and liquid-phase exfoliation 40,41 , it is possible that other four phosphorene allotropes can be obtained in a similar way 1, 2 . In addition, molecular beam epitaxy (MBE) and chemical vapor deposition (CVD) may be alternative approaches to synthesize these two-dimensional materials 1,42 .
With the rapid progress of fabrication techniques, it is reasonable to expect that thermoelectric devices and thermal management in nanoelectronics could be realized in the phosphorene allotropes.