Physical insights on transistors based on lateral heterostructures of monolayer and multilayer PtSe2 via Ab initio modelling of interfaces

Lateral heterostructures (LH) of monolayer-multilayer regions of the same noble transition metal dichalcogenide, such as platinum diselenide (PtSe2), are promising options for the fabrication of efficient two-dimensional field-effect transistors (FETs), by exploiting the dependence of the energy gap on the number of layers and the intrinsically high quality of the heterojunctions. Key for future progress in this direction is understanding the effects of the physics of the lateral interfaces on far-from-equilibrium transport properties. In this work, a multi-scale approach to device simulation, capable to include ab-initio modelling of the interfaces in a computationally efficient way, is presented. As an application, p- and n-type monolayer-multilayer PtSe2 LH-FETs are investigated, considering design parameters such as channel length, number of layers and junction quality. The simulations suggest that such transistors can provide high performance in terms of subthreshold characteristics and switching behavior, and that a single channel device is not capable, even in the ballistic defectless limit, to satisfy the requirements of the semiconductor roadmap for the next decade, and that stacked channel devices would be required. It is shown how ab-initio modelling of interfaces provides a reliable physical description of charge displacements in their proximity, which can be crucial to correctly predict device transport properties, especially in presence of strong dipoles, mixed stoichiometries or imperfections.

A promising route for the next-generation of efficient field-effect transistors (FETs) is to exploit lateral heterostructures (LHs) of two-dimensional (2D) materials, such as transition metal dichalcogenides (TMDs) [1][2][3] . Platinum diselenide (PtSe 2 ) and other noble TMDs are particularly appealing, because their energy gap depends on the number of layers, varying in a broad range from few eVs (semiconductor) to 0 eV (semimetal) [4][5][6][7][8][9][10] . Such tunability enables the fabrication of LH-FETs with high-quality, lattice-matched metal-semiconductor heterojunctions, which is crucial to obtain low-resistance contacts, required for both high performance and low power operation 2,11,12 . Inspired by well-established CMOS processing 13 , it has been recently proposed that by making the TMD channel thin under the gate but thick under the source and drain contacts-e.g. creating a recessed TMD channel 14,15 -it is possible to achieve efficient gate control while keeping contact resistance down to 350-400 Ω µm 14 . Similar improvements have been reported or predicted for LH-FETs based on PtSe 2 7,14,15 , but also other TMDs 8,16 and black phosphorus [17][18][19][20] . Engineering such nanoscale LH-FETs relies on the understanding and control of atomic-scale details of the lateral interfaces, which dominate the nearby charge distribution, thus altering the Schottky barrier shape and affecting the device performance 21 .
From the theoretical perspective, the problem of reliably simulating realistically-sized LH-FETs in far-fromequilibrium conditions is an open issue, especially when it comes to modelling charge displacements or nonidealities at interfaces 21,22 . In most cases, it requires tackling two major challenges. One is to identify a multi-scale procedure which, on the basis of ab-initio atom-by-atom simulations, can yield a small-sized Hamiltonian able to capture the electronic behaviour of the system within an energy window of interest, so that the most relevant transport features can be reproduced with acceptable computational cost. The second one is to include local Scientific Reports | (2021) 11:18482 | https://doi.org/10.1038/s41598-021-98080-y www.nature.com/scientificreports/ dipoles and charge displacements around interfaces from the atomistic ab-initio structure into the mesh used within a device simulator to solve device electrostatics far from equilibrium. Despite its importance, the latter aspect is often overlooked and neglected in device simulations, which rely on determining free-charge redistribution throughout the device according to the solution of the Poisson equation on a grid that is often too coarse to enable the capture of charge redistribution at the interface. In fact, nowadays the ability of including such fine level of detail in a far-from-equilibrium device simulation inevitably comes with the price of a high computational cost or a relatively small device size. State-of-the-art modelling approaches are based on projecting an Hamiltonian obtained with Density Functional Theory (DFT) onto a basis set of Maximally Localized Wannier Functions (MLWF) 23,24 or Pseudo-Atomic Orbitals (PAO) 25 . The former, in particular, allows one to reproduce bands within a relevant energy window with very good accuracy. Both vertical and lateral 2D heterostructures have been studied exploiting such method 12,24,26,27 . For instance, LHs of semimetallic 1 T and semiconducting 2H phases of MoS 2 26 , as well as LHs of multilayer-monolayer noble metal disulphides 12 were investigated by using DFT to estimate the Schottky barrier height, and then using it to rigidly shift the Hamiltonians of the two junction components represented on a MLWF basis. Interactions at the interface were defined therein using a combination of the internal couplings of the bulk components. As a result, the total heterostructure Hamiltonian had a small dimension, as few Wannier Functions (≤ 16) were enough to reproduce the bands of the bulk components over a wide energy window. However, despite its computational advantage, so far this method has not been able to properly capture neither the charge displacements at the interfaces, occurring on a fine spatial scale, nor their finite width.
A direct "wannierization" (i.e. projection onto a MLWFs basis) of the whole heterostructure system has been recently adopted by Szabó et al. 24,27 who used it to study both weakly and strongly coupled 2D heterojunctions, such as TMD-TMD and metal-TMD cases, respectively, directly extracting the Hamiltonian of all the heterostructure sub-regions from the heterostructure MLWF Hamiltonian. Compared to the approach described in Ref. 26 , in this case the interface regions can be reliably represented in real space, with the same basis as that of the individual bulk components. However, the dimension of the device Hamiltonian, and hence the computational burden of the calculations, is much larger, because the number of Wannier Functions needed to model an heterostructure is higher.
In this work we present a very efficient multi-scale procedure that enables the simulation of far-from equilibrium transport in realistically sized LH-FETs, accounting for the effects of local charge displacements at the lateral interfaces at low computational cost. We demonstrate the usefulness and the versatility of the procedure by investigating the performance potential of mono-multilayer PtSe 2 LH-FETs with different channel length, lead thickness and junction smoothness, considering the requirements set out in the 2020 edition of the International Roadmap for Devices and Systems (IRDS) 28 . In particular, we emphasize how ab-initio modelling of interfaces and a reliable physical description of charge distribution in their proximity are crucial to correctly predict charge transport in LH-FETs.
In order to successfully include such fine charge displacements in the coarse electrostatic grid, we propose to directly include their effect on the potential as a fixed contribution. This is equivalent to neglect variations of charge distribution at the interface occurring on a length scale much smaller than the minimum spacing of the electrostatic grid, as it is typical of mean field treatments. Therefore, we extract the on-site energy along the LH from ab-initio DFT simulations, and then map it directly onto the modelled device.
We anticipate that the scope of this work is not to perfectly model an interface, but rather to provide a modular Hamiltonian enabling one to explore several devices, materials or geometries capturing the relevant physics at a reasonably low computational cost. For a deeper chemical understanding of how the junctions atomic-scale details, such as defects, functionalization, or roughness, affect the heterostructure electronic properties, a more accurate orbital-resolved and entirely DFT-based transport code (e.g., Transiesta 29 ) would better be used on a preliminary level, even though on relatively smaller devices and in near-to-equilibrium conditions.
The article is organized in two parts. We first introduce, step by step, our multi-scale procedure in the context of building a model for bilayer-monolayer LH-FET of PtSe 2 . We then apply this method to predict the performance of both bilayer-monolayer and four-layer-monolayer PtSe 2 LH-FETs, with sharp and gradual junctions, in far-from-equilibrium conditions.

Results and discussion
Overview of the multi-scale model. We propose in this work a multi-scale modelling approach consisting in three steps, schematically illustrated in Fig. 1. First, we perform plane-wave DFT simulations of the bulk components of the LH, from which we extract a real-space small-sized Hamiltonian in the MLWF basis. Second, we use DFT to model the LH and calculate the on-site energy profile E on-site (x) along the transport direction x. Third, we construct a larger-scale device Hamiltonian in the MLWF basis by coupling the component Hamiltonians together, taking care of mapping site-by-site the ab-initio variation of E on-site in proximity of the junctions. As will be extensively discussed below, this will ensure that the model reproduces the effects of relevant local charge variations. Finally, the MLWF device Hamiltonian is used within the open-source device simulator NanoTCAD ViDES 30 to simulate ballistic transport in equilibrium and far-from-equilibrium conditions. In the next sections, we go through each of the steps towards the construction of a model for bilayer (2L)-monolayer (1L) PtSe 2 LH-FET.
Ab-initio models of bulk 2L and 1L PtSe 2 . All our DFT calculations are performed using the Quantum Espresso electronic structure package 31 sampling the first Brillouin zone with a 8 × 12 × 1 Monkhorst-Pack grid (increased to 16 × 24 × 1 for computing the DOS). We adopt ultra-soft pseudopotentials, the Perdew-Burke-Ernzerhof exchange-correlation functional 32  www.nature.com/scientificreports/ charge-density cutoff of 400 Ry. We avoid spurious interactions between periodic replicas along the out-of-plane direction z by introducing 3 nm of vacuum above the top-most layer of PtSe 2 in the unit cells. A dipole correction is also used to avoid spurious electric fields along z, using the method described in Ref. 34 and implemented in Quantum Espresso. We describe PtSe 2 using an orthorhombic Bravais lattice. We simulate the 1 T phase and use AA stacking to model multilayer structures, as they were reported to be the most stable configurations 6,35 . We optimize the lattices by relaxing atomic coordinates self-consistently until all residual forces acting on each atom are below 10 -3 Ry Bohr −1 . At equilibrium we find a cell length a x = 0.65 nm, an in-plane nearest Pt-Pt distance of 0.37 nm, an interlayer Pt-Pt distance of around 0.47 nm and around 0.26 nm layer thickness (Se-Se distance), in good agreement with experimental and calculated values available in literature 7,35,36 . The geometries and band-structures for both 1L and 2L systems are shown in Fig. 2. We extract a band-gap of 1.355 eV for the 1L case and a band-gap    7,35,37,38 . We remark that the bandgap energy found by DFT is underestimated, and that a GW correction would most likely provide larger values. This would have the beneficial effect of reducing intra-band tunneling effects when considering short channel devices. Then, we use the open-source code Wannier90 39 to generate a Hamiltonian for the 1L and 2L systems in a basis set of MLWFs, using the same k-points grid of the DFT calculations. For both cases we project the plane-wave Hamiltonian onto N w = 12 Wannier Functions to reproduce the bands shown in red in Fig. 2b,d, which span a range of over 4 eV around the mid-gap level.

Multi-scale construction of MLWF device Hamiltonian.
We build the geometry for the 2L-1L LH of PtSe 2 by embedding 10 cells of 1L PtSe 2 , tiled along x, between 2L regions, as illustrated in Fig. 3a. The obtained supercell, periodic along both x and y directions, contains a 6.4 nm long 1L region. We terminate the top layer at both interfaces with Se atoms, in order to minimize possible in-plane electric fields across the 1L region, and we check that the positions of all edge Se atoms are only slightly modified upon optimization of the geometry close to the interfaces (all residual forces were below 10 -3 Ry Bohr -1 ). From the DFT-calculated ground state of the LH we compute the on-site energy profile E on-site (x) ≡ E midgap, HS (x) along the transport direction x. This is done by following the procedure reported by Cusati et al. 40  (1) The profile of E on-site (x) for the 2L-1L PtSe 2 LH is plotted in Fig. 3b, with reference to the LH Fermi level E F, HS. This is obtained by extracting the electrostatic potential along the black A-B line indicated in Fig. 3a, with z ref fixed as the coordinate of the Pt atoms in the lowest PtSe 2 layer, and y ref fixed at three different values, namely the y coordinate of the Pt atoms labelled "Pt1" and "Pt2" in Fig. 2a, and their mid-point. Their average is taken as final reference on-site energy profile E on-site (x) in the LH and is shown as a black solid line in Fig. 3b. As can be seen, E on-site shows a ~ 0.4 eV variation in the 2L regions next to the interfaces, as a direct consequence of charge redistribution. The slight deviation from zero observed in the 2L regions, far from the junctions, is due to the fact that E F, HS is not exactly localized in the 2L midgap level.
The profile E on-site and the bulk MLWF Hamiltonians for the 2L and 1L systems introduced in the previous section are finally used as input to generate the MLWF Hamiltonian for a LH-FET structure with longer 2L and 1L regions, as schematically illustrated in Fig. 3c,d. The procedure for this Hamiltonian construction consists in two steps. The first step is to connect the MLWF Hamiltonians of the bulk components and is based on the method reported in Refs. 23,26 . This step yields a MLWF Hamiltonian for the LH in block-tridiagonal form, with arbitrary length for the 1L and 2L regions and, most importantly, a very small size (see Supplementary Information). The second step consists in mapping E on-site onto the longer LH. In particular, Fig. 3c shows how E on-site is mapped site by site in proximity of the junctions, and how E on-site is fixed to a constant-equal to its average central (extreme) values-in the sites belonging to the extended 1L (2L) regions, indicated in red in Fig. 3c.
We anticipate here that the variation of E on-site close to an interface is an electrostatic effect due to charge redistribution in its proximity, in turn occurring in response to the junction chemical environment and the band alignment between the interface components. As will be pointed out in the next section, by "freezing" this variation in the static part of the Hamiltonian, and letting the free charge in the system self-consistently adjust boundary conditions of the potential in a device simulator such as NanoTCAD ViDES 30 , we can reproduce with reasonable accuracy the effects of interface atomic details on transport. The main approximation made here is the one typical of mean field approaches, i.e., that we can neglect charge profile variations on a length scale much smaller than the spacing of the grid used for electrostatics in the considered bias voltage range.

Electrostatics in 2L-1L PtSe 2 LH-FET including ab-initio interface modelling.
In this section we present how ab-initio interface modelling ("AbInIM" for short) affects the results of transport simulations in the context of a LH-FET based on the 2L-1L PtSe 2 heterostructure discussed so far.
From the on-site energy E on-site in Fig. 3b one can draw the position of conduction and valence band edges across the LH. This is shown in Fig. 4a, in comparison with the case where mapping from DFT is carried out solely by considering differences in electron affinity (which we will refer to as "no-AbInIM" case for simplicity). From this figure one can immediately observe that the two profiles present some significant differences around the junctions, which reflect the band distortion and consequent charge redistribution naturally occurring when 2L and 1L PtSe 2 regions are in contact. While the energy barrier for electrons is not too sensitive to ab-initio interface modelling, the shape and height of the barrier for holes is significantly different. While in the "no-AbInIM" case the junction is perfectly sharp, in the "AbInIM" case its effective width is found to be roughly 2 nm, with most of the potential variation occurring in the 2L region. For electrons, we find a very close value of ΔE e ~ 0.57 eV and ~ 0.58 eV for the barrier height in the two cases, whereas for holes we obtain ΔE e ~ 0.61 eV in the "no-AbInIM" case and a remarkably higher value of ~ 0.93 eV in the "AbInIM" case, when fully mapping E on-site . This is due to the upward distortion of the valence band in the 2L part of the junction, which leads to the accumulation of holes whenever the chemical potential of the 2L system lies in its valence band. Therefore, important quantitative differences between the two approaches are expected when studying transport in a p-type LH-FET.
Based on these considerations, we studied a 2L-1L-2L PtSe 2 p-type LH-FET in the NanoTCAD ViDES NEGF-Poisson solver, with semi-infinite 2L regions embedding a 6.4 nm long 1L region (same as in DFT). These calculations were performed by fixing the valence band edge of the source and drain 50 meV above the electrochemical potential of source and drain, respectively, which corresponds to ensure charge neutrality for acceptor doping of 4 × 10 12 cm -2 of the 2L source and drain (e.g., see Supplementary Figure S1) 12,23,26,30 .
In order to demonstrate that including the E on-site profile in the NEGF Hamiltonian is reliable, we compare in Fig. 4b the charge density difference (CDD) obtained with DFT for the 2L-1L LH (solid black line) and that obtained with NanoTCAD ViDES (dashed red line) by self-consistently solving the electrostatics at equilibrium without external bias (i.e., no gate or supply voltages are applied), and mapping E on-site as explained in the previous section. In DFT the definition of CDD is the standard , where S is the yz area of the LH supercell, ρ LH is the charge density of the LH and ρ 1L (ρ 2L ) is that of the isolated 1L (2L) region in the LH supercell. In NanoTCAD ViDES the definition is the same, but with ρ LH being the free charge profile per unit volume obtained as a result of the simulation and ρ 1L (ρ 2L ) being the average free charge per unit volume in the 1L (2L) region farthest from the junctions. The DFT CDD clearly shows charge displacements in proximity of the interfaces, stronger around the 1L sites, with slightly different features between left and right ones, which in fact have slightly different atomic structures (see, e.g., Fig. 4a). These occur in a length scale smaller than inter-site distances, which represent the smaller steps of the device simulation grid. Hence, a self-consistent NEGF-electrostatic simulation would not be able to reproduce them. Indeed, it is for this reason that we propose to directly include the effect of charge redistribution on the potential, through the extraction procedure of the E on-site , as a fixed contribution (which is equivalent to "freeze" the charge displacements). Let us stress the fact that, consequently, such displacements do not appear in the NEGF-electrostatic calculation of mobile charge. Indeed, considering Fig. 4b, we see that CDD variations obtained from NEGF-electrostatics simulations are much smaller than those obtained with DFT. This is because all electrostatic contributions, that usually determine the charge density profile, are already contained Charge transport in 2L-1L PtSe 2 LH-FET including ab-initio interface modelling. In order to compare transport with and without ab-initio interface modelling in far-from-equilibrium conditions, we place two gate regions above and below the 1L region, with 0.5 nm equivalent oxide thickness (i.e. 0.5 nm of SiO 2 ), and computed the transfer characteristics for an applied drain-to-source voltage V DS = − 0.2 V, assuming fully coherent ballistic transport, acceptor doping the leads with 6 × 10 12 cm −2 , and setting the valence band edge 0.15 eV above the Fermi level in order to and ensure charge neutrality far from the junctions (see Supplementary Figure S1). The resulting I DS − V GS curves are plotted in Fig. 4c both in linear and semi-logarithmic scale. Here we call V off the value of V GS at which current is I off = 10 -1 A m −1 and we extract I on at V GS = V off + V DS . From the comparison of the two curves we can draw two important conclusions: 1. Both the subthreshold swing (SS; defined as the inverse slope of the I DS − V GS curve in semilogarithmic scale in the subthreshold regime) and the I on /I off ratio are worse in the "AbInIM" case. In particular, we find a SS of 77 mV dec −1 with a I on /I off ratio of 6.8 × 10 2 for the "AbInIM" case, and a SS of 66 mV dec −1 with a I on /I off ratio of 9.3 × 10 2 for the "no-AbInIM" case. This is a direct consequence of the higher barrier for the "AbInIM" case, already discussed in Fig. 4a and still observed in the far-from-equilibrium band diagrams of Fig. 4d. In addition to this, current on the ON state is further suppressed by the holes accumulating at the interfaces and screening other carriers propagating across the device. 2. In the subthreshold regime a higher tunnelling current is present in the "AbInIM" case. This has the important consequence of revealing the unsuitability of this device for low-power applications, which would not be possible with a "no-AbInIM" approach, where currents can reach values as low as 10 -4 A m −1 . The origin of this effect can be related to the valence band edge in the 2L part of the junction being bent upward in the "AbInIM" case, which causes accumulation of carriers and a consequent increase in tunnelling within the bias window. Such observation is confirmed by the spectral currents at V g = 0 V plotted on the right side of Fig. 4d, where the large difference between the two cases is evident. www.nature.com/scientificreports/ Finally, using the "AbInIM" approach presented above, we also assessed the performance of p-type (n-type) 2L-1L PtSe 2 LH-FETs at larger supply voltage V DS = − 0.5 V (+ 0.5 V), which is a more relevant operational condition in the light of the 2020 edition of the IRDS 28 . According to the consensus therein, future logic devices will need to be optimized towards either high-performance (HP) or low-power/high-density (HD) applications. The current in the off state must satisfy I off = 10 -2 A m −1 for HP applications, while a lower I off is required, I off = 10 -4 A m −1 , for HD applications. In each case, we define V off as the gate voltage for which I DS = I off .
We consider two different channel lengths, namely 6.4 nm and 12.8 nm. We fix the valence (conduction) band edge for the p-(n-) type FETs, 0.05 eV above (0.15 eV below) the electrochemical potential in the source and drain, which provides charge neutrality far from the heterointerfaces for a souce and drain acceptor doping of 4 × 10 12 cm −2 (donor doping of 3 × 10 13 cm −2 ). The obtained transfer characteristics for all cases are reported in Supplementary Figure S2, and the SS and I on /I off ratio estimated for HP and HD digital applications are reported in Table 1. We find that the considered LH-FETs can reach very good SS in the range 61-72 mV dec −1 and ratios I on / I off > 10 4 for the HP optimized cases, which increases above 10 5 in the HD ones. The longer channel benefits from lower SS, near the optimum 60 mV dec −1 , enabling for HD applications, contrary to the shorter channel, where tunnelling yields a relatively large current in the sub-threshold regime. We also find no significant differences in performance between p-and n-type FETs for HP applications, whereas the p-type FET clearly outperforms the n-type for HD applications, yielding almost an order of magnitude improvement for both SS and I on /I off ratio.

Charge transport in 4L-1L PtSe 2 LH-FET with sharp and smooth interfaces. A relevant advan-
tage of the method presented in this work is that it enables to model LHs with non-ideal interfaces. In the previous section we have seen the importance of mapping atomic interface details in a situation where small-gap semiconducting leads of 2L PtSe 2 are attached to the larger-gap 1L channel. We now move to analyze a 4L-1L LH of PtSe 2 , where interfaces present the additional degree of freedom of allowing for smooth junctions in the form of slightly tilted 4L-1L interfaces (see Fig. 5a). As a first approximation we decided not to undertake a rigorous study of the actual structural stability of the exposed edges of the 4L regions. As for the 2L-1L structures, we have ensured Se terminations to minimize in-plane electric fields across the 1L region, and we have relaxed the positions of all edge Se atoms until residual forces were below 10 -3 Ry/Bohr. Contrary to other state-of-the-art approaches, atomic-scale details on the channel-lead barrier can now be accounted for in the transport simulations, with the additional advantage of keeping the computational cost down to that required for a perfectly sharp junction.
All details regarding the construction and validation of the MLWF Hamiltonian for the 4L-1L PtSe 2 LH are reported in the Supplementary Information. Hereafter we report on the performance prediction of a doublegated 4L-1L PtSe 2 LH-FETs far from equilibrium, and the variations associated to the ab-initio modelling of the sharp and smooth interfaces. To this end we simulate the transfer characteristics of the device shown in Fig. 5b, with semi-infinite 4L leads and a 6.4 or 12.8-nm long 1L channel. Since the 4L system is metallic, transport is ambipolar. We apply a drain-to-source bias voltage V DS = 0.5 V and scan over both positive and negative top/ bottom gate voltages V GS , corresponding to electron and hole transport, respectively. For all devices, the top/ bottom gates are long as the channel and the dielectric layers (namely SiO 2 ) separating the channel from the gates have an equivalent oxide thickness of 0.5 nm.
The curves obtained for the 6.4 nm channel length are shown in Fig. 5c. The off-state current I off in this case is always larger than 10 -4 A m −1 for any applied V GS , due to intraband tunnelling which dominates transport across such a short channel, therefore we only assessed the performance of this particular device for HP applications. The E on-site profiles and band diagrams at V GS = 0 and V DS = 0.5 V for the LH with 12.8-nm long channel are illustrated in Fig. 5d,e, for both sharp and smooth interfaces. The calculated transfer characteristics, optimized for both HP and HD applications, are shown in Fig. 5f,g. The main figures of merit, namely the SS and the I on / I off ratio, for both HP and HD process optimizations are extracted and reported in Table 2, for all the modelled devices. We find that these do not differ significantly from the values obtained for the 2L-1L LH-FETs, reported in Table 1. Furthermore, we observe that a 3-nm smooth interface between 1 and 4L regions, rather than an atomically sharp one, does not alter the performance of LH-FETs significantly. This is especially true for the n-branch of the transfer characteristics, with only few differences in the p-branch for V GS above a certain threshold, namely V GS − V off = ~ − 0.2 V for HP and ~ − 0.4 V for HD optimizations, where current across the device with sharp junctions becomes slightly larger. We attribute this deviation to the fact that the E on-site profiles in the devices with sharp or smooth interfaces are slightly different (see Supplementary Figure S4), yielding a Schottky barrier for holes which is 0.1 eV lower in the sharp case compared to the smooth case.

Conclusion
We have presented a simulation study of FETs based on LHs of mono-multilayer PtSe 2 , based on a proposed multi-scale procedure which allows accurate ab-initio modelling of device interfaces.
As an application, we have explored various device configurations, namely with two and four layers in the multilayer regions, two different sub-15 nm channel lengths, as well as sharp or 3-nm smooth junction between 4 and 1L regions. All the results were obtained using a multiscale approach which consists in extracting the onsite energy from a DFT model of the LH and including it into a device Hamiltonian expressed in a MLWF basis, which is then used as input for the NanoTCAD ViDES device simulator 30 . This enables far-from-equilibrium transport simulations with ab-initio interface modelling, offering improved accuracy at a convenient computational cost.
We find that almost all devices, especially those with longer channels, yield a nearly optimal SS in the range 61-76 mV dec −1 and average I on /I off ratios larger than 10 4 for the HP case (larger than 10 5 for the HD case). We also find that a modest smoothening of the junction between a 1L and a multilayer region of PtSe 2 does not affect significantly the LH-FET performance. www.nature.com/scientificreports/ Given the degree of ideality of our simulations (no electron-phonon scattering, nor defects or series resistances are considered) these results are to be considered as an upper limit for device performance. The ON currents of planar devices fall short of the requirements of the most recent version of the semiconductor technology Roadmap 28 . However, the roadmap focuses for the next decade on stacked channel devices, whereas here we have considered a single channel transistor. 2D materials are very well suited for stacked channel devices 2 which provide a relatively straightforward way to achieve the ON current requirements, also considering the ON current reduction associated to the non ideality effects (e.g., finite contact resistance and scattering from phonons and defects in the channel) that we have not considered here.
From the methodological point of view, we have shown that correctly capturing charge displacements at interfaces is important when describing transport in p-type 2L-1L PtSe 2 LH-FETs. Namely, holes accumulation in their proximity is reproduced, which has an important impact on the computed current-voltage characteristics, in terms of a significantly lower I on /I off ratio and a higher tunnelling current in the sub-threshold regime. We have also demonstrated that the method allows to simulate heterostructures with non-ideal junctions, such as smooth ones between monolayer and multilayer crystals of the same material.
Despite its computational convenience, chemical details with this approach can only be accounted for indirectly, by means of variations in the ab-initio on-site potential profile. Its accuracy is not comparable to full DFTbased transport codes, which should always be adopted at least to provide an a-priori understanding of interface effects on transport for small device models. For instance, this might provide better insights on how the E onsite profile might distort as a result of doping or applied biases far out of equilibrium. Reproducing such distortions requires a proper ab-initio simulation of the system out-of-equilibrium, and falls outside the scope of this work, for which capturing equilibrium band distortions around interfaces already represents an important improvement with respect to low-cost state-of-the-art transport models, where such effects are not taken into account.
We believe that the mapping procedure presented here will play a critical role in the context of devices based on fully monolayer LHs of different 2D materials [41][42][43] , or in any structure where intense dipoles, mixed stoichiometries or other chemistry-driven effects occur at the interfaces. Table 2. Figures of merit for 4L-1L PtSe 2 FETs. All figures of merit were obtained fixing the supply voltage V DS = 0.5 V, and by mapping E on-site from DFT. The IV curves from which the data were extracted are reported in Fig. 5.