Distinct Correlation between the Vibrational and Thermal Transport Properties of Group \textrm{VA} Monolayer Crystals

TL;DR

This study systematically investigates the vibrational and thermal transport properties of group VA monolayer crystals, revealing correlations and potential for tuning thermal conductivity through alloying for thermoelectric applications.

Contribution

It provides the first comprehensive first-principles analysis of thermal transport in group VA monolayers, highlighting correlations between vibrational and thermal properties and the effects of alloying.

Findings

Phosphorene has the highest thermal conductivity among studied materials.

Anisotropic thermal conductivity observed in all materials.

Alloying can significantly reduce or enhance thermal conductivity.

Abstract

The investigation of thermal transport properties of novel two dimensional materials is crucially important in order to assess their potential to be used in future technological applications, such as thermoelectric power generation. In this respect, lattice thermal transport properties of monolayer structures of the group \textrm{VA} elements (P, As, Sb, Bi, PAs, PSb, PBi, AsSb, AsBi, SbBi, P$_{3}$As$_{1}$, P$_{3}$Sb$_{1}$, P$_{1}$As$_{3}$, As$_{3}$Sb$_{1}$) with black phosphorus like puckered structure were systematically investigated by first principles calculations and an iterative solution of the Phonon Boltzmann transport equation. Phosphorene was found to have the highest lattice thermal conductivity, $\kappa$, due to its low average atomic mass and strong interatomic bonding character. As a matter of course, anisotropic $\kappa$ were obtained for all the considered materials, owing to anisotropy in phonon group velocities and scattering rates (relaxation times) calculated for these structures. However, the determined linear correlation between the anisotropy in $\kappa$ of P, As, and Sb is significant. The results corresponding to the studied compound structures clearly point out that thermal (electronic) conductivity of pristine monolayers might be suppressed (improved) by alloying them with the same group elements. For instance, the room temperature $\kappa$ of PBi along armchair direction was predicted as low as 1.5 Wm$^{-1}$K$^{-1}$, whereas that of P was predicted to be 21 Wm$^{-1}$K$^{-1}$. In spite of the apparent differences in structural and vibrational properties, we peculiarly revealed an intriguing correlation between the $\kappa$ of all the considered materials as $\kappa$=c$_{1}$ + c$_{2}$/$m^{2}$, in particular along zigzag direction.

Tables

Table 1. Table 1: The calculated room temperature κ 𝜅 \kappa of the monolayer structures (bold faces), together with the results reported in the literature. a 𝑎 a , b 𝑏 b , c 𝑐 c , d 𝑑 d , e 𝑒 e , f 𝑓 f , g 𝑔 g , and h ℎ h corresponds to the results extracted from Ref. Zheng et al., 2016 , Cheng et al., 2016 , Qin et al., 2015 , Zeraati et al., 2016 , Liu and Chang, 2015 , Jain and McGaughey, 2015 , Zhu et al., 2014 , and Qin et al., 2016 , respectively. Here, the values are rescaled with the out of plane lattice constant. 5.36 Å (corresponding to the interlayer distance in bulk black phosphorus) was used in our study to make a reliable comparison of monolayer structures with each other. Out of plane lattice constant values are not available for the results represented with bold italic face.

P		As		Sb		Bi
${\kappa }_{ZZ}$	${\kappa }_{AC}$	${\kappa }_{ZZ}$	${\kappa }_{AC}$	${\kappa }_{ZZ}$	${\kappa }_{AC}$
109.6	21.0	20.3	5.6	9.6	5.4	4.5	2.7
109.9a	23.5a	26.8a	5.7a	9.7a	4.7a	7.8b	6.0b
30.2c	13.6c	30.4d	7.8d
48.9e	27.8e
107.7f	35.3f
81.6g	23.8g
15.3h	4.6h

Equations

$$\frac{1}{{\theta }_D^3}=\frac{1}{2}\left(\frac{1}{{\theta }_{LA}^3}+\frac{1}{{\theta }_{TA}^3}\right).$$

$${\theta }_D=\frac{\hslash v_s}{k_B}{\left(\frac{4\pi n}{A}\right)}^{1/2},$$

$${\kappa }_s=A\frac{m{\theta }_D^3\delta n^{1/3}}{{\gamma }^{2}T}$$

$$\frac{m{\theta }_D^3\delta n^{1/3}}{\kappa T},$$

Full text

I Introduction

Triggered by the rediscovery of black phosphorus (BP), group VA materials, specifically monolayer arsenic (As) and antimony (Sb), have started to receive attention due to their peculiar thermal, optical and electronic properties Liu et al. (2014); Luo et al. (2015); Zhang et al. (2015); Çakır et al. (2014); Chaves et al. (2015); Çakır et al. (2015). In fact, this excitement among different researchers is inspired by the fabrication of first phosphorene (P monolayer) based transistor in 2014 Li et al. (2014). The interest not only spreads to the optical and electronic community but also to the research groups studying thermal properties Luo et al. (2015); Zeraati et al. (2016); Zheng et al. (2016); Zhu et al. (2014). From the optical and electronic points of view, phosphorene is appealing because of its high carrier mobility and direct band gap Churchill and Jarillo-Herrero (2014). This band gap is also found to vary from 0.3 eV to 2 eV progressively when its thickness is reduced from its bulk form to single-layer form. On the other hand, the relatively low thermal conductivity and its directional dependency due to the anisotropic lattice structure in particular for phosphorene (P monolayer) and arsenene (As monolayer) is of immense value for thermoelectric applications, providing higher electro-thermal conversion efficiency Luo et al. (2015); Wang et al. (2016a); Zeraati et al. (2016). Hence, the materials with such properties present a high potential to start a new line of research and development leading to completely new technologies.

Group VA elements have the electronic configuration of $s^{2}p^{3}$ in the outermost shell, where this $sp^{3}$ bonding causes a non-planarity, either a puckered or buckled structure Wang et al. (2015); Zheng et al. (2016); Peng et al. (2017a). Van der Waals forces hold the 2D layers where each atom in a layer is 3-fold coordinated with covalent bonds Warschauer (1963). The puckered structure of P, As, and Sb have been reported to be thermodynamically stable with high electron/hole mobility Zhang et al. (2016) and low thermal conductivity Hu et al. (2014) values, which is essential for high-performance innovative thermoelectric research. Therefore, it is very motivating to spend more intense theoretical effort on the nature of phonon transport in these materials in order to shed light on their lattice thermal transport properties. In addition, low thermal conductivity of the VA column inspires controlling the thermal transport properties by material engineering such as using the alloys of these elements. Hence, in this study, it is also an object of interest regarding the change in thermal conductivity when Group VA elements are considered as compound structures.

Therefore, we calculate the intrinsic lattice thermal conductivity of 21 stable single layer two dimensional black phosphorane like structures shown in Fig. 1: P, and As with $Pmna$ space group symmetry, Sb and Bi with $Pmn2_{1}$ space group symmetry, PAs with $Pma2$ space group symmetry, PAs, PSb, PBi, AsSb, AsBi, SbBi with $Pmn2_{1}$ and $P2/m$ space group symmetries, and P1As3, P3As1, P3Sb1, and As3Sb1 with $Pm$ space group symmetry.

II Method

The structural and dynamical properties of all the considered systems were predicted by first-principles calculations based on density functional-theory (DFT) and density-functional perturbation theory (DFPT), as implemented in the Vienna Ab-initio Simulation package (VASP) code Kresse and Hafner (1993); Wu et al. (2005); Nunes and Gonze (2001). The exchange-correlation interactions were treated using the generalized gradient approximation (GGA) within the Perdew-Burke-Ernzerhof (PBE) formulationPerdew et al. (1996a, b). The single electron wave functions were expanded in plane waves with kinetic energy cutoff of at least 450 eV. For the structure optimizations, the Brillouin-zone integrations were performed using a $\Gamma$-centered regular 21$\times$21$\times$1 $k$-point mesh within the Monkhorst-Pack schemeMonkhorst and Pack (1976). The convergence criterion for electronic and ionic relaxations were set as 10*-7* eV and 10*-3* eV/Å , respectively. In order to minimize the periodic interaction along the $z$-direction the vacuum space between the layers was taken at least 15 Å.

Lattice thermal transport properties were determined by the self-consistent solution of Peierls-Boltzmann transport equation Ziman (1960) as implemented in ShengBTE code Li et al. (2014”, 2012). The zeroth iteration solutions corresponding to Relaxation Time Approximation (RTA) were also predicted in order to clarify the influence of the quasi-momentum-conserving normal and the non-quasi-momentum-conserving Umklapp scattering processes on thermal transport results. The second-order inter-atomic force constants (IFCs) required by the ShengBTE were are obtained by PHONOPY code which extracts the proper force-constant file from the results predicted by DFPT as implemented in VASP. The third-order IFCs were derived from the VASP calculations of the structures generated by considering up to 10 next-nearest-neighbor interactions. Here, 8$\times$8$\times$1 and 5$\times$5$\times$1 supercell structures were used in the DFT calculations performed to predict the second- and third-order IFCs, respectively. In lattice thermal transport calculations at least 40$\times$40$\times$1 well converged $q$-grid and 5.36 Å Lindsay et al. (2014); Qin et al. (2014) out of plane lattice constant, $d$ were used for all the considered single layer materials.

III Results

The primitive unit cell of the considered four pristine (P, As, Sb and Bi), thirteen binary compound (PAs, PSb, PBi, AsSb, AsBi, SbBi) and six one-third compound (P3As1, P3Sb1, P1As3, As3Sb1) group VA monolayer crystals has a rectangular lattice with a four-atom basis and distinctive space group symmetries, namely $Pmna$ (53), $Pmn2_{1}$ (31), $Pma2$ (28), $P2/m$ (10), and $Pm$ (6), see Fig. 1. Therefore, except for the three acoustic modes involving the in-phase vibrations of atoms in the unit cell, out-of-phase vibrations of the atoms give rise to nine optical modes for all the structures. The calculated phonon dispersion curves, corresponding to these modes, are depicted in Supplementary Materials (SMs), as Figures 2, 3, and 4, respectively, together with the group velocity of each phonon mode as a line color. As seen in these figures, no negative frequencies, indicating the dynamical instability at T=0 K, were found for the most of the considered structures. However, appreciable negative frequencies for the PSb, PBi, AsSb, AsBi, and SbBi structures with $Pma2$ space group symmetry were obtained, meaning that these structures are not dynamically stable. Also, few THz negative frequencies for P1Sb3 and As1Sb3 with $Pm$ space group symmetry were found out and these structures were not considered in thermal transport calculations. The in-plane longitudinal acoustic (LA) and transverse acoustic (TA) modes of all the materials have a linear dependence as $q$ goes zero (i.e. long wave length limit), whereas the out-of-plane acoustic or flexural modes (ZA) have a quadratic dependence in the long wavelength limit as commonly observed in two-dimensional (2D) materials such as graphene. The considered compounds display similar dispersion curves but with distinctive Debye frequencies, resulting in different phonon group velocity ($v_{qi}=\frac{\partial\omega_{qi}}{\partial q}$) values. Fig. 2(a) shows the calculated average of group velocities of three acoustic modes along zigzag (ZZ) and armchair (AC) directions (corresponding to reduced reciprocal coordinates, $q_{x}$ = 0.05, $q_{y}$ = 0.00 for ZZ, and $q_{x}$ = 0.00, $q_{y}$ = 0.05 for AC). A simple correlation between the average of group velocities of three acoustic modes, $v$ and the average effective mass of elemental configuration of the unit cell, $m$ was established as $v=c_{1}+c_{2}/m$, where $c_{1}$ and $c_{2}$ are constants. Our simple expression is consistent with the fact that the lower the average mass, the higher the group velocity is.

The maximum frequency ($\omega_{m}$) values at the $\Gamma$-point changes as $\omega_{m,\mathrm{P}}$ $>$ $\omega_{m,\mathrm{As}}$ $>$ $\omega_{m,\mathrm{Sb}}$ $>$ $\omega_{m,\mathrm{Bi}}$ and $\omega_{m,\mathrm{PAs}}$ $>$ $\omega_{m,\mathrm{P_{3}As_{1}}}$ $>$ $\omega_{m,\mathrm{PAs}}$, $\ldots$, $>$ $\omega_{m,\mathrm{SbBi}}$ mainly due to the bond stiffness, different atomic masses for pristine monolayers and average atomic mass of the unit cell for the alloy monolayers. When the bond stiffness becomes weaker (i.e. interatomic force constant becomes smaller) and the mass of atom becomes larger, the phonon frequencies shift to lower energies. Thus, the phonon branches become less disperse with exactly the same order, and $v_{qi}$ of phonon modes markedly decrease with increasing the mass of constituent atom(s). In this context, the Debye temperature, calculated from the contributions of individual acoustic modes (i.e. $\theta_{i}=\hbar\omega_{i,max}/k_{B}$, where $i$=$LA$ and $TA$) is given as Peng et al. (2016)

[TABLE]

The calculated Debye temperature ($\theta_{D}$) values for all the considered monolayer materials are listed in Table I in SMs. While $\theta_{D}$ is 206 K for phosphorene, it becomes 50 K for in Bi. For 2D crystals, Debye temperature can also be defined as followsKittel (1986),

[TABLE]

where $v_{s}$ is the sound velocity, $n$ is the number of atoms in the unit-cell, and $A$ is the area of the unit-cell Peng et al. (2017b). Based on this formula we calculated the square root of the unit cell area per number of atoms, $\delta^{2}=A/n$, times Debye temperature (calculated from phonon dispersion via Eq. (1)) with respect to inverse average atomic mass of the crystals (i.e. total atomic mass over number of atoms in the unit cell, $m/n$). Quite surprisingly, the predicted Debye temperatures from our calculations is correlated with lattice constants and inverse of the average atomic mass as expected from Eq. (2), see Fig. 2 (b). Note that, the sound velocity in the plot is considered as inversely proportional with average mass as it is in anharmonic crystals and our results shown in Fig. 2 (a) clearly proves those relations. The obtained perfect $m^{-1}$ correlation in both $v$ and $\theta_{D}$ indicates that these systems might be quite anharmonic.

Phonon dispersion curves shown in SMs are not symmetric along the $\Gamma$-X-S and $\Gamma$-Y-S, giving rise to anisotropic phonon group velocities along the zigzag and armchair directions. In fact, all the considered monolayers resemble the orthogonal symmetry resulting in distinctive interatomic bonding structure along the ZZ and AC directions, that induces anisotropic vibrational properties. Interestingly, the dispersion of acoustic modes along the ZZ and AC directions for Sb and Bi structures is not so different as compared to P and As structures, which would give rise to a smaller difference in thermal conductivity along the ZZ and AC directions. The predicted group velocities of LA modes of phosphorene along the ZZ and AC directions (see Figure 2 in SMs are in a good agreement with the previously reported first principles results (8.6 and 4.5 km/s) Zhu et al. (2014) and experimental measurements (9.6 and 4.6 km/s) Fujii et al. (1982) for bulk black phosphorus. The acoustic modes travel at higher speed along the ZZ direction since all monolayers are stiffer along the ZZ direction as compared to the AC direction. Note that $v_{qi}$ of major heat carrier phonon branches, i.e. acoustic modes, of Sb and Bi both along the ZZ and AC directions near the zone center are about two times smaller than those of black phosphorus as seen in Fig. 2(a). This notable difference in group velocity values of acoustic branches are fairly associated with the difference in $\kappa$ values of these crystals.

On the other hand, the presence of a frequency gap between the low frequency and high frequency optical branches of the most of the considered structure might also be effective on the thermal transport properties of these materials. In fact, in the former studies, the lower thermal conductivity of puckered phosphorus as compared to buckled phosphorus is correlated with the presence of a larger a-o gap in buckled oneZheng et al. (2016). Also, the significant effect of acoustic-optical phonon gaps on the $a+a$ $\leftrightarrow$ $o$ scattering channels was shown for different low dimensional materials Zheng et al. (2016); Lindsay et al. (2013, 2012); Jain and McGaughey (2014). These above-mentioned three phonon processes may also affect the temperature dependence of $\kappa$ since as temperature increases, number of optical phonons increase thereby enhancing three phonon processes involving optical phonons. There is a notable differences in phonon dispersions of compound monolayers with the same concentration but different symmetry, see, for example, Figure 3 in SMs for the phonon dispersion curves for PSb with symmetries $Pmn2_{1}$ and $P2/m$. Therefore, we can expect distinctive phonon scattering mechanisms in the compounds with different symmetries.

Subsequent to phonon dispersion analysis, we first calculated the lattice thermal conductivity of graphene in order to test our methodology and we predicted it as 3290 Wm*-1K-1* which is in quite good agreement with the literature Balandin (2011); Haskins et al. (2011); Nika and Balandin (2017, 2012) (depicted in SMs as Figure 5). Afterwards, we systematically investigated the lattice thermal transport properties of four pristine, thirteen binary and four one-third group VA monolayer structures. The self-consistent solution of Peierls-BTE for the lattice thermal conductivity ($\kappa$) of pristine monolayers as a function of temperature in the range from 100 K to 600 K are depicted with colored symbols in Fig. 3 (a) and (b). $\kappa$ decreases along both ZZ and AC directions when moving from P to Bi. Phosphorene has the highest $\kappa$ due to its highest $\theta_{D}$ temperature. We found that all the monolayers exhibit obvious anisotropic thermal transport, which is attributed to the orientation dependent group velocities spotted in Fig. 2 (a) and phonon relaxation times. For instance, the room temperature $\kappa$ along the ZZ direction is 5.21, 3.63, 1.78, and 2.05 times larger than that along the AC direction for P, As, Sb, and Bi, respectively. Interestingly, this ratio, which is nearly the same for Sb and Bi, changes linearly with atomic mass for P, As, and Sb. This behavior is consistent with the recent theoretical works Zheng et al. (2016); Qin et al. (2015); Liu and Chang (2015); Jain and McGaughey (2015); Zhu et al. (2014) performed using first-principles calculations, see Table 1.

Previously, the thermal conductivity of few layers black phosphorus films with a 15 nm thickness was measured as 40 and 20 Wm*-1K-1* and film with a 9.5 nm thickness was measured Luo et al. (2015) as 20 and 10 Wm*-1K-1* along ZZ and AC directions, respectively. In addition, the thermal conductivity of a single crystal bulk black phosphorus was measured Wang et al. (2016b) as 34 and 17 Wm*-1K-1* along the ZZ and AC directions, respectively. Even though we cannot directly compare the $\kappa$ of a monolayer material with that of its few-layer and bulk form, our results for $\kappa$ of phosphorene along both directions are in fair accordance with the previously reported experimental data for bulk and few-layer systems in terms of order of magnitude. Our calculated $\kappa$ of phosphorene is probably higher than that of phosphorene film with a thickness of 9.5 nm due to the strong interlayer interaction in the measured multilayer samples. The overlap of interlayer wavefunctions rather than a pure van der Waals interactionShulenburger et al. (2015); Hu et al. (2016) leads to an interlayer interaction much stronger than graphene in this crystal. In addition, the thermal conductivity of phosphorene was calculated by several groups. Interestingly, there is a deviation between the calculated $\kappa$ values as seen in Table 1. This can partially be attributed to different computational methods and computational parameters used in these studies. The calculated values for P, As, and Sb structures in this study agree very well with the values reported in a recent studyZheng et al. (2016).

Our calculations revealed that there is a distinct correlation between the $\kappa$ of pristine P, As, Sb, and Bi structures and average atomic mass ($m$). We found a simple correlation between $m$ and $\kappa(m)$ such that $\kappa(m)$ = $c_{1}$+$c_{2}$/$m^{2}$, where $c_{1}$ and $c_{2}$ are constants, see Fig. 4(a). Not surprisingly, the thermal conductivity values predicted by Zheng et al. also fits this function very well. In addition to this surprising correlation, we predicted that $\kappa$ of all the pristine monolayers follow nearly the same trend with temperature as seen in Fig. 3. Therefore, we can certainly claim that the decrease in the $\kappa$ of the pristine structures (i.e. P, As, Sb, and Bi) at temperatures above 100 K follows a $m^{-2}$ behavior except the slight deviation in $\kappa$ of Sb along the AC direction. In fact, the space group symmetry of P and As different from that of Sb and Bi due to the different bond bending angles, $A_{1}$ and $A_{2}$ shown in Fig. 1 and thus one can naturally expect a distinctive thermal transport properties for these two group. Consequently, the marked correlation revealed as a result of our calculations is quite remarkable.

Previously, Slack suggested that the thermal conductivity ($\kappa_{s}$) is determined by four factors, including (1) average atomic mass, (2) interatomic bonding, (3) crystal structure and (4) size of anharmonicitySlack (1973). The first three factors determine harmonic properties. Three phonon scattering dominated $\kappa_{s}$ can be given as

[TABLE]

where $\delta^{2}$ is the unit-cell area per atom, $n$ is the number of atoms in the unit cell (which is equal to four for all systems), $A$ is a constant given in terms of average Grüneisen parameter of modes Peng et al. (2016), $\gamma$, and $m$ is the average atomic mass of the unit-cell. We used this equation to predict size of anharmonicity in these systems. First of all, this simple relation neglects the contribution of the optical modes to the thermal conductivity. This is a reasonable approximation since the contribution of the optical modes to thermal conductivity does not exceeds 10% for the most of the considered crystals as shown in SMs. Inspired from the above relation, we calculated the $\gamma^{2}/A$ as,

[TABLE]

at 300 K. Here, we used the $\kappa$ values shown in Fig. 4. Then, the results were normalized by the calculated $\kappa$ of P. Within the Slack’s approximation, the $\gamma^{2}/A$ term is directly related with the anharmonicity of the crystal (the larger the $\gamma$, the larger the anharmonic effect is). Therefore, our calculations clearly show that there is a notable difference between the anharmonic effects along the ZZ and AC directions. While there is a steep increase in the ZZ direction, the values along the AC direction are very close to each other, and the calculated value for Sb is even smaller than that of P.

Figure 5 shows the calculated lattice thermal transport properties of thirteen binary, (a)-(b), and four one-third, (c)-(d), group VA monolayer structures. Not surprisingly, all the considered compound monolayer crystals exhibit obvious anisotropic thermal transport, which is mainly attributed to the orientation dependent group velocities spotted in Fig. 2(a). The compound monolayers containing P has larger thermal conductivity as compared with the compounds without P atom. The monolayers with Bi have the lowest thermal conductivities for all temperature values considered here. For instance, $\kappa$ of PBi with $P2/m$ space group was found to be 1.5 Wm*-1K-1*. The calculated $\kappa$ of low symmetry compound systems, as low as $\kappa$ of well known 2D thermoelectric materials including single layer SnSe (0.46-0.70 Wm*-1K-1* at 300 K)Zhao et al. (2014), clearly shows that the thermal conductivity of pristine systems can be suppressed up to order of magnitude by alloying in a controllable manner. Indeed, random elemental distribution in alloy systems results significant reduction in phonon mean free path of long wavelength heat carrier phononsZhang and Li (2010); Zhou et al. (2016). Therefore, it can be reasonably proposed that $\kappa$ of black phosphorus can be tuned (reduced to enhance its thermoelectric potential) by alloying it with Sb, Bi, and also As doping even with very low doping concentrations Sevik et al. (2011); Haskins et al. (2011), that will lead an opportunity to engineer these materials. This is because of the fact that substitutional As, Sb or Bi doping improves (suppresses) the electronic (thermal) transport by increasing (decreasing) the density of states around the Fermi level (band gap), see Figure 6, 7, and 8 in SMs. In addition, most of the monolayer structures considered here exhibit multiple extrema in their valence and conduction bands. These extrema in the valence and conduction bands are nearly energetically degenerate with the energy difference less than 0.15 eV. The edges of the valence bands yields several hole pockets. These type of electronic features in electronic band structure may give rise to large Seebeck coefficients in these monolayer structures. For the low band gap systems, bipolar conduction may emerge at high temperatures, which lowers the Seebeck coefficient. Since the considered monolayer structures have band gap values within range of 0.9-0.09 eV, we can maximize thermoelectric efficiency by selecting appropriate monolayers. Combination of the ultra-low $\kappa$ and high electrical conductivity makes these monolayers promising candidates for thermoelectric applications. In our calculations, we did not include spin-orbit coupling (SOC). In fact, SOC becomes more effective for monolayer structures with heavy atoms, such as Bi. It was shown that SOC has a non-negligible effect on the electronic properties of Bi monolayerZhang et al. (2017). However, we do not expect a significant impact of SOC on phonons and thermal conductivity values. For the calculations of conductivity and Seebeck coefficient, SOC should be taken into account to get accurate results.

The correlation between the $\kappa$ of pristine P, As, Sb, and Bi monolayers and their compound structures with average atomic mass was also searched for all the considered materials. Here, the $\kappa$ for a compound structure was calculated as average over the $\kappa$ values corresponding to all the considered symmetries of an individual concentration. As seen in Fig. 6(a), the $\kappa(m)$ = $c_{1}$+$c_{2}$/$m^{2}$ curve fits quite well to the calculated $\kappa$ along the ZZ direction of all the considered materials. However, the calculated $\kappa$ along the AC direction of P3Sb1, PSb, As3Sb1, and Sb do not follow the same relation as good as along the ZZ direction as seen in Fig. 6(b). Nevertheless, when these results considered with the identical temperature dependency of $\kappa$ of all the materials represented in Fig.5, it can be clearly stated that the $\kappa$ of compound group VA monolayer structures can be definitely estimated with $\kappa(m)$ = $c_{1}$+$c_{2}$/$m^{2}$ relation, above 100 K, with less than 10% error, except 36% and 42% underestimation for $\kappa$ of Sb and PSb, and 128% and %89 underestimation for $\kappa$ of P3Sb1 and As3Sb1. Despite the apparent differences, such as group velocities and optical phonon gap values, in phonon dispersion of these materials, the revealed simple trend is quite striking.

The calculated normalized $\gamma^{2}/A$ values are also calculated for the compound structures. As observed for the pristine monolayers, there is no obvious trend with average atomic mass as seen in Fig. 7(a). However, it can be partially asserted that the normalized conductivity along the ZZ direction increases with increasing average mass, meaning that the effect of anharmonicity becomes more and more prominent. However, the values along the AC direction scatter without an explicit trend. Interestingly, the change in the lattice parameters along the ZZ, $a_{0}$ and AC, $b_{0}$ directions, is in line with these results as seen Fig.8. This partially indicates that the change in some structural parameters, such as $A_{1}$, $A_{2}$ and $A_{3}$ angles depicted in Table V in SMs, due to the change in electronic-ionic bonding nature (i.e P-P is different from Sb-Sb bond), leads to a distinctive direction dependent phonon group velocity and scattering rate values in these monolayer structures. This interpretation is also supported by the mean free path, $l$, dependent thermal conductivities shown in Figure 9, 10, 11 in SMs. For instance, the mean free path of Sb and PSb along the AC direction considerably longer than that along the ZZ direction. Conversely, the $l$ of of P, Bi, and PAs ($Pma2$) along the ZZ direction are rather larger.

As mentioned above, Slack suggested that the thermal conductivity is determined by four factorsSlack (1973). In this regard, in order to differentiate the effect of purely interatomic bonding and size of anharmonicity on the thermal conductivity of the crystals, we recalculated the thermal conductivity of all the considered monolayers by equating the mass of all the constituent elements to mass of Phosphor (called as $\kappa_{mp}$). Afterwards, we calculated the ratio of the actual conductivity, $\kappa$, and the one obtained with mass of phosphor, $\kappa_{mp}$ to the thermal conductivity of P, $\kappa_{P}$, as seen in Figure 7(b). The deviation of $\kappa_{mp}$/$\kappa_{P}$ from unity depicts the effect of the interatomic bonding and anharmonic interactions on the thermal conductivity. In other words, it shows difference in the strength of interatomic bonds, which results in distinctive vibrational properties and thus anharmonic nature, of all the considered materials. The results for the ZZ direction shows that these two effects dominate more and more the lattice thermal transport as the average mass of unit cell decreases. However, such correlation for the AC direction is not as strong as the ZZ direction. These results clearly point out the dependence of anharmonicity on the crystal direction.

The zeroth iteration solutions of Peierls-BTE corresponding to RTA were also calculated. The former is able to capture simultaneous interaction of all phonon modes, giving rise to collective phonon excitations being responsible for heat transfer. The ratio of the self-consistent iterative solution to RTA solution for both ZZ and AC directions are depicted in Figure 9, 10, and 11 in SMs, respectively. The difference between the self-consistent and RTA results is greater than 10% and independent of temperature when temperature is slightly greater than the $\theta_{D}$ temperature of the monolayers. This clearly shows that the normal three-phonon processes play an important role in most of the compound crystals investigated in this study. For the 3D crystals with strong Umklapp scattering, the room temperature $\kappa$ calculated with the iterative solution is typically 10% greater than RTA solution Ward et al. (2009) and difference between two methods becomes considerably large for 2D materials. The RTA always severely underestimates thermal conductivity when the normal scattering (being dominant at low T) is strong. $\kappa$ of all the structures is almost inversely proportional to temperature ($\kappa\sim T^{-1}$), indicating that phonon scattering mechanism is dominated by the Umklapp process. We found that the dominant scattering mechanisms along ZZ and AC direction is notably different from each other in particular for the structures with P atom. For instance the calculated ratio for the pristine P and As structures along ZZ direction is larger than 1.8, however, the one obtained for the AC direction is smaller than 1.2. The calculated results for pristine Sb explicitly indicate that Umklapp scattering in this structure is considerably stronger than the other pristine ones in particular along the AC direction. This may explain the slight deviation in $\kappa$ of Sb from $\kappa(m)$ = $c_{1}$+$c_{2}$/$m^{2}$ relation. The iterative solution adds more than 10% to $\kappa$ of most of the crystals considered in our study, in the both ZZ and AC directions.

The ratio of thermal conductivity regarding to the three acoustic and three optical modes to total thermal conductivity, $\kappa_{i}/\kappa$ (where $i$ = ZA, LA, TA, B1u, B1g, and B3g and $\kappa$ is total thermal conductivity) was determined for all the materials and the results are represented in Figure 12, 13, and 14 in SMs. The relative contribution of both acoustic and lowest three optical phonons are almost temperature independent for all materials when the temperature is slightly larger than the $\theta_{D}$ temperature computed using Eq. (1). The total contribution of acoustic branches to $\kappa$ is significantly higher than the contribution of optical modes. However, in some P based materials, the lowest three optical modes, B1u and B1g have a notable contribution as seen in $\kappa_{AC}$ of phosphorene. Also, the contribution of low lying optical mode to conductivity along AC direction is greater than 10% for most of the P and As based systems (P, As, PAs, AsSb, P3As1, P1As3, and As3Sb1). Due to their non-planar structure, monolayers of group VA elements have significantly lower thermal conductivity as compared to graphene, 2000-5000 W/mK Balandin (2011); Haskins et al. (2011). This can be partially attributed to the promoted phonon-phonon scattering of the out-of-plane acoustic (ZA) modes. Due to reflection symmetry perpendicular to in-plane and hexagonal symmetry, three-phonon scattering involving ZA mode (anharmonic scattering of flexural phonons) process is significantly restricted in graphene Nika and Balandin (2017, 2012). However such a symmetry is broken in group VA monolayers and thus the anharmonic scattering is enhanced in these systems and the contribution of ZA mode not as much as in graphene as clearly seen in mode contribution picture of all the considered monolayers.

The ZA mode behaves differently along the ZZ and AC directions. While ZA mode provides a significant contribution to total thermal conductivity along the ZZ direction for the most of the monolayer structures, its contribution along the AC direction is usually smaller than 10% of the total conductivity. This can be explained in terms of distinct structural properties and strength of Umklapp scattering along the ZZ and AC directions. Phosphorene and other group VA monolayers were found to show superior structural flexibility along the AC direction allowing it to have large curvatures. This mechanical flexibility is the origin of less dispersive ZA, LA and TA modes along that direction. The quadratic behavior of the ZA mode is more apparent along this direction with a smaller group velocity. In graphene, the contribution of ZA mode to $\kappa$ is 80% of total conductivity. In general, the contribution of the ZA mode to $\kappa$ is smaller than 50% in present systems. For all the considered materials, thermal conductivity of the acoustic modes is more than 70% of the total thermal conductivity along the ZZ direction. LA and TA modes are always good heat carriers and provide significant contribution to total thermal conductivity along the both ZZ and AC directions.

IV Conclusion

The lattice thermal transport properties of group VA monolayers (P, As, Sb, Bi, PAs, PSb, PBi, AsSb, AsBi, SbBi, P3As1, P3Sb1, P1As3, As3Sb1), found to have no negative frequencies in their phonon dispersions, were systematically investigated by using first principles based calculations. In this manner, the effect of alloying on lattice thermal transport properties of pristine P, As, Sb and Bi were estimated as well. Our calculations revealed that $\kappa$ of all these monolayers can be correlated by their average atomic masses as follows, $\kappa(m)=c_{1}+c_{2}/m^{2}$. In addition, $\kappa$ of all the considered materials are almost inversely proportional to temperature ($\kappa$ $\sim$ T*-1* ) for the T values, $\geq$ 100 K and therefore the $\kappa(m)$ $\sim$ $m^{-2}$ correlation is clearly applicable for the same T range. Their anisotropic crystal structures give rise to significantly different phonon group velocities and peculiarly scattering rates (phonon relaxation times) along the ZZ and AC directions, thereby anisotropic thermal conductivities. Using a simple formula developed by SlackSlack (1973), we revealed that the anharmonic effects, which is becoming more stronger with an increase in the mass of constitute atoms, play an important role in lowering the thermal conductivity of these monolayers. Our investigations clearly point out that thermal transport properties of phosphorene is tunable via substitutional Sb, As and Bi doping. For instance, the thermal conductivity of P along ZZ direction, 109.6 $Wm^{-1}K^{-1}$, can be suppressed as low as 10 $Wm^{-1}K^{-1}$. Efficient thermoelectric applications demands materials with ultra low thermal conductivity and good electrical conductivity. In this respect, we can propose PBi as a good candidate due to its ultra low kappa along AC direction, 1.5 Wm*-1K-1*.

V acknowledgement

We acknowledge the support from Scientific and Technological Research Council of Turkey (TUBITAK-115F024). CS acknowledges the support from Anadolu University (BAP-1407F335, -1705F335). and Turkish Academy of Sciences (TUBA-GEBIP). A part of this work was supported by the BAGEP Award of the Science Academy.