Marco Pala
email: marco.pala@c2n.upsaclay.fr

David Esseni
email: david.esseni@uniud.it

Ab-initio quantum transport with a basis of unit-cell restricted Bloch functions and the NEGF formalism

Keywords: Density functional theory, quantum transport, NEGF, 2D materials

Ab-initio quantum transport with a basis of unit-cell restricted Bloch functions and the NEGF formalism

Marco G. Pala,

I. INTRODUCTION

In recent years many innovations in nanoelectronics devices have been driven by the exploitation of new materials, such as atomically thin 2D materials and their heterostructures [START_REF] Mingda | Single Particle Transport in Two-dimensional Heterojunction Interlayer Tunneling Field Effect Transistor[END_REF]- [START_REF] Marin | Lateral heterostructure field-effect transistors based on two-dimensional material stacks with varying thickness and energy filtering source[END_REF]. Ab-initio methods based on the DFT have been shown to be quite effective to investigate the electronic and optical properties of these new materials, but their direct use for quantum transport calculations has been hindered by the size and complexity of their Hamiltonians. Remarkable exceptions are the contributions employing Hamiltonians based on the linear combination of local orbitals (LCAO) [START_REF] Soler | The SIESTA method for ab initio order-n materials simulation[END_REF], or on the maximally localized Wannnier functions (MLWF) [START_REF] Marzari | Maximally localized generalized Wannier functions for composite energy bands[END_REF], [START_REF] Szabó | Ab initio simulation of singleand few-layer MoS 2 transistors: Effect of electron-phonon scattering[END_REF].

Here we discuss an alternative approach to ab-initio quantum transport [START_REF] Pala | Unit cell restricted Bloch functions basis for first-principle transport models: Theory and application[END_REF], that leverages directly on the Hamiltonian matrix and Bloch functions obtained with plane-wave DFT solvers [START_REF] Kresse | Efficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set[END_REF], [START_REF] Giannozzi | Quantum Espresso: a modular and open-source software project for quantum simulations of materials[END_REF]. We believe that such a methodology, differently from the widely used MLWF approach, can be naturally integrated into the workflow of an ab-initio solver, and thus provide a fairly automatic pathway to ab-initio quantum transport simulations.

II. METHODOLOGY

In order to obtain the DFT Hamiltonian, we start by selfconsistently solving the Kohn-Sham (KS) equations

H KS Ψ n = E n Ψ n , H KS = T + V scf (1)
where the KS Hamiltonian H KS is the sum of the operators corresponding to the kinetic energy T and to the self-consistent potential

V scf = V eI + V H + V xc
, where V eI is the electronion interaction potential, V H the Hartree electrostatic potential, V xc the exchange-correlation potential. To this goal we use a PW basis set and pseudopotentials to represent the valence electron-nuclei interactions. The KS Hamiltonian is then written as

H(k + G, k + G ) = 2 2m (k + G) 2 δ G,G + V L (G -G ) + V N L (k + G, k + G ) (2)
where V L (G -G ) is the Fourier transform of the local part of the total potential and V N L (k + G, k + G ) is the spatially non-local contribution of the pseudopotential, that for most atoms is short ranged and vanishes for |r-r |>r c , with r c being the radius beyond which pseudo and true KS orbitals coincide [START_REF] Giannozzi | Advanced capabilities for materials modelling with Quantum Espresso[END_REF].

For transport calculations it is convenient to use an orthorombic unit cell [START_REF] Pala | Unit cell restricted Bloch functions basis for first-principle transport models: Theory and application[END_REF], where the real space unit vectors can be written as a 1 =(a x , 0, 0), a 2 =(0, a y , 0), a 3 =(0, 0, a z ), with x being the transport direction. Hence the unit vectors of the reciprocal lattice are b 1 =(2π/a x , 0, 0), b 2 =(0, 2π/a y , 0), b 3 =(0, 0, 2π/a z ). The Brillouin zone (BZ) is therefore given by the conditions -π/a s < k s ≤ π/a s (with s = x, y, z).

A. Discussion of the reduced basis set

The DFT Hamiltonian in Eq. ( 2) is a dense matrix with a rank equal to the number, N G , of reciprocal lattice vectors, and thus it cannot be directly employed in transport calculations. Thus our first step is a unitary transformation from PW to real space along the transport direction x, that provides the Hamiltonian in the hybrid xK yz basis given by real-space along x and plane waves in the (y, z) plane [START_REF] Pala | Full-band quantum simulation of electron devices with the pseudopotential method: Theory, implementation, and applications[END_REF]. In this representation the Hamiltonian matrix is a block tri-diagonal matrix, which is crucially important in order to exploit recursive algorithms for non-equilibrium Green's function (NEGF) calculations [START_REF] Anantram | Modeling of Nanoscale Devices[END_REF]. In this hybrid xK yz basis each block describes an a x long spatial region consisting of N dx discretization points, so that the matrix blocks have a rank N G =N dx N Gy N Gz (with N dx =N Gx ). Even if the xK yz basis has been used to perform quantum transport calculations as in [START_REF] Pala | Full-band quantum simulation of electron devices with the pseudopotential method: Theory, implementation, and applications[END_REF], [START_REF] Pala | NEGF based transport modeling with a full-band, pseudopotential Hamiltonian: Theory, Implementation and Full Device Simulations[END_REF], the size N G of the blocks is large and it increases by enlarging the cutoff energy used in DFT calculations. Huge computational advantages can thus be obtained by moving to a basis of unitcell restricted Bloch functions, that preserves the block tridiagonal structure of the Hamiltonian matrix [START_REF] Pala | Unit cell restricted Bloch functions basis for first-principle transport models: Theory and application[END_REF], [START_REF] Van De Put | Scalable atomistic simulations of quantum electron transport using empirical pseudopotentials[END_REF].

To this purpose we first argue that, due to the block tridiagonal form of the Hamiltonian, the unit cell restricted Bloch functions {Ψ 0 k } corresponding to the wave vector k=(k x ,k yz ) satisfy the secular equation

H † 0,1 e -ikxax + H 0,0 + H 0,1 e ikxax {Ψ 0 k } = E(k) {Ψ 0 k } (3
) with H 0,0 , H 0,1 being the Hamiltonian blocks describing respectively one unit cell and its interaction with the nearest neighbor cell. Our reduced basis set consists of a subset of the {Ψ 0 nk } corresponding to a few k x values in the reduced zone -π/a x <k x ≤π/a x along the transport direction and, for each k x , to some tens of bands E n (k) up to the energies of interest for the problem at study.

It is noteworthy that the {Ψ 0 nk } are obtained directly from the PW representation provided by the DFT solver, namely we do not need to solve Eq. ( 3) to find the basis functions.

The size of the basis can be written as

N B = n f n=ni N kB,n (4)
where N kB,n denotes the number of k x values included in the basis for the n-th band, with n ranging between n i and n f . The {Ψ 0 nk } functions are not orthogonal for different k x , but they can be orthogonalized and we denote by Φ the basis set orthonormalized over a unit cell. The Hamiltonian blocks in the Φ basis have a size N B that is typically much smaller than the size N G in the hybrid xK yz basis.

A first validation of the reduced basis set requires that it allows us to calculate the band structure for any k x in the reduced zone, which is discussed below in Sec.II-B.

B. Material system and validation of the reduced basis set

In this paper we exemplify our methodology by focusing on the monolayer PtSe 2 , which is a noble metal dichalcogenure identified as a promising channel material for MOSFETs due to the small effective mass and large density of states (DOS) [START_REF] Almutairi | PtSe 2 field-effect transistors: New opportunities for electronic devices[END_REF]. Figure 1 reports the corresponding band structure along the high-symmetry points of the primitive BZ obtained from the DFT Hamiltonian in the PW basis. The DFT calculation was performed with the Quantum ESPRESSO code [START_REF] Giannozzi | Quantum Espresso: a modular and open-source software project for quantum simulations of materials[END_REF], by using the Perdew-Burke-Ernzerhof [START_REF] Perdew | Generalized gradient approximation made simple[END_REF] approximation to the exchange-correlation functional and a norm-conserving pseudopotential [START_REF] Hamann | Optimized norm-conserving Vanderbilt pseudopotentials[END_REF]. We employed a 12 × 12 × 1 Monkhorst-Pack k-points grid and a cutoff energy of E w = 60 Ry. For transport calculations we first define the orthorombic cell composed of 6 atoms and sketched in Fig. 2. The unit vectors are a 1 =(a x , 0, 0), a 2 =(0, a y , 0), a 3 =(0, 0, a z ) with a x =3.712 Å, a y =6.429372 Å and a z =32 Å, while the vertical distance between the Pt and Se atoms is 1.312 Å. In order to validate the reduced Bloch state basis, we compared the bands computed with the Hamiltonian blocks in the reduced basis with the results of the plane-wave DFT solver. The accuracy of the reduced basis depends on the included number of k x Bloch functions [START_REF] Pala | Unit cell restricted Bloch functions basis for first-principle transport models: Theory and application[END_REF], hence we optimized the basis by using a larger N kB,n for bands with energies in the range of interest, and smaller N kB,n for bands at much lower or higher energies.

In Fig. 3 we show an excellent band-structure reconstruction obtained by using N kB,n =4 for the bands close to the energygap and N kB,n =2 for smaller energy bands, and for an overall number N B =74 of Bloch functions. Such a definition of the reduced basis results in no observable unphysical states [START_REF] Mil'nikov | Equivalent transport models in atomistic quantum wires[END_REF], [START_REF] Shin | Density functional theory based simulations of silicon nanowire field effect transistors[END_REF].

C. Carrier densities and electron current

All relevant physical quantities are defined in terms of the retarded(advanced), [G r(a) ] Φ , and the lesser(greater)-than Green's functions, [G <(>) ] Φ , computed in the reduced Bloch functions basis. In order to compute the carrier concentrations in the real space, we first compute [G <(>) ] Φ , then we reconstruct the Green's function in the hybrid basis according to the expansion

G <(>) (x j K yz , x j K yz ; E) = N B n,m=1 G <(>) (n, m; E) ×Φ n (x j , G yz ) Φ * m (x j , G yz ) (5)
and finally evaluate electron and hole densities as described in Ref. [START_REF] Pala | Unit cell restricted Bloch functions basis for first-principle transport models: Theory and application[END_REF]. Thanks to this procedure we are able to evaluate the 3D electron and hole concentrations, n(r) and p(r), with atomistic accuracy as shown in the example of Fig. 4. However, by leveraging the fact that the electrostatic potential varies over larger spatial scales compared to the discretization used in Eq. ( 5), in the self-consistent solution of the NEGF transport equations and the Poisson equation we evaluate n(r) and p(r) on a coarser mesh. For the 3D Poisson problem we assume Diriclect conditions at the metal boundaries and Neumann conditions elsewhere. The electron current is obtained within the elastic approximation as

I D = ge h T (E) [f D (E) -f S (E)] dE (6)
with g being the degeneracy and f S(D) the Fermi-Dirac function at the source (drain) contact. Figure 5 illustrates the bandstructure and transmission of an infinitely long monolayer PtSe 2 at k y = 0 calculated by using the reduced basis and showing that, as expected, the transmission at a given energy E equals the number of available bands.

The device under investigation is sketched in Fig. 6. It is a single-gate FET with doped source and drain access regions having a length L S,D ≈11.1 nm and an undoped channel region with L G ≈18.5 nm. The donor concentration in the source and drain region is N D =10 20 cm -3 . The equivalent oxide thickness (EOT) of the high-k dielectric layer is 0.65 nm. In the device width direction, y, we assume periodic boundary conditions, described by sampling the BZ along the k y direction with a step ∆k y =0.1 × 2π/a y . Figure 7 shows the transfer characteristics for different V DS values. In the sub-threshold voltage regime we observe an ideal sub-threshold swing of about 60 mV/dec, which stems from the excellent electrostatic integrity of the device due, in turn, to the ultra-thin channel and the aggressive EOT. The device exhibits also a high transconductance resulting in a high on-state current, that confirms the promising properties of this material for high-performance MOSFETs. Of course the transconductance and on-state current in Figure 7 should be considered as upper limit figures, that can be significantly degraded by series resistance and scattering [START_REF] Esseni | Nanoscale MOS Transistors[END_REF].

Quite interestingly, in Figure 7 we can also observe a nonmonotonic I D dependence on V DS , in fact the I D at large V T G increases with V DS for V DS up to 0.3 V, but then it decreases for larger V DS . This behavior has been already observed in simulated FETs with a monolayer MoS 2 channel [START_REF] Chang | Atomistic full-band simulations of monolayer MoS 2 transistors[END_REF], and it has been ascribed to peculiar features of the DOS. We verified that a similar explanation applies to our results, in fact Fig. 8 reports the local density of states (LDOS) at V DS =0.3 V and 0.6 V in the on-state regime (i.e. at V T G =0.8 V), showing that in the drain region the LDOS for energies close to the top of the barrier is larger at V DS =0.3 V than it is at V DS =0.6 V.

Such a non-monotonic I D dependence on V DS , however, is expected to disappear in the presence of significant inelastic scattering mechanisms, such as electron-phonon scattering.

III. CONCLUSIONS AND OUTLOOK

We have shown that a unit-cell restricted Bloch functions basis enables band structure and transport calculations with a significant reduction of the computational burden compared to the original plane-wave DFT formulation. This methodology is a viable approach for ab initio and semi-empirical quantum transport simulations, and it can be seen as an alternative to the MLWF approach. In fact the Bloch functions can be obtained directly from the DFT solvers, thus bypassing the extraction of a MLWF basis and of the corresponding tightbinding Hamiltonian. We foresee that further developments of the methodology will enable us to deal with non homogenous systems (i.e. including hetero-jnctions), and to include inelastic scattering mechanisms.

Fig. 1 .

 1 Fig. 1. Bandstructure of the monolayer PtSe 2 .

Fig. 2 .

 2 Fig. 2. Orthorombic unit cell of the 1T phase monolayer PtSe 2 .

Fig. 3 .

 3 Fig. 3. Bandstructure of the monolayer PtSe 2 for various ky values computed with (circles) the Bloch function basis and (lines) the PW basis.

Fig. 4 .

 4 Fig. 4. Contourplot of the cross section of the electron density at an abscissa x corresponding to the center of the orthorombic unit cell. Colors are in arbitrary units.

Fig. 5 .

 5 Fig. 5. Bandstructure of an infinitely-long monolayer PtSe 2 and the corresponding transmission probability obtained with the Hamiltonian in the reduced basis consisting of unit-cell-restricted Bloch states. ky = 0.

Fig. 6 .

 6 Fig. 6. Sketch of the single-gate monolayer PtSe 2 FET.

Fig. 7 .

 7 Fig. 7. Transfer characteristics of the monolayer PtSe 2 FET for V DS =0.1, 0.3 and 0.6 V.

Fig. 8 .

 8 Fig. 8. LDOS computed for (a) V DS =0.3 V and (b) V DS =0.6 V and for V T G = 0.8 V. The transverse wave-vector is ky = 0.1 × 2π/ay.


M G Pala

P Giannozzi

D Esseni

Unit cell restricted Bloch functions basis for first-principle transport models: theory and application

Unit cell restricted Bloch functions basis for first-principle transport models: Theory and application

I. INTRODUCTION

A quantum transport methodology relying on an ab initio description of the physical system is in several respects the frontier of the transport modelling in nanoscale systems. One reason why such an approach has become necessary to steer the technological developments is that the cross section of many devices has reached truly nanometric dimensions and the transistor length has reached the 10 nm range, where quantum transport effects become important, such as the source-drain tunnelling in MOSFETs 1-3 , or the band-to-band-tunnelling (BTBT) in Tunnel FETs (TFETs) [START_REF] Seabaugh | Proceedings of the IEEE[END_REF][START_REF] Esseni | [END_REF] . Moreover, a quantum transport formalism based on ab initio methods has become indispensable to explore the potentials of new device concepts exploiting the recently discovered atomically thin 2D materials, and their many possible combinations in terms of vertical or lateral hetero-junction options [6][7][8] .

The first-principle electronic-structure calculations are typically based on DFT and on either a plane-wave basis 9,10 , or LCAO 11 . The former basis may be considered the most natural option for periodic crystals, while the latter is closely related to the chemical bonding picture.

The tight-binding method is the most widely used approximated implementation of the Linear Combination of atomic orbitals (LCAO) approach. It is also the most popular method for quantum transport based either on a fixed set of orbitals per atom with empirical coupling parameters (e.g. the sp 3 d 5 s * model 12,13 ), or on maximally localized Wannier functions extracted as a post-processing, sometimes quite delicate and computationally demanding, of first-principle calculations [14][15][16][17][18] .

Plane waves form a complete set of orthogonal functions. They allow for a good control of accuracy and convergence in electronic-structure calculations through a cutoff of the kinetic energy (see also Sec. II). However, frequently they result in a large basis set, particularly for those super-cells that include vacuum regions, where a good description of the exponential wave-function decay demands a large plane-waves set. Consequently, a direct use of first-principle calculations based on plane waves is often considered computationally prohibitive for electronic transport in technologically relevant systems. Some contributions have recently appeared for the empirical pseudopotentials method and using either a quantum-transmitting-boundary approach [19][20][21][START_REF] Fang | Electron Devices[END_REF] , or a Non Equilibrium Green's Function (NEGF) method 23,24 . Until now, DFT-based transport calculations in relatively large systems have been addressed only by using LCAO basis 25 with the adoption of equivalent transport model techniques 26 .

In this work we present a new method for quantum transport in nanoscale devices and physical systems, based on a plane-wave DFT Hamiltonian. The new method employs a basis set of Bloch functions of the underlying system to drastically reduce the size of the transport problem. Our approach does not require the solution of any eigenvalue problem besides those addressed by first-principle calculations 27 , in fact the basis set is obtained directly from the Bloch wave-functions determined by the ab initio solver. An appropriate choice of Bloch functions allows us to effectively avoid the problem of unphysical solutions, whose filtering can be theoretically and computationally challenging 25,26 . We found that the size of the basis set for transport simulations is essentially independent upon the cutoff energy used in first-principle calculations, which is extremely beneficial because it allows one to decouple the size of the transport problem from the computational effort necessary to obtain full convergence and high accuracy in first-principle calculations.

Our results demonstrate that the Bloch functions form an extremely effective basis set, which enables band structure and transport calculations using a basis size that is hundreds times smaller than the plane-wave basis used to calculate the Bloch functions in the ab initio solver. The ability of the Bloch functions to retain most of the physics with a small basis set is not surprising in consideration of the results obtained for band structure calculations with empirical pseudopotential models [28][29][30][31][32] , and it is a precious asset for future developments of quantum transport methods based on a first-principle Hamiltonian. During the writing of this paper we became aware of a recent contribution 33 , where a basis of Bloch functions was used for quantum transport simulations based on an empirical pseudopotential Hamiltonian and a quantum transmitting boundary approach.

The paper is organized as follows. In Sec. II we provide the necessary information about the first-principle methods employed in our calculations, and in particular we clarify the relevant connections to the transport model. In Sec. III we introduce the reduced basis set used in this work and consisting of unit cell restricted Bloch functions. This section also illustrates several tests and comparisons which validate the basis in terms of the reconstruction of the ab initio electronic structure. Sec. IV presents the transport model based on the NEGF formalism, and the procedure to achieve simulations accounting for a self-consistent description of the electrostatics via the Poisson equation. Then in Sec. V we illustrate some examples of complete, self-consistent device simulations for an Mos 2 based nanoscale transistor, and in Sec. VI we finally offer some concluding remarks.

II. AB INITIO HAMILTONIAN

Electronic-structure methods from first principles are typically based on density-functional theory, where one-electron states ("Kohn-Sham orbitals") are obtained by solving selfconsistently the Kohn-Sham equations,

H KS Ψ n = E n Ψ n , H KS = T + V scf (1)
where the Kohn-Sham Hamiltonian H KS is the sum of the kinetic energy T and of the self-consistent potential V scf . In turn, For each atom µ, the nonlocal term v µ N L can be expressed as a sum of N µ projectors, defined via atomic pseudopotential parameters β n (r) and D nn as follows:

V scf = V eI + V H + V xc ,
v µ N L (r, r ) = n,n β µ n (r) D µ nn [β µ n (r ) ] * (2)
and throughout this paper we use a * to denote the complex conjugate of a scalar a, and M † to denote the adjoint of a matrix or a vector M. Eq. (2) corresponds to the "separable" form of pseudopotentials. For the simple norm-conserving pseudopotentials used in this work, the D nn matrix is diagonal: D nn = D n δ nn . The index n is a combined index, running on angular momentum quantum numbers l and m, up to the highest angular momentum values present in the atomic core. Typically, just a few projectors (< 10) per atom need to be taken into account.

The β µ n (r) functions are short-ranged and vanish for r > r c , where r c is the radius beyond which pseudo-atomic and true atomic Kohn-Sham orbitals are the same. For most atoms,

r c ∼ 0.1 ÷ 0.3 nm.
The nonlocal term in the Kohn-Sham Hamiltonian of a crystal is thus

V N L (r, r ) = µ,R v µ N L (r -d µ -R, r -d µ -R) = µ,R nn β µ n (r -d µ -R) D µ nn [β µ n (r -d µ -R)] * (3)
where d µ is the position of atom µ in the unit cell, the R's are lattice vectors, and it is understood that

(V N L ψ)(r) = V N L (r, r )ψ(r )dr . (4)
Kohn-Sham orbitals have the Bloch form and can be expanded into plane waves

P k+G (r) = 1 √ Ω e i(k+G)•r , ( 5
)
where k is the Bloch vector, the G's are reciprocal lattice vectors, Ω is the volume of the crystal. A finite set is obtained by choosing plane waves up to a given kinetic energy value E w (the "cutoff"):

2 2m 0 (k + G) 2 ≤ E w
, where m 0 is the electron mass. The Kohn-Sham Hamiltonian can be expanded into plane waves as well:

k + G|H|k + G ≡ H k (G, G ) = 2 2m (k + G) 2 δ G,G + V L (G -G ) + V N L (k + G, k + G ) (6)
where V L (G -G ) is the Fourier transform of the local part of the total potential (local pseudopotential plus Hartree and exchange potential). The nonlocal contribution, coming from the pseudopotential, is:

V N L (k + G, k + G ) = 1 Ω V N L (r,
r )e -i(k+G)r e i(k+G )r dr dr .

In principle the solution of the secular equation for the Kohn-Sham Hamiltonian matrix

G H k (G, G ) B(G ) = E B(G) (8)
provides the electronic structure E n (k) and the Bloch functions Ψ nk (r), which are completely determined by the eigenvectors B nk (G). In practice H k (G, G ) is a very large matrix and thus it is not stored or directly diagonalized, but rather one resorts to iterative techniques and to on-the-fly computation of H k Ψ products exploiting fast Fourier transform techniques 34 .

The potential and the charge density contain plane waves up to a cutoff energy E ρ =4E w .

In our calculations we used an orthorombic unit cell (see examples in Sec. III B), where the real space unit vectors can be written as a 1 =(a x , 0, 0), a 2 =(0, a y , 0), a 3 =(0, 0, a z ), with x being the transport direction. Hence the unit vectors of the reciprocal lattice are b 1 =(2π/a x , 0, 0), b 2 =(0, 2π/a y , 0), b 3 =(0, 0, 2π/a z ), the reciprocal lattice vectors are

G=n x b 1 +n y b 2 +n z b 3 (with n x , n y , n z =0, ±1, ±2 • • • ).
The Brillouin zone can be taken as the parallelepiped defined by the conditions -π/a s ¡ k s ≤ π/a s (with s=x, y, z), that has a volume Ω RZ =(2π) 3 /Ω cell with Ω cell =a x a y a z being the volume of the unit cell. It is understood that for a 2D crystal in the (x,y) plane, for example, the unit cell includes a relatively large vacuum region in the z direction that makes the extension of the reduced zone along z practically negligible, thus resulting in a 2D electron gas. Ultra-thin films or nanowires consisting of an underlying 3D crystal can be similarly described as a 2D or 1D system by inserting vacuum regions in the unit cell.

III. REDUCED BASIS OF UNIT CELL RESTRICTED BLOCH FUNCTIONS

In our methodology for transport simulations the expansion volume for G vectors is given by the cube inscribed to the sphere used in ab initio calculations, namely the cube set by the condition

2 2m 0 G 2 s ≤ E ρ /3 with s = x, y, z.
The grid of G vectors naturally defines a corresponding grid of points in real space and, if we denote by N Gs the number of G s vectors (with G x , G y , G z lying respectively along the x, y and z axis), the spacing of the grid in real space is d s =a s /N ds , where N ds =N Gs is the number of grid points along s inside the unit cell.

In the remainder of the paper we will often refer to the Hamiltonian matrix in different basis sets. Matrices are denoted by using capital letters into square brackets and, whenever necessary, a subscript indicates the basis. For example [H] K , [H] xKyz and [H] Φ denote the Hamiltonian matrix respectively in the plane-wave basis used in Sec. II, in the hybrid xK yz basis described in Sec. III A, and in the reduced basis of Bloch functions discussed in Sec. III B. The subscript may be omitted to lighten the notation when there is no ambiguity about the basis set. When we refer to the elements of the matrices, instead, we drop the square brackets and the subscript because the symbols used for the elements identify the basis: for example we write H(xK yz , x K yz ) to denote the elements of the Hamiltonian matrix [H] xKyz in the xK yz basis. Finally, we use curly brackets to denote column or row vectors: for example {Ψ} will be used for the column vectors representing the wave-functions.

A. Hybrid xK yz basis

Let us now consider a system having N cx unit cells in the transport direction x and subject to periodic boundary conditions. The Hamiltonian matrices in the plane-wave basis are given by Eq. ( 6) and are identified by k=(k x ,k yz ) with k yz =(k y ,k z ). The corresponding Hamiltonian matrix in the hybrid basis, xK yz consisting of real-space grid points along x and plane waves in the (y,z) directions, can be obtained by using the unitary transformation described below. Before discussing the transformation, however, we here notice that the short-range nature of the non local pseudopotential defined in Eqs. (2, 3) implies that the non local terms practically vanish on distances much smaller than a x , as already mentioned in Sec. II. In particular, we numerically verified that the Hamiltonian in the xK yz basis can be accurately expressed by the block tridiagonal form

[H] xKyz =         H 0,0 H 0,1 0 0 • • • H † 0,1 H † 0,1 H 0,0 H 0,1 0 • • • 0 • • • • • • . . . . . . H 0,1 0 • • • 0 H † 0,1 H 0,0         (9)
where each block describes an a x long region consisting of N dx discretization points

x j =0, d x , 2d x • • • (a x -d x )
, so that blocks H 0,0 and H 0,1 have a rank N G =N dx N Gy N Gz (with

N dx =N Gx ).
The [H] xKyz in Eq. ( 9) corresponds to a given k yz and it has N 2 cx blocks. The matrix is sorted so that, for each discretization point x j , we have all the N Gy N Gz entries corresponding to the spectral components K yz =k yz +G yz (with

G yz =(G y ,G z )).
According to Eq. ( 9) the knowledge of [H] xKyz coincides with the knowledge of H 0,0 , H 0,1 . In Appendix A we discuss in more details the structure of [H] xKyz in Eq. ( 9). In particular we argue that the short range non local pseudopotentials make the off-diagonal blocks other than H 0,1 (i, j) negligible, and result in an H 0,1 block that is a lower triangular matrix, namely we have H 0,1 (i, j) 0 for j ≥ i.

For a system having N cx unit cells along x, the [H] xKyz in Eq. ( 9) can be obtained by transforming from K x to x the Hamiltonian matrices H k in Eq. ( 6) for the corresponding N cx wave-vectors k=(k x , k yz ) with k x values spaced by 2π/(N cx a x ). However we argue that N cx =2 is sufficient to determine H 0,0 , H 0,1 . In fact, for N cx =2 we can rewrite Eq. ( 9) as

[H] xKyz =   H 0,0 H 0,1 + H † 0,1 H † 0,1 + H 0,1 H 0,0   (10)
and Eq. ( 10) allows us to unambiguously determine H 0,0 and H 0,1 in virtue of the above mentioned property H 0,1 (i, j) 0 for j ≥ i. The [H] xKyz in Eq. ( 10) can be obtained by transforming the two matrices H k 0 , H k 1 given by Eq. ( 6) for

k 0 =(k x0 , k yz ) and k 1 =(k x1 , k yz ), with k x0 =0 and k x1 =π/a x .
For any k=(k x , k yz ) the transformation from K x to the home unit cell given by

x j = 0, d x , • • • (a x -d x ) is governed by the N G × N G unitary matrix 35 [U 0 kx ] = 1 √ N dx         I e i(kx+G x,1 )dx I • • • e i(kx+G x,1 )(ax-dx) I I e i(kx+G x,2 )dx I • • • e i(kx+G x,2 )(ax-dx) I • • • • • • . . . . . . I e i(kx+G x,N Gx )dx I • • • e i(kx+G x,N Gx ) (ax-dx) I         ( 11
)
where I is an identity matrix with rank N Gy N Gz =N G /N dx and N dx =N Gx . Hence, the corresponding transformation to the unit cell p extending from x j =pa x to

x j =[(p + 1)a x -d x ] (with p=1, 2, • • • ) is governed by the matrix [U p kx ]=[U 0
kx ] e ikx p ax . Consequently, the transformation from K x to the two unit cells necessary to calculate [H] xKyz in Eq. ( 10) can be written as

[U (2ax) kx ] = 1 √ 2 [U 0 kx ], [U 0 kx ]e ikx ax . (12)
By using Eq. ( 12) we can express the [H] xKyz in Eq. (10) in terms of the two plane-wave

DFT Hamiltonian matrices [H k 0 ], [H k 1 ] as [H] xKyz = k=k 0 ,k 1 [U (2ax) kx ] † [H k ] [U (2ax) kx ] = = 1 2 k=k 0 ,k 1   [U 0 kx ] † [H k ][U 0 kx ] [U 0 kx ] † [H k ][U 0 kx ] e ikx ax [U 0 kx ] † [H k ][U 0 kx ] e -ikx ax [U 0 kx ] † [H k ][U 0 kx ]   (13)
where for k 0 =(0, k yz ) and k 1 =(π/a x , k yz ) the exponential terms simply evaluate to ±1.

Hence, we can finally identify H 0,0 and H 0,1 in the hybrid xK yz basis as 36

[H 0,0 ] = 1 2 [U 0 k x0 ] † [H k 0 ][U 0 k x0 ] + [U 0 k x1 ] † [H k 1 ][U 0 k x1 ] (14a) [H 0,1 ] + [H 0,1 ] † = 1 2 [U 0 k x0 ] † [H k 0 ][U 0 k x0 ] -[U 0 k x1 ] † [H k 1 ][U 0 k x1 ] (14b)
In Appendix A we discuss a generalization of Eqs. ( 12), ( 13), ( 14) that can be used to build the Hamiltonian in the hybrid xK yz basis by using any number N cx of k x Bloch vectors.

B. Reduced basis of unit cell restricted Bloch functions

The blocks H 0,0 , H 0,1 in Eq. ( 9) have an N G rank so that, for high precision ab-initio calculations using a relatively large cutoff energy E ρ , the direct manipulation of such matrices for transport simulations is practically intractable. A drastic reduction of the size of H 0,0 , H 0,1 can be achieved by moving to an appropriate basis set consisting of Bloch functions restricted to a unit cell. To this purpose we first argue that, thanks to the block tridiagonal form of the Hamiltonian in Eq. ( 9), the corresponding Bloch functions in the unit cell p (with p=0, 1, • • • (N cx -1)) take the form 37

Ψ k (x j + p a x , G yz ) = Ψ 0 k (x j , G yz ) e i kxp ax (15)
where the {Ψ 0 k } are Bloch functions restricted to the home unit cell (i.e. for x j =0,

d x , • • • (a x -d x )
), which are in turn the solutions of N cx eigenvalue problems

H † 0,1 e -ikxax + H 0,0 + H 0,1 e ikxax {Ψ 0 k } = E(k) {Ψ 0 k } (16)
with k=(k x ,k yz ). Equation ( 16) requires that the Ψ 0 k (x j , G yz ) functions fulfill the k x dependent boundary condition Ψ 0 k (a x , G yz )=Ψ 0 k (0, G yz ) e i kx ax . The reduced basis employed in this paper is identified as a single basis set suitable for all the eigenvalue problems in Eq. (16). In this respect, we recall that the Bloch functions {Ψ 0 k } are in effect known, because they are determined by the eigenvectors B nk (G) of the secular Eq. ( 8) solved in ab-initio calculations. We can express the column vectors {Ψ 0 nk } in matrix notation as

{Ψ 0 nk } = [U 0 kx ] † {B nk } (17)
where [U 0 kx ] is defined in Eq. ( 11) and k=(k x ,k yz ). Our reduced basis set consists of a subset of the {Ψ 0 nk } corresponding to a few k x values in the reduced zone -π/a x < k x ≤ π/a x ] and, for each k x , to some tens of energies E n (k). There is substantial flexibility in such a definition of the basis set, which can be used to find a good compromise between the size N B of the basis and the accuracy in the reconstruction of the electronic structure. In general the size of the basis can be written as

N B = N kB i=1 N E (k x,i ) (18)
where N kB denotes the number of k x values and N E (k x,i ) the number of energies at k x,i included in the basis.

Appendix B offers more details about the choice of the basis functions. Most of the calculations and simulations reported in this paper were obtained by using either two k x values (i.e. k x =0, π/a x ), or four k x values (i.e. k x =0, ±0.5 π/a x , π/a x ), as exemplified by the results in Fig. 1 discussed below.

Because the {Ψ 0 nk } for different k x are not orthogonal over a x , we apply an orthonormalization procedure. Here, we remark that the use an orthonormal basis is not really necessary, but it is in fact very convenient in the transformations from the reduced basis to real space that are necessary, for example, for the calculation of space charge density in self-consistent simulations (see Sec. IV A). The orthonormalization procedure simply consists in an appropriate linear combination of the Bloch functions {Ψ 0 nk }, but it is prone to numerical instabilities. We overcame this problem by employing the modified Gram-Schmidt algorithm, which is more robust than the standard approach in dealing with rounding errors [START_REF] Watkins | Fundamentals of matrix computations[END_REF] .

Since the {Ψ 0 nk } are column vectors with N G components and we select only a number N B of basis functions much smaller than N G , the orthonormalization procedure can be always successfully completed.

We will denote by Φ m (x j , G yz ) the orthonormalized Bloch functions restricted to a unit cell, with m=1, 2,

• • • N B and x j =0, d x , 2d x • • • (a x -d x )
. The rectangular transformation matrix from the hybrid xK yz basis to the Φ basis is defined as

[U Φ ] = [ {Φ 1 }, {Φ 2 } • • • {Φ N B } ] (19)
and [U Φ ] has N B columns and N G rows. Each {Φ m } is a column vector such that for each

discretization point x j =0, d x , 2d x • • • (a x -d x )
we have all the N Gy N Gz spectral components.

Because the {Φ m } have been orthonormalized, we have

[U Φ ] † [U Φ ]=[I] N B with [I] N B being
the identity matrix with rank N B .

By recalling the expression in Eqs. (14) for the H 0,0 and H 0,1 in xK yz basis, we can readily write H 0,0 and H 0,1 in the reduced basis as

[H 0,0 ] Φ = 1 2 [W k x0 ] † [H k 0 ][W k x0 ] + [W k x1 ] † [H k 1 ][W k x1 ] (20a) [H 0,1 ] Φ + [H 0,1 ] † Φ = 1 2 [W k x0 ] † [H k 0 ][W k x0 ] -[W k x1 ] † [H k 1 ][W k x1 ] ( 20b
)
where we have introduced transformation matrices

[W k x0 ] = [U 0 k x0 ][U Φ ], [W k x1 ] = [U 0 k x1 ][U Φ ]. Eqs. (20) express [H 0,0 ] Φ , [H 0,1 ] Φ directly in terms of the two plane-wave DFT Hamilto- nian matrices [H k 0 ], [H k 1 ]. The expressions for [H 0,0 ] Φ , [H 0,1 ] Φ ,
in turn, allow us to transform the Hamiltonian in Eq. ( 9) and the eigenvalue problems in Eq. ( 16) to the reduced Φ basis, where the size of such blocks is much smaller than N G , as exemplified below. For a given N kB the accuracy of the reduced basis improves by increasing the number N E of Bloch functions at each k x . As it can be seen the error for N B = 160 is small enough for most applications in both Fig. 2(a) and Fig. 2(b), which enables a drastic reduction of the problem size compared to N G =26733. For N B =160 and 240 we also see that, for a given N B , the reconstruction of the electronic structure improves by increasing N kB . This behavior is not unexpected because the {Ψ 0 nk } are continuous functions of k x , so that a linear combination of the {Ψ 0 nk } for a few k x can still approximate well the remaining Bloch functions. Fig. 2 shows, however, that very few k x values are sufficient to achieve a close agreement with the reference results.

We finally notice in Fig. 2(b) that, for N B =120, the energy error in the conduction band steeply increases for energies above about 2 eV. This is because in this system we have 26 valence bands, all of which are included in the reduced basis set. Consequently, for N E =30 the basis set includes the Bloch functions for only the four lowest conduction bands, which results in relatively large errors for higher conduction bands. In this respect we verified that, by including in the basis set the Bloch functions for all the valence bands at each k x , we can effectively suppress the problem of unphysical solutions 25,26 . Hence the number of valence bands in the system sets a lower bound for the size N B of the reduced basis.

As for the unphysical states that have been sometimes observed upon the introduction of a reduced basis, we have extended our analysis by inspecting the transmission across a single layer MoS 2 in the flat band condition, namely with neither built in nor externally applied potential. The transmission calculated by using the reduced basis (not shown) is very steeply, exponentially suppressed in the energy gap. In other words we do not observe any feature of the transmission indicating the presence of evanescent, spurious states in the energy gap. Moreover, for energies belonging to the valence or the conduction band the transmission at a given lateral wave-vector k y equals, as expected, the number of available bands. The inspection of the transmission is a good complement to the analysis of the electronic states, and it reinforced our confidence in an effective suppression of unphysical effects in our reduced basis calculations.

IV. NEGF BASED TRANSPORT MODEL

In this section we describe the procedure based on the NEGF formalism to achieve a self-consistent simulation of a nanodevice or a mesoscopic system subject to external bias conditions. All the relevant physical quantities, such as density of states, carrier concentration and currents, were computed in terms of retarded and Green's function matrices, which are calculated in the reduced Bloch function basis.

A. Charge, current and self-consistent calculations

In the Bloch function basis the retarded(advanced), [G r(a) ] Φ , and lesser(greater)-than Green's functions, [G <(>) ] Φ at a given energy E are defined as

[G r(a) ] Φ = (E + iη)[I] Φ -[H] Φ -[Σ r(a) ] Φ -1 (21) and [G <(>) ] Φ = [G r ] Φ [Σ <(>) ] Φ [G a ] Φ (22)
where η is a positive(negative

) infinitesimal, [Σ r(a) ] Φ = [Σ r(a) L ] Φ + [Σ r(a) R ] Φ + [Σ r(a) ph ] Φ and [Σ <(>) ] Φ = [Σ <(>) L ] Φ +[Σ <(>) R ] Φ +[Σ <(>)
ph ] Φ are the retarded(advanced) and the lesser(greater)than self-energies describing the connection to contacts (i.e. left lead, L, and right lead, R), or possible interaction with photons or phonons [START_REF] Mahan | Many-Particle Physics[END_REF] . Thanks to the block-tridiagonal structure of the Hamiltonian matrix, the sub-matrices of the retarded(advanced) and lesser(greater)than Green's functions that are needed to calculate carrier concentrations and current density can be efficiently computed with well-known recursive algorithms [START_REF] Anantram | Proceedings of the IEEE[END_REF] , and by manipulating matrix blocks of rank N B .

More precisely, in order to calculate the 3D real space concentration of mobile carriers, we need to compute the diagonal terms of the real space Green's functions starting from the Green's functions in the Bloch functions basis. To this end we first transform the Green's functions from the Bloch functions to the xK yz basis, and then compute the charge in real space. Inside each unit cell the Green's functions in the xK yz basis can be concisely written

in matrix notation as [G <(>) ] xKyz =[U Φ ] [G <(>) ] Φ [U Φ ] † , with [U Φ ]
given by Eq. ( 19). An explicit expression is given by

G <(>) (x j K yz , x j K yz ; E) = N B n,m=1 G <(>) (n, m; E) Φ n (x j , G yz ) Φ * m (x j , G yz ) ( 23
)
Then we evaluated the free electron concentration on the fine mesh grid with discretization

steps (d x , d y , d z ) as n(x j , r yz ) = -i d x d y d z ∞ E 0 (x j ) dE 2πN Gy N Gz Gyz,G yz G < (x j K yz , x j K yz ; E) e i(Gyz-G yz )•ryz (24)
and similarly for the free hole concentration

p(x j , r yz ) = i d x d y d z E 0 (x j ) -∞ dE 2πN Gy N Gz Gyz,G yz G > (x j K yz , x j K yz ; E) e i(Gyz-G yz )•ryz (25)
with K yz =(k yz +G yz ) and with E 0 (x j ) being the neutrality point that we assumed to be at the center of the energy bandgap.

It is understood that all equations in this section refer to a given k yz and that, if the system is periodic along either y or z, a sum over k yz is necessary to calculate all physical quantities.

In order to simulate the transport properties of realistic devices, it is necessary to evaluate the electrostatic potential induced by external biases, ionized dopants and mobile carriers.

Such an electrostatic potential φ(r) can be accurately described within the Hartree approximation, namely by self-consistently solving the equations for the Green's functions (that in turn give the carrier concentrations via Eqs. (24,25)), with the 3D Poisson equation

∇ • [ (r)∇φ(r)] = -e [p(r) -n(r) + N D (r) -N A (r)] (26)
where (r) is the material-dependent permittivity, and n(r), p(r), N A (r), N D (r) are the electron, hole, acceptor and donor concentration, respectively.

Here, we assume that the electrostatic potential φ(r) varies over a relatively large spatial scale compared to the discretizations (d x , d y , d z ) used to calculate the atomistic Green's functions. Consequently, in order to reduce the size of the Poisson equation problem in devices with a technologically relevant size, we transfer the free electron and hole concentrations computed in Eqs. (24,25) on a coarser mesh with a discretization (∆ x , ∆ y , ∆ z ) ∼ 0.1 ÷ 0.2 nm.

The conversion from the finer to the coarser mesh is performed so as to conserve the integral of the carrier concentration. The effect of external biases was imposed by setting Dirichlet boundary conditions at the contacts and Neumann boundary conditions to non contacted boundary regions.

Finally, once the self-consistent solution has been obtained, we express the electron current as a function of the Green's functions and self-energies computed at the contact L(R)

as

I L(R) = -i e dE tr [Γ L(R) ] [G r L(R) ] -[G r L(R)) ] † f L(R) + [G < L(R) ] (27)
where

[Γ L(R) ] = i [Σ L(R) ] -[Σ L(R) ] † , tr{• • • } is for the trace operation and f L(R) is the
Fermi-Dirac distribution. The calculation of the current can be carried out by using the Green's functions in the reduced Bloch functions basis, where the size of the matrices is the smallest.

B. Implementation and computational burden

Our simulation procedure takes full advantage of the first principle calculations carried out by the ab initio solver 10 . In fact, after a duly converged ab initio simulation has been achieved for the unit cell of the physical system, our approach can be summarized in the following steps.

1. Assuming N kB = 2, for example, select a subset of the Bloch states for k 0 =(0,k yz ) and k 1 =(π/a x ,k yz ) obtained by DFT calculations and transform them to the xK yz basis using Eq. (17). Then orthogonalize the basis functions so as to obtain {Φ k 0 }, {Φ k 1 with n=1, 2, • • • N B and assemble the U Φ in Eq. (19).

2. Calculate [H 0,0 ] Φ , [H 0,1 ] Φ in the reduced Bloch function basis by using Eq. (20).

3. Solve for the Green's functions in the reduced basis by using an initial guess of the electrostatic potential φ(r).

4. Calculate carrier concentrations from Eqs. (24, 25) and then solve the Poisson equation for a new guess of φ(r). Loop between steps 3 and 4 until a specified convergence is reached and finally calculate the current using Eq. ( 27).

Here it should be mentioned that, for any physical system or nanoscale device, steps 1 and 2 above have to be performed only once. In fact they correspond to material properties, so that the relevant quantities calculated in these steps can be stored and re-used in subsequent simulations corresponding to different bias conditions. The number N B of Bloch basis functions is the most important parameter affecting the computational load of the Green's functions equations.

Several important optimizations are possible and have been introduced in the implementation of our methodology. For example, it is apparent from Eqs. (24, 25) that the sums over G yz , G yz for any couple of basis functions Φ n , Φ m can be carried out only once for a given physical system. In other words we can introduce the new quantity K n,m (x j , r yz ) defined as

K n,m (x j , r yz ) = 1 N Gy N Gz Gyz,G yz Φ n (x j , G yz )Φ * m (x j , G yz ) e i(Gyz-G yz )•ryz (28)
and then notice that the free electron concentration can be written in terms of K n,m (x j , r yz )

as

n(x j , r yz ) = -i d x d y d z ∞ E 0 (x j ) N B n,m=1 G < (n, m; E) K n,m (x j , r yz ) dE 2π . (29)
A similar expression holds for the free hole concentration p(x j , r yz ).

As it can be seen, once K n,m (x j , r yz ) has been calculated for a given physical system, then the carrier concentrations during self-consistent simulations can be obtained by using Eq. ( 29), that is by skipping the sums over G yz , G yz .

V. SIMULATION RESULTS

We present in this section an example of a nanoscale transistor that we could efficiently simulate by using the reduced basis of Bloch functions computed directly from plane-wave DFT calculations.

The device under investigation is composed by the monolayer MoS 2 n-MOSFET sketched in Fig. 3(a), where the 2D semiconductor is sitting on a 10 nm thick back-oxide having a dielectric constant ε ox = 3.9ε 0 (with ε 0 being the vacuum permittivity). Source and drain regions are considered chemically doped with a donor concentration of 10 20 cm -3 and the top gate oxide has an equivalent oxide thickness of about 1 nm. The lateral direction was assumed to be periodic and was described by including a discrete sampling of the wave-vector k y with a constant ∆k y = 0.1 π/a y .

Figure 3 shows the drain current, I DS , versus the top gate voltage, V T G , characteristic for a gate length of L G = 64a x 20 nm, 32 a x 10 nm and 16 a x 5 nm. The large I DS values are due to the fact that neither scattering nor series resistance are included in simulations. Thanks to the sub-nanometer thickness of the MoS 2 layer the I DS versus V T G characteristics of the transistor are still well behaved for a channel length of about 5 nm, even if a degradation of the sub-threshold swing is observed with respect to longer FETs.

The onset of short channel effects also manifests itself in a significant left-shift of the I DS versus V T G characteristic in the shortest gate lengths.

The sub-threshold swing degradation in the shortest device is mainly due to the onset of a sizeable source-to-drain quantum tunneling, as it is illustrated in Fig. 4 reporting the current spectra and the profile of the lowest conduction band along the transport direction for the MOSFET with either with L G = 5 nm or 10 nm. For the shortest gate length the spectral current is spread over energies well below the top of the barrier, thus confirming a significant source-to-drain tunneling contribution to the off current in this specific bias condition. In the longer transistor, instead, the off current is dominated by thermionic emission above the top of the barrier.

VI. CONCLUSIONS

We presented new theoretical developments and a sound implementation for a firstprinciple transport model based on the NEGF formalism, and on a basis set obtained directly from the ab initio Bloch functions. Differently from previous papers proposing a similar approach for DFT calculations based on an LCAO basis, we used plane-wave ab initio calculations and we argue that, thanks to an appropriate choice of the basis functions, we could effectively suppress the problem of unphysical solutions, whose treatment is delicate and computationally demanding 25,26 .

We found that the unit cell restricted Bloch functions basis enables band structure and transport calculations with hundreds times reductions of the size of the problem compared to the original DFT formulation. Moreover, while we have here reported results only for a homogeneous system consisting of a single layer MoS 2 , we envision that our approach can also be used for hetero-structures, thus paving the way for a number of technologically and physically important applications, such as contacts between metals and 2D materials, as well as vertical or horizontal hetero-junctions between 2D semiconductors.

The methods of this work can be applied also to first-principle calculations based on an LCAO basis, however we think that the herein reported demonstration for plane-wave DFT

Hamiltonians is particularly promising for future developments concerning the electronphonon interaction. In fact the NEGF formalism can naturally include electron-phonon scattering [START_REF] Mahan | Many-Particle Physics[END_REF] and, moreover, the plane-wave DFT approach is especially suitable for the calculation of phonon spectra [START_REF] Baroni | [END_REF] and electron-phonon coupling coefficients 45 .

The benefits of a Bloch functions basis are not confined to first-principle methods, on the contrary they directly apply also to empirical pseudopotential Hamiltonians in both their local and non-local formulation 46,[START_REF] Cohen | Electron structure and optical properties of semiconductors[END_REF] , and promise large computational advantages compared to the methods recently proposed by some of the present authors 23,24 .

We believe that the results of this work qualify the methodology based on unit cell restricted Bloch functions as a viable approach for ab initio and semi-empirical quantum transport simulations and, in particular, as an alternative to maximally localized Wannier functions. In this respect, while the direct use of Bloch functions is attractive in that it circumvents the a posteriori determination of the Wannier functions basis, further work is admittedly needed to demonstrate the general applicability of the methods of this paper, and the feasibility of the above mentioned extensions to hetero-structures and dissipative transport.

where [U 0 kx ] has been defined in Eq. ( 11). Then we can use Eq. (A1) and reformulate Eq. ( 13) as

[H] xKyz = k [U (Ncx) kx ] † [H k ] [U (Ncx) kx ] (A2)
where

[H 0,0 ] = 1 N cx k [U 0 kx ] † [H k ][U 0 kx ] (A3a) [H 0,1 ] = 1 N cx k [U 0 kx ] † [H k ][U 0 kx ] e ikxax (A3b)
Eq. ( A3) is a generalization of Eq. ( 14) for N cx ≥ 3, that immediately leads to the corresponding extension of Eq. ( 20) for the Hamiltonian blocks in the reduced basis. In the remainder of this Appendix we discuss in more details the assumption in Eq. ( 9)

about the structure of the Hamiltonian in the xK yz basis.

As it is discussed in Secs. II and III A, the non local pseudopotential V N L (r, r ) defined in Eq. ( 3) is the dominant non local term of the Hamiltonian matrix in the xK yz basis.

However [H] xKyz has to be calculated numerically by using Eq. ( 13) for N cx =2 or Eq. (A2) for larger N cx values. Consequently, while it is expected that the short-range nature of the non local pseudopotential results in a single off-diagonal block of the Hamiltonian matrix as assumed by Eq. ( 9), we have no analytical expression for V N L (x j , G yz ) ensuring a priori that, for example, a second off-diagonal block H 0,2 is indeed negligible. Hence we carried out a numerical analysis of this aspect.

In this respect we firstly notice that, if we hypothesize that H 0,2 is not negligible then, in order to identify H 0,0 , H 0,1 and H 0,2 , it is necessary to write the [H] xKyz across at least four unit cells and thus by using four k x values (i.e. k x =0, ±0.5 π/a x , π/a x ). In fact, for N cx =4

we can extend Eq. ( 10) and write

[H] xKyz =        H 0,0 H 0,1 H 0,2 + H † 0,2 H † 0,1 H † 0,1 H 0,0 H 0,1 H 0,2 + H † 0,2 H † 0,2 + H 0,2 H † 0,1 H 0,0 H 0,1 H 0,1 H † 0,2 + H 0,2 H † 0,1 H 0,0        (A4)
which allows us to identify H 0,2 (besides H 0,0 and H 0,1 ), by assuming that the H 0,2 block is a lower triangular matrix (namely H 0,2 (i, j) 0 for j ≥ i). In the presence of H 0,2 , the secular Eq. ( 16) can be rewritten as

H † 0,2 e -ikx(2ax) + H † 0,1 e -ikxax + H 0,0 + H 0,1 e ikxax + H 0,2 e ikx(2ax) {Ψ 0 nk } = E n (k){Ψ 0 nk } (A5)
and Eq. (A5) allows one to determine the electronic structure in the xK yz basis and duly accounting for H 0,2 .

For the single layer MoS 2 studied in this work, we calculated [H] xKyz for N cx =4 by using Eq. (A2) and we verified that the electronic structure obtained with Eq. (A5) and accounting for H 0,2 is practically identical to the results of Eq. ( 16) that instead neglects H 0,2 . This means that H 0,2 can be neglected in the case at study, which legitimates the tridiagonal form of [H] xKyz in Eq. ( 9) and implies that H 0,1 block is a lower triangular matrix. This latter observation, in turn, allows us to determine H 0,0 , H 0,1 from Eq. (10).

We reiterate that the negligibility of H 0,2 is an expected result from a physical perspective, in virtue of the relatively short range nature of the non-local pseudopotential discussed in Sec. II.

Quite interestingly the Ψ 0 nk (x j , G yz ) for two different k x values appear to be sufficient to obtain solutions of Eq. ( 16) fulfilling the boundary condition Ψ 0 (a x , G yz )=Ψ 0 (0, G yz ) e i kx ax for any k x value. In order to show this we first recall the standard notation for Bloch functions Ψ 0 nk (x j , G yz )=u n,k (x j , G yz ) e i kxx j (where u n,k (x j , G yz ) is the periodic part of Ψ 0 nk ), and then we expand the unknown Ψ 0 for a generic k x value in terms of the Ψ 0 nk 0 and Ψ 0 nk 1 for k 0 =(0, k yz ) and k 1 =(π/a x , k yz ). By recalling Eq. ( 17) we write

Ψ 0 (x j , G yz ) = M 0 n=1 C n,k 0 u n,k 0 (x j , G yz ) + M 1 m=1 C m,k 1 u m,k 1 (x j , G yz ) e i π ax x j (B1)
where x j =0, d x , 2d x • • • (a x -d x ) and M 0 , M 1 denote the number of respectively Ψ 0 nk 0 and Ψ 0 nk 1 functions included in the expansion. Because u n,k 0 and u n,k 1 are periodic over a x , it is straightforward to see that the boundary condition Ψ 0 (a x , G yz )=Ψ 0 (0, G yz ) e i kx ax becomes 1 -e i kx ax M 0 n=1 C n,k 0 u n,k 0 (0, G yz ) = 1 + e i kx ax M 1 m=1 C m,k 1 u m,k 1 (0, G yz ) .

(B2)

For k x =0 the Eq. (B2) can be fulfilled by taking all C m,k 1 = 0 so that, as expected, Ψ 0 can be obtained by using only the Ψ 0 nk 0 functions. Likewise for k x =π/a x we can take all C m,k 0 = 0 and build the solution by using only the Ψ 0 k 1 functions. For any other k x value, Eq. (B2) provides the relation between the C n,k 0 and C m,k 1 coefficients ensuring that Ψ 0 fulfills the boundary condition Ψ 0 (a x , G yz )=Ψ 0 (0, G yz ) e i kx ax , hence it can be a solution of Eq. ( 16).

While two k x values appear to be sufficient to build a basis, we found that sampling the reduced zone with more than two k x values can improve the accuracy in the reconstruction of the electronic structure for a given overall number of basis functions, as it is exemplified in Fig. 2. Moreover we found that, by including in the basis set the Bloch functions for all the valence bands at each k x , we can effectively suppress the unphysical solutions sometimes observed upon the introduction of a reduced basis 25,26 . 16) with [H 0,0 ], [H 0,1 ] blocks), or in the reduced basis (i.e. Eq. ( 16)

  where V eI is the electron-ion interaction potential, V H the Hartree electrostatic potential, V XC the "exchange-correlation" potential. The two latter terms depend upon the charge density, which can be written as the sum of the squares of all occupied Kohn-Sham orbitals. Let us use a plane-wave bases set and pseudopotentials to represent the valence electronnuclei interactions. The solution of the Kohn-Sham equations reduces to a secular problem, in which the potential is computed self-consistently. Leaving apart the problem of how to compute the charge density and the self-consistent potential, the only difference between Kohn-Sham and empirical-pseudopotential Hamiltonians is the presence of a non-local term in the Kohn-Sham Hamiltonian. Both the Hartree and the exchange-correlation potentials, the latter in the typical GGA (generalized gradient approximation) form, are local functions v(r) of the position r. The nonlocal term stems from atomic norm-conserving pseudopotentials, which contain two types of contributions: (a) a local v L (r) part with the expected asymptotic v L (r) ∼ -Z v e 2 /r behavior at large r (e is the electron charge, Z v the number of valence electrons of the atom); (b) a nonlocal, short-ranged v N L (r, r ) part.

Figure 1

 1 Figure 1(a) illustrates the unit cell of the single layer MoS 2 that we used in this work as a baseline material for electronic-structure calculations and device simulations; the unit vectors are a 1 =(a x , 0, 0), a 2 =(0, a y , 0), a 3 =(0, 0, a z ) with a x = 3.18818 Å, a y = 5.52208 Å and a z = 20.2 Å, while the vertical distance between the Mo and S atoms is 1.564 Å. This cell was built by expanding the relaxed primitive unit cell of the monolayer MoS 2 , whose bandstructure along the high-symmetry points of the primitive Brillouin zone matches very well the results reported in Ref. 39 (not shown). Figure 1(b) reports the corresponding band structure obtained from the DFT Hamiltonian in the plane-wave basis (solid line). The DFT calculation was performed by means of the Quantum ESPRESSO code 10 , using a normconserving pseudopotential 40 , and the Perdew-Burke-Ernzerhof 41 (PBE) approximation to the exchange-correlation functional. The self-consistent solution was obtained by employing a 15 × 12 × 1 Monkhorst-Pack k-points grid and a cutoff energy of E w = 90 Ryd, resulting in a number on plane waves N G =26733. Fig. 1(b) also shows the band structure obtained with Eq. (16) in the hybrid xK yz basis (symbols), still having a size N G =26733.

Figure 2

 2 Figure2addresses the reconstruction of the electronic structure in the reduced basis Φ, and in particular it reports the absolute energy difference ∆E (x-axis) between the energies obtained by using Eq. (16) either in the reduced basis (i.e. for [H 0,0 ] Φ , [H 0,1 ] Φ with an N B size), or in the complete xK yz basis (i.e. for [H 0,0 ], [H 0,1 ]) with an N G size). The results of Fig.2(a) were obtained by using two k x values to build the reduced basis (i.e. N kB =2 and k x =0, π/a x ), whilst those in Fig.2(b) correspond to N kB =4 (i.e. k x =0, ±0.5π/a x , π/a x ).

Figure 5

 5 Figure 5 illustrates, for the device in Fig. 3(a) subject to a specific external bias, the dependence of the self-consistently calculated I DS on the number N B of Bloch states in the reduced basis. The corresponding CPU time versus N B is also reported. It can be observed that, for the case at study, the self-consistent solution of the Poisson-NEGF equations rapidly converges to a stable value for N B ≥ 90, which in turn enables the simulation of one bias point with a wall clock computation time slightly longer than one hour by using only ten CPU cores.

  By using Eq. (A3) we verified that the reconstruction of [H 0,0 ], [H 0,1 ] from the plane-wave DFT Hamiltonian matrices [H k ] is practically independent of N cx for N cx > 2. Namely any further increase of N cx has a negligible effect on the electronic structure calculated in the xK yz or in the Φ basis that has been analyzed, for example, in Figs. 1, 2.

FIG. 1

 1 FIG. 1. (Left) Cross section and the top view of the unit cell for the single layer MoS 2 employed in this work. (Right) Electronic structure for the single layer MoS 2 versus k x and for k y =0. Ab initio calculations (solid lines) are compared to results obtained by using Eq. (16) in the hybrid xK yz basis (diamonds).

FIG. 2 .

 2 FIG.2. Absolute energy difference, ∆E, between the electronic structure calculated either in the hybrid xK yz basis (i.e. Eq. (16) with [H 0,0 ], [H 0,1 ] blocks), or in the reduced basis (i.e. Eq.(16)

  with [H 0,0 ] Φ , [H 0,1 ] Φ blocks). (a) ∆E for reduced basis calculations obtained with two k x values (i.e. k x =0, π/a x ) and different number N E of basis functions at each k x . (b) Same as (a) but for reduced basis calculations obtained with four k x values (i.e. k x =0, ±0.5π/a x , π/a x ).

FIG. 3 .FIG. 4 .

 34 FIG. 3. (Left) Sketch of the single-gate MOSFET simulated in this work and consisting of a single layer MoS 2 channel material. The length of the source and drain extensions is L S = L D 9 nm and L G indicates the top gate length. (Right) Simulated drain current, I DS , versus top gate voltage, V T G , characteristics at V DS =0.6 V for n-type MoS 2 MOSFETs featuring different gate lengths L G 20, 10 and 5 nm. The I DS curves reported in semi-logarithmic scales have been V T G shifted so as to have the same I DS =1µA/µm at V T G = 0 V for all gate lengths, whereas the curves in linear scales have not.

  the sum runs over N cx Bloch vectors k=(k x ,k yz ), with k x taking N cx values in the range -π/a x < k x ≤ π/a x with a spacing 2π/(N cx a x ) and including k x =0. Eq. (A2) allows us to calculate [H] xKyz for any value of N cx and the matrix [H] xKyz has a rank N cx N G .Eq. (A2) is useful in several respects. Firstly we see that, according to the assumption for [H] xKyz in Eq. (9), the blocks [H 0,0 ] and [H 0,1 ] can be determined from the element (1,1) and (1,2) of the matrix [H] xKyz , so that Eq. (A2) allows us to write

R. Kim, U. E. Avci, and I. A. Young, IEEE Trans. on Electron Devices 62, 713 (2015).

D. Esseni, M. Pala, and T. Rollo, IEEE Trans. on Electron Devices 62,

(2015).

ACKNOWLEDGMENTS P. G. acknowledges support from the EU via the MaX Centre of Excellence (Project No. 824143).

In Sec. III the Hamiltonian in the hybrid basis was built by using two DFT Hamiltonian matrices [H k 0 ], [H k 1 ], but the extension of the methodology to more than two k values is quite natural. To this purpose we first generalize Eq. ( 12) and define the transformation

Appendix B: Selection of Bloch functions for the reduced basis

As already mentioned in Sec. III B, there exists a significant flexibility in the definition of the Bloch functions basis set, so that we add a few conceptual and practical remarks about the choice of the basis. The first is that the basis set cannot be formed by using the Ψ 0 nk (x j , G yz ) corresponding to a single k x value (with k=(k x ,k yz )), because all such functions fulfill the same boundary condition Ψ 0 nk (a x , G yz )=Ψ 0 nk (0, G yz ) e i kx ax for their specific k x value, so that it is impossible to build solutions of Eq. ( 16) for a different k x value.