Strong electron–phonon coupling influences carrier transport and thermoelectric performances in group-IV/V elemental monolayers

The interactions between electrons and phonons play the key role in determining the carrier transport properties in semiconductors. In this work, comprehensive investigations on full electron–phonon (el–ph) couplings and their influences on carrier mobility and thermoelectric (TE) performances of 2D group IV and V elemental monolayers are performed, and we also analyze the selection rules on el–ph couplings using group theory. For shallow n/p-dopings in Si, Ge, and Sn, ZA/TA/LO phonon modes dominate the intervalley scatterings. Similarly strong intervalley scatterings via ZA/TO phonon modes can be identified for CBM electrons in P, As, and Sb, and for VBM holes, ZA/TA phonon modes dominate intervalley scatterings in P while LA phonons dominate intravalley scatterings in As and Sb. By considering full el–ph couplings, the TE performance for these two series of monolayers are predicted, which seriously downgrades the thermoelectric figures of merits compared with those predicted by the constant relaxation time approximation.


INTRODUCTION
Two-dimensional (2D) materials offer challenging opportunities for a wide range of applications in nanoscale devices. Among the extraordinary properties in 2D materials, carrier transport properties play the key role in determining the performance in microelectronic, optoelectronic and thermoelectric devices. Until now, extensive studies have been devoted on carrier transport properties in 2D materials, which are dominantly determined by the interaction between electrons and phonons 1,2 . Recently, it has been reported that the electron-phonon couplings can in turn efficiently limit the phonon transport properties and reduce the lattice thermal conductivities in graphene 3 , MoS 2 , and PtSSe 4 . In the simplest way, the electron-phonon (el-ph) interaction can be understood as originated from the electrostatic potential (known as deformation potential) generated by exciting a phonon in the crystal, which, in turn, affects all carriers directly. Subsequently, the deformation potential approximation (DPA) method based on the intravalley coupling between electrons and long-wavelength longitudial acoustic (LA) phonon modes 5,6 , was widely used to estimate or understand the carrier mobilities in many non-polar 2D semiconductors including graphene 7 , phosphorene 8 and transition metal dichalcogenide (TMD) 9 . However, in some polar or highly anisotropic systems, the DPA method usually misestimates the intrinsic carrier mobility, in which the coupling from the longitudinal optical (LO) phonon modes described as the Fr€ ohlich interaction, and the direction-dependent contributions can not be neglected. For example, as previously reported 2,10,11 , the scatterings from flexural ZA phonons contribute dominantly to the total scattering rates in group IV elemental monolayers (silicene, germanene, and stanene) due to the lack of σ h -symmetry in buckled honeycomb structure belonging to D 3d point group, which fails the DPA method. Therefore, depending on the specific materials and the situation of the underlying mechanism, systematically computing energy-and momemtum-dependent el-ph interactions is necessary to accurately predict phononassisted processes.
Inspired by the extraordinary properties of graphene, the group IV elemental materials of silicene (Si), germanene (Ge), and stanene (Sn) as the promising alternatives, have attracted tremendous interests. Studies show that such monolayers all possess Dirac fermions similar to graphene, high mechanical flexibility and high electron mobility, leading to the potential applications in batteries and topological devices 12,13 . Monolayer Sn was predicted to possess an intrinsic spin-orbit coupling (SOC) gap~100 meV, which is ideal to realize the quantum spin Hall (QSH) effect at room temperature 14,15 . Recently, Jian et al. successfully synthesized large-scale and high-quality stanene on Sb (111), providing a good platform to study interesting phenomena such as QSH effect 16 . In contrast to 2D group IV materials, 2D group V elemental materials of phosphorene (P), arsenene (As), and antimonene (Sb) with the same buckled honeycomb structure are semiconductors with wide band gaps, probably leading to the applications in transistors with high on/ off ratios 17 or optoelectric devices operating in the visible range 18 . Recent works reported that, the surface plasmon resonance sensor based on β-phosphorene/MoS 2 heterojunction possesses significantly higher sensitivity compared with conventional or graphene-based SPR sensors 19 , and the monolayer arsenene has the potential to serve as anode materials in highperformance Magnesium-ion batteries 20 . Interestingly, the βantimonene undergoes a transformation from a topological semimetal to a topological insulator at 22 bilayer, then to QSH phase at 8 bilayers and finally to a topological trivial semiconductor at 3 or thinner bilayers 21 .
Herein, based on the first principles method, we systematically investigate the full el-ph coupling effects in 2D group IV and V materials. These two series of 2D materials indeed possess a broad range of carrier mobilities, which may lead to various applications. Compared with group IV counterparts, group V materials are more buckled, strengthening the overlap of the p z orbitals, thus may lead to the enhancement of the interaction between the electrons and ZA phonons. The conduction band and valence band for group V materials are flatter than group IV materials, which is beneficial to satisfy the energy conservation condition and further increases the el-ph coupling near the Fermi level. To further investigate the roles of different phonon modes, we also derive in details the selection rule for the full el-ph scatterings in these two series of materials. Based on the group-theory analysis, we find that, the LA phonons play a dominant role in the hole scattering via acoustic phonon modes in arsenene and antimonene, which is in contrast to the previous report that group V materials have the same scattering mechanism as group IV materials 22 . For monolayer phosphorene, arsenene, and antimonene, despite the different bands and locations of valence-band maximum (VBM), LA phonon modes with A irreducible representation (irreps) of C 2 group near Γ point dominate the intravalley scatterings due to the identical symmetry and irreps of the intial and final electronic states. However, phosphorene is the only material with degenerate VBM where the intervalley scattering from one VBM to another is dominated by ZA and TA phonon modes. The VBM of monolayer arsenene and antimonene locates at Γ point without degeneracy, hence the intravalley scattering becomes dominant. Finally, by considering the full el-ph scatterings, an accurate prediction of the carrier transport properties and thermoelectric performances of these two series of materials are achieved.

RESULTS AND DISCUSSION Optimized configurations and band structures
Compared with the planar geometry of graphene, 2D group IV and V crystals have buckled structure with two sublayers. The top and side views of configurations of the 2D group IV materials (silicene, germanene, stanene) and group V materials (phosphorene, arsenene, antimonene) are shown in Fig. 1a. The optimized lattice constants a, the distance between two sublayers h, bond angle θ and bond length d are shown in Supplementary Table 1 which are consistent with previous reports [23][24][25][26] .
In graphene, the 2s orbital is hybridized with p x and p y orbitals to form three planar σ bonds via sp 2 hybridization, and the remaining p z orbital forms π bonding between adjacent C atoms to guarantee the planarity 27 . However, when the bond length increases, the p z -p z overlapping decreases and the weaker π bonding can not ensure the stability of the planar geometry. By buckling, the planner sp 2 bonding is dehybridized to form sp 3 -like hybridization to enhance the overlapping between p z orbitals 28 . Figure 1b shows the calculated electronic localization function (ELF) for silicene and phosphorene. The ELF is a positiondependent function with values ranging from 0 to 1. ELF = 1 reflects the maximum probability to find the localized electrons while ELF = 0.5 means the electron-gas-like behavior 29 . The strength of the π bond of 2D group V materials is much higher than their group IV counterparts. The one more electron in P, As and Sb increases the degeneracy of electronic states, leading to the unstable state according to the Jahn-Teller theorem, and by a high degree of buckling the degeneracy can be relieved, resulted from the broken symmetry 30 . Hence, compared with group IV materials, the group V counterparts possess larger buckling heights, smaller lattice constants and bond angles, which together increase the overlap of the p z orbitals.
In order to further investigate the strength of chemical bonding, we also calculate the bulk (B H ) and shear moduli (G H ) for group IV and group V materials, which determine the ductile-brittle nature of materials by Pugh's ratio B H /G H 31 . According to Pugh's rule, a small Pugh's ratio indicates brittle nature of a material. As shown in Supplementary Table 1, the Pugh's ratios for phosphorene, arsenene, and antimonene are all lower than their buckled group IV counterparts, indicating that the high degree of buckling in group V materials strengthens the π bonding. The trend of Pugh's ratios B H /G H (Si) < B H /G H (Ge) < B H /G H (Sn) and B H /G H (P) < B H /G H (As) < B H /G H (Sb) indicates that increasing bond length d lowers the bonding strength by decreasing orbital overlap for the same group materials, compensating the effect from the increase of buckling.
The calculated band structures of silicene, germanene, and stanene are shown in Supplementary Fig. 1a-c. All these group IV materials are semimetals with a Dirac-cone located at the K point, and the SOC effect opens the band gaps of 1.48 meV, 23 meV, 72 meV, respectively. Based on the band structures projected by the atomic orbitals, the band structure near the conduction band minimum CBM (C 1 ) and valence band maximum VBM (V 1 ) are dominantly contributed from p z orbitals. However, as shown in Supplementary Fig. 1d-f, as a comparison, the group V materials are all semiconductors with sizable indirect band gaps. The calculated band gaps without (with) SOC are 1.94 eV (1.94 eV), 1.59 eV (1.47 eV), 1.26 eV (0.99 eV) for phosphorene, arsenene, and antimonene respectively. The CBM locates close to the middle of the Γ-M, which is mainly contributed from p x and p z orbitals, and the VBM locates at Γ point for arsenene and antimonene which is dominantly contributed from p x and p y orbitals, while for phosphorene the VBM locates along Γ-K direction mainly  contributed by p z orbitals. The multiple local valleys and peaks, which may play important roles in intervalley el-ph scattering, are labeled as C 2 , C 3 and V 2 , V 3 , and their difference from CBM or VBM are shown as E C2 , E C3 , E V2 , E V3 in Table 1. For instance, the energy difference between V 1 and V 2 valleys (E V2 ) is 50 meV in phosphorene.
Full el-ph scatterings: selection rules and carrier mobilities The calculated carrier mobilities at 300 K along armchair (a) and zigzag (b) directions for these two series of 2D materials, based on the modified DPA theory 32 , are shown in Table 2, accompanied with effective mass m * , deformation potential constant E l , elastic modulus C 2D , and the relaxation time τ defined as τ = μm * /e, which is determined by the interaction between carriers and LA phonons. The deformation potential constant E l reflects the coupling strength of electron and long-wavelength LA phonons, defined as ΔE CBM(VBM) /(δl/l 0 ) 33 , in which ΔE CBM(VBM) is the shift of the VBM and CBM energy level with respect to the lattice deformation δl. The elastic modulus C 2D is related to the interatomic forces. Hence, both E l and C 2D are influenced by the overlap of atomic orbitals. The values of E l and C 2D of 2D group V materials are larger than their group IV counterparts due to the enhanced overlap of atomic orbitals by buckling. The obtained mobilities for group IV materials based on the DPA method are comparable to graphene in the order of 10 5 -10 6 cm 2 V −1 s −1 , and those of group V materials are in the range of about 10-1000 cm 2 V −1 s −1 . However, as mentioned above, the DPA method may fail for the band-convergence systems such as the group IV and V elemental crystals shown in Table 2, therefore the full el-ph scatterings should be investigated to accurately predict the carrier transport properties. According to Eq. (1), the carrier mobility is mainly determined by the electron velocities and the relaxation times τ.
The calculated el-ph scatterings for CBM electrons via acoustic and optical phonon modes for these two series of 2D materials are shown in Fig. 2 and Supplementary Fig. 4, respectively, and the calculated el-ph scatterings for VBM holes via acoustic and optical phonon modes for these two series of 2D materials are shown in Fig. 3 and Supplementary Fig. 5, respectively. All the scattering rates are calculated at 300 K within~0.4 eV. Obviously, the calculated scattering rates for group V materials are much larger than their group IV counterparts, in consistency with the behaviors of deformation potential constant E l and elastic modulus C 2D mentioned above, which is probably resulted from two mechanisms: (i) The enhanced p z -p z orbitals overlapping enlarges the el-ph matrix element g mnν k; q ð Þ described in Eq. (2) due to the stronger sensitivity to the atomic positions. (ii) The band structures of group V materials near Fermi level are much flatter than the those of group IV materials, which is more beneficial to satisfy the energy conservation conditions in the scattering processes.  Table 2. Calculated effective mass m*, deformation potential constant E l , 2D elastic modulus C, carrier mobility μ DPA , relaxation time τ and μ el−ph along a and b directions, and roles of dominant contributions to total acoustic electron-phonon couplings.
Direction Carrier type Effective mass(m 0 ) E l (eV) Starting from the CBM, with the increase of electronic energy, the scattering rates increase gradually for group IV materials while dramatically for group V materials especially for phosphorene, which possesses a wide and nearly flat region near the C 1 valley as shown in Supplementary Fig. 1d, leading to about 4 times larger in the scattering rate at E = 0.07 eV than those in arsenene and antimonene.
Due to the small phonon energies, for monolayer Sn, As, and Sb, since the energy of CBM electrons is well smaller than E C2 , only the intravalley scattering and intervalley scatterings between degenerate C 1 valleys are allowed for shallow dopings. When the energy of C 1 -pocket electrons is higher than E C2 , C 1 -C 2 intervalley scattering is allowed, and the scattering rate increases abruptly as a result. In particular, for monolayer Sb, when the energy of C 1 -pocket electrons is higher than E C3 , both C 1 -C 2 and C 1 -C 3 intervalley scatterings are allowed, leading to a further increase of the scattering rates, as shown in Fig. 2f and Supplementary Fig. 4f.
The mode-and energy-resolved el--ph scattering rates of V 1peak holes for these two series of materials at 300 K within~0.4 eV are shown in Fig. 3 and Supplementary Fig. 5. The behaviors of the scattering rates of holes for silicene and germanene are similar to their electron case, due to the Dirac-cone dispersions for both CBM and VBM carriers. For stanene, as the hole energy descreases to E V2 , the scattering rate experiences a sharp rise, implying the activation of V 1 -V 2 intervalley scatterings. The intravalley and V 2 -V 1 intervalley scatterings are also allowed for V 2 -peak holes with energy lower than E V2 , as shown in Fig. 3c and Supplementary Fig. 5c. For the hole scatterings in phosphorene, due to the flat valence band near Fermi level, also reflected in the effective mass as listed in Table 2, the scattering rates of V 1 -peak holes with energy near Fermi level are very high compared with those for monolayer As and Sb. When the hole energy is lower than E V2 , the V 1 -V 2 intervalley scattering are activated and the scattering rates further increase up to~3.7 × 10 15 s −1 near −0.05 eV. Figure 4a shows the contour map of the valence band in Brillouin zone for phosphorene. It is clear that the energy difference between V 1 and V 2 peaks is rather small, which is beneficial for satisfying the conservation of energy. Figure 4b shows the scattering rates of holes via ZA phonons in the valence band in the Brillouin zone, which reveals that, holes at the junction of the two peaks undergo strong scatterings, due to the typical quadratic dispersion of ZA mode and flat electronic dispersion.
The respective mode-resolved scatterings rates for electrons at CBM and holes at VBM for these two elemental crystals are shown in Supplementary Figs. 6 and 7. For group IV elemental monolayers (Si, Ge, Sn) with shallow n-and p-dopings, the intravalley scattering rates for CBM electrons or VBM holes via TO and LO phonon modes are significantly larger than those via acoustic phonon modes, and ZO modes also contribute significantly in monolayer Ge. Overall, the intervalley scatterings rates significantly surpass the intravalley scatterings, and the contributions from ZA, TA, and LO phonon modes dominate the intervalley scatterings. However, the situations for group V elemental monolayers (P, As, Sb) are different, especially for intravalley scatterings. For intravalley scatterings in monolayer P, the contributions from ZO, LA, and TA phonon modes dominate. In monolayer As, the contributions from TO, LO, LA, and TA phonon modes dominate, and in monolayer Sb, the contributions from LA, ZO, and LO phonon modes dominate. Similarly, intervalley scatterings for CBM electrons are much stronger than intravalley scatterings in group V monolayers, and the contributions from ZA and TO modes dominate. For el-ph scatterings for VBM holes in group V elemental crystals, ZA and TA modes dominate interpeak scatterings in monolayer P, as shown in Supplementary Fig. 7. As shown in Supplementary Figs. 1 and 7, since the global VBM located at Γ point is well above the other local VBMs, interpeak scatterings are not considered in monolayer As and Sb for shallow p-dopings.
In order to further clarify the underlying mechanisms behind electron-phonon interactions, based on the group-theory analysis, the selection rules in terms of the symmetry of systems for the interaction between phonons and electrons are derived. The electron-phonon interaction depends on the electronic wave functions of the initial states ϕ i (k), final states ϕ f (k + q) and the coupled phonon with the eigenvector q described by the phononinduced deformation potential δ q,ν U 34-36 . Based on the Fermi's golden rules, in order to calculate the scattering probability from initial state ϕ i (k) to final states ϕ f (k + q) of electrons or holes via a specific phonon mode, we only need to obtain the core matrix, which could be expressed as 37 where the perturbation potential δU ¼ P l;a q l;a ∂Uðri Þ ∂q l;a , in which q l,a is the ion displacement of the a atom in the l-th unit cell, and r i is the electron coordinate. Therefore, the expression for the selection rules for such scattering process M based on the group theory is written as, where Γ i/f represents the irreps for electron initial/final states, and Γ ph denotes the symmetry of the corresponding phonon involved in this process. The involved electrons and phonons must obey the rules of momentum conservation and the compatibility relation between little groups, since they possibly locate at different points in the Brillouin Zone. non null/null means the allowed/forbidden transition channels. Due to the lack of in-plane symmetry σ h in monolayer Si, Ge, Sn of group IV and monolayer P, As, Sn of group V, their electronic and phonon symmetry is lowered considerably compared to plane graphene. For simplification without loss of generality, we only discuss on the cases for el-ph scatterings via acoustic phonon modes, and the cases for those via optical phonon modes can be discussed in a similar way. As shown in Supplementary Fig. 1 and Supplementary Table 2, for the intrapeak scatterings of holes in monolayer Si, Ge, and Sn, the hole states in the V 1 peak along K-M and K-Γ involved in the intrapeak scatterings transform with the A and B irreps respectively, both belonging to the little group C 2 . Mathematically, it is allowed for the initial states and final states to belong to the same irreps in the el-ph scattering, i.e. A Ð A or B Ð B, but the contribution of such scattering channel is negligible compared to the scatterings with different irreps for initial and final states in the real systems we studied here. Because, as shown in Supplementary Fig. 1a-c, the energy difference of the scattering channel A Ð B is smaller than that for the scattering channel of A Ð A or B Ð B, leading to the larger density of states of the former scattering channel, which is beneficial for the el-ph scattering in real systems. Therefore, for intrapeak scatterings in monolayer Si, Ge, and Sn, we only need to consider the initial and final states transforming with the A and B irreps, i.e. A Ð B processes.  Fig. 3a-c. It should be noticed that, some phonon modes locating near Γ point transform with the irreps belonging to the little group of C S , but it is reasonable to exclude them since the initial and final states of carriers are with C 2 symmetry, which are forbidden to be coupled to the phonon modes with C S group considering the incompatibility relation between C 2 and C S 38 . Furthermore, as shown in Fig. 5, for the interpeak scatterings of holes in group IV monolayer materials, since the degenerate VBMs are located at six corners in the Brillouin Zone, there are at least three different interpeak-scattering channels, including the scattering to the nearest, next-to-nearest and opposite-side K points. However, the interpeak scattering to the next-to-nearest VBMs is actually forbidden since the phonon modes involved are with C S symmetry, which is incompatible to the C 2 symmetry. Therefore, only TA modes are allowed restricted by the selection rules and conservation of momentums, as shown in Supplementary Table 4, which is consistent with previous work 2 . Moreover, as shown in Supplementary Tables 3 and 5, since the CBM and VBM points touch at the K points forming the so-called Dirac points in group IV elemental crystals, therefore, in the same way, the intravalley scatterings for CBM electrons scattered via the ZA and TA phonon modes are allowed, and the intervalley scatterings for CBM electrons between degenerate CBMs via TA modes are allowed. It should be noted that, although ZA phonon mode is allowed in intravalley/intrapeak el-ph scatterings for group IV elemental crystals obeying the selection rule, the calculated scattering rates for intravalley/intrapeak el-ph scatterings via ZA phonon mode as shown in Supplementary Figs. 6a-c and 7a-c are indeed trivial. This disagreement can be understood by the fact that, CBM/VBM carriers with linear dispersion as shown in Supplementary Fig. 1 can not scatter to each other via ZA phonon modes with quadratic dispersion near Γ point, restricted by the energy conservation.  Table 4, for the interpeak scatterings of the VBM holes for monolayer P, As and Sb, the situation is quite different due to the different locations of VBMs. For monolayer As and Sb, since the VBM locates at Γ point, which lacks of degeneracy in the Brillouin Zone, the interpeak scatterings between VBMs no longer exist for the VBM holes at shallow p-dopings. As a result, only LA mode dominates the el-ph scatterings of VBM holes for monolayer As and Sb, since the initial and final states of VBM holes for monolayer Sb and As near the Γ point transform with A irreps belonging to the little group C 2 , which is in good agreement with the mode-resolved scattering rates of VBM holes for monolayer As and Sb as shown in Fig. 3e, f, respectively. However, special care should be taken for interpeak scatterings of VBM holes for monolayer P, since the degeneracy of the VBM points in the Brillouin Zone for monolayer P is six, as shown in Fig. 5d. Three interpeak-scattering channels to the nearest neighboring VBMs are possible for VBM holes, and the little groups for the initial and final hole states are C 2 . However, the phonon involved in the interpeak scattering to the secondly nearest neighboring K point transforms with the little group of C S , incompatible to C 2 , therefore this scattering channel is forbidden. Thus, in monolayer P, the scatterings of VBM holes to the nearest and the thirdly nearest neighboring VBMs are allowed, and respectively dominated by the TA and ZA phonon modes. The selectrion-rule analysis is in good agreement with the numerical results shown in Fig. 3d.

As shown in Supplementary
Similarly, as shown in Supplementary Fig. 1 and Supplementary Table 3, the CBMs of group V materials (monolayer P, As, and Sb) locate at the points along Γ-M with the little group C S , and thus the corresponding irreps of the initial states and final states of CBM are A′. The acoustic phonon modes of ZA, TA, and LA near Γ points and along Γ-M transform with the A′, A″, and A′ irreps, respectively as shown in Supplementary Fig. 3a-c. Therefore, restricted by the selection rules, both ZA and LA modes are allowed in the intravalley scatterings of CBM electrons for monolayer P, As, and Sb. For the intervalley scatterings of electrons between degenerate CBMs, the ZA and LA modes are allowed, as shown in Supplementary Table 5. And Fig. 6 gives the transitional paths of electrons scattered by phonons in conduction bands. Based on the calculated full el-ph materix elements, the temperature-dependent and mode-resolved relaxation time and the subsequent carrier mobilites along special directions can be calculated according to Eqs. (1) and (2), and the calculated mobilities for these two series of 2D elemental materials are shown in Fig. 7a, b, respectively. The calculated mobilities at 300 K, by a comparison with those calculated by the DPA method are shown in Table 2 as well. By comparison, three distinctions about the carrier mobilities of 2D group IV and group V materials can be identified: (1) Generally, the results by the DPA method are overestimated by 2-3 orders in magnitude due to the only consideration of LA phonons and neglecting the intervalley scattering process; (2) The hole mobilities of group V materials calculated by these two methods are in roughly good agreement except phosphorene owing to the dominant LA phonons-holes coupling in monolayer β-Sb and As; (3) The obvious anisotropy in carrier transport for group V materials obtained by the DPA method no longer exists in the results considering full el-ph coupling.

Thermoelectric performance
The thermoelectric (TE) performance of a material is quantified by a dimensionless figure of merit zT(zT = S 2 σT/(κ e + κ L )), where S represents the Seebeck coefficient, σ is the electron conductivity, and κ e/L is the electronic/lattice thermal conductivity. High TE performance needs high zT value of materials with excellent electrical transport properties and poor thermal transport properties simultaneously. However, optimization of zT values in the TE materials is always challenging due to the intercorrelations of the parameters, and optimizing one leads to deteriorating the other. In recent years, several optimization strategies have been proposed to enhance successfully zT values in TE materials, and among them, the strategies of the dimension confinement and band convergence attract broad interests owing to the successful discovery of layered TE materials with extremely high zT values [39][40][41] . Furthermore, considering the appropriate band gaps and heavy atoms in arsenene, stanene, and antimonene, high zT values can be expected in these two series of 2D elemental materials. As reported previously 42 , antimonene has the zT value of 2.15 at room temperature and can be further enhanced to 2.90 with strain engineering, predicted by the DPA and constant relaxation time approximation (CRTA) methods. For arsenene, different constant relaxation times lead to ten times deviation in zT values 26 , implying the importance of the accurate calculations of relaxation time in the zT values for a material.
Here, we also calculate the TE and the related transport properties of these two series of 2D elemental crystals by considering full el-ph coupling. Figure 8a shows the calculated Seebeck coefficient S for the group IV materials as a function of chemical potential E f at 300 K using band-and momentumdependent relaxation times considering full el-ph couplings. For comparison, the results using the CRTA method are also presented. As we know, the behavior of the Seebeck coefficient can be understood by the Mott relation [43][44][45] , Where N(E) and τ(E) are the energy-dependent density of states (DOS) and electronic relaxation time, respectively. The Seebeck coefficient can be separated into the band term and the scattering term. A large enhancement of the DOS near the Fermi level leads to a high S. The scattering term indicates that S is also related to the logarithmic energy derivatives of τ(E), which is inversely proportional to the scattering rates and usually ignored in TE community. In the CRTA method, the scatterings term is equal to zero considering the definition as shown in Eq. (3). Considering the n-type systems (E f > 0), the scattering rates of electrons gradually increase as shown in Fig. 2a-c, therefore the scattering term is positive and considerably large, compared with the zero scattering term in the CRTA method. Moreover, when the scattering term exceeds the band term, S even undergoes an abnormal sign-reversal. Hence, regarding the relaxation time as a constant for group IV materials brings significant errors to S. However, for the group V materials as shown in Fig. 9a, the energydependent relaxation time influences little on the Seebeck coefficient resulted from the overwhelming band term in the Seebeck coefficient. The maximum Seebeck coefficient of group V materials is~1500 μVK −1 , which is nearly 8 times larger than the group IV materials and many traditional bulk TE materials including Bi 2 Te 3 (215 μVK −1 ) 46 , PbTe (185 μVK −1 ) 47 , and SnSe (~510 μVK −1 ) 48 .
For the electrical conductivity σ as shown in Figs. 8b and 9b, near E f = 0, σ reaches~1 × 10 5 Ω −1 m −1 ,~1 × 10 3 Ω −1 m −1 ,~1 × 10 4 Ω −1 m −1 for silicene, germanene, and stanene respectively, nearly 1, 4, and 3 orders smaller in magnitude than those calculated using the CRTA method respectively. The σ for n-type group V materials are nearly overestimated by an order of magnitude. For a p-type system, the DPA method gives a better estimation of σ for arsenene, and antimonene compared with phosphorene. As mentioned above, in phosphorene the intervalley scatterings via ZA phonons are strong which fails the DPA method.
The total thermal conductivities is the sum of the electronic contribution κ e and the lattice contribution κ L , i.e. κ = κ e + κ L . The electronic thermal conductivity κ e exhibits similar trend with σ for both DPA and the full el-ph scattering methods. According to the Wiedemann-Franz law, κ e = LσT, where L is the Lorenz number, increasing σ leads to the increase of κ e . The lattice thermal conductivities κ L for these two series of 2D elemental crystals are chosen from our previous reports 49 .
The results of dimensionless figure of merit zT are shown in Figs. 8c and 9c. For group IV materials, due to the highly overestimated σ by the CRTA method, despite the relatively small S, the maximum zT value by the CRTA method are in the range of 1.4-1.7, which is much larger than many traditional thermoelectric materials, e.g., Bi 2 Te 3 (1.2) 46 , PbTe (0.30) 50 , SnSe (0.70) 51 . However, the zT values for group V materials predicted by considering full el-ph coupling are only 0.02-0.2, much smaller than those predicted by using the DPA method. For example, the p-type phosphorene was reported as a good room-temperature thermoelectric material by the CRTA method with zT = 0.48, but the zT value decreases to 0.008 when considering full el-ph coupling. The maximum zT value at 300 K in antimonene reaches~0.4, which is 5 times smaller than that predicted by the CRTA method, and is the highest zT among these materials. As shown in Supplementary Fig. 8, when the temperature rises, the maximum zT value for antimonene at 700 K reaches 1.3, which is comparable to conventional TE materials such as PbSe (1.1 at 900 K) 52 , Cu 2 Se (1.5 at 1000 K) 53 , and PbTe (1.64 at 770 K) 54 .
In summary, using first principles calculations, we systematically investigate the effects of the electron-phonon couplings in 2D honeycomb group IV and group V materials with D 3d symmetry. High buckling increases the overlap between p z orbitals and enhances the interaction between carriers and ZA phonons. The less dispersive band structure of group V materials further enhance the electron-phonon coupling due to easier satisfaction of energy conservation condition. We find that the D 3d symmetry is not a sufficient condition for ZA phonon scattering to dominate. Based on our results of selection rules, the symmetry of transition process only permits LA phonon to be involved in intravalley scattering within VBM in monolayer phosphorus, arsenene and antimonene. However, TA and ZA modes from intervalley scattering between degenerate VBM is dominant in phosphorus and far beyond LA modes. Yet only LA modes dominate in monolayer arsenene and antimonene since no intervalley scattering between VBM at Γ point due to lack of degeneracy. By comparison to the carrier mobilities calculated by the DPA and the full el-ph couplings, respectively, we find the DPA method always overestimates carrier mobilities except hole mobilities in monolayer β-Sb and As. Furthermore, we evaluate the TE performance of all these materials with the constant relaxation time obtained from the DPA method, and band-and momentum-resolved relaxation rates respectively, which indicate that the TE properties of materials may be seriously misestimate with constant relaxation time approximation due to neglecting the intervalley scatterings or intravalley scatterings via phonons other than LA modes.

Electron-phonon coupling
We carry out density functional theory (DFT) using the QUANTUM ESPRESSO code 55 with the local density approximation (LDA) 56 . Normconserving pseudopotentials (NCPP) method with a kinetic energy cutoff of 90 Ry are used to perform the self-consistent DFT calculation. A vacuum space of 28 Å; is set along the perpendicular direction to eliminate the interlayer interacions due to the periodic boundary conditions. For the phonon dispersion, the coarse Monkhorst-Pack k-mesh and q-mesh for all the materials discussed are taken as 16 × 16 × 1 and 8 × 8 × 1 respectively. The Wannier interpolation method is used to generate the ultradense fine grid to describe the processes of electron-phonon scattering accurately in the Wannier90 57 and EPW code 58,59 . Accordingly, much fine a b Fig. 7 Carrier transport properties. The calculated carrier mobility for electrons and holes for a group IV and b group V materials considering full electron-phonon coupling with temperature ranging from 100 K to 500 K.
k and q meshes, 240 × 240 × 1 and 120 × 120 × 1 can be used to calculate the el-ph coupling matrix guaranteeing the numerical convergence of the intrinsic carrier mobility. In the Boltzmann transport theory, the electron mobility can be defined as 60 Where n e is the number of electrons, Ω and Ω BZ denote the volume of the unit cell and the first Brillouin zone, respectively, f 0 nk is the Fermi-Dirac distribution, v nk,α = h −1 ∂ϵ nk /∂k α is the velocity of the single-particle electron eigenvalue ϵ nk . The interaction between electron and phonon lies in the parameter τ nk which is expressed as 60 1 τnk ¼ 2ImΣ FM nk ðωÞ ¼ 2π _ P mν R dq ΩBZ g mnν k; q ð Þ j j 2 ð1 À f 0 mkþq þ n qν Þδðϵ nk À ϵ mkþq À _ω qν Þ h þ ðf 0 mkþq þ n qν Þδðϵ nk À ϵ mkþq þ _ω qν Þ i where the sum is over all the final band index m for electrons and all the phonon mode index ν and wavevector q. ω qν and n qν are the frequency and Bose-Einstein distribution of phonons. ϵ mk+q and f 0 mkþq are the electron eigenvalue and the Fermi-Dirac distribution of the final state with band index m and wavevector k 0 ¼ k þ q. The electron-phonon matrix element g mnν k; q ð Þ is calculated using 61 g mnν ðk; qÞ ¼ ϕ mkþq Δ qν V KS ϕ nk j i with ϕ nk and ϕ mk+q being the initial and final electronic Bloch state, respectively. Δ qν V KS is the variation of the self-consistent Kohn-Sham (KS) potential experienced by the electrons.

Electronic transport properties
We use BoltzTraP2 code based on the rigid band approach to calculate the transport parameters for arsenene and antimonene monolayers 62 . The temperature-and doping-dependent electrical transport properties, including carrier concentrations for holes and electrons n h/e , electronic conductivity σ, electronic thermal conductivities κ el and Seebeck coefficient S are computed by solving the semiclassical Boltzmann transport equation (BTE), which can be written as 63 where Ω is the unit cell volume, f 0 is the Fermi-Dirac distribution, μ is the chemical potential, D(ε) is the density of states, σ αβ ðεÞ is the energy-dependent conductivity tensor, which can be obtained by σ αβ ðεÞ ¼ 1 N P n;k σ αβ ðn; kÞ δðεÀεn;k Þ dε , in which N is the number of sampled k points and σ αβ ðn; kÞ can be calculated by the formulae based on the Fig. 8 Thermoelectric properties of group-IV materials. The a-c seebeck coefficients d-f electrical conductivity g-i zT value as a function of chemical potential at 300 K for silicene, germanene, and stanene. kinetic theory, i.e., σ αβ ðn; kÞ ¼ e 2 τ n;k v α ði; kÞv β ðn; kÞ. The velocity of carrier v α,β (n, k) is defined by v i ¼ 1 _ ∂εn;k ∂ki ði ¼ α; βÞ. In this code, the electron-phonon relaxation time is treated as both energy and direction dependent from Eq. (2).

DATA AVAILABILITY
The data that support the findings of this study and the code for the first principles methods proposed in this work are available from the corresponding author (Hao Zhang) upon reasonable request. Fig. 9 Thermoelectric properties of group-V materials. The a-c seebeck coefficients d-f electrical conductivity g-i zT value as a function of chemical potential at 300 K for phosphorene, arsenene, and antimonene.