Introduction

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 recently3,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 suggested8,9,10. For example, Liao et al.10 showed that the power factor of phosphorene can reach as high as 70 μWcm−1K−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 conductivity11. 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 formulation12,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-Pack17 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 code18,19,20. In this approach, the thermal conductivity along the α direction can be calculated by:

$${\kappa }_{\alpha }=\frac{1}{{N}_{q}V}\sum _{\overrightarrow{q},j}{C}_{\overrightarrow{q},j}{v}_{\overrightarrow{q},j,\alpha }^{2}{\tau }_{\overrightarrow{q},j}.$$
(1)

Here \({C}_{\overrightarrow{q},j}\) is the specific heat of the phonon mode with the wave vector \(\overrightarrow{q}\) and polarization j, \({v}_{\overrightarrow{q},j,\alpha }\) is the corresponding phonon group velocity, and \({\tau }_{\overrightarrow{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 second- and 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 package21 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 black (α-) and blue (β-) phosphorene. The Gaussian function with a scale parameter of 1.0 for broadening is used to enforce the conservation of energy in the three-phonon process. Well-converged 60 × 60 × 1, 99 × 99 × 1, 60 × 60 × 1, 50 × 50 × 1, and 50 × 50 × 1 q-meshes are respectively used for the α-, β-, γ-, δ- and ζ-phosphorene. 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 measurements23,24,25.

Results and Discussion

The crystal structures of phosphorene allotropes are schematically depicted in Fig. 1. Following the well-known notation for phosphorene allotropes2, 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 calculations2, 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 (\({\rm{\Delta }}E\)) with respect to the most stable α-phosphorene by using the formula:

$${\rm{\Delta }}E=E-E(\alpha -{\rm{phase}}),$$
(2)

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 polymorphs27. 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 1
figure 1

Top- and side-views of (a) α-phosphorene, (b) β-phosphorene, (c) γ-phosphorene, (d) δ-phosphorene, and (e) ζ-phosphorene. The coordinate axes (x, y, z) and lattice vectors (\(\overrightarrow{{a}_{1}}\), \(\overrightarrow{{a}_{2}}\)) are indicated.

Table 1 The optimized lattice constants (a 1, a 2), total energy (E), and energy relative to the α-phase (ΔE) of the five phosphorene allotropes.
Figure 2
figure 2

Structural snap shots of the ζ-phosphorene at (a) 300 K, (b) 500 K, and (c) 700 K during the molecular dynamics simulations. The corresponding side-views are shown in (d), (e), and (f), respectively.

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 method1, 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 structure22, 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 3
figure 3

The phonon dispersion relations of (a) α-phosphorene, (b) β-phosphorene, (c) γ-phosphorene, (d) δ-phosphorene, and (e) ζ-phosphorene. The highlighted areas are calculated phonon gaps with values of 0.92, 3.41, 0.35, 0.22, and 0.47 THz, respectively.

Table 2 The Γ point group velocity of longitudinal acoustic phonons along the x- and y-direction for the five phosphorene allotropes.

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 calculations22, 28 and molecular dynamics simulations31, 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 and Peng et al.33, but a bit different from that reported by Jain and McGaughey22. 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 MoS2 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/mK2 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.

Figure 4
figure 4

Calculated thermal conductivity of five phosphorene allotropes as a function of temperature.

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) mode22, 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 phonons35. In this regard, neglecting the optical-phonon contributions will result in an underestimated thermal conductivity9. 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 CoSb3 skutterudite36, 37.

Table 3 Percentage contribution of different phonon branches to the room temperature thermal conductivity for the five phosphorene allotropes.

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 103 nm for the α-, β-, and ζ-phase, and from 10 to 104 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 conductivity38. 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.

Figure 5
figure 5

The normalized accumulative thermal conductivity at room temperature as a function of cutoff phonon MFP for: (a) α-phosphorene, (b) β-phosphorene, (c) γ-phosphorene, (d) δ-phosphorene, and (e) ζ-phosphorene. The dashed line denotes the MFP value corresponding to 50% κ P accumulation.

Table 4 The MFP value (in unit of nm) corresponding to 50% κ P accumulation for the five phosphorene allotropes.

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. 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.

Figure 6
figure 6

Room temperature group velocities of different phonon modes as a function of frequency for: (a) α-phosphorene, (b) β-phosphorene, (c) γ-phosphorene, (d) δ-phosphorene, and (e) ζ-phosphorene.

Figure 7
figure 7

Room temperature relaxation time of different phonon modes as a function of frequency for: (a) α-phosphorene, (b) β-phosphorene, (c) γ-phosphorene, (d) δ-phosphorene, and (e) ζ-phosphorene.

Figure 8
figure 8

Total scattering phase space for: (a) α-phosphorene, (b) β-phosphorene, (c) γ-phosphorene, (d) δ-phosphorene, and (e) ζ-phosphorene in three-phonon scattering processes.

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 cleavage3, 4, 39 and liquid-phase exfoliation40, 41, it is possible that other four phosphorene allotropes can be obtained in a similar way1, 2. In addition, molecular beam epitaxy (MBE) and chemical vapor deposition (CVD) may be alternative approaches to synthesize these two-dimensional materials1, 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.