Landauer‐QFLPS Model for Mixed Schottky‐Ohmic Contact Two‐Dimensional Transistors

Abstract Two‐dimensional material‐based field‐effect transistors (2DM‐FETs) are playing a revolutionary role in electronic devices. However, before electronic design automation (EDA) for 2DM‐FETs can be achieved, it remains necessary to determine how to incorporate contact transports into model. Reported methods compromise between physical intelligibility and model compactness due to the heterojunction nature. To address this, quasi‐Fermi‐level phase space theory (QFLPS) is generalized to incorporate contact transports using the Landauer formula. It turns out that the Landauer‐QFLPS model effectively overcomes the issue of concern. The proposed new formula can describe 2DM‐FETs with Schottky or Ohmic contacts with superior accuracy and efficiency over previous methods, especially when describing non‐monotonic drain conductance characteristics. A three‐bit threshold inverter quantizer (TIQ) circuit is fabricated using ambipolar black phosphorus and it is demonstrated that the model accurately predicts circuit performance. The model could be very effective and valuable in the development of 2DM‐FET‐based integrated circuits.


Introduction
Two-dimensional material-based field-effect transistors (2DM-FETs) have attracted significant interest for their potential to continue Moore's law [1,2].Their atomic thickness and dangling-bond-free interface with gate oxide enable high tunability when applied in FET devices [3].For example, the ambipolar 2DM-FETs, which can transport the electrons and holes simultaneously, are extensively reported for a variety of channel materials, such as graphene (Gr) [4], black phosphorus (BP) [5], tungsten diselenide (WSe2) [6], and molybdenum ditelluride (MoTe2) [7].Its conduction modes (hole's mode and electron's mode) have opened a brand new way of fabricating highly efficient computational components, which could bring bonus for broad applications, including signal processing [8][9][10], logic operation [11,12], communication [13], hardware-security [14], 2DM-based memory [15] and in-memory computing [16,17].To explore the system-level applications mentioned above for 2DM-FETs, an experiment-accessible and circuit-deployable physical model that can link the lab data and industrial applications is desired by the electronic design automation (EDA) software.
To this end, computational efficiency and accuracy are equally important.It has been reported that an abnormal nonlinearity that is distinct from the traditional models' prediction will arise from the output curves if significant Schottky barriers are formed at the source/drain contact [5,[18][19][20].Unlike silicon, which can form Ohmic contacts through heavy doping, contacts between 2DMs and metals mainly belong to the Schottky category [1].Due to the contact issue, the models currently developed generally make a trade-off between physical accuracy and computational efficiency, which has become a key bottleneck in developing 2DM-FETs' models applicable for EDA tools.Since carriers pass through a series-connection of "source contact/channel/drain contact," an appropriate modeling would not describe the contact and the channel flows separately, where the former is essentially quantum and the latter remains in the semiclassical or even classical domain in most cases.Therefore, to self-consistently describe the channel current, including contact effects, it should be accomplished by seamlessly connecting two physical pictures.This is a challenging task, and according to traditional views, it would lead to a very complex model.
The previously reported models, represented by the virtual-source model [21][22][23] and the Landauer model [24][25][26][27], completely ignore the drift-diffusion process in the channel.The virtual-source model is an empirical model that affords little information regarding transport mechanisms.The Landauer model is considered accurate in describing coherent transport at the quantum regime, but it completely ignores the possible carrier scattering inside the channel.Recently reported models that simultaneously handle the contacts and channel are still based on the scheme of equivalent circuits.
One is the equivalent parasitic resistance (EPR) method, which regards the contacts as parasitic nonlinear resistors and substitutes them into the circuit to solve the resulting circuit equation [21][22][23][28][29][30][31][32][33].These equations usually have to be iteratively solved, sacrificing the efficiency.The other is the non-self-consistent model, which regards the contact as a device in parallel with the channel, calculates their currents separately, and then takes a weighted average as the equivalent current [34][35][36].Despite its remarkable efficiency, this method is an ad hoc approach that lacks a sound physical foundation.
Among them, although Ref. [34] attempted to introduce the electrostatic relation of the channel when calculating the junction energy bands of the source and drain contacts, a theoretical problem persists that the impact of the contacts on the boundary values of the quasi-Fermi levels (QFLs, or "electrochemical potential") was ignored.In other words, the contact and the channel are still treated independently in that model, and a "post-processing" had to be performed to splice the results, failing to consider the coupling between the contacts and the channel.Other works [37,38] did incorporate the contact current effect into the FET models, but the time consumption during simulation is a hindrance to their application in practical circuit design.
It turns out that an ideal model solution should comprehensively reflect the transport characteristics of the contacts and channel, while maintaining the computational efficiency that is feasible for circuit level simulation.This is an extremely difficult task and fundamentally requires innovative bottom-up theoretical methods.A breakthrough may be expected through our recently developed quasi-Fermi level phase space (QFLPS) theory [39], that can be used to describe the intrinsic channel very efficiently.This theory rigorously considers the influence of QFL splitting on transport characteristics and self-consistently includes channel electrostatics.The present work, signified by "Landauer-QFLPS modeling", further proposes that the Landauer formula describing contact transport can be cogently connected to the QFLPS model describing channel transport.A rather surprising conclusion of the work is that the Landauer formula can be written in QFLPS-model form augmented by a barrier attenuation factor based on 2D carrier state density, intimately connecting the physics of the two regions.With this theoretical result, an efficient expression of the current that does not depend on intermediate variables can be derived.
Based on the established model, three verification sections follow thereby.In "Model verification: BP-FET", we use BP, a typical ambipolar 2DM, as an example, to demonstrate the model's capability of mastering experimental data.Based on standard parameter extraction methods, complete I-V data collected from a dozen of devices have been tested.In "Tape-out verification: BP-ATIQ circuit", we tape out a three-bit ambipolar threshold inverter quantizer (ATIQ) circuit with the ambipolar-BP process to demonstrate the model's circuit simulation capability.ATIQ circuit can be used in flash analogto-digital converters to reduce chip area and power consumption.Here, the study focuses on predicting and optimizing the circuit's performance to demonstrate its value for EDA software.In "Unipolar version: MoS2-FET verification", we also verify the model against an MoS2-based FET, which typically exhibits n-type transport, to show that the coverage of the model in unipolar transport as well.
Recent studies have proposed novel contact technologies that can achieve Ohmic contacts [40][41][42], but their feasibility for industrial production is still an open question due to process-compatibility issues.For example, Bi can suppress the generation of metal-induced gap states (MIGS) due to its semi-metallic band dispersion characteristics, resulting in Ohmic contacts [42].However, both Bi and Sn are low-melting-point metals, which are incompatible with back-end-of-line (BEOL) processing, hindering the practicality of this technology.2DM transistors prepared by conventional contact processes with large-scale production potential generally exhibit Schottky characteristics.Hence, it is generally eager to model devices by combining carrier injection and channel transport [2], and this work would like to fill this gap.

Landauer-QFLPS modeling
The Landauer-QFLPS model is constructed as follows.The channel region is divided into three parts: (i) the source-contact region [  ,   ], (ii) the intrinsic channel [  ,   ], and (iii) the drain-contact region [  ,   ].With the gate-source voltage   and drain-source voltage   (assuming   > 0 for clarity) applied, the current   flows from the drain to the source (as shown in Fig. 1a) and includes both electron and hole components, i.e.,   =   +  ℎ .It is assumed that the Schottky barrier at the source blocks the injection of electrons into the channel, but a finite electron flow is present due to the thermal effect or tunneling effect, leading to the electron-injection current   on the [  ,   ] interval.At the same time, electron-QFL   drops from   at   to   at   (as shown by the blue dotted line in Fig. 1b).Similarly, holes overcoming the Schottky barrier at the drain generates a hole-injection current  ℎ on the [  ,   ], interval and the hole-QFL   's rising from   at   to   at   (as shown by the red dotted line in Fig. 1b).The energy band profile assumed here can be generalized when necessary.The internal QFLs,   and   , are implicitly determined through the conservation laws of electron and hole currents, respectively, i.e.,   =   and  ℎ =  ℎ .Once   and   are determined, the conservation quantities   and  ℎ (or   and  ℎ ) can be calculated, and thus   can be obtained.However, solving the internal QFLs is cumbersome.A central message of this work is that   can be determined without explicitly finding   and   .Here, the intrinsic channel currents   and  ℎ are already described by the QFLPS model, i.e., the integrals of carrier densities with their QFLs' paths.Therefore, the key proposal is developing QFLPS-model-like forms for the contact currents   and  ℎ , which is achieved through some intriguing derivation below.For simplicity, we here focus on the case for electrons, and the derivation for the holes is similar.As the injection of electrons from the source to the channel occurs in an atomically thin space, the relation between the amount of QFL lowering and the resulting current   should be described by the Landauer formula [24] where  is the channel width,  is the elementary charge, and ℏ is reduced Planck constant.Γ  () is the transmission energy spectrum of injected electrons at the source, while   () is the corresponding density-of-mode (DOM) energy  The transmission energy spectrum Γ  () is modeled considering the following fact: for electron energy higher than the source-contact barrier ( >   , where   represents the conduction band edge   at   and the Schottky barrier for source electrons Φ , ≔   −   , as shown in Fig. 1b), the thermal emission mechanism dominates, while for electron energy lower than the barrier ( <   ), quantum tunneling dominates.This physical picture (Fig. 1c) can be summarized by the following equation: where  = (2ℏ where   is the valley degeneracy,   * is the relative effective mass of source electrons, and   () is the collective energy spectrum of electrons, taking into account the energy-level occupancy described by the Maxwell-Boltzmann distribution in thermal equilibrium.However, what is of interest is not the thermal equilibrium but the situation where the system is driven away from thermal equilibrium by applied   .Therefore, the   () spectrum should be modified accordingly to lift the collective momentum to a higher level, as shown in Fig. 1d.With this acceleration effect,   () is written as where  , represents the elemental kinetic energy when   = 0, and Φ , represents the acceleration barrier for the source electrons, while   represents the temperature of the source contact.
The significance of Eq. (5) lies in that we connect the incoherent and coherent transports with the BAF factor.On the one hand, the integral of  in Eq. ( 5) represents the QFLPS model describing an incoherent transport and arises by integrating the drift-diffusion current along the spatial dimension.On the other hand,   originates from the coherent transport Eq. (1), where the current is given an integral over energy dimension.Therefore, the united form indicates that although the two current mechanisms are different, they can be unified in the QFL dimension.During this transformation, the unique carrier DOS of 2D materials and the abrupt transmission spectrum near the band edge play a crucial role.
After the electrons are injected from the source into the channel, the source-injection current   is converted as the intrinsicchannel current   , which is described by the QFLPS model [39], i.e., the integral of (  ⁄ )   over the interval [  ,   ].
Using the electron current conservation condition   =   , the internal QFL   can be implicitly eliminated from Eq. ( 5), which leads to the Landauer-QFLPS formula (Supplementary Note 1) The pre-factor in Eq. ( 7) is always less than 1.Hence, it reflects the contact effect included in this model.Apparently, the effect becomes more pronounced with greater CCL index   or the acceleration barrier Φ , .Therefore, the model continuously describes the transition from Schottky to Ohmic contact by varying   or Φ , parameters.Since   comprehensively incorporates nearly all the model parameters, the following analysis is mainly focused on it.
1 eV Basically,   measures the relative transport capabilities between the intrinsic channel and source contact according to its expression Eq. ( 6), where the numerator actually gives the carrier lifetime   ≔   *    ⁄ due to the intra-band relaxation in the channel, while the denominator defines an effective source-contact-electron lifetime  , for the source contact.
Significant contact effect occurs if  , is remarkably shorter than   .With channel materials fixed, properly choosing the contact metal to reduce the barriers Φ , and Φ , or increase  , can extend  , .The channel length  in the denominator actually reminds that an appropriate contact is much important in the short-channel device.We can have a rough idea for   's orders of magnitude by estimating the parameters given in Table I.
With parameters in Table I, we can estimate that   ~Φ, 2  ⁄ , which means that   will approach 0 if the Schottky barrier is lowered to a negligible level.In this case, the pre-factor in Eq. (7) will tend to 1, making the contact effect negligible and consistent with our physical intuition of device behavior.Thus, the above analysis confirms the model's validity to a certain degree.
For holes, the physical pictures drawn above need inverting on the energy dimension.For example, for drain-injection holes' the transmission spectrum Γ ℎ () (Fig. 1e), hole's thermal emission mechanism is activated when the energy is lower than the valence band maximum, while its quantum tunneling dominates when the energy is higher than it.For  ℎ (), the adaptions are similar (Fig. 1f).Hence, the Landauer-QFLPS model formula for holes can be derived as (Supplementary Note 2) where   is the hole mobility,  is the 2D-hole density, Φ ,ℎ is the acceleration barrier for the drain holes, and  ℎ is the hole's CCL index at the drain contact defined as where  ℎ * is the relative effective mass of holes, Φ ,ℎ is the thermal emission barrier for drain holes,  ℎ * is the drain holes' relative effective mass,  ,ℎ is the elemental thermal equilibrium kinetic energy of holes in the drain contact,   is the equivalent temperature of the drain contact, and Φ ,ℎ is the corresponding Schottky barrier.
In a word, the Landauer-QFLPS model gives the total current as   =   +  ℎ , where   and  ℎ are described by Eq. ( 7) and (8), respectively.They distinguish with the intrinsic QFLPS model by the respective pre-factors to describe the contact effects at the source and drain.It is worth mentioning that this unique pre-factor eliminates the conductivity divergence that occurs in the intrinsic QFLPS model for short-channel limits ( → 0) and provides a finite value.
The modeling equations for 2D-carrier densities  and  have been described in detail in the previous work [39,43].For the self-consistency of this work, these formulae are included in the Method for reference.The input parameters required there include the interface trap densities of electrons and holes,  , and  ,ℎ , and the thermal activation barrier of the channel carrier, Φ  .Other physical parameters, such as the relative effective masses of carriers,   * and  ℎ * , can be found in published literature.Therefore, a nine-parameters list {  ,  , ,   ,  ,ℎ , Φ  , Φ , ,   , Φ ,ℎ ,  ℎ } needs specified.For practical ambipolar devices, electrons and holes can share a set of model parameters that characterize the contact effect (because ambipolar transport phenomena can hardly be observed if the CCL indices of electrons and holes differ too much), and {Φ , ,   , Φ ,ℎ ,  ℎ } can be simplified as {Φ  , }, so a seven-parameters list {  ,  , ,   ,  ,ℎ , Φ  , Φ  , } can be used instead.Furthermore, for unipolar devices, such as n-type (p-type) devices,   (  ) can be set to 0, and the remaining model variables can be used to describe the electrical characteristics.
In the ideal case of the model assumptions, all model parameters should be constants.However, physical processes in practical devices, such as channel carrier scattering with lattice phonons, carrier-carrier scattering, and Coulomb scattering induced by ionized impurities and defects, result that the mobility is not constant.The local electrical environment should influence the ionization degree of interface traps for electrons and holes.As for the conductivity of the contact, it should also be regulated by the gate-source voltage, as the gate electric field can affect the barrier height.Based on the above analysis, the model parameters should vary with the external bias.Due to the physical processes' complexity, which is constrained by experimental conditions and fluctuations in environmental and material properties, it is not easy to systematically model them.However, thanks to careful consideration along the   variation dimension, a simple but efficient model parameterization strategy is proposed here, that is, the model parameters are assumed to only depend on   and can be approximated by a universal low-rank function.The coefficients of the low-rank function serve as parameters to be calibrated by experimental optimization.It has been found that linearized Gaussian wavelet functions [44] are well-suited for approximating the low-rank function.Technical details about this part are described in the Method section.

Model verification: BP-FET
BP is a classic 2DM platform demonstrating ambipolar current transport [45].Unlike the transition-metal dichalcogenides (TMDs) family, BP can maintain a direct bandgap from the single-layer to the multi-layer range.The bandgap can be tuned from 2 eV in the single layer to 0.3 eV in the bulk material, which perfectly bridges the bandgap gap between Gr and TMDs.
BP with a few-nanometer thickness obtained by mechanical exfoliation typically has a few hundred milli-electron volts bandgap and a switch ratio of about 10 3 at a drain-source bias of 0.1 V [45].Taking advantage of the ambipolar transport characteristics of BP, numerous promising applications have been reported in the literature.
Here, a batch of BP transistors was experimentally prepared by mechanical exfoliation.The device's optical microscope image is shown in Fig. 2a.The device adopts a back-gate structure (BG).The channel part is thermally evaporated with electrodes, which can be used as the gate, source, and drain of the transistor, respectively, to apply gate-source voltage   and drain-source voltage   .The fabrication process is described in the Method section.Mechanical exfoliation was chosen as the film preparation method mainly for two reasons: (i) the preparation process of mechanical exfoliation is simpler and faster, which is suitable for BP, a material that is readily oxidized by water vapor and facilitates the preservation of its intrinsic properties; (ii) the material obtained by mechanical exfoliation has inherent feature of random dispersion due to the uneven external mechanical stress introduced during the tearing and transfer process.Thus, it can provide device data with rich performance, which can be used for more harsh model testing.As shown in Fig. 2b, the BP film in the channel region was characterized by atomic force microscopy (AFM), and the thickness was determined to be approximately 9 nm (±0.7 nm).Transmission electron microscopy (TEM) combined with energy dispersive spectroscopy (EDS) characterizes ~ 2 nm natural phosphorus oxide (POx) layers existing at the interfaces of BP with Ti-and Al2O3-layers (Supplementary Note 4).Hence, the intrinsic BP thickness is around 5-7 nm and is less than the thickness characterized by AFM, which is normal for BP [46].After the device fabrication, Raman spectrum of the channel region was measured (Fig. 2c), which displayed three distinct and sharp characteristic peaks, indicating that the BP sample did not experience significant degradation during the preparation process.
We sweep the drain-source voltage   from 0 to 3 V and the gate-source voltage   from −3 V to 3 V to test the BP-FETs.
As the result, 15-sets BP-FETs' I-V data are collected.The selected voltage's range covers all possible voltage inputs within the power supply voltage   = 3 V to make a full comparison, i.e., −  <   <   and 0 <   <   .This test setting should be standard for ambipolar-device's experimental calibrations.However, it was rarely obeyed in practice [47], which limits the persuasiveness of the model's accuracy.A relevant cause is that the threshold voltage of the prepared device usually falls in a non-standard range, so the test voltage has to be shifted.In contrast, the BP-FETs in this article allows performance evaluation in the standard voltage range, and is advantageous for demonstrating the application of the model in circuit simulation, as shown in the next section.
The measured I-V data, their model simulations, and extracted parameters are given in Supplementary Note 5 due to their large data-volume.We note that the extracted η parameter exhibits a bell-shaped distribution with the lowest value about 0 and the highest value approaching 4, while most of devices'  concentrate on 1-3 interval (Fig. 2d).According to our theory, the  factor measures the degree of contact effect.The extracted distribution means that devices with nearly ideal Ohmic The drain conductance   ≔     ⁄ shows significantly different patterns among the devices, especially where the electron dominates (roughly when   > 1.2 V).For example,   for the typical case of   = 3 V is studied (Fig. 2k-m).
The experimental data for the drain conductance curve is obtained from the finite difference approximation (FDA) of the discrete data.Due to the inevitable noise signals contained in the experimental signals, the results obtained by FDA method will be affected by local fluctuations.Varying FDA's range can, to some extent, offset the noise's influence, which leads to the error bars in Fig. 2k-m.The device with a smaller η, Dev#5-6, exhibits a decreasing   -trend (Fig. 2k), which indicates the intrinsic saturation of the electron-dominated drain current for an Ohmic-contact device (Fig. 2k).In contrast,   of Dev#1-2 and Dev#2-3, undergoes a non-monotonic change (Fig. 2l and Fig. 2m), which is an important feature caused by Schottky contact.It can be seen that the model can reproduce all these conductivities well within the error range.
Accurate drain conductance is essential for circuit simulation because, as is well known, the gain of a transistor amplifier is strongly influenced by the drain conductance (whose inverse is the output impedance).Even the predicted trend would be incorrect if the contact effect is not considered correctly.As a comparison, the output and drain conductance characteristics of the η device Dev #2-3 with the most pronounced contact effect at   = 3 V were simulated using the traditional EPR method, as shown in Fig. 3.Although the equivalent parasitic resistance   was scanned over a wide range of values (10 kΩ ~ 10 MΩ), the trend obtained by EPR was still inconsistent with the experimental results, significantly inferior to the Landauer-QFLPS model.As for the drain conductance, the difference between the two is even more apparent:   predicted by EPR, which monotonically decreases with   , is entirely wrong, while Landauer-QFLPS can correctly reproduce the single-peak of the   curve rather than a monotonically changed pattern.The significant improvement in simulation quality should be attributed to the non-zero BAF.It is worth comparing the simulation results of the Landauer-QFLPS model and the EPR model with   = 10 kΩ in Fig. 3a to illustrate its role.The difference between them is most prominent near   = 0, and as   increases, the two tend to overlap gradually.In the EPR model with such a small resistance, the BAF term is equivalently set to 0, and the electron density excited by the gate field gradually decreases as the drain potential increases due to the depletion of the gate-drain field, while the current gradually shows a concavely saturated trend due to the increased lateral drain-source field.In the Landauer-QFLPS model with a non-zero BAF, a significant contact effect represented by a large BAF near   =0 strongly suppresses the output curve.When   increases to a significant level, the scattering barrier factor is almost no longer affected by   , and the trend of the corrected model is consistent with that of the intrinsic model.
In addition to accuracy, it should be noted that the model's computing efficiency is significantly superior to the EPR method.This is crucial since the Landauer-QFLPS model does not require the solution of equations in order to characterize the contact effect.By comparing the CPU running time used by a single calculation obtained by averaging cumulative timings, as shown in Fig. 3c, it was determined that the former is roughly two orders of magnitude faster than the latter for the identical task (calculating the data in Fig. 3a).To guarantee a fair comparison, non-elementary calculations required by EPR and Landauer-QFLPS computations employed the identical numerical methodologies.The speed improvement is therefore solely attributable to the time saved by not needing to solve the KCL equation.
Overall, the proposed model can satisfactorily match all the experimental measurement data obtained over the entire voltage range.In addition, this example revealed that the  factor recovered by the model quantitatively reflects varying contact effect, emphasizing the model's benefit in considering the intrinsic channel and contact holistically.

Tape-out verification: BP-ATIQ circuit
The previous section demonstrated that the developed Landauer-QFLPS model could achieve customized modeling of transistor performance for specific processes.This technique enables the prediction of device performance after circuit deployment, which is the basis for chip design using emerging 2DM transistors.Here, an ambipolar threshold inverter quantizer (ATIQ), designed theoretically in previous work [39], is chosen as a chip example to test the circuit simulation capability of the model.ATIQ can be used as a quantization element in a flash analog-to-digital converter (ADC), which has a compact structure and low power consumption advantages by fully utilizing the ambipolar characteristics of 2DMs [39].A simple summary of its principal characteristics is shown in Fig. 4a.The conventional threshold inverter quantizer (TIQ) uses a series of inverters with different flip thresholds to quantize the analog input into digital levels [58,59].Because the quantization process is completed in a single step in parallel, it is the fastest ADC architecture [60].However, the full parallel structure also rapidly increases chip area and dynamic power consumption with increasing quantized bits, thus limiting its application.ATIQ employs the physical properties of the BP-FET's ambipolar transport to strike a balance between high speed, compact size, and low power consumption.As shown in Fig. 4b, BP-FETs are stacked in series and the output voltage signal is still yielded in parallel, but most of BP-FETs (except the lowest and highest ones) in the series is multiplexed during operation since they can work as a p-type or n-type transistor, saving half of the devices compared with CMOS-based TIQ.The area savings are even more noteworthy.Previous study [39] has shown that the scaling properties of chip area with the number of bits have been compressed from the square law ( 2 ) to the quasi-linear law ( ln ) by this design, where  represents the number of quantization levels.Additionally, dynamic power consumption is significantly reduced because all the transistors share the only one current branch drawn from the power supply   .Therefore, the ATIQ circuit is a circuit design example suitable for verifying models and has practical application significance.
The chip fabricated here includes the minimum circuit configuration that can demonstrate the operations of ATIQ, which is the two-output-port structure with three BP-FETs, as shown in Fig. 4b.Its ideal voltage transfer curves are shown in Fig. 4e.Due to the different equivalent networks seen by the two ports, two different code boundary voltages  cb,1 and  cb,2 are successively generated, thereby realizing a three-bit one-hot-coded quantization.Without a doubt, if the large-area thin film preparation method is adopted, the quantization level that can be verified will be higher.Although the mechanical exfoliation method is more difficult in preparing the circuit, the advantage is that it can provide a more comprehensive check of the model, so this preparation method is still adopted here.
Based on the ambipolar-BP process, the core circuit unit of the 3-bit ATIQ chip was fabricated in the laboratory.The optical microscope photo of the chip is shown in Fig. 4c.A local back-gate (BG) was fabricated as the input port, and three flakes of the appropriate thickness (Flake I-III) were transferred to the local gate region.Among them, Flake I has a larger area, which allows the fabrication of three BP-FETs, T1a, T1b, and T1c, by sharing a common drain electrode (Port 3) so that different BP-FETs can be accessed by choosing the corresponding ports for   in subsequent circuit tests.The other two smaller flakes, Flake II and Flake III are used to construct T2 and T3, respectively, where the drain of T2 is connected to the source of T3 through Port 5, the source of T2 is connected to Port 3, and the drain of T3 is led out through Port 6 for power supply   .The high two-bit devices of the 3-bit ATIQ are fixed to T2 and T3.In contrast, the low-bit device is selected among T1a, T1b, and T1c, which are denoted as ATIQ#6-1, ATIQ#6-2, and ATIQ#6-4, respectively.The mapping relationship among circuits, devices, and ports on the chip is summarized in the embedded table in Fig. 4c.
First, the I-V data of the BP-FETs on the chip were measured to calibrate the model.The bias conditions were consistent with the previous section.The comparison between the model simulations and the experimental data is shown in Fig. 4f-j (parameters are given in Supplementary Note 6).Again, the model can accurately reproduce the experimental results for both T1a/b/c and T3, which show negligible contact effect, and for T2, which shows significant contact effects.
Next, the three 3-bit ATIQ circuits were simulated based on the calibrated model.The test scheme is summarized in Fig. 4d, and is illustrated as follows.It is assigned that the input analog signal  in was fed from the local back gate BG, the power supply voltage   was connected to Port 6, and the output voltages  o,1 and  ,2 were measured from Port 3 and Port  4k-m, demonstrating that the model predicts the three circuits all exhibiting the expected three-value quantization.Based on the simulations, the key circuit parameters, the code boundary voltages  cb,1 and  cb,2 , were extracted, as summarized in the table beside Fig. 4m.Note that the code boundary voltage of the ATIQ#6-1 configuration is higher than the other two configurations.This can be explained by the fact that the lowest transistor T1a of the ATIQ#6-1 has a higher impedance compared to T1b and T1c used in the other two circuits, consistent with the fact that T1a has the longest channel (see Fig. 4c).Optimizing the circuit performance is an important function of EDA software.Typically, the ideal code distribution should be uniform.That is, the overall quantization space is evenly divided.Hence, we should have  cb,1 =1 V and  cb,2 =2 V.However, the initial boundary codes shown above are lower than these ideal values.Although many model parameters can be adjusted to achieve the optimization goal, the most realistic optimization strategy in practical circuit design is still to adjust the channel widths.
The following calculations show that this is feasible.
Note that the drain current remains unchanged under overall proportional scaling, so we can fix the scaling factor for one device and modify the channel width design of the remaining devices.Since we need to increase the boundary codes, we need to increase the conductivity of the pull-up network, so it can be expected that the correction factors for T2 and T3 will be greater than 1.Furthermore, due to the physical half-pitch of the process, the allowed width design is not continuous on the real axis but can only take some discrete values.Obviously, the larger the device size, the more quasi-continuous the values.However, too large dimension will impair the chip area.Now we consider an extreme case, assuming that the original device width is already twice the half-pitch, so the acceptable width-amplification factor can only be a multiple of 0.5.
Under the constraints, the device widths were redesigned using the calibrated model.The results are shown in Fig. 5a-c, which indicate the magnification ratio of the new widths relative to the original values.With the new designs, a sufficiently uniform code distribution was achieved (Fig. 5d-f), with fluctuations relative to the ideal distribution not exceeding 40 mV.layered semiconductor that has attracted much attention due to its suitable mobility and bandgap.Here, the MoS2 transistors [61] were prepared based on back-gate technology for research purposes (Fig. 6a-b).
The power supply voltage was selected as 4 V. Unlike BP, the MoS2 transistor completely turns off at a negatively higher gate-source voltage, so   here is only measured till −1V.The experimental data for the transfer, output, and drain conductance characteristics of the prepared MoS2 transistor are shown in the circular plots in Fig. 6c-e.It can be seen that the output characteristics exhibit a significant contact effect, while the corresponding drain conductance curve shows a clear single-peak feature.The model extraction results (parameters given in Supplementary Note 7) show that, for unipolar devices, the model can still describe the device well and achieve first-order derivative accuracy.

Conclusion
This article reports a theoretical model, Landauer-QFLPS model, which combines the contact current with the channel current to describe a 2DM-FET.The unique advantage of the Landauer-QFLPS model is manifested in its simulation ability to reproduce the non-monotonic drain conductance observed in the experiments.The universality and practicality of the model is examined in detail to achieve a full-stack coverage from underlying physics to top-level system applications.The Landauer-QFLPS model provides critical support for ending the long-standing situation where no standard model is available for 2DM-FETs and will actively promote 2DM transistors as a new electronic device for mature large-scale system-level applications.

Method
The carrier density model, low-rank approximation formula, and experimental preparation processes are described below.

Carrier density model
The electron () and hole () densities' functional relations with their QFLs (  and   , respectively) are implicitly determined by the electrostatic-statistic relations (ESRs), which comprise three equations given as follows [39].
(i) Gauss's law for the gate-oxide-2DM channel system where Ψ is the electrostatic energy, and   =  ,  0   ⁄ is the specific gate-oxide capacitance with  , as the relative permittivity,  0 as the vacuum permittivity, and   as the gate-oxide thickness.
Efficient algorithms to solve the ESRs defined by Eqs.(10)-( 12) and the integrals in Eqs.(7) and ( 8) have been reported by our previous numerical work [43].Hence, there is no obstacle to planting the model into a standard circuit simulator.

Low-rank approximation formula
An  -th (  ≥ 2 ) order linearized Gaussian Wavelet (LGW) with  -extension on the bounded interval ℬ = [ , ,  , ] is constructed as follows [44]  , ( where ℎ{  } 1≤≤ () represents the linear interpolation function of the control points {  ,   } 1≤≤ located on interval ℬ, and the interpolation function value outside the interpolation interval ℬ is set to 0. Although the integral form is employed by Eq. (15), the value can be evaluated analytically with the help of the error function (given in Supplementary Note 3), which possesses a standard fast algorithm by rational Chebyshev approximations [62].

Control points setting
In principle, both the {  } 1≤≤ and {  } 1≤≤ coordinates can be selected as optimization variables, but this would significantly increase the difficulty of the optimization convergence.Therefore, we evenly distribute the {  } 1≤≤ coordinates on the interpolation interval ℬ and leave only the {  } 1≤≤ coordinates as optimization variables.

Rank optimization
In principle, the higher the value of , the more accurate the approximation.However, due to considerations of convergence cost, we limit  to be less than or equal to 3 in this case.

Boundary effect
We employed the odd-extension technique in the computation here to suppress the well-known boundary effects [63].

BP-FET
First, deposit a 35 nm Pd/3.5 nm Ti metal stack as a buried gate electrode on a 300 nm SiO2 substrate.Then, by atomic layer deposit (ALD), a 20 nm Al2O3 dielectric layer was prepared.Next, use a mechanical exfoliation process to peel BP from the bulk material into thin layers and transfer it onto the prepared buried gate substrate using a dry transfer method.Next, a 3.5 nm Ti/35 nm drain-source electrode stack was defined using electron beam lithography.Finally, a 20 nm Al2O3 passivation layer was encapsulated to protect the BP channel.

MoS2-FET
First, we sputtered a 10 nm TiN electrode onto a 300 nm SiO2 substrate and left a predefined pattern with a lift-off process.
Next, a 15 nm HfO2 dielectric layer was deposited with ALD as the gate insulator, followed by a spin-coating process to deposit a uniform layer of MoS2 nanosheet [61].Heat the sample in an N2 atmosphere for 1 hour at 300 ℃.Finally, we used electron beam evaporation to deposit source/drain Ti/Pd electrodes.
where  , is the elemental thermal kinectic energy for electron in the source metal,   is the effectively local temperature for the source junction, and Φ , is the acceleration barrier for the electrons injected from the source junction.

C. Formulization of 𝐼 𝑒𝑠
Given the definitions for Γ  () (Eq.(3) and (7)) and   () (Eq.(8)), one can introduce an index function () to simply the calculation so that one can obtain the final result for  ℎ by substituting the second integral above with Eq. (48).V DS = 0.3 V : 0.27 V : 3.0 V

Fig. 1 |
Fig. 1 | The schematic diagram of the Landauer-QFLPS model principle.a. 2DM-FET device structure and electrical testing schematic; b.Energy band diagram of current transport mechanism; c.Source electron transmission energy spectrum model; d.Source electron collective kinetic energy spectrum model; e. Drain conductance hole transmission energy spectrum model; f.Drain conductance hole collective kinetic energy spectrum model.
contact and devices with significant Schottky contact are rare.Most of devices are in an intermediate stage.Actually, both Ohmic and Schottky contacts for 2DM FETs have been experimentally reported[5, 13, 18-20, 38, 48-53].Experiments within similar materials systems can even exhibit both types of devices[6,49,[54][55][56][57], too.Hence, it is essential to consider the coexistence of Ohmic and Schottky contacts in modeling practical devices.Specifically, the experimental data and their model simulations of three devices with representative η values, Dev#5-6 (with the lowest η), Dev#1-2 (with medium η), and Dev#2-3 (with the highest η), are compared and presented.As shown in Fig.2e-m, they are the transfer (Fig.2e-g), output (Fig.2h-j), and drain conductance (Fig.2k-m) curves of the three devices, respectively.The simulated transfer curves in Fig.2e-g reflect the matching on the global current range, while the output curves (Fig.2h-j) focus on the matching for the on-state current.Their transfer and output curves show a typical ambipolar feature, i.e., the hole current dominates and decreases with increasing   when   < 0 , then, the electron current dominates and increases with increasing   when   > 0. The simulation results show that in the complete voltage range, the model simulation results can highly match the experimental test data.

Fig. 2 |
Fig. 2 | Experimental Verification ---BP-FETs.a. Optical microscopy images of the prepared devices, with white dashed boxes highlighting the channel region used in the device; b and c show the AFM and Raman characterizations, respectively.AFM measurements were done for the colored dashed lines in a. Raman was performed for a region pointed by a yellow circle in a; d.The distribution of the model parameter  in the device, where the color of the block represents the exact value for ; With the classification in d, three typical- devices, Dev#5-6, Dev#1-2, and Dev#2-3, were selected, and their transfer (e-g), output (h-j), and   = 3 V drain conductance (k-m) curves are shown in the figure.

Fig. 4 |
Fig. 4 | Circuit chip design verification: ATIQ circuit.a. Threshold inverter quantizer (TIQ) circuit structure based on conventional CMOS process; b.ATIQ circuit structure based on ambipolar-BP process; c. optical microscope image of 3-bit ATIQ chip die, where the surface contact electrode is colored with a red-yellow gradient to visually distinguish it from the back gate (pale yellow), and the light green color is the BP thin layer material; d. Circuit test signal code table; e.The ideal electrical characteristic curve of the circuit; f-j are the model calibration results of the on-chip BP transistor's device output characteristics; k-m are the predicted circuit operating characteristics based on the calibrated model, and the table on the right summarizes the predicted boundary code element voltages; n-p and the corresponding actual test results are shown in the table on the right.
5, respectively.To test ATIQ#6-1, one should connected   to Port 1 while floating Port 2 and Port 4; To test ATIQ#6-2, one should connected   to Port 2 while floating Port 4 and Port 1; And, to test ATIQ#6-4, one should connected   to Port 4 while floating Port 1 and Port 2; The circuit simulation results are shown in Fig.

Fig. 5 |
Fig. 5 | Optimizing 3-bit ATIQ circuits based on the Landauer-QFLPS model.a-c are schematics for the width optimizations of the three circuits ATIQ#6-1, ATIQ#6-2, and ATIQ#6-4, respectively.d-f are the optimized circuit operating curves.The table on the right summarizes the optimized boundary codes.The measured data of the circuits are shown in Fig. 4n-p.Comparing with the simulations in Fig. 4k-m, it is evident that the model's simulations successfully reproduce the experimental data.The extracted code boundary voltages from the experiment show an error within 50 mV compared to the model's extracted results, thus quantitatively verifying the model's circuit simulation capability.

Fig. 6 |
Fig. 6 | Verification of MoS2 devices.a, Schematic diagram; b, Raman spectrum; (c-e) Comparison of simulated and measured data for the device's transfer, output, and drain conductance characteristics.
FET device is studied in the following to demonstrate this model for unipolar devices.MoS2 is a typical n-type 2D

Fig
Fig. S3 EDS mapping images

Table I
Parameter estimations for With the aid of (), one can go further to simplify the   in Eq. (1) as Re-write the Fermi window in the square bracket as the integral of the derivative as of states for conduction band electrons in the channel semiconductor with the spin valley degeneracy set as 1, distribution function  can be expressed by the partial derivative of  as Since a leap of the derivative arises for () (c.f.Eq.(11)),  ′′ () yields a Dirac delta function as With the definition of () and noting that   −   = Φ , , one further obtains where approximation (  ,   ) ≈ (  ,   ) has been made, and the index   is defined as one can obtain the final result for   by substituting the second integral in Eq. (25) with Eq. (23).With the definition of () and noting that   −   = −Φ ,ℎ , one further obtains Simulated output/transfer curves in linear/logarithm scale benchmarked with experimental data are shown below