Electrical-equivalent van der Waals gap for 2D bilayers
Dhirendra Vaidya Saurabh Lodha Swaroop Ganguly


























[]dhirendra22121987@gmail.com



Department of Electrical Engineering, IIT Bombay, Mumbai 400076, India

Department of Electrical Engineering, IIT Bombay, Mumbai 400076, India

[]swaroop.ganguly@gmail.com
Department of Electrical Engineering, IIT Bombay, Mumbai 400076, India











 Vertical stacks of two-dimensional (2D) materials, separated by the van der Waals gap and held together by the van der Waals forces, are immensely promising for a plethora of nanotechnological applications. Charge control in these stacks may be modeled using either a simple electrostatics approach or a detailed atomistic one. In this paper, we compare these approaches for a gated 2D transition metal dichalcogenide bilayer and show that recently reported electrostatics-based models of this system give large errors in band energy compared to atomistic (Density Functional Theory) simulations. These errors are due to the tails of the ionic potentials that reduce the electrical-equivalent van der Waals gap between the 2D layers, and can be corrected by using the reduced gap in the electrostatic model. For a physical van der Waals gap (defined as the chalcogen to chalcogen distance) of 3 in a 2D bilayer, the electrical-equivalent gap is less than 1 Å. For the example of band-to-band tunneling based ultra low-power transistors, this is seen to lead to errors of several hundred millivolts and more in the threshold voltage estimated from electrostatics.





Electrical-equivalent van der Waals gap for 2D bilayers
    Swaroop Ganguly
    December 30, 2023
=======================================================



Recently there has been burgeoning interest in 2D vertical stacks for a wide variety of electronic applications. When a 2D monolayer (ML) of transition metal dichalcogenides (TMD) or graphene is placed on top of another 2D ML, such a bilayer (BL) is held together by the weak van der Waals force between the two MLs. These homogeneous or heterogeneous BLs or multilayers may be called van der Waals (vdW) stacks; they have been imagined as the building blocks of many novel electronic <cit.>, optoelectronic <cit.>, spintronic <cit.> and power devices <cit.>. 

Possibly the most significant use of these 2D BLs would be in next-generation tunnel-FETs (TFETs), to fulfil the extreme energy-efficiency requirement for internet of things (IoT) applications projected by the International Roadmap for Devices and Systems  <cit.>. The workhorse of conventional electronics, viz. the CMOS switch, works by controlling the height of a thermal barrier; and therefore it needs at least (kT/q)ln10=60 mV to  change the current by an order of magnitude; in other words, it has a subthreshold swing (SS) greater than 60 mV/decade. TFETs, on the other hand, work by controlling the thickness of a tunneling barrier; so they can, in principle, have an SS < 60 mV/decade<cit.> and are therefore of immense interest for extreme low-energy IoT applications. TFETs have been fabricated and studied using group IV <cit.> as well as III-V <cit.> semiconductors. However, because of the usually high density of interface states <cit.> - leading to trap-assisted tunneling, achieving an SS < 60 mV/decade in these systems has proven to be difficult. This issue might be alleviated by using a vertically stacked 2D BL and controlling the tunnel current between the top and bottom MLs through gate voltages, leading to super-steep sub-threshold characteristics. Even though realizing a vertical 2D stack <cit.> of two MLs is a challenging task, several novel device ideas have been explored in some depth in the literature <cit.>. 

As in most electronic devices, electrostatics plays the leading role in the functionality of 2D devices, including BL devices. In prior work, e.g. <cit.>, the band alignments of the vertical 2D stacks have been modeled using simple electrostatics -  i.e. assuming charge sheets corresponding to the MLs, separated by the their physical distance. This ignores the spatial distribution of the ionic potentials, in particular their tails falling off from the atomic sites, and thereby underestimates (overestimates) the effective electrical thickness (separation) of the MLs. We will show that this can lead to significant errors in modeling the gate control when the two layers are separated by a distance comparable to the vdW gap (defined as chalcogen to chalcogen distance) viz. ≈ 3-3.5 <cit.>.



The aforementioned effect may be illustrated by considering the cross-section of a vdW BL and the ionic potentials therein. Figure <ref>a shows the specific example of a homojunction WS2 BL TFET. A vdW BL can in general be realized with vertical stacking of any suitable combination of 2D MLs. We choose a homojunction BL here for simplicity of illustration and intuitive explanation. Figure <ref>b shows the laterally averaged electrostatic potential energy of WS2 BL system in equilibrium. We can clearly see that the electrical thickness of each ML is larger than the physical thickness (defined by the lattice sites of the boundary chalcogen atoms) on account of the tails of the ionic potentials. A comparison of atomistic simulations using density functional theory (DFT) versus electrostatic simulations reveals that large errors are possible in the regular electrostatic simulations which ignore the effect of the ionic potential tails.



We will begin with a discussion of the intuitive picture of vertical 2D MLs under an electric field. There we also explain the band alignment in these vertical 2D stacks, and present the popular purely electrostatic approach that has been used in modeling these vdW stacks. Following that, we discuss the tails of the ionic potentials and the necessity of including these for predictive modeling. Thereafter, we compare DFT and electrostatic simulations of the vdW stacks under an electric field, and propose a simple correction to the interlayer separation used in the electrostatic calculation in order to model the effect of the ionic potential tails. Lastly, the importance of the potential tail correction on device design is illustrated by comparing the threshold voltages for band to band tunneling (BTBT) obtained from the simple (uncorrected) electrostatic model versus our corrected one.

 Figure <ref>a shows a four terminal electronic device. Such a device can be used as a p-n junction diode<cit.>, a tunnel diode<cit.> or a TFET<cit.>. This flexibility in application is possible because the top and bottom gate can modulate the band offsets of the top and bottom MLs.Figure <ref>a shows a gated TMD bilayer system. CB and VB offsets in two MLs of the system can be modified by applying a gate bias. Figure <ref>b shows the band alignment of the 2D MLs when no bias is applied. Under this condition, the band offsets are obviously zero as may be expected. Figure <ref>c shows the band alignments when a positive voltage is applied to the top gate and a negative voltage is applied to the bottom gate. The positive potential on the top gate makes the top ML n-type and the conduction band minimum (CBM) is near the Fermi level; whereas the negative potential on the bottom gate makes the bottom ML p-type and the valence band maximum (VBM) moves closer to the Fermi level. The condition reverses on reversing the bias polarities (Figure <ref>d). It is this gate control of band offsets in the vdW BL systems that opens up many interesting device application possibilities.

Since the scope of this paper is gate-controlled band alignment, we ignore the transport completely in what follows, and drop the source and drain terminals from our discussion (Figure <ref>a).
The charge density (ρ_2D,T and ρ_2D,B) in the 2D MLs (TMD_Top and TMD_Bottom) can be related to the electric fields (ℰ_Tox, ℰ_IL and ℰ_Box) using Gauss' law.


    D_IL-D_Tox   = ρ_2D,T
    
    D_Box-D_IL   = ρ_2D,B
    ϵ_ILℰ_IL - ϵ_Toxℰ_Tox   = ρ_2D,T
    ϵ_Boxℰ_Box - ϵ_ILℰ_IL   = ρ_2D,B

and


    ℰ_Tox   = (χ_T+qϕ_n,T)-(qϕ_M,T-qV_TG)/t_Tox
    ℰ_IL   = (χ_B+qϕ_n,B)-(χ_T+qϕ_n,T)/t_IL
    ℰ_Box   = (qϕ_M,B-qV_BG)-(χ_B+qϕ_n,B)/t_Box

where, χ_T, χ_B are the electron affinities of top and bottom MLs. t_Tox,  t_IL and  t_Box are the physical thicknesses of the top, bottom and inter layer dielectrics. ϕ_M,T and ϕ_M,B are the metal workfunctions of the top and bottom gates and qϕ_n,T = E_c,T-E_F  and  qϕ_n,B = E_c,B-E_F. The electron and hole densities in each monolayer are then related to ϕ_n,T and ϕ_n,B as,


    n_T    = g_valley,Tm_c,Tk_BT/πħ^2ln[exp(-qϕ_n,T/k_BT)+1] 
    
    n_B    = g_valley,Bm_c,Bk_BT/πħ^2ln[exp(-qϕ_n,B/k_BT)+1] 
    
    p_T    = g_valley,Tm_v,Tk_BT/πħ^2ln[exp(-qϕ_p,T/k_BT)+1] 
    
    p_B    = g_valley,Bm_v,Bk_BT/πħ^2ln[exp(-qϕ_p,B/k_BT)+1]

where, qϕ_p,T=E_g,T-qϕ_n,T and qϕ_p,B=E_g,B-qϕ_n,B. E_g,T and  E_g,B are the bandgaps of top and bottom 2D MLs respectively. The total charge density in each ML is given by,


    ρ_2D,T   = q(p_T - n_T + N_D,T - N_A,T) 
    ρ_2D,B   = q(p_B - n_B + N_D,B - N_A,B)

We have solved equations (3), (4), (12) and (13) self-consistently, using a Newton-Raphson method.


 Figure <ref>b shows the laterally averaged electrostatic potential of a WS2 homojunction BL system when the interlayer separation is 4 Å. It can be seen that the ionic potentials extend a few angstroms beyond the lattice sites of the chalcogens which are at the edges of the MLs. The extent of these tails into the gap between the MLs is comparable to the vdW gap (≈ 3-3.5 Å). 

Assuming 2D charge sheets corresponding to the TMD MLs, separated by the physical distance between them, would seem to be a good approximation when the interlayer separation is large compared to the extent of the ionic potential tail. It is obviously suspect when these become comparable. However, this latter case is of great practical interest since the interlayer separation in vertical stacks of 2D MLs is expected to be the vdW gap. Figure <ref> shows the laterally averaged electrostatic potential energy, obtained from DFT simulations, for BL homojunctions of MoS2, MoSe2, MoTe2, WS2, WSe2 and WTe2 as a function of position, starting from the lattice site of the outermost chalcogens. Since the chalcogen atoms are on the edges of the TMD MLs, it is the tails of the chalcogen ion potentials that extend a few angstroms into the space between the two MLs. Now, in order to quantify the extent of the tail, we observe that the potential energy for MoS2 and WS2 decays to 10% of the value that it has at the boundary chalcogen lattice site at a distance of 1.05 Å. Thereafter we consider this 10% value for the disulfides (viz. -1.93 eV) as the reference. We then observe that the potential energy for both MoSe2 and WSe2 reach that value at a distance of 1.12 from the boundary chalcogen atom; while for both MoTe2 and WTe2, this distance is 1.25 Å.



Before proceeding further, we discuss another possible explanation for the error in the pure electrostatic model that we considered: namely, the spread in the electron/hole wavefuctions. We modeled each 2D MLs as a delta-function potential well. Then the wavefunction of electrons and holes are of the form,


    ψ(z)    = √(κ)e^κ (z-z_0) for  z<z_0
        = √(κ)e^-κ (z-z_0) for  z>z_0

where z_0 is the position of the ML. Now, the spatial distribution of electrons/holes can be captured by multiplying equations (8-11) by |ψ(z)|^2 from equations (14-15). This would also suggest a finite thickness of the MLs. However, using a 1D finite difference Poisson solver<cit.> we found that the band alignments in vdW stacks are insensitive to the wave nature of electrons/holes for reasonable values of κ. Therefore we concluded that it is not the spread in the wavefunction that causes the error in the electrostatic model.



§ RESULTS AND DISCUSSION


 A fair comparison between atomistic DFT and phenomenological electrostatic models requires, firstly, a consistent set of TMD ML parameters to be used therein. Thus, for the electrostatic simulations, we use parameters that are extracted from the DFT simulations; these are shown in Table <ref>. The extracted values of the lattice constant, bandgap and effective masses are in reasonable agreement with previous abinitio studies and experiments<cit.>.





 



As depicted in figure <ref>, the applied gate electric field modulates the band alignments in the two MLs. The band offsets can be inferred from the conduction band edges at K-point in the Brillouin zone. The bandstructure of the overall BL system such as the one shown in inset of figure <ref> is just the superposition of the ML bandstructures, which are degenerate at zero bias. This degeneracy is broken by an applied gate electric field. Then, one ML becomes n-type while the other ML becomes p-type (figure <ref>). In the bandstructure of the BL system, the CBM can be seen shifted towards (away from) the Fermi level for the n-type (p-type) ML. The difference between the CBM of the two MLs gives the band offset (Δ E_C). At higher biases, the CBM of the p-type ML moves far away from the Fermi level and mixes with the higher order bands of the other ML. This makes the direct measurement of the CBM of the p-type ML difficult but measuring the VBM of this ML at the K-point is straightforward. 

Figure <ref> shows the bandstructures of BL systems of WS2 having 7 (a and b) and 4 (c and d) interlayer separations. Figure <ref>a is a bandstructure of a WS2 BL system with interlayer separation of 7 under the zero-bias condition. Since the interlayer separation is 7 the two MLs in the BL system are practically isolated from each other. The bandstructure of the entire BL system is a doubly degenerate ML bandstructure. When the BL system is subjected to an electric field, the degeneracy breaks; this is illustrated in figure <ref>b. Figure <ref>b shows the bandstructure of the overall BL system for V_TG = +V and V_BG = -V. Due to the applied gate field in this case, the top ML becomes n-type and the CBM of the top ML moves towards the Fermi level. The VBM of the bottom ML moves towards the Fermi level indicating that the bottom ML has become p-type. It can be seen in figure <ref>b that the bandgap of the overall BL system at K-point is smaller than the ML bandgap. This is not true bandgap modulation with electric field as suggested by <cit.> but rather an effect of the potential drop across the BL system, with each ML retaining its ML bandgap value (as shown in the Supplementary Material). Figure <ref>c (zero bias) and <ref>d (V_TG = +V and V_BG = -V) show the bandstructure of the overall BL system of WS2 when the interlayer separation is 4 Å. It can be seen that this is not exactly the double degenerate ML bandstructure at zeros bias. CB near Q-point (which is halfway between the Γ-point and K-point) and VB near Γ-point are not doubly degenerate. This is because of significant interlayer interaction at the smaller separation, which modifies the Hamiltonian of each layer. 

The interlayer interaction <cit.> of course increases as the separation is reduced. This necessitates modification of the BL Hamiltonian to incorporate off-diagonal interaction terms between the MLs, which would in turn modulate the solutions to the Schrödinger equation, viz. the bandstructure. However, the methodology for solving the Poisson equation to capture electrostatic effects due to space charge would remain the same. To elucidate this methodology, while excluding the complication due to the interacting BL Hamiltonian, we constrain our discussion to interlayer separation of 4Å and more. However, we emphasize that the electrostatics treatment presented here would continue to hold for smaller separation.   

Figure <ref> shows the comparison between electrostatic and DFT simulations for a bilayer system of WS2. The comparison for other TMDs is presented in supplementary material. The top gate is biased from 0 V to 4 V and the opposite polarity bias is applied to the bottom gate. This makes the top ML n-type and the bottom ML p-type.  The DFT simulations are subject to the constraint of overall charge neutrality; because the total number of electrons are fixed by the number of atoms. For a reasonable comparison, the electrostatic simulations must also ensure overall charge neutrality. Applying opposite polarity bias to the top and bottom gates imposes overall charge neutrality in the BL system (ρ_2D,T=-ρ_2D,B). We compare DFT and electrostatic simulations for interlayer separation of 4, 5, 6 and 7 Å. It can be seen that the electrostatic simulations and the DFT simulations are not in quantitative agreement. At a given top gate bias (+V) the electrostatic model suggests that the Fermi level is closer to the CB (VB) edge of the top (bottom) ML compared to the DFT simulations. The error shown in figure <ref> shows that the electrostatic and atomistic simulations agree (disagree) for large (small) interlayer separation. Figure <ref> shows that the error is maximum at the bias point where each ML is close to becoming degenerate n-type or p-type. Comparing the nearly intrinsic bias region (V_TG = 0 to 0.8 V) to the extreme degenerate region (V_TG > 3.0 V), we see that the error is smaller in the former. The additional error in the latter bias region is attributed to the fact that the Q-valley (Γ-valley) in CB (VB) starts to populate with the electrons (holes) as the Fermi level is moved closer to the CB (VB) of the top (bottom) monolayer. 



Figure figure <ref>c (<ref>d) shows the maximum error in CBM of the top ML (VBM of bottom ML) divided by the ML bandgap. The maximum error in CBM and VBM is presented for all TMD BL systems considered in this study with interlayer separation of 4 and 5 Å. The bandgap normalized errors are obviously large enough (1/4^th of the bandgap) to affect the device design quite substantially.



Since the extent of the ionic potential tails beyond the chalcogen lattice sites is comparable to the vdW gap, the electrical thickness of the thin MLs is a few angstroms more than the physical thickness determined by the lattice atom positions. Now it is verified using a finite difference 1D Poisson solver<cit.> that the thickness of the MLs per se do not impact the electrostatic simulations. However,  their increased electrical thickness reduces the effective (electrical) interlayer separation. This introduces significant error in the estimation of the CBM (VBM) of the top (bottom) ML in both the non-degenerate and degenerate regions. To account for the reduction in the effective separation and achieve agreement between the electrostatic and approaches, we introduce a simple empirical correction to the interlayer separation. The starting guess for the correction, in the illustrative case of the homogeneous BL, is twice the tail length estimated in the figure <ref>; it is then fine-tuned to match the electrostatic simulation results with that of the . Figure <ref> presents the values of the final tail length correction for MoS2, MoSe2, MoTe2, WS2, WSe2 and WTe2. Clearly, its primary dependence is on the chalcogen element in the TMD; more specifically, the atomic number thereof. This is as expected, given that potential tail in question belongs to the chalcogen atom. It may be observed from figure <ref> that the tail length correction leads to a near-perfect agreement between the electrostatic model and DFT simulations in the non-degenerate region. In the degenerate region, there is some discrepancy in the CBM of the top ML and VBM of the bottom ML even after applying the correction. This is more apparent from the equivalent of figure <ref> for the other TMDs, shown in the supplementary material. The reason behind it is the nonzero occupancy of the Q-valley and Γ-valley which are not included at all in the electrostatic simulations. The error for the largest bias in the CBM (VBM) of the top (bottom) ML depends inversely on the energy separation between the K-valley and the Q-valley (Γ-valley) as shown in table <ref>. The larger the separation, the smaller would be the occupation in the higher energy valleys, and hence the smaller would be the error. The occupation of higher valleys may be included within the electrostatic model by increasing the density of states effective mass; however the detailed treatment of this correction is outside the scope of this paper.

To illustrate the impact of the tails of the ionic potential for the specific example of a vdW-TFET performance, we apply the potential tail length correction to a previously studied vdW-TFET designed using a MoS2-WTe2 vertical heterojunction <cit.>. We use the same transistor design and simulation parameters as suggested by <cit.> and calculate the tunneling window and threshold voltage for BTBT.
The interlayer dielectric in the above design is intended to be two MLs (6 Å) of the 2D material boron nitride (BN), which has a dielectric constant of 4.2. Their electrostatic model assumes 6 of a homogeneous material with this dielectric constant. In reality, such a MoS2-BN-WTe2 vdW-TFET would be realized by stacking the MLs of MoS2, BN and WTe2, wherein each layer would be separated from its neighboring layer/s by the vdW gap. This picture suggests that while the effective dielectric constant for the BN may indeed be reasonably taken to be 4.2, the potential tail correction for the MoS2 and WTe2 should consider the vacuum values calculated above. The thickness of BN is then assumed to be the chalcogen-chalcogen separation minus the tail correction.
Figure <ref>a shows the CBM, VBM and Fermi levels for the MoS2-WTe2 vdW-TFET when drain voltage V_D=0.3  V and source is grounded, V_S=0  V. The gate field is varied by varying the top gate voltage while keeping the bottom gate at 0 V. It can be seen that at V_TG=0 V the band alignment is of type-II, which is illustrated in figure <ref>b. In this situation, the CBM (VBM) of the MoS2 (WTe2) aligns with the bandgap of WTe2 (MoS2); this is not favorable for elastic tunneling of electrons from the VB of WTe2 since no states are available in the MoS2. Now, the band alignment can be modified by applying a gate field. At about V_TG=0.27  V the band alignment changes to type-III, which is illustrated in figure <ref>c. It can be seen that the CBM of MoS2 is now below the VBM of WTe2 which enables BTBT. The applied bias on the top gate at which the band alignment changes from type-II to type-III may be defined as the threshold voltage V_T of the vdW-TFET. The V_T determined by the tail length corrected electrostatic model is 0.66 V. The pure electrostatic model underestimates the V_T by 380 mV, rendering it utterly useless.

The situation is spectacularly aggravated when we replace BN as the interlayer dielectric by vacuum of thickness equal to vdW gap. Note that this choice of interlayer dielectric has been proposed for vdW-TFETs in the literature <cit.> in order to maximize BTBT. We consider the same MoS2-WTe2 BL as above, but this time with the vdW gap separation between the MLs. Since the precise value of the vdW gap for this stack was not available, we consider it to be 3-3.5 for the purpose of V_T estimation. The V_T estimated by the pure electrostatic model shows no variation with the interlayer separation (note that the threshold voltage coming out to be nearly 0 V in this case is coincidental - it is a function of the gate metal workfunction chosen by <cit.>). The potential tail length corrected electrostatic model, on the other hand, exhibits strong dependence as shown in figure <ref>. The error in the uncorrected V_T is seen to be extremely large in this scenario. This follows in turn from the much-larger relative error in the interlayer separation (viz. the difference between the physical and electrical interlayer separation, divided by the former) for the very-small vdW gap separation: note that the electrical equivalent interlayer separation is just 1.1Å (0.6Å) when the physical chalcogen-chalcogen interlayer separation is 3.5Å (3Å).





§ CONCLUSION

Angstrom-level insights drawn from atomistic DFT simulations suggest that pure electrostatics-based design of 2D bilayer devices may be prone to large errors. This is because the tails of the ionic potentials protrude into the interlayer gap, greatly reducing the electrical-equivalent interlayer separation with respect to the physical separation. The errors in the band energies within each ML could be as large as a fourth of the bandgap. In the example of 2D bilayer TFETs, the threshold voltage estimated from pure electrostatics may be off by hundreds of millivolts to over a volt; the smaller the physical separation, the larger the error. It is shown that employing an electrical-equivalent interlayer separation - obtained by deducting the ionic potential tail length from the physical interlayer separation - suffices to correct the electrostatic model.



§ METHODS


 For the atomistic DFT simulations here, we have used the Atomistix Simulation Toolkit (ATK) from Quantumwise <cit.>. Figure <ref> shows the simulated atomic geometry of BL WS2 under an electric field. The BLs are seen to have an anti-symmetric arrangement, with each ML having a 2H geometry. We study the BL system by varying the interlayer separation between the two MLs, between 4-10 for MoS2, MoSe2, MoTe2, WS2, WSe2 and WTe2 BLs. The interlayer separation is defined as the distance between the lattice sites of boundary chalcogen atoms (i.e. bottom chalcogen of top ML, and top chalcogen of bottom ML). We use a 21x21x1 K-mesh for the DFT calculation, with a 110 Hartree density mesh cutoff. The convergence criterion used for the simulations is 10^-6 Hartree; for the electron occupancy, we use Fermi-Dirac statistics with broadening equivalent to 300 K. We use the Open MX basis sets<cit.> which are shipped along with the ATK package. We do not relax the bilayer system intentionally to focus on the interlayer separation and its impact on the band offsets. However the individual MLs are relaxed with the criteria of 0.01 eV/for force and 0.001 GPa for stress. The two relaxed MLs are then stacked to form a BL. The dielectric constant between the metallic contacts and the BL is set to give a 1 nm effective oxide thickness (EOT) i.e. 1 nm SiO2. On the other hand, the space between the two MLs is considered to be vacuum - with the dielectric constant of free space.

















