Dissipative transport and phonon scattering suppression via valley engineering in single-layer antimonene and arsenene field-effect transistors

Abstract

Two-dimensional (2D) semiconductors are promising channel materials for next-generation field-effect transistors (FETs) thanks to their unique mechanical properties and enhanced electrostatic control. However, the performance of these devices can be strongly limited by the scattering processes between carriers and phonons, usually occurring at high rates in 2D materials. Here, we use quantum transport simulations calibrated on first-principle computations to report on dissipative transport in antimonene and arsenene n-type FETs at the scaling limit. We show that the widely-used approximations of either ballistic transport or simple acoustic deformation potential scattering result in large overestimation of the ON current, due to neglecting the dominant intervalley and optical phonon scattering processes. We additionally investigate a recently proposed valley engineering strategy to improve the device performance by removing the valley degeneracy and suppressing most of the intervalley scattering channels via an uniaxial strain along the zigzag direction. The method is applicable to other similar 2D semiconductors characterized by multivalley transport.

Introduction

Recently, two-dimensional (2D) materials have shown a wealth of interesting properties and possible technological applications^{1,2,3,4,5,6,7,8}. Their mechanical properties, ultimate thinness, and absence of surface dangling bonds are particularly suitable for conventional and flexible electronics^{1,9,10,11,12,13,14}. The 2D materials with a finite energy bandgap have the potential to replace Si as the channel material in ultra-scaled metal-oxide-semiconductor field-effect transistors (MOSFETs) for future nanoelectronics^1,9. However, two main issues hinder the adoption of 2D materials for MOSFET applications, namely the high resistance between metallic contacts and semiconducting 2D materials¹⁵ and their relatively low mobility.

Experimental measurements of room temperature electron mobilities in 2D materials have so far reported values lower than 400 cm² V⁻¹ s⁻¹,^{9,16,17,18,19} much smaller than the value for bulk Si (about 1400 cm² V⁻¹ s⁻¹)^20,21. When employed in electron devices in place of 3D semiconductors such as Si, Ge, and III–V compounds, 2D-material-based MOSFETs do not suffer from surface roughness scattering^22,23, but have a room temperature mobility limited by defects and electron–phonon coupling (EPC). The former, just as the contact resistance, can in principle be reduced by developing dedicated engineering processes²⁴. On the contrary, the electron–phonon scattering represents an intrinsic mechanism and cannot be easily suppressed at room temperature.

Unfortunately, 2D semiconductors usually have larger electron–phonon scattering rates with respect to their 3D counterparts, due to a higher density of electronic and phononic states at energies close to the band extrema²⁰. The low phonon-limited mobility in monolayers has been also confirmed by several first-principle calculations based on density functional theory (DFT) and on the Boltzmann transport equation^{6,19,25,26,27,28,29}.

In light of these considerations, it appears compulsory to take into account the electron–phonon interactions in the simulation of 2D material-based electron devices, in order to obtain reliable quantitative predictions of the device performance at room temperature. Nevertheless, only a few recent simulation studies have fully included the effect of EPC^30,31. In other works, the mobility has been evaluated through the Takagi formula³², which takes into account only the Bardeen-Shockley acoustic-deformation-potential (ADP) scattering³³ and neglects all the optical and intervalley phonon scatterings^13,34. However, many 2D semiconductors have multiple degenerate conduction band valleys, which results in a large phase space for intervalley scatterings and in large scattering rates^19,35. The inelastic nature of these scattering processes has been already shown to deeply affect the transport in electronic devices^30,36.

More often, only the ballistic regime is considered in simulations. This approximation is usually justified by invoking the shortness of the considered device, though no experimental evidence of ballistic transport at room temperature for 2D materials other than graphene has been provided yet. Therefore, it is still unclear to which extent the electron–phonon scattering affects the room temperature transport even in devices with characteristic lengths of less than ten nanometers.

In this paper, we report the results of room temperature numerical simulations of arsenene- and antimonene-based short-channel MOSFETs with double-gate architecture. Arsenene and antimonene are among the most promising candidates as channel materials in emerging MOSFET devices, due to their small effective masses and correspondingly high theoretical mobility^13,34,37. However, previous application-oriented studies have been based either on purely ballistic simulations or have considered only acoustic phonon scattering. A full treatment of EPC in these materials is therefore still missing.

To this purpose, we use a non-equilibrium Green’s function (NEGF) transport solver including acoustic and optical intra- and intervalley EPCs with parameters fully determined from first-principles calculations. The NEGF transport solver is coupled with a Poisson solver for self-consistently computing the electrostatic potential in the whole device domain.

The main findings of this work are not limited to the specific case of arsenene and antimonene, but have more general implications and suggest that a more careful treatment of the electron–phonon coupling is necessary in simulating MOSFETs based on 2D materials in order to avoid severe overestimations of performance. In view of our results, we investigate the approach proposed by Sohier et al.³⁵ to mitigate the detrimental intervalley scattering. This technique is applicable not only to the case of antimonene and arsenene, but also to many similar 2D materials characterized by multivalley transport.

Results

Conduction band valleys

For the purpose of electronic transport simulation, the conduction band of arsenene and antimonene (a single-layer of arsenic/antimony atoms forming a buckled hexagonal honeycomb structure, see Fig. 1a) can be described as six energetically degenerate valleys, located at the midpoint of the Γ-M line. We note that the next low-energy conduction band valley is 0.23 eV (0.27 eV) above the conduction band minimum for antimonene (arsenene), thus playing negligible roles in the electronic transport at 300 K. In our simulations, we restrict ourselves to consider the six low-lying valleys, modeled within an anisotropic effective-mass approximation, as described in the section “The anisotropic effective mass model”. We have verified that this is enough to include the relevant states for transport at the doping level and temperature (300 K) considered in this work. In this regard, we compare in Fig. 2 the transfer characteristics computed in the ballistic regime within our model for an antimonene-based MOSFET with the corresponding data obtained by using an atomistic full-band model in ref. ¹³. The good agreement between our results (dashed line) and those obtained within the atomistic full-band model (dots) confirms the validity and accuracy of our approach. The spin-orbit coupling is neglected in this work since it does not affect appreciably the conduction band, as also shown in ref. ¹³.

Fig. 1: Atomic, band structure, and device schematic.

a Top and side views of the single-layer antimonene and arsenene atomic structure. The unit cell and the direction of the uniaxial strain is also represented. b Energy contour plots of the lowest conduction band (CB) in single-layer antimonene. Six-fold degenerate CB valleys are indicated. Three types of intervalley scattering channels are represented by red arrows. c Schematic view of the double-gate field-effect transistor investigated in this work. The gate length L_G varies between 4 and 9 nm. The source and drain extensions are doped with a donor concentration N_D = 3.23 × 10¹³ cm⁻². The doping region is marked by blue color and stops exactly at the channel-source and channel-drain interfaces. Two t_OX = 2 nm thick dielectric layers (HfO₂) sandwich the 2D material channel and separate it from the gate electrodes. Transport occurs along the x-axis, z is vertical confinement direction, and the device is assumed to be periodic along y. d local k-space reference frame (k_∥, k_⊥) of a CB valley and the global one (k_x, k_y), where k_x and k_y are the wavevector components along the transport and the transverse (periodic) direction, respectively.

Fig. 2: Calculated transfer characteristic including different electron–phonon scattering for an Sb-based double-gate transistor with L_G =5 nm.

I_DS–V_GS in a linear and b logarithmic scale. The considered electron–phonon scattering processes are: intravalley scattering mediated by acoustic phonons (blue line), intravalley scattering mediated by acoustic and optical phonons (black line), and both intervalley and intravalley scattering mediated by acoustic and optical phonons (red line). The results obtained in the ballistic transport regime (dashed lines) are compared to the values (gray dots) reported in ref. ¹³.

Electron–phonon couplings from first-principles

We compute the relevant electron–phonon couplings from first-principles by using the density functional perturbation theory (DFPT). The values of the deformation potentials and phonon mode frequencies of intravalley and intervalley EPCs are listed in Table 1. The details about the DFPT calculations can be found in the section “First-principles calculation details”. As shown in Fig. 1b, three types of intervalley EPCs connecting two electronic states located in different valleys (V1 ↔ V2, V1 ↔ V3, and V1 ↔ V4) are considered. In the transport calculations, all the intravalley EPCs are included, while only the intervalley EPCs with a significant deformation potential (>4 × 10⁷ eV cm⁻¹) are taken into account. Ignoring the smaller EPCs has very little effect on the results, since the electron–phonon self-energies depend quadratically on the deformation potential (see Eq. (6) in section “Quantum transport”). The intravalley electron-ZA phonon interactions (not included in Table 1) are neglected in our simulations. The validity of this approximation can be checked from Supplementary Tables I and II in the Supplementary Information, that report the scattering rates at the conduction band minima associated to each phonon mode. The data show that the intravalley scattering rates between electrons and ZA phonons are extremely small compared to all other scattering processes. Recently, it has been shown that in antimonene the electron-transverse acoustic (TA) phonon couplings are forbidden by symmetry for both intra- and intervalley scattering³⁸, in good agreement with our numerical calculation here.

Table 1 Conduction band effective masses (m_∥ and m_⊥), phonon frequencies (ω), and electron–phonon deformation potentials Ξ calculated from DFT and DFPT for all the phonon modes in antimonene and arsenene.

From Table 1, we can notice that some intervalley EPCs have deformation potential values almost an order of magnitude larger than those of the intravalley EPCs. Due to the relatively low phonon frequencies in antimonene and arsenene, those phonon modes can be thermally populated at room temperature, which leads to increased scattering of electrons through both emission and absorption of phonons. Therefore, the intervalley EPCs are expected to play a non-negligible role in the electronic transport and sensibly affect the MOSFET performance. This is confirmed by our quantum transport simulations, presented in the next section.

Quantum transport simulations with dissipation

Our aim is to provide quantitative and precise predictions of the optimum performance of these devices, when all the extrinsic sources of scattering are eliminated. Under these hypotheses, the electron–phonon interaction represents the dominant scattering mechanism at room temperature.

Based on the electronic model and the EPC models presented in the previous sections, we have built a dissipative quantum transport solver including all the elastic and inelastic electron–phonon scattering mechanisms, while still keeping the computational burden at the minimum level. More details on the quantum transport model are given in section “Quantum transport”.

We focus on n-type MOSFET devices. A sketch of the simulated device structure is drawn in Fig. 1c. We consider a double-gate MOSFET architecture with source and drain extensions chemically doped at a donor concentration N_D = 3.23 × 10¹³ cm⁻² and gate lengths L_G ranging from 4 to 9 nm. The channel is undoped and has the same length of the gates. The effective oxide thickness is chosen to be 0.42 nm (corresponding to a 2 nm thick HfO₂ high-κ dielectric layer), and the supply voltage is V_DD = 0.55 V, according to the ITRS specifications³⁹.

We show in Fig. 2 the transfer characteristics of armchair-antimonene-based MOSFET with L_G = 5 nm. To illustrate the role played by the different electron–phonon scattering channels we carry out the simulations in three cases: (1) including only the intravalley acoustic scattering (elastic), (2) including only the intravalley scattering (elastic and inelastic), and (3) including all the intravalley and intervalley scattering (elastic and inelastic).

For a clearer comparison of different calculations and following the common practice, we shift all the transfer characteristics in order to set the OFF state, chosen to be the device state in which I_DS = I_OFF ≡ 0.1 A m⁻¹ at V_GS = 0, following the ITRS³⁹. Accordingly, the ON state is at V_GS = V_DD = 0.55 V, and the corresponding current defines the ON current (I_ON). The subthreshold swing (SS) is defined as the inverse slope of the I_DS–V_GS curve in semi-logarithmic scale in the subthreshold regime (I_DS < 1 A m⁻¹), SS = ${[\partial {{\rm{log}}}_{10}({I}_{{\rm{DS}}})/\partial {V}_{{\rm{GS}}}]}^{-1}$.

As shown in Fig. 2 for an antimonene-based MOSFET with L_G = 5 nm, the SS is unaffected by the EPC. However, the ON state performance is strongly degraded, especially by the intervalley scatterings, consistently with the high amplitude of the associated deformation potentials. In quantitative terms, when all the significant scattering processes are included, the ON current drops by a factor larger than 4 with respect to the ballistic or ADP scattering approximation. Similar results are obtained in the case of the MOSFET based on arsenene and can be found in the Supplementary Information. More in general, these data strongly suggest that modeling the transport in 2D-material-based transistors as elastic, either in the ballistic or in the ADP scattering approximation, is not sufficient to reliably predict the device performance. On the contrary, inelastic intravalley and intervalley scattering mechanisms dominate and significantly affect the transport, even at gate lengths as short as 5 nm.

To elucidate the effects of EPC on the charge transport all along the device, we plot in Fig. 3 the current spectrum in the OFF and in the ON states for L_G = 5 nm. We consider both the case in which all the relevant scattering processes are included (panels (a) and (b)) and the case in which only the elastic scattering between electrons and acoustic phonons is taken into account (panels (c) and (d)). We first compare the current spectra in the OFF state. In Fig. 3a, it can be seen that the injection energies associated to the inelastic current spectrum gather around the source Fermi level (set to 0 eV). However, close to the barrier the spectrum broadens over a wider energy window, similar to that over which the elastic current concentrates (Fig. 3c). This effect is driven by the significant optical phonon absorption, which rises electrons at energies for which the tunneling probability through the barrier is higher. The net result is that the barycenter of the inelastic current spectrum crossing the barrier is close to that observed in elastic transport conditions. Moreover, due to the smallness of the charge density in the channel in the subthreshold region, the conduction band profile is unaffected by the transport and only determined by the device geometry and by the dielectric properties of the materials (see the Supplementary Information for related data). Particularly, the efficiency of the gate in modulating the barrier turns out to be independent of the electron–phonon scattering processes taken into account. Since the SS is completely determined by the energy distribution of carriers and by the gate modulation efficiency, the previous considerations explain the closeness between the SS values in the elastic and inelastic transport cases. This argument holds respectless of the considered gate length (see the Supplementary Information for the data at L_G = 9 nm).

Fig. 3: Current density spectrum in an Sb-based MOSFET with L_G = 5 nm.

a, c in the OFF state and b, d in the ON state, considering a, b inter- and intravalley acoustic and optical scattering altogether, c, d only considering intravalley elastic acoustic scattering. The black dashed lines indicate the conduction band profile. The Fermi level at the source contact is set at 0 eV, while the Fermi level at the drain contact is set at −0.55 eV. The current density spectrum is in arbitrary unit.

We now focus on the ON state. The current spectrum in Fig. 3b and d, shows that, contrarily to the behavior in the OFF state, electrons are now mainly injected at higher energies with respect to the top of the channel barrier. The current spectrum in Fig. 3b, obtained by including the inelastic scattering processes, demonstrates that electrons gradually lose energy by emitting optical phonons. The strong emission processes entail a significant backscattering and explains the considerable decrease of the current with respect to the case in which only elastic electron–phonon interactions are considered.

Figure 4 summarizes the simulated transfer characteristics of the antimonene- and arsenene-based MOSFET for L_G varying from 4 to 9 nm, and the transport axis along the armchair and zigzag directions. All the relevant electron–phonon scattering processes are included. In all the considered cases the antimonene- and arsenene-based devices exhibit very similar transfer characteristics. The extracted SS and the ON current for the antimonene-based MOSFET are reported in Fig. 5. The very similar data obtained for the arsenene-based transistor can be found in the Supplementary Information. We observe that the armchair and zigzag transport directions are not equivalent, since the rotation by π/4 needed to switch between them is not a symmetry of the crystal. Particularly, while the density of states and the scattering rates associated to the six degenerate valleys are invariant under arbitrary rotations, the anisotropy of the valleys makes the average effective mass seen in the armchair and zigzag directions different. The higher ON current in the zigzag direction arises from a smaller average transport effective mass $\langle {m}_{t}^{* }\rangle$, and thus a higher group velocity. On the other hand, the smaller $\langle {m}_{t}^{* }\rangle$ entails a larger source-to-drain tunneling (STDT), which results in higher SS values.

Fig. 4: Calculated transfer characteristic of Sb and As based MOSFETs.

I_DS–V_GS in linear and logarithmic scale along armchair and zigzag directions, and for L_G = 4, 5, 7, and 9 nm.

Fig. 5: Figures-of-merit of the Sb-based MOSFETs.

ON current and SS as a function of the gate length L_G extracted from the transfer characteristics in Fig. 4 of the antimonene-based MOSFET.

Reducing phonon scattering by valley engineering

Since the six-fold degenerate conduction band valleys are inherited from the crystal rotational symmetry, applying an external symmetry-breaking uniaxial strain along the zigzag direction (see Fig. 1) can split the six valleys into four higher-energy and two lower-energy valleys on the armchair axis. On the contrary, if the uniaxial strain is along the armchair direction, only two valleys will be pushed up to higher energy³⁵. In the following, we will consider the uniaxial strain along the zigzag direction. A sizable splitting of ≈ 0.16 eV can be obtained for a small uniaxial strain of 2% in the case of arsenene, while the splitting is ≈ 0.1 eV in antimonene³⁵. The effect on EPC and effective masses by such small strain has been shown to be negligible³⁵. Instead, since the strain shifts the energy of the valleys, it effectively turns off the intervalley scattering channels between the lower-energy and the higher-energy valleys. In addition, the transport becomes anisotropic, since the two lower-lying valleys both have the large effective mass in the armchair direction and the small one in the zigzag direction.

Figure 6 compares the transfer characteristics of the strained-arsenene-based MOSFET along the armchair and zigzag directions, for L_G = 4 and 9 nm. The transfer characteristics of the unstrained device are also plotted for reference. The simulations clearly shows an improvement of the ON current for all the considered configurations. Such an improvement is in excess of a factor 2.5 in the zigzag direction for L_G = 9 nm and amounts to approximately a factor 2 in the armchair direction for L_G = 4 nm. Similar improvements (a factor of 2.4 in the zigzag direction for L_G = 9 nm and 2 in the armchair direction for L_G = 4 nm) are obtained for MOSFETs based on strained antimonene. The corresponding data are reported in the Supplementary Information.

Fig. 6: Calculated transfer characteristics for MOSFETs based on As with an uniaxial strain of 2% along the zigzag direction.

I_DS–V_GS in linear and logarithmic scale. Two different gate lengths (L_G = 4 and 9 nm) and transport directions (armchair and zigzag) are considered.

The SS and I_ON for both arsenene- and antimonene-based MOSFETs are quantified in Fig. 7 for L_G = 4, 5, 7, and 9 nm. The plots shows that, while the transistors fabricated by using the two monolayers exhibit almost identical SSs as a function of L_G, that based on arsenene can achieve higher values of I_ON. In both cases, at short channel lengths (L_G < 5 nm) the best performance are obtained when the transport occurs along the armchair direction. Indeed, due to the higher effective mass along this direction, the I_ON is less degraded by the STDT, which results in a better scaling behavior. On the contrary, at longer gate lengths, the transport along the zigzag direction is more advantageous, due to the higher injection velocity associated to the small effective mass. We remark that also the low m_⊥/m_∥ ratio plays a significant role in enhancing the carrier velocity, since at each given energy it entails a larger kinetic component in the transport direction.

Fig. 7: Figures-of-merit of the MOSFETs with strained Sb and As.

ON current and SS as a function of the gate length L_G extracted from the transfer characteristics in Fig. 6 of the MOSFETs based on a strained antimonene (Sb) and b arsenene (As).

Discussion

We have investigated the dissipative transport in single-layer arsenene and antimonene MOSFETs by performing first-principle-based NEGF quantum simulations including all the relevant intra- and intervalley electron–phonon scattering mechanisms. We have shown that even for a gate length of 5 nm the ON current can be strongly decreased by the optical intravalley and intervalley phonon scatterings, and that, as a consequence, the simulations in the ballistic regime and/or only including the acoustic phonon scattering can largely overestimate the performance of real devices. Due to the combination of the six-fold degeneracy of the conduction band valleys, the large magnitude of the electron–phonon couplings, and the low optical phonon frequencies, the intervalley electron–phonon scattering is the dominant scattering mechanism in the arsenene- and antimonene-based MOSFETs at room temperature.

Moreover, we have investigated a viable approach to mitigate the impact of intervalley phonon scattering by removing the valley degeneracy via a small uniaxial strain along the zigzag direction. Our calculations indicate that for gate lengths larger than 5 nm, the ON current in the zigzag direction can be increased by a factor 2.5. For gate lengths shorter than 5 nm, selecting the armchair direction as the transport direction minimizes the source-to-drain tunneling effect and results in a two times ON-current improvement with respect to the case with unstrained material.

Overall, our work can provide useful guidelines to the simulation of transistors based on 2D materials, and suggests that there are room and opportunities to overcome the obstacles on the way towards the development of a future 2D-based electronics.

Methods

The anisotropic effective mass model

The conduction band of Sb and As monolayers has six-fold degenerate valleys located at the midpoint of the Γ-M paths. These valleys have a larger effective mass along the high-symmetry Γ-M direction, denoted by m_∥. The effective mass along the direction perpendicular to Γ-M is denoted by m_⊥. By using local k-coordinates for each valley, we can write the conduction band energy in the anisotropic effective mass model as

$$\frac{{\hslash }^{2}}{2{m}_{0}}\left(\frac{{k}_{\perp }^{2}}{{m}_{\perp }}+\frac{{k}_{\parallel }^{2}}{{m}_{\parallel }}\right)=E,$$

(1)

where ${{\hslash}}$ is the reduced Planck constant, m₀ the electron mass, and k_⊥, k_∥ the wavevector components along the Γ-M and the orthogonal direction, respectively. For transport calculations, we need to express k_⊥ and k_∥ in terms of k_x and k_y, namely the wave vector components in the transport and transverse directions, respectively (see Fig. 1d):

$$\left\{\begin{array}{lll}{k}_{\parallel }&=&\cos \theta {k}_{x}+\sin \theta {k}_{y}\\ {k}_{\perp }&=&-\sin \theta {k}_{x}+\cos \theta {k}_{y}\end{array}\right.$$

(2)

where θ denotes the angle between the k_∥-axis and the transport direction x.

The kinetic energy of electrons in terms of k_x and k_y reads:

$$\begin{array}{ll}\frac{{\hslash }^{2}}{2{m}_{0}}&\left[{k}_{x}^{2}\left(\frac{\cos {\theta }^{2}}{{m}_{\parallel }}+\frac{\sin {\theta }^{2}}{{m}_{\perp }}\right)+{k}_{y}^{2}\left(\frac{\cos {\theta }^{2}}{{m}_{\perp }}+\frac{\sin {\theta }^{2}}{{m}_{\parallel }}\right)\right.\\ &\left.+2{k}_{x}{k}_{y}\sin \theta \cos \theta \left(\frac{1}{{m}_{\parallel }}-\frac{1}{{m}_{\perp }}\right)\right]=E.\\ \end{array}$$

(3)

The Hamiltonian in the form employed in the simulations is obtained by replacing k_x with the differential operator − i∂/∂x. We enforce periodic boundary conditions in the transverse direction, thus the Hamiltonian maintains a parametric dependence on k_y. The value of θ to be considered for each valley depends on the chosen transport direction (armchair or zigzag).

Quantum transport

To simulate the electron transport, we adopt the Keldysh-Green’s function formalism⁴⁰, which allows us to include the electron–phonon coupling by means of local self-energies. We self-consistently solve the following equations for every valley

$$\left\{\begin{array}{l}(E{\bf{I}}-{\bf{H}}-{{\boldsymbol{\Sigma }}}_{{\rm{SD}}}-\left.{{\boldsymbol{\Sigma }}}_{{\rm{ph}}}\right]){\bf{G}}={\bf{I}}\\ {{\bf{G}}}^{\lessgtr }={\bf{G}}({{\boldsymbol{\Sigma }}}_{{\rm{SD}}}^{\lessgtr }+{{\boldsymbol{\Sigma }}}_{{\rm{e-ph}}}^{\lessgtr }){{\bf{G}}}^{\dagger }\end{array}\right.$$

(4)

where E is the electron energy, I is the identity matrix, H is the anisotropic effective mass Hamiltonian matrix, Σ_SD and Σ_e-ph are the retarded self-energies corresponding to source and drain contacts and to the electron–phonon coupling, G is the retarded Green’s function, and ${{\boldsymbol{\Sigma }}}_{{\rm{SD}}}^{\lessgtr }$, ${{\boldsymbol{\Sigma }}}_{{\rm{e-ph}}}^{\lessgtr }$ and G^≶ are the corresponding lesser-than and greater-than quantities. Equations (4) are spatially discretized over a uniform mesh with stepsize equal to 0.5 Å.

For the intravalley acoustic phonons, the local self-energies at the i-th discrete space site along the transport direction at room temperature can be expressed in the quasi-elastic approximation as⁴¹

$${{\boldsymbol{\Sigma }}}_{{\rm{ac}},{\rm{\nu }}}^{\lessgtr }(i,i;{k}_{y},E)=\frac{{{{\Xi }}}_{{\rm{ac}}}^{2}{k}_{{\rm{B}}}T}{{C}_{{\rm{2D}}}}\mathop{\sum}\limits_{{q}_{y}}{{\bf{G}}}_{\nu }^{\lessgtr }(i,i;{k}_{y}-{q}_{y},E),$$

(5)

where C_2D is the elastic modulus, Ξ_ac the acoustic deformation potential, k_B the Boltzmann constant, T the temperature, q_y the transverse phonon wavevector and ν is the valley index.

For the intravalley optical phonons, the self-energies related to a phonon with branch λ can be expressed as

$$\begin{array}{l}{{\boldsymbol{\Sigma }}}_{{\rm{\lambda }}{{\Gamma }},\nu }^{\lessgtr }(i,i;{k}_{y},E)=\frac{\hslash {{{\Xi }}}_{\lambda {{\Gamma }}}^{2}}{2\rho {\omega }_{\lambda }^{{{\Gamma }}}}\mathop{\sum}\limits_{{q}_{y}}\left\{{{\bf{G}}}_{\nu }^{\lessgtr }(i,i;{k}_{y}-{q}_{y},E\mp \hslash {\omega }_{\lambda }^{{{\Gamma }}})N(T,{\omega }_{\lambda }^{{{\Gamma }}})\right.\\ \qquad\qquad\qquad\qquad+\,\left.{{\bf{G}}}_{\nu }^{\lessgtr }(i,i;{k}_{y}-{q}_{y},E\pm \hslash {\omega }_{\lambda }^{{{\Gamma }}})\left[N(T,{\omega }_{\lambda }^{{{\Gamma }}})+1\right]\right\},\end{array}$$

(6)

where ρ is the mass density, $N(T,{\omega }_{\lambda }^{{{\Gamma }}})$ is the Bose–Einstein distribution at temperature T of a phonon with frequency ${\omega }_{\lambda }^{{{\Gamma }}}$ at Γ point, and the upper/lower sign of the term $\hslash {\omega }_{\lambda }^{{{\Gamma }}}$ corresponds to lesser/greater-than self-energies.

The self-energies associated to the intervalley scattering are analogously obtained as:

$$\begin{array}{l}{{\boldsymbol{\Sigma }}}_{{\lambda}{{\boldsymbol{q}},\nu }}^{\lessgtr } (i,i;{k}_{y},E)=\frac{\hslash {{{\Xi }}}_{\lambda {\boldsymbol{q}}}^{2}}{2\rho {\omega }_{\lambda }^{{\boldsymbol{q}}}}\mathop{\sum}\limits_{\nu ^{\prime} \ne \nu }\mathop{\sum}\limits_{{q}_{y}}\left\{{{\bf{G}}}_{\nu ^{\prime} }^{\lessgtr }(i,i;{k}_{y}-{q}_{y},E\mp \hslash {\omega }_{\lambda }^{{\boldsymbol{q}}})N(T,{\omega }_{\lambda }^{{\boldsymbol{q}}})\right.\\ \left.\qquad\qquad\qquad\qquad\,+\,{{\bf{G}}}_{\nu ^{\prime} }^{\lessgtr }(i,i;{k}_{y}-{q}_{y},E\pm \hslash {\omega }_{\lambda }^{{\boldsymbol{q}}})\left[N(T,{\omega }_{\lambda }^{{\boldsymbol{q}}})+1\right]\right\},\end{array}$$

(7)

where q is the phonon wavevectors connecting two degenerate valley minima, and the three scattering processes shown in Fig. 1 are considered.

The total lesser-than and greater-than self-energies ${{\boldsymbol{\Sigma }}}_{{\rm{e-ph}}}^{\lessgtr }$ are obtained by summing the self-energies associated to all the different scattering processes. The retarded self-energy is then computed by

$${{\boldsymbol{\Sigma }}}_{{\rm{e-ph}}}^{}(i,i;{k}_{y},E)=\frac{1}{2}[{{\boldsymbol{\Sigma }}}_{{\rm{e-ph}}}^{ \,{>}\,}(i,i;{k}_{y},E)-{{\boldsymbol{\Sigma }}}_{{\rm{e-ph}}}^{\,{<}\,}(i,i;{k}_{y},E)],$$

(8)

The real part of the retarded self-energy, which only contributes to a shift of the energy levels, is neglected.

Equations ((4)–(8)) are iteratively solved within the so-called self-consistent Born approximation (SCBA) until the convergence is achieved. The SCBA convergence ensures that the electronic current is conserved along the device structure. Once the SCBA loop achieves the convergence, the electron and transport current densities are determined from the Green’s functions of all the valleys

$$n(i)=\frac{-i}{2\pi }\mathop{\sum}\limits_{\nu ,\ {k}_{y}}\mathop{\int}\nolimits_{{E}_{{{\min}}}}^{{E}_{{{\max}}}}{\rm{d}}E{{\bf{G}}}_{\nu }^{\,{<}\,}(i,i;{k}_{y},E),$$

(9)

and

$$\begin{array}{l}{J}_{x}(i)=\frac{e}{h}\mathop{\sum}\limits_{\nu ,\ {k}_{y}}\mathop{\int}\nolimits_{{E}_{{{\min}}}}^{{E}_{{{\max}}}}{\rm{d}}E\left[{{\bf{H}}}_{i,i+1}{{\bf{G}}}_{\nu }^{<}(i+1,i;{k}_{y},E)\right.\\ \qquad\quad\,\,-\left.{{\bf{H}}}_{i+1,i}{{\bf{G}}}_{\nu }^{<}(i,i+1;{k}_{y},E)\right],\end{array}$$

(10)

where the integration range from E_min to E_max includes all the energy states contributing to transport. We discretize the energy into a uniform grid with a stepsize ΔE = 4 meV for antimonene and ΔE = 5 meV for arsenene. Each energy interval of width ΔE is further discretized by defining 3 quadrature nodes. The charge and current densities are computed by first performing an integration over each interval by using the Gauss-Legendre quadrature rule and then adding the contributions of all the intervals.

Due to the extreme thinness of the As and Sb monolayers, the electrostatics of the MOSFETs is expected to be strongly dominated by the dielectric environment. Based on this consideration, in the Poisson equation we modeled the monolayers as ideal surfaces and did not account for their specific dielectric properties. Accordingly, the mobile charge density is treated as a 2D charge, distributed at the interface between the top and bottom gate oxide layers. This 2D charge density was used to set up the 2D nonlinear Poisson equation

$$\nabla \cdot \varepsilon ({\bf{r}})\nabla \phi ({\bf{r}})=-e[n({\bf{r}})-{N}_{D}({\bf{r}})]$$

(11)

where r = (x, z), ε(r) is the oxide dielectric constant, ϕ(r) is the electrostatic potential and N_D(r) is the donor concentration. Equation (11) was discretized in the x-z cross-section with a finite-element method by assuming Dirichlet boundary conditions at the metal-gate nodes and von Neumann’s conditions otherwise, then solved with a Newton-Raphson algorithm. The simulation runs over a loop between the Green’s functions and the Poisson equations until an overall self-consistent solution is attained.

First-principles calculation details

The DFT simulations were performed in a plane wave basis by means of the Quantum Espresso suite^42,43. The antimonene and arsenene monolayers were simulated by using scalar relativistic norm-conserving pseudopotentials based on the local density approximation exchange-correlation functional. In both cases, a 16 × 16 × 1 Monkhorst-Pack (M-P) k-point grid was used, with a cutoff energy of 90 Ry. A vacuum layer of 28 Å was considered along the perpendicular direction to eliminate the interlayer interactions due to the periodic boundary conditions.

The phonon dispersion relations were computed within the DFPT framework by using a 8 × 8 × 1 M-P q-point grid and a 16 × 16 × 1 M-P k-point grid.

The electron–phonon matrix elements between two electronic states $\left|n{\boldsymbol{k}}\right\rangle$ and $\left|m{\boldsymbol{k}}+{\boldsymbol{q}}\right\rangle$ connected by a phonon mode λ of wavevector q were computed within the DFTP framework as

$${g}_{{\boldsymbol{k}}{\boldsymbol{q}},mn\lambda }=\sqrt{\hslash /(2M{\omega }_{\lambda {\boldsymbol{q}}})}\left\langle m{\boldsymbol{k}}+{\boldsymbol{q}}\right|{{\Delta }}{V}_{{\boldsymbol{q}}\lambda }^{{\rm{KS}}}\left| n{\boldsymbol{k}}\right\rangle,$$

(12)

where M is the unit cell mass, m and n are band indices, and ${{\Delta }}{V}_{{\boldsymbol{q}}\lambda }^{{\rm{KS}}}$ is the lattice periodic part of the perturbed Kohn–Sham potential expanded to first order in the atomic displacement. The g_kq,mnλ were computed on a 16 × 16 × 1 k-point grid and a 8 × 8 × 1 q-point grid, and then interpolated on a denser grid by projection on a basis of maximally localized Wannier functions by using the wannier90 code⁴⁴ and the EPW code⁴⁵. Six Wannier orbitals were employed to reproduce accurately the conduction and valence bands close to the bandgap and achieve good spatial localization of the Wannier orbitals. The scattering rates for each phonon mode were computed through the EPW code, by interpolating the electron–phonon matrix elements on a 120 × 120 × 1 q-point grid. The intravalley and intervalley electron–phonon coupling were described by the matrix element at q = 0 and ${\boldsymbol{q}}={\boldsymbol{k}}^{\prime} -{\boldsymbol{k}}$, respectively, where ${\boldsymbol{k}}^{\prime}$ and k correspond to two degenerate valley minima. The deformation potentials were finally computed as

$${{{\Xi }}}_{\lambda {\boldsymbol{q}}}=| {g}_{{\boldsymbol{k}}{\boldsymbol{q}},mn\lambda }| \times \sqrt{2M{\omega }_{\lambda {\boldsymbol{q}}}/\hslash }.$$

(13)