Fa
Zhang
ab,
Xiong
Zheng
b,
Huimin
Wang
c,
Liang
Ding
*a and
Guangzhao
Qin
*b
aState Key Laboratory of Robotics and System, Harbin Institute of Technology, Harbin, 150001, People's Republic of China. E-mail: liangding@hit.edu.cn
bState Key Laboratory of Advanced Design and Manufacturing for Vehicle Body, College of Mechanical and Vehicle Engineering, Hunan University, 410082, Changsha, P. R. China. E-mail: gzqin@hnu.edu.cn
cHunan Key Laboratory for Micro-Nano Energy Materials & Device and School of Physics and Optoelectronics, Xiangtan University, Xiangtan, 411105, Hunan, China
First published on 10th May 2022
Phosphorene, as a new type of two-dimensional (2D) semiconductor material, possesses unique physical and chemical properties, and thus has attracted widespread attention in recent years. With its increasing applications in nano/optoelectronics and thermoelectrics, a comprehensive study of its thermal transport properties is necessary. It has been concluded from previous studies that there exist vast differences and uncertainties in the theoretically predicted thermal conductivity, which is generally attributed to the selection of an XC functional. However, there is no comprehensive investigation on this issue. In this study, based on first-principles calculations using 12 different exchange–correlation (XC) functionals, the phonon transport properties of phosphorene are systematically studied by solving the Boltzmann transport equation (BTE). The results show noticeable differences in the phonon transport properties of phosphorene under different XC functionals. For instance, the thermal conductivity of phosphorene along the zigzag direction ranges from 0.51 to 30.48 W m−1 K−1, and that along the armchair direction ranges from 0.15 to 5.24 W m−1 K−1. Moreover, the anisotropy of thermal conductivity is between 3.4 and 6.9, which fundamentally originates from the special hinge-like structure. This study conducts an in-depth analysis of phonon transport to reveal the effect and mechanism of different XC functionals in predicting the thermal transport properties, which would provide a reference for future research on phosphorene and other novel materials.
The thermal transport properties play a vital role in lots of applications as mentioned above. Although the thermal conductivity of phosphorene and bulk phosphorus has been experimentally characterized well, the thickness of the phosphorene samples for measuring thermal conductivity is too large to reach the ideal state due to the limitation of synthesis technology.21,22 Thus, in addition to experimental measurements, researchers tend towards the theoretical prediction of the thermal conductivity of phosphorene and the exploration of phonon transport properties.23–26 The thermal conductivity of phosphorene is still ambiguous even if there are a number of studies in the literature. There are a few possible reasons for this ambiguity. First, different empirical potentials or force fields are fitted following different procedures and are used to describe the interatomic interaction in classical molecular dynamics (MD) simulations, which have a significant effect on the results. Second, lots of pseudopotentials and exchange–correlation (XC) functionals are employed in the first-principles calculations, which are constructed in different frames, such as the local density approximation (LDA) and the generalized gradient approximation (GGA). The different pseudopotentials and XC functionals have a significant effect on the evaluation of thermal conductivity. For instance, XC functional groups affect the optimization of the crystal structure and the lattice constant, which have a decisive influence on the lattice vibration and heat transfer performance.27 In previous studies, researchers have found a vast difference in the thermal conductivity of single-layer phosphorene. Specifically, the thermal conductivity of phosphorene along the zigzag direction is between 15.33 and 152.7 W m−1 K−1, and the thermal conductivity in the armchair direction is between 4.59 and 63.9 W m−1 K−1.28,29 Previous studies attributed the large uncertainty of thermal conductivity of phosphorene to the selection of XC functionals. Even though the selection of functional groups in the first-principles calculations is particularly essential for accurately predicting the thermal conductivity of phosphorene, there is no comprehensive investigation on this issue. It is well known that black phosphorene is entirely different from planar graphene and the flexural silicene crystal structure.27,30–32 2D phosphorene has a hinge-like structure and is folded along the armchair direction. The hinge-like structure leads to different atomic arrangement periods and different degrees of density in different directions, which results in the fantastic anisotropy in the thermal transport properties of phosphorene. Thus, the goal of this study is to systematically investigate the effects of different XC functionals on the thermal conductivity of phosphorene.
In this study, the thermal transport properties of phosphorene are quantitatively studied by solving the phonon Boltzmann transport equation (BTE) based on first-principles calculations using different exchange–correlation (XC) functionals. The results show that the anisotropy in the thermal transport process of phosphorene exhibits significant differences with different XC functionals. Moreover, the thermal conductivity in the zigzag direction ranges from 0.51 to 30.48 W m−1 K−1 and the thermal conductivity in the armchair direction is between 0.15 and 5.24 W m−1 K−1. An in-depth study on phonon behavior is conducted to analyze the mechanism underlying the anisotropy of thermal transport in phosphorene. The results achieved in this study provide a valuable reference for the accurate quantitative calculation of the thermal conductivity of phosphorene and future research on other novel materials, which is of considerable significance to nano-engineering applications.
Based on the full convergence test, the cut-off values of the kinetic energy of the wave function are set to 700 and 450 eV for PAW and ultrasoft (USPP_GGA and USPP_LDA) pseudopotentials, respectively. A Monkhorst–Pack k-mesh of 15 × 11 × 1 is used to sample the Brillouin Zone (BZ), and the energy convergence threshold is set to 10−8 eV. In order to prevent the interaction caused by periodic boundary conditions, a vacuum layer of 20 Å is used along the out-of-plane direction. To completely optimize the shape and volume of the unit cell, all the atoms have been relaxed until the maximum Hellmann–Feynman force acting on each atom is not higher than 10−8 eV Å−1.
Before calculating the second-order and third-order interatomic force constants, the unit cell needs to be expanded to predict the dynamics of the lattice accurately. By judging the convergence of phonon dispersion, the supercell size of phosphorene is chosen as 6 × 6 × 1. For the calculation of all IFCs in phosphorene, space group symmetry is used to reduce the calculation cost and the numerical noise of IFCs.45 The cut-off distance (rcutoff) is thoroughly tested in the anharmonic IFC calculations, and the specific results can be obtained from the analysis and discussion. A detailed discussion on convergence will be presented later based on Fig. 2 and the ESI.† The use of the Lagrangian multiplier method46,47 strengthens the translation and rotation invariance of IFCs. Based on the density functional perturbation theory (DFPT), the dielectric constant (given in Table 1) is obtained. κ is obtained by using an iterative process to solve the linearized phonon BTE based on force constants, which is implemented in the ShengBTE package.27,46
XC | Direction | Lattice constant | κ | Dielectric constant | MFP | Grüneisen parameter | Heat capacity (105 J m−3 K−1) | Phase space (10−3 a.u.) | A1 = A2 | A3 | B1 = B2 |
---|---|---|---|---|---|---|---|---|---|---|---|
(Å) | (W m−1 K−1) | (nm) | (Å) | ||||||||
LDA | Zigzag | 3.267 | 4.93 | 4.373 | 2.943 | 1.281 | 18.094 | 3.299 | 2.199 | 2.225 | 3.324 |
Armchair | 4.367 | 1.02 | 8.864 | 1.756 | |||||||
PBE | Zigzag | 3.298 | 14.89 | 3.971 | 11.538 | 0.716 | 17.044 | 3.642 | 2.220 | 2.259 | 3.545 |
Armchair | 4.625 | 2.42 | 5.261 | 4.419 | |||||||
rPBE | Zigzag | 3.313 | 11.81 | 3.834 | 7.976 | 0.391 | 16.569 | 3.763 | 2.228 | 2.273 | 3.66 |
Armchair | 4.753 | 2.31 | 4.652 | 3.955 | |||||||
PBEsol | Zigzag | 3.282 | 5.91 | 4.267 | 4.104 | 1.603 | 17.769 | 3.369 | 2.228 | 2.272 | 3.66 |
Armchair | 4.436 | 1.06 | 7.025 | 2.035 | |||||||
PW91 | Zigzag | 3.304 | 7.22 | 3.965 | 4.585 | 0.502 | 17.039 | 3.668 | 2.234 | 2.262 | 3.545 |
Armchair | 4.625 | 1.62 | 5.223 | 2.448 | |||||||
USPP_GGA | Zigzag | 3.301 | 0.51 | 3.951 | 0.298 | 2.251 | 17.084 | 3.677 | 2.220 | 2.260 | 3.548 |
Armchair | 4.626 | 0.15 | 5.249 | 0.192 | |||||||
USPP_LDA | Zigzag | 3.264 | 5.81 | 4.346 | 3.541 | 0.209 | 18.123 | 3.338 | 2.196 | 2.224 | 3.329 |
Armchair | 4.371 | 1.59 | 8.828 | 2.274 | |||||||
optB88 | Zigzag | 3.322 | 30.48 | 3.961 | 21.611 | 0.332 | 17.185 | 3.612 | 2.235 | 2.275 | 3.503 |
Armchair | 4.580 | 5.11 | 5.440 | 8.276 | |||||||
optB86b | Zigzag | 3.304 | 15.93 | 4.077 | 10.717 | 0.179 | 17.476 | 3.458 | 2.263 | 2.259 | 3.438 |
Armchair | 4.506 | 3.26 | 6.063 | 3.314 | |||||||
optPBE | Zigzag | 3.323 | 23.19 | 3.904 | 17.968 | 0.938 | 17.017 | 3.722 | 2.236 | 2.278 | 3.545 |
Armchair | 4.627 | 5.24 | 5.164 | 7.409 | |||||||
vdW_DF2 | Zigzag | 3.378 | 14.56 | 3.771 | 10.328 | 0.074 | 16.443 | 4.156 | 2.263 | 2.321 | 3.685 |
Armchair | 4.782 | 2.41 | 4.589 | 4.104 | |||||||
vdW_DF22 | Zigzag | 3.378 | 13.67 | 3.771 | 9.246 | 0.868 | 16.445 | 4.189 | 2.263 | 2.321 | 3.685 |
Armchair | 4.782 | 1.98 | 4.589 | 3.673 |
Fig. 1 Comparison of the phonon dispersions of 2D black phosphorene calculated using different XC functionals. The experimentally measured data of bulk black phosphorus26,27,48 are also plotted. The drawn phonon density of states (DOS) is calculated using the optB86b functional. Inset: the IBZ and the side view of phosphorene where arrows mark the vibration mode of the FA phonon branch. |
According to phonon spectrum analysis from Fig. 1, the phonon modes below the band gap play a decisive role in the thermal transport. Moreover, these modes dominate the anisotropy of heat transfer. Readers can see the detailed discussion presented in Fig. 3.
In the past few years, first-principles calculations combined with BTE methods have been recognized as good methods to predict the heat transport properties of materials. However, there are still some factors that can affect the final calculation results. One of the important factors is the cut-off radius of the third-order force constant and the convergence of the Q-mesh division; in order to obtain accurate thermal conductivity prediction results, they should be carefully judged. As shown in Fig. 2 and Fig. S1 (ESI†), under the condition that the phonon–phonon scattering caused by anharmonicity is fully respected, we have performed the close truncation and Q-grid convergence tests on all 12 XC functionals.
The results show that for all XC functionals, the thermal conductivity of black phosphorene has a good convergence. According to the judgment of convergence, the cut-off radius for different XC functionals is as follows: USPP_LDA (7.3 Å), optB86/LDA (7.7 Å), USPP_GGA/optPBE (8.0 Å), and other XC functionals (8.2 Å). Besides, for the convenience of comparison, all XC functional Q grids are selected as 100 × 100 × 1.
It can be seen from Fig. 3(a and b) that there are differences in the predicted values of thermal conductivity using different functionals. However, common conclusions on the properties of thermal transport in phosphorene can be achieved. For instance, the thermal conductivity of single-layer phosphorene is lower in both the zigzag direction and the armchair direction, both of which are two orders of magnitude lower than the thermal conductivity of graphene (3000–5000 W m−1 K−1).33,47,48 The reason lies in the fact that the FA phonon mode contributes less to the two thermal conductivities in black phosphorene. Thus, the scattering rate of the FA phonon branch in phosphorene is very high, resulting in a small contribution of FA to the thermal conductivity.
From Fig. 3(a), it can be seen that the transport of phosphorene has a strong anisotropy. Previous results showed that the thermal conductivity of phosphorene along the zigzag direction is larger than that along the armchair direction, and the ratio ranges from 2.2 to 5.5.27 The anisotropy ratio calculated using vdW_DF2, vdW_DF22, optB88 (considering the vdW interaction), and PBE (without considering the vdW interaction) is higher than the value interval previously reported. It can be seen that although the prediction of thermal conductivity is larger when considering the vdW interactions, it has little effect on anisotropy.
In addition to the anisotropy, Fig. 3(b) shows that the thermal conductivity calculated using the “opt” functionals (optB88, optB86b, and optPBE) and the “vdW” functionals (vdW_DF2 and vdW_DF22) is higher. The reason may lie in the inclusion of van der Waals (vdW) interactions in the XC functionals. Due to the special hinge structure of the monolayer phosphorene, lots of adjacent P atoms do not form bonds, but there are still some interactions. These effects can be captured using the XC functionals including vdW interactions, which have a significant impact on the properties of bulk BP and phosphorene. Thus, the vdW interaction will disrupt the formation of resonance bonds by changing the lattice constant, especially the distance between the non-covalent bond P atoms, leading to long-range interactions. With the above discussions, it is understood that the thermal conductivity calculated from the “opt” functionals and “vdW” functionals is higher under the same cut-off radius.
The anisotropy of thermal conductivity is attributed to the anisotropy of phonon group velocity determined by phonon dispersion.23,49,50 As shown in Fig. 1, the thermal conductivity along the armchair direction of black phosphorene is mainly contributed by the LA phonon mode. For the thermal conductivity along the zigzag direction, TA/FA participates in addition to the contribution of the LA phonon modes, and the contribution of FA is more significant than that along the armchair direction. The reason is that the unique hinge structure of phosphorene causes the “uneven” effect in the armchair direction, which leads to the difference in the phonon mode in the two directions, and finally causes the thermal conductivity anisotropy.
To observe the anisotropy in the heat transport process of phosphorene more intuitively, the modal contribution is further examined as shown in Fig. 4(a and b). The XC functional used here is the representative optB86b. Fig. 4(a) clearly shows that the acoustic phonon branch is the main contributor to the thermal conductivity of phosphorene. Fig. 4(b) shows that the phonon group velocity in the zigzag direction in the acoustic phonon branch is significantly higher than in the armchair direction. In the ESI,† we have given the relationship among the phonon frequency of the other 11 functionals, the thermal conductivity, and group velocity.
Fig. 4(c and d) show the influence of different XC functionals on thermal conductivity and group velocity. In Fig. 4(c), the thermal conductivity of USPP_GGA in both zigzag and armchair directions is abnormally low. This is consistent with the result reflected in Fig. 3(a). It is worth noting that although Fig. 1 shows that the phonon spectrum using USPP_GGA has no imaginary frequency indicating the thermodynamic stability, this does not mean that USPP_GGA can accurately predict the thermal conductivity of the materials. The reason may be that the features of USPP and GGA are not compatible. Therefore, we suggest that the use of USPP_GGA functional groups should be avoided when calculating the thermal transport properties of black phosphorene. From the group velocities shown in Fig. 4(d), it can be observed that the group velocities of the acoustic phonon branches of vdW_DF2 and vdW_DF22 in the armchair direction are smaller than those of all other functionals, and the zigzag direction is more significant compared to that in other functionals. This makes the thermal conductivity calculated from vdW_DF2 and vdW_DF22 in the armchair direction small, and the thermal conductivity in the zigzag direction is relatively large. Therefore, the thermal conductivity anisotropy ratio of these two XC functionals is higher than those of all other functionals.
The understanding of thermal conductivity can be deepened through the analysis of phonon lifetime. The phonon lifetime predicted by optPBE and optB88 is higher than those predicted using the rest of the XC functionals. This is the reason for the high thermal conductivity obtained using optPBE and optB88. As we all know, the thermal conductivity of micro–nano materials is affected by three factors, which are specific heat capacity (Cv), group velocity (v), and phonon lifetime (t).51 Among them, the specific heat capacity and group velocity are wholly affected by the harmonic characteristics, while the phonon lifetime is affected by the scattering phase space and scattering intensity. In Table 1 and Fig. 5(c), we can see that the scattering phase space is similar. Now it is imperative to consider the scattering intensity, that is, the phonon anharmonicity quantified by the Grüneisen parameter.
Fig. 5(b) shows the comparison of Grüneisen parameters after using different XC functionals. As can be seen in Fig. 5(b), the Grüneisen parameters from optB88 are relatively small at 300 K, which results in the thermal conductivity of black phosphorene from optB88 being significantly higher than those from other functionals. LDA, PBEsol, and USPP_GGA have larger parameter values, which cause a smaller thermal conductivity using these XC functionals. One thing that should be noted is that the thermal conductivity is derived from the coupling of three factors, that is to say, the smaller the Grüneisen parameter, the longer the phonon lifetime obtained using different functionals in phosphorene, but it has not necessarily resulted in higher thermal conductivity. Table 1 shows that the total Grüneisen parameter of optB86b is the smallest, but the thermal conductivity obtained at 300 K is not the largest.
Phosphorene has a hinge-like structure, as shown in Fig. 6(a), and its wrinkled structure can actually be seen as a deformation of a two-dimensional planar structure.27 Each P atom has five bonds formed. However, since the sp-hybridization in phosphorene is weak,53 only three p electrons can be used for bonding, which leads to unsaturated covalent bonding in phosphorene. Therefore, the true bond state in phosphorene is the hybridization or resonance between the different electronic configurations where the electrons occupy the p orbital. In fact, compared with the resonance bonding of the conventional rock salt structure, the resonance bonding in phosphorene exists in a weakened form.54 This is due to the structural deformation of the perfect rock salt structure of the phosphorene hinge-like structure. However, the weakened resonance bond still has a non-negligible effect on the phonon transport properties of phosphorene. The stronger the resonance bond is, the lower the thermal conductivity is.53 As shown in Fig. 6(a), due to the deformation, the first nearest neighbor distances in phosphorene are A1 and A2, and the second nearest neighbor distance is A3. As shown in Table 1, the optimized structure of USPP_LDA has the smallest first and second nearest neighbor distance (A1 = A2 = 2.196 Å, and A3 = 2.224 Å). Therefore, the resonant bonding is further enhanced in USPP_LDA. This structural deformation further enhances long-distance interaction, which makes the thermal conductivity calculated by the USPP_LDA XC functional become lower.
The armchair and zigzag directions show significantly different bonding characteristics. This is another explanation for the anisotropy of thermal conductivity. The electron (s, px, py, and pz) orbital projection density of states (pDOS) is shown in Fig. 6(b). From the calculated results, it can be seen that the bonding state near the VBM is mainly controlled by the pz orbitals hybridized with the px and py orbitals. The s orbital is mainly restricted to around 9 eV below the VBM, showing a weak hybridization relationship with the px, py and pz orbitals.
Fig. 6(c) shows the in-plane electron density of phosphorene. Anisotropic behavior can also be understood from the basic point of view of atomic bonds. In order to depict the probability of occurrence of an electron pair, we plot the ELF of phosphorene in Fig. 6(c).10,52 The ELF shows the position and size of the bond and lone pair electrons. The ELF ranges from 0 to 1, where 0 means no electrons, and 1 means the position with the greatest degree of electronic localization. The ELF of phosphorene along the zigzag direction is greater than 0.5, which means that the electrons are localized. For the armchair direction, the ELF is less than 0.5, which means that the electrons are delocalized. The ELFs along the two directions show great differences for phosphorene, which is the physical reason for the anisotropy of a series of elements with a hinge-like structure.
Fig. 7 The cumulative thermal conductivity of phosphorene along (a) zigzag and (b) armchair directions with respect to the phonon MFP for different XC functionals at 300 K. |
The size of the nanostructure is smaller than the feature size, which can effectively control the thermal conductivity. So, the measurement of the MFP is an essential consideration in the thermal design of the nanostructure. Fig. 7 plots the cumulative thermal conductivity (k) of the phonon free path in the zigzag direction and armchair direction at 300 K. Comparing Fig. 7(a and b), when the size is 1 × 104 nm, all the phosphorene κ values calculated using the XC functionals converged. On the whole, due to the anisotropy of black phosphorene, the MFP in the zigzag direction is more substantial, which causes the zigzag direction κ to increase more drastically with the MFP. In order to understand ballistics and phonon transport better, Table 1 lists the MFP when k is 50%. For most XC functionals, the MFP in the zigzag direction is between 2 and 10 nm, and the MFP in the armchair direction is between 2 and 4 nm. After comparison, the MFP calculated using optPBE, optB88 and PBE is higher than the average value. The MFP obtained by LDA functional calculation is much smaller than the average value.
Footnote |
† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d2ma00407k |
This journal is © The Royal Society of Chemistry 2022 |