Resonant bonding driven giant phonon anharmonicity and low thermal conductivity of phosphorene

Two-dimensional (2D) phosphorene, which possesses fascinating physical and chemical properties distinctively different from other 2D materials, calls for a fundamental understanding of thermal transport properties for its rapidly growing applications in nano- and optoelectronics and thermoelectrics. However, even the basic phonon property, for example, the exact value of the lattice thermal conductivity ($\kappa$) of phosphorene reported in the literature, can differ unacceptably by one order of magnitude. More importantly, the fundamental physics underlying its unique properties such as strong phonon anharmonicity and unusual anisotropy remains largely unknown. In this paper, based on the analysis of electronic structure and lattice dynamics from first principles, we report that the giant phonon anharmonicity in phosphorene is associated with the soft transverse optical (TO) phonon modes and arises from the long-range interactions driven by the orbital governed resonant bonding. We also provide a microscopic picture connecting the anisotropic and low $\kappa$ of phosphorene to the giant directional phonon anharmonicity and long-range interactions, which are further traced back to the asymmetric resonant orbital occupations of electrons and characteristics of the hinge-like structure. The unambiguously low $\kappa$ of phosphorene obtained consistently by three independent ab initio methods confirms the phonon anharmonicity to a large extent and is expected to end the confusing huge deviations in previous studies. This work further pinpoints the necessity of including van der Waals interactions to accurately describe the interatomic interactions in phosphorene. We propose in 2D material that resonant bonding leads to low thermal conductivity, despite that it is originally found in three-dimensional (3D) thermoelectric and phase-change materials. Our study offers insights into phonon transport from the view of orbital states, which would be of great significance to the design of emerging phosphorene-based nanodevices.

Figure 1

Charge density difference (total density minus the sum of atomic densities) and three-dimensional (3D) phonon dispersion of phosphorene. (a) Top view and (b) side view of the change density difference (yellow, positive; cyan, negative) with zigzag ($x$) and armchair ($y$) directions indicated (isosurface level, 0.007). (c) Phonon dispersion and density of states (DOS) of phosphorene. The soft out-of-plane transverse optical (${\mathrm{TO}}_z$) phonon branch and the three acoustic phonon branches are labeled and plotted in different colors. Inset: The Brillouin zone (BZ) with high-symmetry $k$-points indicated. The 3D phonon dispersion in the whole BZ with color indicating the phonon mode contribution (percentage) to the total thermal conductivity along (d) zigzag and (e) armchair directions. Thermal conductivity and its anisotropy are mainly contributed (85%) by phonon modes below the gap (see Sec. 1 of the Supplemental Material [46]).

Figure 2

Analysis of phonon anharmonicity in phosphorene by first-principles calculations. (a–c) Potential energy change per atom (together with harmonic and anharmonic fitting) with respect to displacement of the P atom along the zigzag ($x$), armchair ($y$), and out-of-plane ($z$) directions. The anharmonic fitting is performed with the polynomial function up to the cubic term. Inset: Schematic of atomic displacement. (d) Convergence of Grüneisen parameters ($\gamma$) with respect to $r_{\mathrm{cutoff}}$ for zone-center ($\mathrm{\Gamma }$) optical phonon modes polarized along the $x, y$, and $z$ directions, respectively. (e–k) The orbital projected density of states (pDOS) of electrons ($s, p_x, p_y$, and $p_z$) for structures at equilibrium $[d=0$, (h)] and modulated by negative (e–g) and positive (i–k) displacements of the P atom along the $x, y$, and $z$ directions.

Figure 3

Characteristics of bonding and behavior of electrons in phosphorene. (a) Side view of the highest occupied molecular orbital (HOMO) ($p$-type $\sigma$ bond; blue, positive; red, negative), corresponding to the valence-band maximum (isosurface level, 0.03). (b) Absolute charge density change due to the displacement (0.02 Å along the $z$ direction) of the central nether P atom within a $5\times 5\times 1$ supercell, showing long-range interactions along the [110] direction. The plot is top view and atoms are marked on site (upside atoms are in dark purple). The unit is ${10}^{-5} e$ Å${}^{-3}$. (c) Top view of electron localization function (ELF) within a $2\times 2\times 1$ supercell.

Figure 4

Phonon transport properties of phosphorene calculated by three independent ab initio methods. (a) Lattice thermal conductivity ($\kappa$) of phosphorene at 300 K vs the interaction cutoff distance obtained by solving the phonon BTE based on the ALD method. Inset: $\kappa$ vs $Q$-grid ($N\times N$). The green (gray) ellipse indicates the cutoff distance that is finally chosen for the reported $\kappa$ values. (b) Extrapolated $\kappa$ by linearly fitting the thermal conductivity obtained from SED/EAIMD method with respect to system size. The error of the converged $\kappa$ is estimated by the standard error in the linear fitting. (c) $\kappa$ and the anisotropy at 300 K as a function of system size. The anisotropy is defined as ${\kappa }_{\mathrm{zigzag}}/{\kappa }_{\mathrm{armchair}}$. The $\kappa$ along the armchair direction calculated by NEAIMD is plotted for comparison. The error bar is determined based on the calculation of the temperature gradient ($\mathbf{\nabla }T$) and heat flux ($J$). (d) The comparison of mode level phonon lifetime between ALD/BTE and SED/EAIMD methods. (e) The comparison of $\kappa$ from our calculation with experiments and other reported results. Left: $\kappa$ of phosphorene from different methods by others and the current work [11, 16, 20, 21, 22, 23, 24, 25, 26]. Middle: Experimental $\kappa$ values of phosphorene films with different thickness [17, 18, 19]. Right: $\kappa$ of bulk BP from both experiment (with error bar) and theoretical calculation [16].