The lattice thermal conductivity in monolayers group-VA: from elements to binary compounds

The lattice thermal conductivity ( κL ) of the monolayers of partial group-VA elements and binary compounds are systemically investigated by the first-principles calculations and phonon Boltzmann transport equation (PBTE), including aW-antimonene, α-arsenene, black phosphorus, α-SbAs, α-SbP and α-AsP. The κL values decrease with the increasing of atomic mass for these materials with similar geometry and valence structures. It is ascribed to phonon branches softening, low phonon group velocity, and large Grüneisen parameters. Due to the neutralization of phonon group velocity and phonon lifetime, κL of binary compounds is between their corresponding elements. As the atomic radius and mass increase, the bond strength and the phonon group velocity decreases. Furthermore, the dimensionless parameter γ2/A , which comes from the Slack equation and only has the dependence of Grüneisen parameter, grows up with the atomic mass rising, which indicates that a larger anharmonicity is present in the heavier V-V monolayers. For SbAs and SbP compounds, the thermal conductivity anisotropy mainly results from the anisotropy of elastic coefficients along armchair and zigzag directions. Our results highlight the impact of atomic arrangement on the thermal conductivity of group VA binary compounds. This work paves a way to modulate the thermal conductivity of 2D VA elements by incorporation atoms with suitable mass and may guide to improve thermoelectrical performance via the alloying method.


Introduction
Because of the overconsumption of energy and increasing environmental problems in recent years, the urgency to seek more efficient energy conversion devices and clean and renewable energy resources is becoming stronger. Thermoelectric (TE) materials, which can utilize waste heat to convert to electrical energy directly during energy conversion, are useful to improve efficiency of energy conversion. The TE conversion efficiency is characterized by TE figure of merit zT, which is a dimensionless figure and calculated by zT are Seebeck coefficient, electrical conductivity, electric thermal conductivity, lattice thermal conductivity and absolute temperature, respectively. To achieve a high zT value, TE materials should have a high power factor (S 2 s) and a low total thermal conductivity ( E L k k + ). Because of the coupling relationship of S, σ and E k [1,2], carrier optimization is inefficient in improving zT [3][4][5]. Instead, it is efficient to take strategies on reducing L k or seeking low L k materials for high-performance TE materials.
A lot of studies have witnessed that group VA 2D materials can be as promising TE candidates. Three monoelemental phosphorene, α-arsenene, and aW-antimonene are identified as candidates for high TE energy conversion efficiency because of their extremely high carrier mobilities and large Seebeck coefficients, but the main challenge comes from their high thermal conductivities [6]. First-principles simulations reveal that zT of Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. phosphorene may reach the criterion for commercial development at 500 K and is close to 1 at room temperature after moderate doping [7,8]. Puckered arsenene has a low lattice thermal conductivity at roomtemperature and its preferential thermal transport direction and electrical transport direction are perpendicular to each other, which makes α-As a promising TE material [9]. Buckled antimoene exhibits a high TE performance with maximal room-temperature zT value of 2.15 [10,11].
In order to further increasing group VA 2D TE performance, reducing lattice thermal conductivity is a significant issue. Alloying is an effective method to decrease lattice thermal conductivity, which has been verified by many experiments [12][13][14]. Therein only a stoichiometric composition is considered in this work for the simplicity. β-SbAs holds a lower lattice thermal conductivity relative to monolayer antimonene and arsenene [15], which is attributed to a lower phonon lifetime of β-SbAs compared to Sb and As. Puckered AsP with three different chemical bonding styles all have a lower thermal conductivity than black phosphorus [16]. The κ of β-AsP monolayer approximately equals to the As monolayers, due to neutralization between the reduction of phonon lifetimes and the enhancement of phonon group velocities from As to AsP monolayers [17]. The roomtemperature thermal conductivities of β-NX (X=P, As, and Sb) binary compounds are several times higher than those of their monoelemental P, As, and Sb monolayers, respectively. The higher phonon group velocities, as well as the longer phonon lifetimes, are responsible for such an enhancement [18]. The κ of pristine P, As, Sb, and Bi monolayer and their compound structures with average atomic mass meet C C m 1 2 2 k = + along the zigzag direction, while the calculated κ values do not follow the same relationship along the armchair direction [19].
A method using elastic properties has been established and applied to quickly screen the potential TE materials with low lattice thermal conductivity [20]. Low elastic constants of α-SbAs and α-SbP has been found [21] among the two-dimensional (2D) group V-V binary compounds. It means that both of them may have low lattice thermal conductivity, in accordance with the above method. In this work, we calculate the lattice thermal conductivity of stable monolayer group-VA binary compounds, α-SbAs, α-SbP and α-AsP, and their constituents materials. The thermal conductivities in these materials are evaluated on the base of phonon Boltzmann transport equation. The results reveal that phonon transport properties of binary compounds, including lattice thermal conductivity, phonon group velocity and phonon lifetime, are in the middle compared with the monoelemental monolayer. Owing to similar geometry and valence configuration, the L k of 2D V-V group binary compounds and elemental monolayers is strongly dependent on the atomic mass.

Structure optimization
The structure are optimized by the first-principles theory based on density functional theory (DFT), which is used in Vienna ab initio simulation package (VASP) [22,23]. A plane-wave basis set is employed with kinetic energy cut off of 800 eV. We choose the generalized gradient approximation (GGA) within Perdew-Burke-Ernzerhof (PBE) formulation. A 13 13 1´Γ-centered k-points mesh is used during structural optimization for the unit cell until the energy differences are converged within 10 −8 eV, with a Hallman-Feynman force convergence threshold of 10 −6 eV/Å. In order to minimize interaction between adjacent layers, the vacuum space between layers along z-direction is set to 20 Å at least.

Lattice thermal conduction
The lattice thermal conductivity L k is based on a full iterative solution to the Boltzmann transport equation (BTE) [24]. The L k can be calculated using the ShengBTE package [24], which only needs harmonic and anharmonic interatomic force constants. The harmonic force constants are calculated from density functional perturbation theory (DFPT) [23,25], while for the cubic force constants the finite displacement difference method is employed. Based on the convergence test of L k dependence on the cutoff radius, 0.6 nm (0.7 nm) is used as the cutoff distance for α-SbP, α-As, black phosphorus (α-SbAs, aW-Sb), respectively. The q grid mesh is 110 110 1´at least. These parameters are balanced between convergent conditions and computational sources. Figure 1 shows the optimized structure of α-AsP. A similar puckered structure is also found in α-SbAs, α-SbP, asymmetric washboard structure of antimonene (aW-Sb), puckered arsenene (α-As), and black phosphorus (BP). The optimized parameters are listed in table 1, which is quite consistent with previous results [9,21,[26][27][28][29][30][31][32][33]. The out-of-plane bond distance (d 1 ), in-plane bond distance (d 2 ), and puckering height (h) increase with their relative atomic mass per unit (M). This variation is attributed to the the enhancement of atomic radius resulting in the weakness of bond strength. The dynamic stabilities of these six materials can be verified by their phonon spectra, as shown in figure 2. No imaginary frequency indicates these materials are dynamically stable [21,32,34,35]. With the average atomic mass increasing, the optical phonon branches markedly shift down and soften. Meanwhile, the slopes of acoustic phonon branches near Γ point decrease, which implies that lower phonon group velocities are present in heavier elemental monolayer. The frequency gaps for optical braches  between the three branches with lower frequency and left six branches clearly narrow, which means that these phonon scattering channels grow up. All these features in phonon spectra are beneficial to suppress the phonon transport, giving rise to a low thermal conductivity.     [36]. In addition, the κ values of As and Sb monolayer slightly deviates from the results reported in [35]. As it is known that different computational methods and computational parameters used can lead to the L k variation. Taking an SbAs monolayer as an example, the L k decreases from 4.57 to 1.56 W/mK, when the lattice constant extends by 3%. Here these parameters and methods in this work, including pseudopotential, the cut-off radius of 3rd force constants, have been verified in accordance with previous studies. The lattice thermal conductivity shows temperature dependence, as shown in figure 3, which approximately follow T 1 -. The lattice thermal conductivity decreases with total atomic mass in unit cell increasing. Three groups can be classified in terms of the constituents of binary compounds, including BP, α-AsP, and α-As, BP, α-SbP and aW-Sb, α-As, α-SbAs, and aW-Sb. Lattice thermal conductivity of binary compounds is between its counterpart elements, i.e., the lattice thermal conductivity is neutralized. It is similar to the L k of β-AsP and its counterpart elements [37]. In the light of kinetic theory of thermal conductivity, phonon group velocity and its relaxation time are closely relevant to the phonon transportation. Moreover, the phonon relaxation time is strongly dependent on the population of phonon scattering channels P 3 and scattering strength Grüneisen parameter γ. The contribution of acoustic phonon to thermal conductivity is usually predominant in contrast to optical phonon branches. Figure 4 shows v g , τ, P 3 , and γ of acoustic phonon branches of BP, α-AsP, and α-As. The phonon group velocity and relaxation time become smaller and smaller from the BP to AsP and to As monolayer, which leads to the L k decrease. While P 3 and γ have opposite trends from the BP to AsP and to As monolayer, which means that scattering channels and scattering strength increase with the increase of average atomic mass. For other two groups, the group of BP, α-SbP and aW-Sb, the group of α-As, α-SbAs, and aW-Sb, there are similar variation trends for v g , τ, P 3 , and γ. It is noteworthy that the lattice thermal conductivity of monolayer β-SbAs is lower than those of both monolayer As and Sb [15]. It is different from our work. The structural difference between β-SbAs and α-SbAs may be responsible for the L k difference.

Results and discussion
The phonon transport anisotropy is characterized by the ratio of lattice thermal conductivity along zigzag and armchair directions, Lzz Lac k k . The thermal conductivity is strongly dependent on phonon group velocity and phonon relaxation time, therefore its anisotropy is also dependent on the anisotropy of both factors. Similar to the monatomic linear chain model, the phonon group velocity is determined by bond strength and average atomic mass. Thus the anisotropy of group velocity is determined by that of its elastic constant in the same material. The ratio 2.36 (2.35) of thermal conductivity for SbAs (SbP) is almost equal to the ratio 2.31 (2.30) of square root of elastic constants along zigzag and armchair directions. It reflects that the anisotropy of phonon group velocity plays a predominant role in the anisotropy of heat transfer for the SbAs (SbP) compounds.
3.1. Phonon group velocity: v g One determinant factor of the phonon group velocity is an interatomic bond strength, which can be reflected by the largest second order force constants. The largest second order force constants of these six materials, α-SbAs, α-SbP, α-AsP, aW-Sb, α-As and BP, reach 7.72, 8.84, 9.86, 6.82, 8.32 and 12.83 eV/Å 2 , respectively. The larger value means a strong interatomic bonding and thus enhances phonon group velocity, increasing the lattice thermal conductivity. The variation of group velocity with force constant square root is shown in figure 5(a), in which the relationship between the both almost satisfies a linear equation. Other determinant factor of a phonon group velocity is an atomic mass. The dependance of phonon group velocity with the inverse of average atomic mass square root m 1 2 is shown in figure 5(b). The phonon group velocity is inversely proportional with m 1 2 , which is accordance with the simple harmonic model. Surprisingly, it is different with m −1 correlation in [19], which is considered due to the anharmonicity in these materials. In this work, the anharmonicity in this system is described by the usual Grüneisen parameter γ.
In terms of the dependance of atomic mass and bond strength with phonon group velocity, shown in figure 5(c), the group velocities of α-AsP is higher than that of α-SbP and α-SbAs because of stronger bonds and lighter atomic mass. Additionally, the group velocities of α-AsP, α-SbP and α-SbAs lie between their corresponding monoelemental monolayers. Binary compounds have neutralized phonon group velocity in comparison with their constituent monolayers, due to middle bond strength and average atomic mass.

Phonon anharmonicity in the Slack frame:
The lattice thermal conductivity is mainly governed by an anharmonic phonon scattering, because of harmonic phonons with no chance to scatter from each other. At room temperature, three-phonon scattering processes, which corresponds to lowest order anharmonic terms, are stronger than the other higher order anharmonic scattering processes and thus dominate the performance of the lattice thermal conductivity [38]. In the Slack model, the lattice thermal conductivity is determined by four factors, including (i) average atomic mass, (ii) interatomic bonding, (iii) crystal structure, and (iv) anharmonicity [39], waves, and v L and v S can be calculated by using the bulk modulus B and shear modulus G, According to the Voigt-Reuss-Hill (V-R-H) theory [40], the corresponding elastic properties, including B and G, can be calculated from the elastic constants. The elastic constants of α-SbAs and α-SbP are in well agreement with previous calculated values [21], which shows the reliability on using DFPT. According to above equation, the calculated Debye temperature of monolayer P, As, Sb, AsP, SbP and SbAs are 201, 99, 87, 137, 100, and 82, and they are almost same with Kocabas' work [19]. While Debye temperatures of those materials in Kocabas' work [19] are derived from the phonon dispersion relation. As a result, it is reliable to calculate the average sound speed and Debye temperature through elastic properties. is used to predict the size of anharmonicity in these systems [19]. A 2 g is directly related to the anharmonicity of the crystal. The almost linear relationship between A 2 g and average atomic mass is present in figure 6. The anharmonicity for these systems enhances with the average atomic mass increasing.

Summary
In this work, we calculate the lattice thermal conductivity of α-AsP, α-SbAs, α-SbP, aW-Sb, α-As and black phosphorus using first-principles calculations and an iterative solution of the PBTE. The κ decreases with the atomic mass increasing for these materials with similar geometry and valence structures. Neutralized κ in binary compounds are found in comparison with corresponding monoelemental monolayers, which is owing to neutralized phonon group velocity and phonon scattering anharmonicity. The thermal conductivity of group VA monoelemental monolayers can be reduced by replacement with heavier atoms and vice versa. This work paves a way to modulate thermal conductivity by incorporation atoms with suitable mass.