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

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 15 and their relatively low mobility.
Experimental measurements of room temperature electron mobilities in 2D materials have so far reported values lower than 400 cm 2 V −1 s −1 , 9,[16][17][18][19] much smaller than the value for bulk Si (about 1400 cm 2 V −1 s −1 ) 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 24 . 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 20 . The low phononlimited 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 32 , which takes into account only the Bardeen-Shockley acousticdeformation-potential (ADP) scattering 33 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 shortchannel 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 firstprinciples 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. 35 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.

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. 13 . 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. 13 .

Electron-phonon couplings from first-principles
We compute the relevant electron-phonon couplings from firstprinciples 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 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 singlelayer 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 13 cm −2 . 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 2 ) 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. 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 7 eV cm −1 ) 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 electrontransverse acoustic (TA) phonon couplings are forbidden by symmetry for both intra-and intervalley scattering 38 , in good agreement with our numerical calculation here.
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.   The values of C 2D and Ξ Γ AC are taken from ref. 13 The intervalley EPCs with a deformation potential >4 × 10 7 are marked in bold.

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 13 cm −2 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 2 high-κ dielectric layer), and the supply voltage is V DD = 0.55 V, according to the ITRS specifications 39 .
We show in Fig. 2 the transfer characteristics of armchairantimonene-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 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).
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.   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 hm Ã t i, and thus a higher group velocity. On the other hand, the smaller hm Ã t i entails a larger source-to-drain tunneling (STDT), which results in higher SS values.

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 higherenergy 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 35 . 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 35 . The effect on EPC and effective masses by such small strain has been shown to be negligible 35 . 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 strainedarsenene-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.
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.

DISCUSSION
We have investigated the dissipative transport in single-layer arsenene and antimonene MOSFETs by performing first-principlebased NEGF quantum simulations including all the relevant intraand 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 ONcurrent 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.

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 where _ is the reduced Planck constant, m 0 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): 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: 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 40 , 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 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 selfenergies corresponding to source and drain contacts and to the  electron-phonon coupling, G is the retarded Green's function, and Σ ≶ SD , Σ ≶ eÀph 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 41 Σ ≶ ac;ν ði; i; k y ; EÞ ¼ 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 where ρ is the mass density, NðT; ω Γ λ Þ is the Bose-Einstein distribution at temperature T of a phonon with frequency ω Γ λ at Γ point, and the upper/ lower sign of the term _ω Γ λ corresponds to lesser/greater-than selfenergies.
The self-energies associated to the intervalley scattering are analogously obtained as: 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 Σ ≶ eÀph are obtained by summing the self-energies associated to all the different scattering processes. The retarded self-energy is then computed by Σ eÀph ði; i; k y ; EÞ ¼ 1 2 ½Σ > eÀph ði; i; k y ; EÞ À Σ < eÀph ði; i; k y ; EÞ; 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 selfconsistent 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Þ ¼ Ài 2π X ν; ky Z Emax Emin dEG < ν ði; i; k y ; EÞ; and 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 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 exchangecorrelation functional. In both cases, a 16 × 16 × 1 Monkhorst-Pack (M-P) kpoint 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 kpoint grid.
The electron-phonon matrix elements between two electronic states nk j i and mk þ q j iconnected by a phonon mode λ of wavevector q were computed within the DFTP framework as where M is the unit cell mass, m and n are band indices, and ΔV KS qλ 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 44 and the EPW code 45 . 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 q ¼ k 0 À k, respectively, where k 0 and k correspond to two degenerate valley minima. The deformation potentials were finally computed as Ξ λq ¼ jg kq;mnλ j ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2Mω λq =_ q :

DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author upon reasonable request.