Bilayer TMDs for Future FETs: Carrier Dynamics and Device Implications

1. Introduction

The semiconductor industry has experienced extraordinary growth, initially driven by silicon, which can form large, single-crystal structures ideally suited to modern electronic manufacturing. Since the 1970s, silicon’s advantageous properties and scalability have enabled continuous advancements, allowing transistor densities to increase exponentially—a trend famously captured by Moore’s Law [1]. The miniaturization of silicon-based devices has reduced costs, boosted speed, and lowered power consumption, making it the foundation of digital technology. However, as silicon transistor dimensions have approached their physical limits, it has become challenging to further downscale without encountering fundamental limitations [2,3]. These issues include short-channel effects, increased leakage currents, and degradation of device reliability, which hinder the performance of ever-smaller silicon transistors [4,5].

Motivated by the need for high-performance devices beyond silicon’s capabilities, researchers have turned to alternative materials [6,7,8,9,10]. Since the early 2000s, the search for novel materials to sustain and surpass silicon’s legacy in both classical and unconventional computing has intensified. These new materials, especially those with lower dimensionality, offer unique physical properties that address the challenges faced by silicon at nanoscale dimensions. One class of promising candidates is two-dimensional (2D) materials [11,12,13], which, unlike bulk materials, are formed from layers just one atom thick. These materials exhibit unique electrostatic control due to their thin nature, allowing for reduced short-channel effects and exceptional gate control.

Among two dimentional materials, Transition Metal Dichalcogenides (TMDs) [14,15,16] have garnered significant interest. Structurally TMDs consist of a transition metal layer sandwiched between two layers of chalcogen atoms, with each M-X-M atomic layer forming strong covalent bonds within the layers but weak van der Waals bonds between them. This unique structure enables TMDs to be scaled to atomic thicknesses while maintaining their semiconducting properties, providing a band gap that graphene lacks and allowing for effective switching in field-effect transistors (FETs). The versatility and tunability of TMDs, along with their stable semiconducting properties, make them ideal for future applications in both energy-efficient and high-performance devices.

Bilayer transition metal dichalcogenides (TMDs) offer unique advantages over monolayers for applications in next-generation transistors and electronic devices. They exhibit a higher density of states, enabling greater charge storage, as well as improved carrier mobility for efficient current flow [17,18].Bilayer structures enhance fabrication yields [19] by protecting the bottom layer during processing, minimizing defects, and providing better stability under environmental stress. In particular, WSe₂ and WS₂ have attracted significant attention for pFET applications due to their intrinsically high hole mobility. With their potential to support further device miniaturization and improved performance, bilayer TMDs represent a promising path for advancing both transistor technology and the broader field of nanoelectronics.

To evaluate the transport properties of bilayer TMDs, it is crucial to consider the dielectric environment, as real-world devices often feature these materials supported by substrates and surrounded by gate dielectrics. This work employs a combination of first-principles and semiclassical computational methods to investigate carrier transport in bilayer TMDs. Density Functional Theory (DFT) and Density Functional Perturbation Theory (DFPT) as implemented in Quantum ESPRESSO [20,21] are used to calculate the full-band electronic structures and phonon dispersions, while electron–phonon scattering rates are obtained using the Electron–Phonon Wannier (EPW) [22,23,24] approach. These first-principles results are incorporated into semiclassical full-band Monte Carlo simulations, solving the Boltzmann transport equation and accounting for intrinsic electron–phonon scattering, screening by the dielectric environment, scattering with the hybrid plasmon/phonon interface excitations, and with ionized impurities. By addressing these critical factors, this study provides a robust framework for understanding and optimizing the electronic transport properties of TMDs, helping their integration into next-generation nanoelectronics.

Our discussion is organized as follows: Section 2 outlines the geometry of the double-gate transistor considered and the theoretical framework underlying our study. Section 3 and Section 4 describe the physical models used for the band structure, phonon dispersion, and electron–phonon scattering rates, as well as the treatment of IPP and impurity scattering. Finally, Section 5, Section 6, and Section 7 present and analyze our transport and FET results.

2. Device Geometry and Monte Carlo Simulation Framework

In our work we have considered a double-gate field-effect transistor (FET) structure with a bilayer TMD as the channel material, shown in Figure 1. The device has a $L_{\mathrm{ch}}=10$ nm channel, with source and drain extensions of equal length $L_{\mathrm{S}/\mathrm{D}}=10$ nm. The source/drain extensions are doped, while the channel is left undoped. A self-aligned $L_{\mathrm{g}}=10$ nm gate is centered on the channel.

The top and bottom gate insulators each result in an equivalent oxide thickness (EOT) of 0.7 nm. This double-gate structure captures the impact of the dielectric environment (substrate and top oxide) on carrier transport, including interface phonons, interface plasmon–phonon (IPP) modes (often improperly called “remote phonons”), and dielectric screening induced by free carriers and by the environment.

Electronic transport is treated by solving the two-dimensional Boltzmann Transport Equation (BTE) self-consistently with the Poisson equation using an ensemble Monte Carlo (MC) method. The BTE describes how the carrier distribution evolves over time in real and momentum space under electric fields and scattering. We assume the device to be infinitely wide, such that the gate is much wider than the channel length; this allows us to solve the 2D Poisson equation and neglect edge effects.

The electrostatic potential is obtained by solving the 2D Poisson equation, which accounts for charge distribution and material permittivity. This solution updates the electric field profile used for carrier trajectories in the Monte Carlo simulation. The spatial domain is discretized using non-uniform meshes: one along the transport direction (x) and another along the out-of-plane direction (z), which is refined near the TMD and oxide interfaces. The Poisson equation is solved using the finite element method adapted from [25] with the following form:

$$\nabla ·\mathbf{D}\left(\mathbf{r}\right)=-\nabla ·[\epsilon \left(\mathbf{r}\right)·\nabla V\left(\mathbf{r}\right)]=\rho \left(\mathbf{r}\right)=e[N_D\left(\mathbf{r}\right)-n_{\mathrm{el}}\left(\mathbf{r}\right)],$$

(1)

where $V\left(\mathbf{r}\right)$ is the electrostatic potential, $\rho \left(\mathbf{r}\right)$ is the charge density, $N_D$ is the doping profile, and $n_{\mathrm{el}}$ is the electron density. The dielectric tensor is defined piece-wise for each element. In the TMD, we use an anisotropic dielectric tensor:

$$\epsilon =\left[\begin{array}{cc}{\epsilon }_{\Vert }& 0\\ 0& {\epsilon }_{\perp }\end{array}\right],$$

(2)

where ${\epsilon }_{\Vert }$ and ${\epsilon }_{\perp }$ are the in-plane and out-of-plane dielectric constants, respectively. In contrast, the oxide layers use isotropic dielectric constants. The TMD is discretized along z using two mesh elements, each of thickness $\Delta z=d/2$, where d is the total bilayer thickness. The cloud-in-cell method [26] is employed to map the charge of Monte Carlo particles onto the grid nodes at the center of the TMD layer. Dirichlet boundary conditions are assumed at the source and drain terminals, while Neumann conditions are applied elsewhere.

The MC algorithm alternates between two steps: (i) carriers drift through the device under the electric field and (ii) scattering events are considered based on pre-tabulated rates, which alter energy, momentum, or direction. Carrier positions are then updated, the charge density is recalculated, and the Poisson equation is solved again. Convergence is reached once steady-state transport is established, typically by tracking 2,000–5,000 carriers over several picoseconds to obtain statistically reliable current–voltage and field–energy characteristics. After each MC iteration and before the next free flight, we record carrier energy, velocity, current, charge density, and electrostatic potential. These are ensemble-averaged over the total simulation time:

$$\langle A\rangle ={\displaystyle \frac{1}{{\tau }_{\mathrm{sim}}}}\sum _i{\tau }_i{\langle A\rangle }_{{\tau }_i},$$

(3)

where ${\tau }_{\mathrm{sim}}$ is the total simulation time and ${\tau }_i$ the i-th free flight duration. In the low-field regime, mobility is extracted from the diffusion constant $D_{\theta }$ using the Einstein relation in the non-degenerate limit:

$$D_{\theta }={\displaystyle \frac{1}{2}}{\displaystyle \frac{d}{dt}}〈{\left(x_{\theta }-\langle x_{\theta }\rangle \right)}^2〉,{\mu }_{\theta }={\displaystyle \frac{e{D}_{\theta }}{k_{B}T}},$$

(4)

where $x_{\theta }$ is the carrier position along the transport direction $\theta$. This approach minimizes stochastic noise when the drift velocity is much smaller than the thermal velocity. At high fields, the steady-state drift velocity, carrier energy, and field-dependent mobility are computed as

$$v_{\mathrm{drift}}={\displaystyle \frac{1}{N_p}}\sum _i{v}_i,\langle E\rangle ={\displaystyle \frac{1}{N_p}}\sum _i{E}_i,$$

(5)

where $N_p$ is the total number of simulated carriers, $v_i$ is the component of the velocity along transport direction (x-axis), and $E_i$ the kinetic energy. Device-level simulations also yield drain current and spatial profiles of charge density, potential, electric field, and carrier energy for each gate and drain bias point. These data allow us to construct I–V characteristics and analyze how scattering mechanisms affect carriers across the device.

We perform Density Functional Theory (DFT) calculations to obtain the bilayer TMD band structure and phonon dispersion. The band structure reveals key electronic properties, such as effective masses and valley configurations that govern carrier kinematics, while the phonon dispersion characterizes the lattice vibrations responsible for intrinsic scattering processes. From these fundamental quantities, we tabulate the rates for electron scattering with the TMD phonons, with IPPs, and with ionized impurities (dopant atoms). The following sections provide details on the scattering mechanisms considered in these simulations.

3. Band Structure and Phonon Spectrum

Band structures of freestanding bilayer TMDs were calculated using density functional theory (DFT) as implemented in Quantum ESPRESSO [20,21]. We employed the Optimized Norm Conserving Vanderbilt (ONCV) pseudopotentials [27] and GGA-PBE exchange-correlation functional [28,29]. The structures were relaxed until atomic forces were below ${10}^{-4}$ eV/nm, ensuring energy minimization (Figure 2). A vacuum spacing of approximately 1.4 nm was used between supercells to prevent unphysical interlayer interactions.

The phonon dispersions were computed using density functional perturbation theory (DFPT) [30] with a coarse $\mathbf{k}$-point mesh of $12\times 12\times 1$ and $\mathbf{q}$-point mesh of $6\times 6\times 1$. To obtain accurate phonon energies, especially for low-energy acoustic modes, Wannier interpolation was performed on a finer $30\times 30\times 1$ mesh using the EPW package [22,23,24]. Spin–orbit coupling (SOC) was included in all calculations as it influences the electronic structure and carrier mobility. The DFT input parameters are summarized in Table 1.

Figure 3 shows the calculated band structure and phonon dispersion of bilayer WS₂ along the high-symmetry path M-$\Gamma$-K-M. Bilayer WS₂ exhibits a direct bandgap of 1.56 eV with the conduction band minimum at the K-point and a satellite conduction valley at Q, 58 meV higher in energy. The electron effective masses are isotropic: 0.3 $m_e$ at K and 0.74 $m_e$ at Q. The valence band maximum lies at K with a local maximum at $\Gamma$ 87 meV lower, with hole effective masses 0.27 $m_h$ at K and 2.1 $m_h$ at $\Gamma$.

Previous literature reports bilayer WS₂ as an indirect-bandgap material with the conduction band minimum (CBM) at $\Gamma$. We observed, however, that the bandgap character strongly depends on the interlayer spacing (Figure 4). At an interlayer distance of 0.68 nm, the CBM shifts to $\Gamma$, producing an indirect gap. Reducing the spacing to 0.64 nm lowers the Q valley relative to K, restoring a direct gap between Q and $\Gamma$. This sensitivity of valley ordering to interlayer distance highlights the delicate balance between competing conduction-band valleys in bilayer WS₂. Even small structural changes such as strain alter the nature of the bandgap.

Experimentally, angle-resolved photoemission spectroscopy (ARPES) can probe the valence bands, but direct identification of the CBM remains challenging. This valley sensitivity has direct consequences for transport, as different conduction-band minima (K, Q, or $\Gamma$) result in varying effective masses, density of states, and intervalley scattering, which can markedly influence carrier mobility and device performance.

4. Scattering Mechanisms

4.1. Electron–Phonon Scattering

Intrinsic electron–phonon scattering is included using first-order perturbation theory. The matrix element for a transition from band m at wavevector $\mathbf{k}$ to band n at $\mathbf{k}+\mathbf{q}$ via phonon mode $\nu$ is

$$g_{mn}^{\nu }(\mathbf{k},\mathbf{q})=\sqrt{{\displaystyle \frac{\hslash }{2{\omega }_{\mathbf{q}\nu }}}}\langle u_{n,\mathbf{k}+\mathbf{q}}\left|\Delta V_{\mathbf{q}\nu }^{\mathrm{SCF}}\right|u_{m,\mathbf{k}}\rangle ,$$

(6)

where ℏ is the reduced Planck constant, ${\omega }_{\mathbf{q}\nu }$ is the frequency of phonon mode $\nu$ at wavevector $\mathbf{q}$, and $\Delta V_{\mathbf{q}\nu }^{\mathrm{SCF}}$ is the change in the self-consistent potential due to the phonon. The electronic eigenstates $|u_{n,\mathbf{k}}\rangle$ are used in evaluating the matrix elements. Matrix elements are computed using density functional perturbation theory (DFPT) and interpolated using maximally localized Wannier functions (MLWFs) via EPW [22,23,24]. The scattering rate is obtained using Fermi’s golden rule:

$${\displaystyle \frac{1}{{\tau }^{\left(\nu \right)}(\mathbf{k},n)}}={\displaystyle \frac{2\pi }{\hslash }}\sum _{{\mathbf{k}}^{\prime },n^{\prime }}{\left|g_{mn}^{\nu }(\mathbf{k},\mathbf{q})\right|}^2\left(N_{\mathbf{q}\nu }+\frac{1}{2}\pm \frac{1}{2}\right)\delta \left(E_{n^{\prime }}\left({\mathbf{k}}^{\prime }\right)-E_n\left(\mathbf{k}\right)\pm \hslash {\omega }_{\mathbf{q}\nu }\right),$$

(7)

where ${\tau }^{\left(\nu \right)}(\mathbf{k},n)$ is the carrier lifetime for phonon mode $\nu$, $g_{mn}^{\nu }(\mathbf{k},\mathbf{q})$ is the electron–phonon coupling matrix element, $N_{\mathbf{q}\nu }$ is the Bose–Einstein occupation number representing phonon population, and $E_n\left(\mathbf{k}\right)$ denotes the energy of the electronic state in band n at wavevector $\mathbf{k}$. In all our calculations scattering with out-of-plane acoustic (ZA) phonons is neglected because 2H-phase TMDs possess horizontal mirror (${\sigma }_h$) symmetry, which forbids first-order carrier–ZA interactions [31]. Higher-order (two-phonon) processes are much weaker in these materials and therefore contribute negligibly to carrier scattering [32]. The sum over final states $({\mathbf{k}}^{\prime },n^{\prime })$ accounts for all possible scattering processes to different bands and wavevectors. The final-state density of states is evaluated using a 2D Gilat–Raubenheimer scheme [33,34]. Energy-resolved scattering rates are obtained by summing the contributions from all electronic states within each energy bin and normalizing by the density of states, yielding average rates that are further separated into acoustic and optical phonon contributions. Figure 5 presents phonon-limited scattering rates for electrons (a) and holes (b) in bilayer WS₂, showing the contributions of acoustic and optical phonon modes plotted against carrier energy.

4.2. Dielectric Screening

The top and bottom dielectrics (dielectric environment) screen the electron–phonon interaction [35] and a full DFT calculation of the electron–phonon matrix elements for the entire double-gate structure would require an impractically large supercell. To circumvent this, we employ a semi-empirical approach [35]: the matrix elements are first computed for a free-standing 2D layer using Quantum ESPRESSO and EPW. Since the Hartree term dominates the scattering potential, the free-standing Green’s function (implicitly used in DFT) is replaced with one that satisfies the boundary conditions of the double-gate structure, including both dielectrics. This modification rescales the DFT scattering rates by the squared ratio of the two Green’s functions:

$${\left|{\displaystyle \frac{{\mathcal{G}}_{Q,{\omega }_Q^{\left(\eta \right)}}^{\left(\mathrm{env}\right)}(d,d)}{{\mathcal{G}}_Q^{\left(\mathrm{vac}\right)}(d,d)}}\right|}^2,$$

(8)

where d is the position at the center of the layer (along the z-direction normal to the plane). Here, ${\mathcal{G}}_{Q,{\omega }_Q^{\left(\eta \right)}}^{\left(\mathrm{env}\right)}(d,d)$ is the Green’s function of the Poisson equation evaluated at $z=d$ for the dielectric environment of the double-gate geometry, at wavevector Q and phonon frequency ${\omega }_Q^{\left(\eta \right)}$. ${\mathcal{G}}_Q^{\left(\mathrm{vac}\right)}(d,d)$ is the corresponding Green’s function for the free-standing layer. The Green’s function ${\mathcal{G}}_{Q,{\omega }_Q^{\left(\eta \right)}}^{\left(\mathrm{env}\right)}(d,d)$ incorporates the effect of the 2D electron gas (2DEG) in the TMD layer, so the procedure described also captures free-carrier screening. A thorough discussion of the assumptions behind the scaling factor in Eq. (8) is provided in [35]. Implemented in the Monte Carlo simulations via the rejection method, this rescaling reduces the scattering rates relative to the free-standing case, leading to higher carrier mobility compared to the freestanding layer.

4.3. Interface Plasmon–Phonon (IPP) Scattering

The plasma oscillations of electrons in the bilayer couple with the optical phonons of the top-gate dielectric, the substrate, and the polar 2D layer, giving rise to hybrid interface plasmon–phonon (IPP) modes. This hybridization is particularly important as it significantly affects the scattering rates, especially at small momentum transfer Q. The full theoretical framework is presented in [35] and [36], which provide all necessary details; here, we adopt the notation from [35] and summarize the main results for convenience. Because electron–IPP interactions transfer momentum only through the phonon-like content of each hybrid mode, we must evaluate the phonon content ${\Phi }^{\left(\alpha \right)}\left({\omega }_Q^{\left(i\right)}\right)$ of phonon $\alpha$ (TO1, TO2, TO3, TO4, or ZO) for every mode i [37,38]. In terms of this quantity, the scattering potential amplitude ${\mathcal{A}}_{Q,{\omega }_Q^{\left(i\right)}}^{\left(\alpha \right)}$ for phonon $\alpha$ of mode i is given by [35]:

$$\begin{array}{c}{\left|{\mathcal{A}}_{Q,{\omega }_Q^{\left(i\right)}}^{\left(\alpha \right)}\right|}^2={\Phi }^{\left(\alpha \right)}\left({\omega }_Q^{\left(i\right)}\right){\left(1-e^{-2Q{t}_{\mathrm{b}}}\right)}^2\hfill \\ \hfill {\displaystyle \frac{e^2\hslash {\omega }_Q^{\left(i\right)}}{2Q}}\left|{\displaystyle \frac{1}{{ϵ}_{\mathrm{TOT}}^{(\alpha ,\mathrm{high})}(Q,{\omega }_Q^{\left(i\right)})}}-{\displaystyle \frac{1}{{ϵ}_{\mathrm{TOT}}^{(\alpha ,\mathrm{low})}(Q,{\omega }_Q^{\left(i\right)})}}\right|.\end{array}$$

(9)

In Eq. ( 9), ${ϵ}_{\mathrm{TOT}}(Q,\omega )$ is the dielectric function of the system. The superscripts high and low correspond to the absence and presence, respectively, of the response from phonon $\alpha$, allowing us to isolate its contribution to the Fröhlich-like potential. A complete derivation of Eq. ( 9) and the resulting scattering rates can be found in Ref. [35].

As discussed in Ref. [35], scattering with interface hybrid excitations controls carrier transport in high-$\kappa$ environments. In such dielectrics, this mechanism counteracts the benefits of screening and suppresses the mobility of the 2D channel, often below that of the free-standing case. The reduction becomes more severe as the dielectric constant of the top-gate insulator increases, e.g., when replacing hBN with ZrO₂.

4.4. Landau Damping

Landau damping occurs when plasmons fall inside the single-particle excitation continuum, i.e., $E(k_F-Q)-E\left(k_F\right)\le \hslash {\omega }_P\left(Q\right)\le E(k_F+Q)-E\left(k_F\right)$ ($k_F$ being the Fermi wave vector), so that they are no longer well-defined eigenmodes of the system and decay into incoherent electron–hole excitations. The dielectric response of the 2D TMD monolayer is described by

$${\epsilon }_{2D\perp }(Q,\omega )={\epsilon }_{2D\perp }^{(\infty )}\left[1-e^2{G}_Q(h/2,h/2){\Pi }_{2D}(Q,\omega )\right]+\left({\epsilon }_{2D\perp }^{\left(0\right)}-{\epsilon }_{2D\perp }^{(\infty )}\right){\displaystyle \frac{{\omega }_{ZO}^2}{{\omega }_{ZO}^2-{\omega }^2}},$$

(10)

where ${\epsilon }_{2D\perp }^{\left(0\right)}$ and ${\epsilon }_{2D\perp }^{(\infty )}$ are the static and optical out-of-plane dielectric constants of the monolayer, h is the effective thickness, $G_Q(z,z^{\prime })$ is the in-plane Fourier transform of the Poisson Green’s function of the device geometry (evaluated at $z=z^{\prime }=h/2$), ${\Pi }_{2D}(Q,\omega )$ is the 2D carrier polarizability given by Stern [39], and ${\omega }_{ZO}$ is the frequency of the out-of-plane optical (ZO) phonon. In mirror-symmetric (2H) monolayers, the ZO-induced potential is antisymmetric and the electron–phonon matrix element vanishes to first order, so electrons do not couple to this component of the hybrid modes. In the long-wavelength limit ($Q\to 0$), the carrier contribution in Eq. (10) can be approximated by

$${\epsilon }_{2D}^{\mathrm{el}}(Q,\omega )\approx {\epsilon }_{2D}^{(\infty )}\left(1-{\displaystyle \frac{{\omega }_P{\left(Q\right)}^2}{{\omega }^2}}\right),$$

(11)

where the 2D plasmon frequency is defined as ${\omega }_P{\left(Q\right)}^2=e^{2}nQ/2{\epsilon }_{2D}^{\infty }{m}^{*}$, with n the carrier density and $m^{*}$ the effective mass. The onset of Landau damping is set by the cutoff wave vector ($Q_{LD}$):

$$Q_{LD}={\left[k_F^2+{\displaystyle \frac{2m^{*}}{{\hslash }^2}}{\omega }_P\left(Q_{LD}\right)\right]}^{1/2}-k_F,$$

(12)

where $k_F$ is the Fermi wave vector. For $Q<Q_{LD}$, plasmons remain coherent and dynamically screen the phonon-like hybrid modes; for $Q>Q_{LD}$, we replace the dynamic long-wavelength term $1-{\omega }_P{\left(Q\right)}^2/{\omega }^2$ in Eq. 11 with the static Thomas-Fermi expression ${\epsilon }_{TF}\left(Q\right)=1+Q_{TF}^2/Q^2$ where $Q_{TF}$ is the screening wave vector give by,

$$Q_{TF}={\displaystyle \frac{e^{2}n}{2{\epsilon }_{2D}{k}_{B}T}}\left(\mathrm{non}\text{-}\mathrm{degenerate}\right),Q_{TF}={\displaystyle \frac{e^2{m}^{*}g}{2\pi {\hslash }^2{\epsilon }_{2D}}}\left(\mathrm{degenerate}\right),$$

(13)

with g the valley degeneracy. This framework captures the crossover correctly: Below $Q_{LD}$, plasmons dynamically hybridize with phonon-like modes, whereas above $Q_{LD}$, only statically screened phonon-like branches survive. We consider two double-gate dielectric configurations — hBN/bilayer TMD/SiO₂ (low-$\kappa$) and HfO₂/bilayer TMD/SiO₂ (high-$\kappa$) — to capture the range of remote phonon scattering effects expected in realistic device environments. Figure 6 shows the dispersion of fully hybridized IPP modes in the SiO₂/WS₂/HfO₂ stack, revealing strong phonon-plasmon coupling at small wave vectors Q. Above the Landau damping cutoff $Q_{LD}$, plasmons are suppressed and only surface optical (SO) phonons remain.

Figure 7 shows the hole-IPP scattering rates separated by wave vector regimes ($Q\le Q_{LD}$ for hybridized modes and $Q>Q_{LD}$ for decoupled SO modes). Scattering is dominated by long-wavelength hybridized modes, especially in the HfO₂ configuration due to stronger coupling. In contrast, the hBN stack has weak IPP scattering due to its weak coupling.

Overall, these results underscore the fundamental mobility trade-off when employing high-$\kappa$ dielectrics, such as HfO₂: While they provide stronger electrostatic gate control, this advantage is offset by increased remote phonon scattering that degrades carrier mobility. Since IPP scattering rates depend on the local carrier density—which varies spatially within a FET—we pre-tabulate these rates for five carrier densities ranging from $1\times {10}^{10}$ to $1\times {10}^{14}$ cm⁻². During the Monte Carlo simulations, the local carrier density at each point is used to interpolate the appropriate IPP scattering rate on a logarithmic scale, ensuring accurate representation of spatially varying IPP rates.

4.5. Impurity Scattering and Dielectric Impact

The impact of charged impurities in the TMD channel is evaluated by solving the Poisson Green’s function for the double-gate geometry considered in this work, following Ref. [40]. The corresponding screened Coulomb potential for a charge located at the bilayer center is expressed as,

$${\varphi }_Q^{\mathrm{sc}}(z=d)={\displaystyle \frac{e^2{G}_Q^{\left(22\right)}(d,d)}{1-e^2{G}_Q^{\left(22\right)}(d,d){\Pi }_{2D}(Q,\omega =0)}},$$

(14)

where $\Pi 2D(Q,\omega =0)$ represents the static electronic polarizability of the 2D system, originally formulated by Stern [41] at zero temperature and later generalized to finite temperatures by Maldague [42], and $G^{\left(22\right)}Q(d,d)$ denotes the Fourier transform of the Poisson Green’s function evaluated with the source charge located at the center of the 2D layer [40]:

$${\Pi }_{2D}(Q,\omega ;T,E_{\mathrm{F}})={\int }_0^{\infty }\mathrm{d}\mu {\displaystyle \frac{{\Pi }_{2D}(Q,\omega ;0,\mu )}{4k_{\mathrm{B}}Tcosh\left(\frac{E_{\mathrm{F}}-\mu }{2k_{\mathrm{B}}T}\right)}}$$

(15)

where T the temperature, ${\Pi }_{2D}(Q,\omega )$ is the zero temperature polarizability, $E_{\mathrm{F}}$ the Fermi level and $k_{\mathrm{B}}$ is the Boltzmann constant. Lastly, the impurity scattering rate corresponding to an impurity concentration $N_I$ can be expressed as

$${\displaystyle \frac{1}{{\tau }^{\mathrm{imp}}(\mathbf{k},n)}}={\displaystyle \frac{2\pi }{\hslash }}{N}_I\sum _{{\mathbf{k}}^{\prime },n^{\prime }}|\langle {\mathbf{k}}^{\prime },n^{\prime }|{\varphi }_Q^{\mathrm{sc}}(\mathbf{k}-{\mathbf{k}}^{\prime }){|\mathbf{k},n\rangle |}^2\delta (E_{{\mathbf{k}}^{\prime }}-E_{\mathbf{k}}).$$

(16)

The impact of impurity scattering is strongly modulated by the choice of gate dielectric. High-$\kappa$ top oxides, such as HfO₂, enhance dielectric screening because their large dielectric constant reduces the electric field generated by impurities. This weaker Coulomb field lowers the scattering strength experienced by carriers, thereby improving impurity-limited mobility. However, as disccused in the previous section, a high-$\kappa$ dielectric environment also strengthens interface plasmon–phonon (IPP) scattering in the channel, which can significantly degrade carrier mobility. In contrast, low-$\kappa$ dielectrics like hBN do not amplify IPP scattering as strongly, and in such cases impurity scattering may be the dominant mechanism depending on the carrier density. This interplay between IPP and impurity scattering has been studied studied previously [43]. Figure 8(a), adapted from [43], shows that with a HfO₂ top-gate insulator, IPP scattering dominates transport across the entire range of impurity concentrations, whereas in the case of a low-$\kappa$ hBN top insulator, its influence becomes significant only when the impurity density is below the mid-10¹¹ cm⁻² level.

5. Electronic Transport in Bilayer Transition Metal Dichalcogenides

This section presents the carrier transport properties of bilayer transition metal dichalcogenides (TMDs). Table 2 and Table 3 report the calculated hole and electron mobilities at a carrier density of $5\times {10}^{12}$ cm⁻² and $T=300$ K for freestanding bilayers and for double-gated structures with hBN or HfO₂ as the top dielectric. Previous theoretical reports [44] have predicted a high hole mobility in WSe₂ and WS₂. Similary, in our study we found particularly high hole mobilities for bilayer WS₂ and WSe₂ underscoring their promise for high-performance devices. At the same time, the introduction of a high-$\kappa$ dielectric such as HfO₂ significantly reduces the mobility due to enhanced IPP scattering, whereas low-$\kappa$ dielectrics like hBN preserve much of the intrinsic transport advantage. A more detailed discussion of these trends is provided below.

Freestanding bilayers exhibit the highest mobilities because carriers scatter only with intrinsic acoustic and optical phonons. For example, we calculate hole and electron mobilities of 2300 cm²/Vs and 161 cm²/Vs for bilayer WS₂, and 1300 cm²/Vs and 201 cm²/Vs for bilayer WSe₂, respectively.

When hBN is introduced as the top-gate dielectric, mobilities remain close to or in some cases even exceed those of the freestanding case. For instance, the hole mobility in bilayer WSe₂ increases from 1300 cm²/Vs to 1500 cm²/Vs, while the electron mobility in WS₂ improves slightly from 161 cm²/Vs to 172 cm²/Vs. This beneficial behavior stems from two factors. First, hBN exhibits inherently weak interface plasmon phonon (IPP) scattering due to its low ionic polarization and the high optical phonon frequencies associated with the light B and N ions. Second, hBN provides additional dielectric screening of Coulomb interactions, which further suppresses scattering from charged impurities and enhances transport characteristics. The combined effect of weak IPP coupling and efficient dielectric screening makes hBN an exceptionally favorable dielectric environment for TMD devices [38].

In contrast, HfO₂ leads to strong degradation in mobility across all TMDs studied. For WS₂, the hole mobility drops from 2300 cm²/Vs in the freestanding configuration to only 100 cm²/Vs, and for WSe₂, the electron mobility falls from 201 cm²/Vs to 36 cm²/Vs. This sharp reduction is due to the high dielectric constant and strong ionic polarization of HfO₂, which significantly enhance remote phonon and IPP scattering. Although the high-$\kappa$ oxide provides strong electrostatic control the additional scattering channels dominate, severely limiting mobility. These findings are consistent with earlier reports [37,45] for TMDs, showing that mobility decreases nearly monotonically as the gate insulator’s dielectric constant increases. While hBN offers the best transport performance and can even outperform the freestanding case in some configurations, it comes with trade-offs such as a low dielectric constant that limits gate capacitance and a relatively small barrier height that could increase gate-leakage currents if used alone as a gate dielectric.

6. Velocity-Field Characteristics

Next, we studied high-field carrier transport for freestanding and double-gated bilayer TMDs (Figure 9 and Figure 10). The drift velocity initially increases linearly with field (Ohmic regime), saturates, and then exhibits negative differential mobility (NDM) at high fields due to carriers populating satellite valleys with heavier effective masses.

For freestanding WS₂, NDM occurs near ${10}^5$ V/cm as holes gain energy to populate the $\Gamma$ valley with larger effective mass, reducing drift velocity. In HfO₂/bilayer WS₂/SiO₂, increased scattering suppresses energy gain, shifting NDM onset to fields above $5\times {10}^6$ V/cm. The hBN configuration behaves similarly to freestanding but with slightly delayed NDM.

Bilayer WSe₂ exhibits a similar behavior with NDM onset around ${10}^4$ V/cm in freestanding conditions but delayed under HfO₂. Bilayer MoS₂ has lower mobility and delayed NDM onset due to stronger intrinsic scattering and dielectric effects, with NDM absent in the HfO₂ environment over the field range studied.

In summary, bilayer TMDs exhibit transport strongly modulated by dielectric environment. Freestanding layers achieve the highest mobilities, limited only by intrinsic phonon scattering. Low-$\kappa$ dielectrics like hBN moderately reduce mobility by remote phonon scattering, whereas high-$\kappa$ materials such as HfO₂ drastically degrade mobility due to enhanced interface plasmon-phonon scattering. High-field transport shows NDM related to satellite valley population, which is suppressed by increased scattering in high-$\kappa$ stacks. These findings emphasize the critical trade-offs in device design between gate control and mobility retention in bilayer TMD-based transistors.

From the mobility results, it is evident that bilayer WSe₂ and WS₂ are strong candidates for future FET technologies. Compared to their monolayer counterparts, the bilayers generally exhibit higher hole mobilities (e.g., bilayer WS₂: 2300 cm²/Vs vs. monolayer WS₂: 750 cm²/Vs) but similar electron mobilities.

7. Bilayer TMD-Based Double-Gate MOSFETs

In this section, we analyze the performance of double-gate MOSFET architectures employing TMD bilayers as channels, focusing on how design choices affect electrostatic control, scattering, and series resistance. As device dimensions scale aggressively into the sub-10 nm regime, maintaining strong gate control is critical. To address this, the semiconductor industry has adopted high-$\kappa$ dielectrics such as HfO₂, which enable reduced equivalent oxide thickness (EOT) while mitigating gate leakage. While hBN provides low scattering and yields high mobilities, its relatively low dielectric constant makes it less attractive for scaled CMOS technologies. Motivated by this tradeoff, we concentrate on evaluating device performance with HfO₂ as the top-gate dielectric, assessing its impact on transport and electrostatics in bilayer TMD channels.

Our first goal is to identify the most promising device architecture. We begin with a 30 nm long double-gate nFET based on bilayer WSe₂ and HfO₂ as the top oxide. In our discussions, we will first study the electrostatics of this baseline device using phonon-limited transport. For the final device design, we include the complete scattering set (phonon, interface plasmon phonon, impurity) and evaluate multiple configurations involving HfO₂ and hBN as top oxides, and bilayer WSe₂ and WS₂ as channels for both nFET and pFET.

The first design is a double-gate n-type MOSFET with bilayer WSe₂ as the channel and HfO₂ as the top gate oxide, shown in Figure 1. The device has a channel length of $L_{\mathrm{ch}}=10$ nm and source/drain extensions of $L_{\mathrm{S}/\mathrm{D}}=10$ nm. The channel is intrinsic, and the source/drain regions are degenerately doped. A self-aligned top gate of $L_{\mathrm{g}}=10$ nm is centered over the channel. The EOT is given by:

$${ϵ}_{{\mathrm{SiO}}_2}·t_{{\mathrm{SiO}}_2}={ϵ}_{\mathrm{ox}}·t_{\mathrm{ox}}$$

Assuming a 0.7 nm EOT and using HfO₂ (${ϵ}_r\approx 25$), the physical oxide thickness is $t_{{\mathrm{HfO}}_2}\approx 4.2$ nm. Figure 11 shows the transfer characteristics at $V_{\mathrm{DS}}=0.2$ V and $0.5$ V. The device achieves an on-current $I_{\mathrm{on}}\approx 1000$ A/m at $V_{\mathrm{DS}}=0.5$ V, a threshold voltage $V_{\mathrm{T}}\approx -0.05$ V, and a subthreshold slope (SS) of ≈ 89 mV/dec.

In scaled FETs with high-$\kappa$ dielectrics, electrostatics can be affected by strong fringing fields at the gate edges. These are electric field lines that do not pass vertically through the channel but instead fringe laterally into the source and drain regions. Fringing fields degrade short-channel performance by distorting the potential profile and weakening gate control. TMDs like bilayer WSe₂ have high out-of-plane permittivity (${\epsilon }_{\perp }=12$–15) and low in-plane permittivity (${\epsilon }_{\Vert }=5$–7), which amplifies fringing field penetration through the high-$\kappa$ gate dielectric. As noted in [46,47], this results in stronger short-channel degradation especially increased DIBL and worse SS.

Figure 12 (a) shows the potential energy profile along the channel. The steep potential profile near the drain end at higher $V_{\mathrm{DS}}$, corresponds to electric fields that exceed the Ohmic region. In Figure 12 (b), the average electron kinetic energy increases near the drain for higher biases, exceeding 250 meV at $V_{\mathrm{DS}}=0.4$ V. This corresponds to substantial carrier heating in high-field regions. As dicussed in Section 4.3, IPP scattering is very strong in high-k dielectric environment. But, we need high-$\kappa$ to maintain the gate control. A solution to mitigate this issue is using low-$\kappa$ spcaer dielectric in the source and drain region but keeping the high-$\kappa$ HfO₂ above the channel region. So, we introduced low-$\kappa$ SiO₂ spacers in the source/drain extensions and extended the drain length to 30 nm Figure 13 to also allow the carriers to cool before reaching the contact.

Figure 14 shows the transfer characteristics of the large-drain device with spacers and only phonon scattering. We observe marginal improvement in subthreshold swing (SS ≈ 84 mV/dec). However, the extended drain region introduces parasitic resistance that reduces both the ON-current ($I_{\mathrm{on}}$) and transconductance.

While these results account only for intrinsic phonon scattering, additional scattering mechanisms such as interface plasmon phonons are expected to further degrade $I_{\mathrm{on}}$ by increasing the total series resistance. Figure 15 shows the spacers comparison transfer characteristics of the large-drain device with phonon and IPP scattering included. As expected we see improvement in $I_{\mathrm{on}}$ with the low-k spacers. At low bias we see higher current for the device without spacers owing to the fringe-induced barrier lowering (FIBL) as discussed in [47].

High $I_{\mathrm{on}}$ remains a critical performance metric for logic transistors, as it directly influences switching speed, drive strength, and energy efficiency. A higher $I_{\mathrm{on}}$ enables faster charging and discharging of load capacitances, reducing gate delay and improving overall circuit throughput. Maintaining high drive current is critical for ensuring fast switching and high performance in scaled logic devices, while also improving power efficiency by reducing switching delays [48,49,50].

To improve $I_{\mathrm{on}}$, reducing source/drain extension length is necessary to suppress parasitic resistance, as shown in Figure 16, which illustrates the final device design with shorter source/drain extensions. This design choice, however, comes at the cost of increased hot carrier generation. Hot-carrier degradation in aggressively scaled devices is a well-characterized reliability challenge. However, comprehensive modeling of HCD under various bias conditions enables designers to predict degradation and implement mitigations early in the design process [51,52].

To quantify the impact of various scattering mechanisms on device performance, we simulate pFETs with increasing scattering complexity: starting with only bulk phonons, then adding interface plasmon phonon (IPP) scattering, and finally including impurity scattering. We use pFETs for this analysis because bilayer WSe2 exhibits higher hole mobility, leading to higher baseline $I\mathrm{on}$ values—making the impact of added scattering mechanisms more evident. As shown in Table 4, the long-channel (30 nm drain) pFET suffers a sharp drop in $I_{\mathrm{on}}$ from 950 A/m to 283 A/m as scattering increases. In contrast, a short-channel (5 nm drain) pFET with all scattering mechanisms still delivers a high $I_{\mathrm{on}}$ of 820 A/m, demonstrating that aggressive scaling can recover performance even under significant scattering.

Finally, we present all the results for a ultra short device with 5nm source and drain extension. Figure 17 and Figure 18 show the I–V characteristics of bilayer WSe₂ nFETs and pFETs, comparing devices with hBN and HfO₂ as top gate dielectrics and SiO₂ spacers.

Table 5 and Table 6 summarize the ON-current, transconductance, and SS for short-drain nFET and pFET designs based on bilayer WSe₂ and WS₂. Using HfO₂ as the top gate insulator and including all scattering mechanisms (bulk phonons, IPP, and impurity), the best performance was achieved with the SiO₂–HfO₂–SiO₂/bilayer-WSe₂/SiO₂ pFET, which reached $I_{\mathrm{on}}=820$ A/m at $V_{\mathrm{GS}}=V_{\mathrm{DS}}=0.4$ V. Replacing HfO₂ with hBN as the top oxide further improved $I_{\mathrm{on}}$ by $\sim 10\%$ to 890 A/m as shown in Figure 19.

8. Conclusions

We have presented a comprehensive investigation of carrier transport and device performance in bilayer transition metal dichalcogenides (TMDs), using full-band Monte Carlo simulations. Our study first examined the low-field mobility and high-field velocity-field characteristics for several bilayer TMDs—including WS₂, WSe₂, and MoS₂—in freestanding and double-gated configurations with hBN and HfO₂ dielectrics.

The results clearly show that the dielectric environment plays a pivotal role in modulating transport. Freestanding bilayers exhibit the highest mobilities, limited only by intrinsic phonon scattering. hBN, due to its weak ionic polarization and high phonon energy, retains high mobility and in some cases even improves upon the freestanding case by screening out-of-plane fields while minimally enhancing interface plasmon-phonon (IPP) scattering. In contrast, high-$\kappa$ materials like HfO₂ drastically degrade mobility due to strong IPP and remote phonon scattering, despite offering better electrostatic control. These trends were observed consistently across all bilayer TMDs studied, reinforcing the importance of dielectric selection in device design. Our FET simulations indicate that scaling down the source/drain extensions significantly improves the ON current ($I_{\mathrm{on}}$) by reducing series resistance, though care must be taken to manage potential reliability issues such as hot-carrier degradation. High-field transport characteristics reveal strong velocity saturation and negative differential velocity (NDV), especially in freestanding and hBN configurations. These effects are suppressed in high-$\kappa$ environments due to increased scattering.

By considering HfO₂ as the top gate dielectric and accounting for bulk phonon, interface plasmon phonon (IPP), and impurity scattering mechanisms, the optimal device configuration was identified as the SiO₂-HfO₂-SiO₂/bilayer-WSe₂/SiO₂ pFET with 5 nm source/drain extensions. This setup achieves an ON current ($I_{\mathrm{ON}}$) of 820 A/m at $V_{\mathrm{GS}}=V_{\mathrm{DS}}=0.4$ V, meeting ITRS performance benchmarks. Replacing the top gate dielectric with hBN leads to a modest increase in $I_{\mathrm{ON}}$ to approximately 890 A/m under the same bias conditions, marking a 10% improvement. These findings highlight the critical role of channel-dielectric engineering and device architecture optimization in enhancing the performance of bilayer TMD-based transistors.

In summary, bilayer WSe₂ and WS₂ could, in principle, enable sub-10 nm FETs that meet the latest ITRS high-performance targets, offering favorable transport and electrostatic characteristics. However, realizing high-performance devices requires careful trade-offs between electrostatic control, mobility degradation, and reliability, as well as attention to fabrication challenges such as high-quality channel growth, gate insulator deposition, and low-resistance contacts. Our results underscore the critical role of accurate physical modeling, including scattering mechanisms and dielectric effects, in guiding the optimization of bilayer TMD-based transistors.