Electronic Transport and Thermopower in 2D and 3D Heterostructures--A Theory Perspective

TL;DR

This review summarizes recent theoretical and experimental advances in understanding electronic and thermoelectric transport in 2D and 3D heterostructures, highlighting progress in simulation techniques and device applications.

Contribution

It provides a comprehensive overview of recent progress in modeling and understanding transport phenomena in heterostructures, emphasizing new theoretical approaches and their implications for device development.

Findings

Progress in understanding electronic transport in 2D materials and heterojunctions.

Advances in quantum transport simulation techniques like NEGF and Wigner functions.

Potential for improved energy harvesting and quantum electronic devices.

Abstract

In this review, we discuss the impact of interfaces and heterojuctions on the electronic and thermoelectric transport properties of materials. We review recent progress in understanding electronic transport in two-dimensional (2D) materials ranging from graphene to transition metal dichalcogenides (TMDs), their homojunctions (grain boundaries), lateral heterojunctions (such as graphene/MoS$_2$ lateral interfaces), and vertical van der Waals (vdW) heterostructures. We also review work in thermoelectric properties of 2D heterojunctions, as well as their applications in creating devices such as resonant tunneling diodes (RTDs). Lastly, we turn our focus to work in three-dimensional (3D) heterostructures. While transport in 3D heterostructures has been researched for several decades, here we review recent progress in theory and simulation of quantum effects on transport via the Wigner and non-equilibrium Green's functions (NEGF) approaches. These simulation techniques have been successfully applied toward understanding the impact of heterojunctions on the thermoelectric properties, with applications in energy harvesting, and electron resonant tunneling, with applications in RTDs. We conclude that tremendous progress has been made in both simulation and experiments toward the goal of understanding transport in heterostructures and this progress will soon be parlayed into improved energy converters and quantum nanoelectronic devices.

Equations

$${\tau }_b(k_i)=\frac{\mid 4k_{i\perp }{k}_{t\perp }\mid }{\mid k_{i\perp }+k_{t\perp }{\mid }^2}$$

$${\Gamma }_b(E)=\frac{1/4{\pi }^2\int {\tau }_b(k_i)\delta (E-E_b(k_i))d{k}_i}{1/4{\pi }^2\int \delta (E-E_b(k_i))d{k}_i}$$

$$\Xi (E)=\sum _b{v}_b(E){\Gamma }_b(E)D_b(E)$$

$$G_{GB/int.}=\frac{e^2}{2}{\int }_{E_C}^{E_{max}}\Xi (E)\left(-\frac{\partial f(E-E_F,T)}{\partial E}\right)dE$$

$${S_{Mott}=-\frac{{\pi }^2{k}_B^{2}T}{3|e|\sigma }\frac{d\sigma (E)}{dE}|}_{E=E_F}.$$

$$\left(\frac{\partial }{\partial t}+v_r{\nabla }_r+\frac{eF}{\hslash }{\nabla }_k\right)f_w\left(r,k,t\right)=Q{f}_w\left(r,k,t\right)+{\left(\frac{\partial f_w}{\partial t}\right)}_{coll}$$

$$f_w\left(r,k,t\right)=\frac{1}{2\pi }\int d{r}^{\prime }{e}^{-i{r}^{\prime }k}\rho \left(r+\frac{r^{\prime }}{2},r-\frac{r^{\prime }}{2}\right)$$

$$Q{f}_w\left(r,k,t\right)=\int d{k}^{\prime }{V}_w\left(r,k-k^{\prime }\right)f_w\left(r^{\prime },k,t\right)$$

$$V_w\left(r,k\right)=\frac{1}{i\hslash {\left(2\pi \right)}^d}\int d{r}^{\prime }{e}^{-i{r}^{\prime }k}\left(V\left(r+\frac{r^{\prime }}{2}\right)-V\left(r-\frac{r^{\prime }}{2}\right)\right)$$

$$V_w\left(r,k\right)=\frac{2}{\pi \hslash }Im\{e^{2ikr}\hat{V}(2k)\}$$

$$\hat{V}(k)=\underset{-\infty }{\stackrel{\infty }{\int }}V(r)e^{-ikr}dr.$$

$$V(x)=V_{cl}(x)+V_{qm}(x).$$

$$\left(v_r{\nabla }_r+\frac{eF}{\hslash }{\nabla }_k\right)f_w\left(r,k,t\right)=Q{f}_w\left(r,k,t\right)+{\left(\frac{\partial f_w}{\partial t}\right)}_{coll}.$$

$$g_{i+1}(k)=\left(I(k)+\frac{eF}{\hslash }\frac{\partial f_0}{\partial k}-v(k)\frac{\partial f}{\partial r}\right)/S_0(k).$$

$$g_{i+1}(k)=\left(I(k)+\frac{eF}{\hslash }\frac{\partial f_0}{\partial k}-v(k)\frac{\partial f}{\partial r}+Q{f}_w\right)/S_0(k).$$

$$V(r)=\sum _{n=-\infty }^{\infty }{V}_p(r-n{L}_p).$$

$$Q{f}_w=\sum _{m=1}^{\infty }{W}_m(r)\left[f_w\left(r,k-\frac{m\pi }{L_p}\right)-f_w\left(r,k+\frac{m\pi }{L_p}\right)\right].$$

$$V_p(r)=\frac{V_0}{2}\left\{-erf[\beta (r-a)]+erf[\beta (r+a)]\right\}$$

$$W_m(r)=\frac{2V_0}{\pi \hslash m}{e}^{\frac{-m^2{\pi }^2}{{\beta }^2{L_p}^2}}\sin \left(\frac{2\pi ma}{L_p}\right)\sin \left(\frac{2\pi mr}{L_p}\right).$$

$$V_q(r)=V_0(1+\cos (K_{0}r))/2,\text{mbox}where\text{leavevmode}{K}_0=2\pi /L_p$$

$$Q{f}_w=W_m(r)\left[f_w\left(r,k-\frac{K_0}{2}\right)-f_w\left(r,k+\frac{K_0}{2}\right)\right]$$

$$\text{mbox}withquantumweight{W}_m(r)=\frac{A\sin (K_{0}r)}{\pi \hslash }.$$

$${\hat{G}}^R(E)=[E\hat{I}-\widehat{H_0}-\widehat{V_c}-{\hat{\Sigma }}^R(E){]}^{-1}$$

$${\hat{G}}^{<}(E)={\hat{G}}^R(E){\hat{\Sigma }}^{<}(E){\hat{G}}^A(E)$$

$$\hat{\Gamma }=i({\hat{\Sigma }}^R-{\hat{\Sigma }}^A)=i({\hat{\Sigma }}^{<}-{\hat{\Sigma }}^{>}),$$

$${\hat{\Sigma }}^R(t)=\Theta (t)({\hat{\Sigma }}^{<}-{\hat{\Sigma }}^{>}).$$

Full text

1 Introduction and overview

Progress in understanding charge transport in two-dimensional (2D) materials is forging along at a furious pace. It follows an arc paralleling that of three-dimensional (3D) materials a few decades prior: after their intrinsic (bulk) properties were well studied, researchers turned to enhancing the intrinsic characteristics by combining materials together into heterostructures. This approach brings both tremendous complexity and opens up vistas to emerging properties not available in individual materials or atomic monolayers. Such heterostructures have device applications ranging from CMOS transistors with few-nanometer channel lengths that continue to track the scaling predictions of Moore’s law to beyond-CMOS devices such as thermoelectric energy converters and resonant tunneling diodes based on quantum tunneling. The common thread across all the aforementioned materials and applications is that their functionality revolves around electron transport across interfaces between dissimilar materials, called heterojunctions. The resulting heterostructures alter the transmission of electrons and lead to increased contributions from tunneling, resonant quantum states, or energetic electrons with improved thermoelectric coefficients.

In order to cover the recent progress, while at the same time drawing parallels between 2D and 3D heterostructures, we organize the paper as follows: first we cover 2D materials, starting with a brief history before delving into graphene and MoS2 grain boundaries. Next, we review lateral heterostructures, such as those formed between graphene and MoS2 in Sec. 3, and include band alignment theory and resonant tunneling diode (RTD) applications. We round out the coverage of 2D materials with vertical (stacked) heterostructures in Sec. 4.1, including progress on applications in thermoelectric energy conversion. Then we switch over to 3D heterostructures and periodic superlattices in Sec. 5, covering the history of tunneling diodes, the theory and numerical simulation methods for quantum transport in heterostructures, and their applications in enhancing thermoelectric energy conversion efficiency of semiconductor materials. We end with conclusions and future outlook.

2 Charge transport in two-dimensional materials

Although the idea of graphene (a single layer of graphite) existed theoretically for several decades, it was believed that crystalline two-dimensional (2D) materials, including graphene, could not exist in nature because of thermodynamical instability [1]. It was in 2004, when a group of scientists from University of Manchester in United Kingdom and Institute for Microelectronics Technology in Russia successfully separated a single layer from bulk graphite for the first time [2]. Their groundbreaking work resulted in the 2010 Nobel Prize in Physics being awarded to Novoselov and Geim. They followed it up by preparing single-layer exfoliated forms of other van der Waals 2D atomic crystals, including transition metal dichalcogenides (TMDs) such as MoS2 and NbSe2, and boron nitride [3]. This cohort of 2D materials includes graphene, hexagonal boron nitride (hBN), and transition metal dichalcogenides (TMDs), which have covalent in-plane bonds and weak van der Waals (vdW) bonds across atomic planes. The weak interlayer coupling allows these vdW monolayers to be separated, manipulated, and assembled into vertical and lateral heterostructures. These materials have a host of unique and superlative electronic, optical, thermal, and thermoelectric properties, thus forming a basis for future low-power logic devices, high-efficiency energy materials, and high-performance optoelectronics.

The unique electronic structure of graphene, where the low energy charged carriers mimic relativistic particles with zero rest mass having an ’effective speed of light’ of 106 m s*-1* [4], along with experimental demonstrations of integer quantum Hall states [5] and Klein tunneling [6], not only opened up avenues for rich, fundamental science but also led to a possibility of realizing new electronic and magneto-electronic device applications. Earlier reports suggested that carrier mobilities, which were found to be about 15,000 cm2 V*-1* s*-1* for both electron and holes, were independent of temperature over the range between 10 and 1000 K [4], suggesting a very strong influence from the interaction with the underlying substrate. This has led to a series of studies trying to establish the upper bound of the intrinsic phonon-limited mobility by fabricating and measuring highly crystalline, suspended graphene [7, 8, 9].

An electron mobility of 230,000 $\text{cm}^{\text{2}}\leavevmode\nobreak\ \text{V}^{\text{-1}}\leavevmode\nobreak\ \text{s}^{\text{-1}}$ at a carrier concentration of $2\times 10^{11}$ cm*-2* is among the highest ever measured mobility in suspended graphene [7]. This value can be considered as the phonon-limited upper bound [10]. Despite such high mobilities, the weak electrostatic gating and absence of band gap makes it an undesirable candidate for 2D field-effect transistor (FET). However, it did not deter scientific community from continuing research on graphene because of its rich and unique fundamental physical properties. On the other hand, the inherent band gap of TMDs made them promising candidates for making FETs. So, research turned to this new family of 2D materials based on TMDs. Early studies of MoS2 found significantly lower mobility but with a much better electrostatic control. The subsequent introduction of a substrate [11, 12, 13, 14] for better electrostatic gating provided an opportunity to measure large-area graphene sheets, although substrate impurity [15], interface, and remote phonon [16] scattering significantly reduce carrier mobility.

In the subsequent years, researchers developed improved techniques for growing large-area samples of single crystalline graphene. Hao et al. [18] demonstrated reproducible growth of centimeter-sized single-crystalline graphene grains. This paved the way for wafer-scale growth of single crystalline domains of graphene [19] and MoS2 [20] in a lab environment. However, for commercial purposes, cost-effectiveness and high-throughput synthesis of wafer-scale single-crystalline 2D materials are the major requirements to compete with the existing silicon technology. To date, chemical vapor deposition (CVD) is the cheapest and fastest method of growing such large-area sheets of 2D materials,  [21, 22, 23, 24] but CVD-grown sheets are inherently polycrystalline in nature. The grains grow from each metal seed point until the point where they meet neighboring grains, where they form grain boundaries (GBs). Each grain is randomly oriented in space, and thus has different crystallographic orientation with angles varying between 0∘ to 60∘. Different crystallographic orientation of adjacent grains causes a mismatch in the crystal structure at their boundary, leading to the formation of pentagon-heptagon pairs [25, 17]. Macroscopically, the relative mismatch between the crystallographic orientations of adjacent grains is referred to misorientation angle. Figure 1 shows atomic-resolution ADF-STEM images of graphene crystals and a GB with misorientation angle of 27∘[17]. To realize practical applications of large-area 2D devices made up of graphene or TMDs, we believe that it is crucial to develop a detailed understanding of the impact of GB misorientation angles on electrical transport in polycrystalline sheets. To this end, numerous experimental  [26, 27, 17] and theoretical studies  [28, 29] have been carried out.

2.1 Charge transport across graphene grain boundaries

The earliest studies focused electrical resistance across graphene grain boundaries. Like GBs in bulk semiconductors, GBs in graphene were also found to impede electrical conduction. Experimental studies show that the resistivity of graphene GBs varies over a broad range from few $\Omega$ $\mu\text{m}$  [31] to several tens of k$\Omega$ $\mu\text{m}$  [32, 33]. There are many factors that are responsible for this large variation of resistivity. These include grain misorientation angles  [34], GB roughness, wrinkles and transition width [35], impurities and structural defects at the boundaries [36, 37, 38, 39], position of the Fermi level, and other differences arising due to processing and measurement techniques. Most of the time, it is challenging to isolate the effect of one of these factors from others. To minimize the large spread in resistivity, there has been a constant effort on improving device quality  [40, 41, 42, 43, 44, 45] by minimizing impurities and structural defects at the boundaries. Even in high quality large-area samples, the effect of misorientation angles on charge transport is quite significant.

Figure 2(a) shows a schematic of the Brillouin zones (BZs) of two adjacent grains. BZ of the left grain is rotated by $\Theta_{L}$ in the anti-clockwise direction, whereas the one on the right by $\Theta_{R}$ in the clockwise direction such that the relative misorientation is given by $\Theta_{M}=\Theta_{L}+\Theta_{R}$. $\Theta_{B}$ is the angle that boundary (blue-dashed line) makes with respect to the reference, shown in the Figure 2(b), (the hexagons are omitted for clarity and the red-dashed lines represent the crystallographic orientation of each grain). For a given misorientation angle $\Theta_{M}$, the type of GB is decided based on $\Theta_{B}$. For example, when $\Theta_{B}=0^{\circ}$ GB is symmetric, and is called twin GB. On increasing $\Theta_{B}$, that is by rotating the boundary towards the crystallographic direction of one of the grains, the GB becomes asymmetric, and is called tilt GB. Finally, when the boundary completely aligns with the crystallographic direction of the left/right grain, that is $\Theta_{B}=\Theta_{M}/2$, it is the most asymmetric GB that can be obtained for a given $\Theta_{M}$.

The transport of carriers across any GB requires quantum-mechanical wave continuity [46], based on which the energy as well as transverse momentum (momentum parallel to the GB/interface) of the incident carrier should be conserved on transmission. Mathematically, for every incident carrier $k_{i}$ we find a transmitted wavevector $k_{t}$ such that $k_{i_{\parallel}}=k_{t_{\parallel}}$ and $E(k_{i})=E(k_{t_{\parallel}}+k_{t_{\perp}})$, where $k_{i_{\parallel}}$ and $k_{t_{\parallel}}$ are the components of incident and transmitted wavevectors respectively along the GB direction. $k_{i_{\perp}}$ and $k_{t_{\perp}}$ are the components of incident and transmitted wavevectors respectively perpendicular to the GB direction. This is analogous to the Acoustic Mismatch model in phonon transport where the transmission across a GB is considered to be specular, that is the effect of GB roughness is negligible. Then the mode-dependent transmission coefficient $\tau_{b}(k_{i})$ for every $k_{i}$ per branch $b$ is calculated by {ceqn}

[TABLE]

Next, $\tau_{b}(k_{i})$ is converted to energy-resolved transmission coefficient $\Gamma_{b}(E)$ by averaging over the contour with constant energy, described by $\delta(E-E_{b}(k_{i}))$. This approach uses the 2D version of the linear extrapolation described by Gilat and Raubenheimer[47] as {ceqn}

[TABLE]

The denominator in Eq. (2) is the density of states per branch $D_{b}(E)$. The transport distribution function TDF can then be calculated as {ceqn}

[TABLE]

where mode-dependent velocity in the direction of transport $v_{b_{\parallel}}(k_{i})$ is converted to energy-resolved velocity $v_{b_{\parallel}}(E)$ in the same way as transmission coefficient. The TDF is then used to numerically compute GB conductance in the Landauer Formalism as follows {ceqn}

[TABLE]

$f(E)=[1+exp((E-E_{F})/k_{B}T)]^{-1}$ is the equilibrium Fermi-Dirac distribution function. In the past, models have also been developed to study the microscopic effects—atomic structure[48] and buckling [49]—on GB resistance, where periodic GBs are constructed using a dislocation model described by [25]. An asymmetric GB with $\Theta_{L}$ and $\Theta_{R}$ equal to 8.2∘ and 30∘ respectively, constructed from the dislocation model is shown in Figure 3. The quantum transmission and electronic transport properties were calculated from first principle Density Functional Theory (DFT) using a non-equilibrium Green’s function (NEGF) formalism implemented in the TRANSIESTA code[50].

Clark et al.  [51] measured GB resistivity $R_{GB}$ and found that there exists a positive correlation between $R_{GB}$ and misorientation angles $\Theta_{M}$, but the effect on resistivity due to a GB transition width could not be separated. In a recent study, we demonstrated that, although $\Theta_{M}$ plays an important role in determining $R_{GB}$, it is the relative angle of adjacent grains with respect to the GB that has more influence on $R_{GB}$ than $\Theta_{M}$ alone  [30]. A GB formed between a pair of grains which can be represented by same translational vectors (component of an unit vector representing crystal orientation of each grain along GB direction) are found to be highly conductive (about 80$\%$), whereas the ones with different translational vectors behave as reflectors of carriers  [48]. In other words, symmetric GBs are highly transparent to carriers in intrinsic (defect-free) graphene  [34].

Figure 4(a) and 4(b) show the dependence of GB resistance with carrier density for various misorientation angles in twin and tilt GBs, respectively. Several approaches have been employed to dope 2D materials which include diffusion or implantation of metals[52] and physi-/chemisorption of gases[53], but direct doping has been found to be challenging in 2D materials. Consequently, modulation of carrier density in 2D FETs is mostly achieved by electrostatic gating. The typical carrier densities obtained from gating vary between 1012 to 1013 cm*-2*. We found that GB resistance in twin GBs show a sharp dependence on carrier density, reaching the ballistic resistance value for nC equals 1013 cm*-2*. For tilt GBs, the resistance shows a decreasing trend even beyond 1013 cm*-2* due to the conical bandstructure of graphene.

The effect of misorientation angles on GB resistance is illustrated in Figure 4. Figure 4(a) shows a weak dependence of $R_{GB}$ on misorientation angle in symmetric GBs, whereas in a completely asymmetric GB, the resistivity varies over 20 orders of magnitude with misorientation angles shown in Figure 4(b). Symmetric rotation of BZ in twin GBs still conserve parallel momentum, and hence misorientation angles have negligible effect on GB resistivity. For tilt GBs, a fraction of the carriers cannot conserve parallel momentum, and hence get reflected which results in increase in GB resistivity. As graphene has a very steep and conical bandstructure near Dirac point (group velocity of about $10^{5}$ ms*-1*), even a small rotation in BZs results in a large fraction of wavevectors which cannot conserve energy across GBs. This appears as an opening in the energy spectrum where there is no transmission, often referred to as transport gap. Figure 4(c) clearly shows that for a given misorientation angle there is a wide spread of GB resistivity depending on boundary angle $\Theta_{B}$. Despite numerous studies on graphene GBs, we believe that a systematic experimental demonstration of the impact of misorientation angle on electrical transport in clean (defect-free) samples is still missing.

2.2 MoS2 Grain Boundaries

Owing to their intrinsic bandgap, TMDs and especially MoS2, have received significant research attention in recent years for electronic applications  [54, 55]. Although intrinsic room-temperature (RT) carrier mobilities in bulk and monolayer MoS2 have been theoretically predicted to be about 400 and 1000 cm2 V*-1* s*-1* respectively [56, 57], earlier studies reported a very low RT-carrier mobility of about 3 cm2 V*-1* s*-1* [3]. This was attributed to the dominant charged impurity scattering and trapped charges at the substrate interface. Over the years, there have been significant efforts on fabricating ultra-clean samples, which led to a gradual improvement in RT-carrier mobilities [58] upto more than 500 cm2 V*-1* s*-1* in bulk samples of MoS2  [59]. Beside the effect of surface impurities and trapped charges in single crystalline MoS2, there had been a significant research dedicated to understand the effect of GBs in CVD-grown polycrystalline sheets. The field-effect mobility has been found to depend weakly on the channel length [60, 61] for lengths up to several tens of microns in CVD-grown MoS2, shown in Figure 5. Moreover, electrical properties of CVD-grown monolayer MoS2 has been reported to be comparable to that of their exfoliated counterparts [62]. These studies indicate that GBs do not play a significant role in determining electrical properties in polycrystalline MoS2 samples [63]. Despite several studies on electrical properties of CVD-grown MoS2, there is little known about the impact of misorientation angle on charge transport in CVD-grown MoS2 sheets.

In our recent study [30], we explained why GB resistivity is weakly dependent on misorientation angle in both twin and tilt MoS2 GBs. In Figure 6(a), it can be seen that GB resistivity of twin GBs vs. misorientation angle is comparable to that of twin graphene GBs. Figure 6(b) shows that the resistivity of the most asymmetric-tilt GBs in MoS2 varies over two orders of magnitude as compared to 20 orders of magnitude variation of resistivity in graphene tilt GBs with misorientation angle. We attribute this to the flatter parabolic conduction band (CB) of MoS2. Even for large misorientation angles caused by the rotation of BZs of grains on either side of GB, there is a significant overlap in their CBs. As a result, a larger fraction of wavevectors conserve energy as compared to that in graphene. Consequently, the boundary angle parameterized by $\Theta_{B}$ does not largely influence GB resistivity in MoS2 GBs for a given misorientation angle $\Theta_{M}$, as can be seen in Figure 6(c). It has been attributed to the flat parabolic bandstructure of MoS2, in stark contrast to the steep cone-shaped bandstructure near the $K$-point in graphene. The density of states for MoS2 is step-like, which makes it difficult to move the Fermi level deeper inside the CB by gating. Thus, GB resistance has been plotted for the typical range of carrier densities—1012 to 1013 cm*-2*.

3 Lateral heterojunctions in two-dimensional materials

Owing to the potential of realizing all-2D devices including transistors, p-n junctions  [64, 65], superlattices  [66], tunneling devices  [67, 68], and lateral heterojunctions have drawn significant research interest in recent years. Few popular 2D lateral heterojunctions include graphene/TMD  [30, 69, 70, 71] TMD/TMD [64, 72, 73], and graphene/non-TMD heterojunctions  [74, 75, 76, 77, 78]. A recent review [79] has tabulated different applications of various types of 2D heterojunctions. Here we are restricting ourselves to graphene/TMD and TMD/TMD junctions. Metal-semiconductor junctions such as graphene/TMD heterostructures give rise to Schottky barrier height (SBH), whereas junctions between two semiconductors, such as TMD/TMD heterostructures, are categorized by their band offsets. Based on the band offsets, semiconductor-semiconductor interfaces can be either straddling gap (type-I alignment) or staggered-gap type (type-II alignment). MoSe2/ WS2, MoTe2/ WSe2, and WSe2/ WTe2 are examples of the former heterojunction-type, whereas MoS2/ WS2 and MoSe2/ WSe2 are the examples of the latter category [80].

3.1 Band alignment

There are different approaches to computing as well as aligning the electronic structure at the junction. One such approach is atomistic, based on first principles described by the widely-used Density Functional theory (DFT). In practice, most popular DFT implementations use a periodic basis, which makes it challenging to simulate non-periodic structures such as interfaces. A commonly used work-around is to make the structure artificially periodic within a larger expanded unit cell, called supercell shown in Figure 7. One of the major drawbacks of this method is that it requires large supercells to effectively isolate the interfaces in adjacent supercells from each other, thus magnifying the computational cost. Therefore, simpler and more efficient methods, validated by DFT calculations, are desirable. The most widely used such empirical model is the Schottky-Mott rule [81, 30], which says that the vacuum energy levels of the semiconductors on either side of the junction must be aligned, and thus allows band bending at the interface to account for different work functions and electron affinities. However, Schottky-Mott rule does not capture the charge redistribution that might happen when two materials are brought in contact. In an interface comprising of two three-dimensional materials, a local charge transfer leads to the formation of a plate capacitor. In contrast to this plate capacitor, a 2D heterojunction forms a pair of line dipoles [82]. In 2D heterojunctions, it has been shown that the effect of this interface dipole becomes negligible when the overall device dimensions are much larger than the characteristic junction-width, typically about 10 nm [83]. Consequently, the band alignment in such heterojunctions was found to be less sensitive to the details of interfacial structure. Furthermore, the alignment was found to closely follow the Schottky-Mott rule [30] that is widely used to describe 3D interfaces. Figure 8(a) shows an schematic of band-bending at the interface for unbiased (Vbg=0) and biased (Vbg¿0) gate conditions. The bandstructure of graphene with respect to MoS2 at the interface is shown in Figure 8(b) for different biasing conditions. Besides intrinsic band-bending, there could be additional band-bending due to extrinsic factors such as local strain arising from lattice mismatch [84], confinement effects [85], and the presence of impurities [55] at the junction. All these additional effects could change the transmission coefficient of carriers across the barrier, which in turn, will affect GB resistivity (see Equations 2 and  4).

3.2 Transport in graphene/MoS2 heterojunctions

Earlier studies used metals for contacts in MoS2 FETs, whose performance was largely undermined by the contact resistance at their interface. In general, it is difficult to achieve zero energy barrier between a semiconductor and metal, which is an important criterion that influences device performance. Since graphene is a semi-metal whose work function can be tuned, it was found to form Ohmic contact with MoS2 [81, 86]. This led researchers to investigate electrical transport in lateral heterojunctions with a particular focus on graphene/MoS2 heterostructures [82]. The interface resistance between graphene and MoS2 has been quantified as a function of gate voltage and found to vary with the height of the potential barrier between them, which was found to depend on the band alignment and vary with carrier concentration [69].

Despite several studies, both theoretical [87] as well as experimental [88], there is little known about the impact of misorientation angle on graphene/MoS2 lateral heterojunctions. Sun et al. [71] carried out first-principles calculations to study the relation between contact geometries and electrical properties of graphene/MoS2 heterostructures in four configurations—zig-zag/zig-zag, armchair/armachair, armchair/zig-zag, and zig-zag/armchair. They found that MoS2 is metalized due to the occurrence of band gap states at the contact interface which also reduces the Schottky barrier height (SBH) in all the four configurations. Recently, we systematically investigated the effect of misorientation angles on interface resistivity at graphene/MoS2 interface [30] and found that the interface resistivity in symmetric graphene/MoS2 interfaces (Class-I heterojunctions) is less sensitive to mismatch angle as shown in Figure 9(a), but the resistivity is much larger than those of symmetric homojunctions. On the other hand, Figure 9(b) shows that the resistivity of asymmetric graphene/MoS2 interface (Class-II heterojunctions) strongly depends on mismatch angle. Like in graphene tilt GBs, this is also due to the underlap between CBs of graphene and MoS2 on rotation of their BZs.

3.3 Resonant tunneling diodes in 2D lateral heterostructures

Tuning the band alignment at lateral heterojunction interface can open up interesting applications by optimizing their electronic and optical properties. Using one-step tube-moving chemeical vapor deposition method [89], a lateral WS2/WS2(1-x)Se2x hetrostructure was synthesized. Controlling the ratio of WSe2 and WS2 in evaporaton sources enables the tuning of the band alignment in heterostructures. Another approach to tuning the bandgap [90] in MoS2 nanoribbons is introducing the sulfur vacancies, which enabled the creation of double barrier quantum well structure (DBQW). Following this approach [91] and using NEGF based on tight-binding calculations, a DBQW of A-MoS2 nanoribbons with two sulphur line vacancies was simulated. The resultant analysis showed that the double barrier structure has a PVR (peak-and valley-current ratio) of about 78 at room temperature.

One of the first demonstrations [92] created pristine multi-junction heterostructures by direct synthesis based on graphene, MoS2, MoSe2 and WSe2. This lead to the creation of resonant tunneling diode (RTD) in an atomically-thick stack with negative differential resistance (NDR). The spectrally-narrow NDR created in this structure, as shown in Figure 10, is superior to manually stacked heterostructures that were studied previously. TMDs in the symFET architecture  [93] have been shown to achieve PVR up to $10^{9}$ compared to few hundred observed in graphene or III-V RTDs. Also, the ability to achieve such high PVR even in nanoscale devices makes them promising candidates for digital logic applications. A vertical heterostructure based on phospherene/rhenium disulfide (BP/ReS2) demonstrated NDR with a PVR up to 6.9 at 180 K  [94]. A ternary inverter was also created from this heterostructure that exhibited three distinct logic values, demonstrating a proof of concept for future multi-valued logic devices.

4 Vertical 2D heterostructures

For several device applications, the 2D material has to be placed on a substrate, typically SiO2. Impurity scattering from the substrate significantly deteriorates carrier mobility as compared to their intrinsic phonon-limited mobility [95, 11, 96]. To minimize the reduction of mobility of graphene on highly disordered SiO2 [97, 98], heterojunction made up of graphene and hexagonal-boron nitride had been proposed which demonstrated improved carrier mobility [13]. Moreover, the possibility of fabricating low power devices using band-to-band tunneling (BTBT) and achieving steep subthreshold slope below 60 mV/dec paved the way for devices made up of vertical heterostructures. Vertical heterostructures, which are typically referred to as van der Waals heterostructures, have been recently found to be of research interest because of the large contact area at the interface between the two materials and absence of dangling bonds. As the name suggests, the layers are held together in the vertical direction with weak van der Waals forces as compared to the strong covalent bonding in lateral heterojunctions. Exciting physics and potential device applications in stretchable/flexible electronics have been demonstrated by making field-effect transistors [99, 100], memory devices  [101, 102, 103], and photodetectors [104, 105] by vertically stacking two-dimensional materials.

Earlier studies used mechanical exfoliation to peel off layers of various materials from their bulk counterparts and then used a micro-manipulator to stack and align the layers in a vertical configuration. The mechanical transfer techniques make the interfaces between layers highly susceptible to impurities and contaminants [106, 107], which made it difficult to fabricate vertical heterostructures with sharp and clean interfaces. In subsequent years, researchers devised novel fabrication techniques by employing direct-synthesis of such vertical stacks using single-step vapor-phase growth process, which yields much cleaner and sharper interfaces [108, 109]. Recently, CVD growth of centimeter-scale 2D vertical heterostructure has also been demonstrated [110]. All the aforementioned fabrication techniques and potential applications of 2D van der Waals heterostructures have been discussed in great details in the following reviews [111, 112].

4.1 Charge transport in vertical 2D heterojunctions

The charge conduction across different materials in a vertical stack takes places either through BTBT or modulation of SBH at a 2D-3D interface. BTBT, also called interlayer tunneling, was first demonstrated in 2D vertical heterostructures by Britnell et al. [99] where an ultrathin layer of hexagonal boron nitride (hBN) was sandwiched between two graphene monolayers. The tunneling current was regulated by varying the gate voltage, which altered both carrier densities in graphene layers as well as internal barrier voltage. Georgiou et al. [100] replaced hBN with WSe2 to improve the ON-current and switching ratio from $10^{5}$ to $10^{6}$. Roy et al. [113] fabricated a MoS2/WSe2 vertical heterostructure in a dual-gate architecture, where the carrier densities of each material were modulated independently using a pair of symmetric gates. The same device exhibited large reverse-biased current caused by BTBT, an Esaki diode behavior with negative differential resistance (NDR) or a forward-biased rectifying behavior depending on different gating conditions. Nourbakhsh et al. [114] attributed the presence of NDR in such devices to lateral BTBT. Their results also confirmed that lateral BTBT is more dominant than the tunneling through the overlapped vertical section.

Recently, Sarkar et al. [115] demonstrated a 2D/3D semiconductor planar transistor based on BTBT which achieved the record-low subthreshold swing with a minimum value of 3.9 mV/dec and an average of 31.1 mV/dec (5.5, 12.8, 22 mV/dec) over four (one/two/three, respectively) decades of drain current. They named it ATLAS, which stands for atomically thin and layered semiconducting-channel tunnel FET. They strategically chose germanium as the 3D material among all other Group IV and III-V compounds and MoS2 as the 2D material among other TMDs to form a staggered heterojunction with a low tunneling barrier height. Another interesting vertical three-terminal device is called barristor, a portmanteau of barrier and transistor, in which the gate voltage is used to regulate tunneling through the SBH formed at their interface instead of BTBT. Yang et al. [116] demonstrated the first graphene barristor using a tri-layer stack (metal/graphene/silicon) with a large on-off ratio of $10^{5}$. Yu et al. [117] used the same idea to fabricate a 2D-2D barristor using metal/MoS2/graphene stack to improve the current density by about two to five orders of magnitude while maintaining a high on-off ratio ($>10^{3}$). Recently, twist angle (misorientation angle between layers in vertical direction) has been used in bilayer graphene to demonstrate superconductivity at few magic angles, the first being at $1.1^{\circ}$. At such magic angles, twisted bilayer-graphene exhibits flat bands near the Dirac point resulting in increase in density of states [118]. In another recent work, Liao et al.[119] found that twist angles play a crucial role in determining vertical conductivity in a graphene/MoS2 stack. Like our recent study on lateral graphene/MoS2 heterojunction[30], the conductivity has been found to be highest for twist angle equal to $0^{\circ}$ and minimum for $30^{\circ}$ in graphene/MoS2 vertical heterojunction. The reason of this has been attributed to the reduction in transmission coefficient with twist angles.

4.2 Electronic thermopower in 2D heterostructures

In addition to having large electron mobilities, some 2D materials have exhibited interesting thermoelectric (TE) properties, including possessing a large Seebeck coefficient, also called thermopower. It is defined as the ratio of the voltage produced by thermally diffusing electrons and the temperature gradient that is driving their diffusion from the hot side to the cold. Research in this area started with the measurements of TE transport in graphene [120], which exhibited anomalous TE transport of Dirac particles [121]. More recently, the enhanced TE Seebeck coefficient in graphene [122] was attributed to the role of hydrodynamic transport through inelastic scattering, which also leads to a violation of the Mott relation for the Seebeck coefficient (not to be confused with the Schottky-Mott rule for band alignment)

[TABLE]

Proposals to enhance the TE properties of graphene used heterostructures [123] and functionalization [124]. Graphene also exhibited a significant phonon drag component of the Seebeck coefficient [125], which is defined as the additional thermopower caused by the exchange of momentum between the heat-carrying lattice vibrations (phonons) and the electrons. This effect is particularly prominent at low temperatures where the diffusion thermopower is typically low.

Using both carrier and phononic engineering on a ZrSe2/HfSe2 single-layer superlattice [126], a ZT of 5.3 in n-type and 3.2 in p-type device was achieved. First-principles calculations along with Boltzmann transport equation showed that this is achieved due to the high degenerate nature of the conduction bands in n-type device. A partially overlapped graphene/graphene vertical heterostruture [127] studied using atomistic tight-binding Hamiltonian demonstrated that the ZT can reach unity at room temperature. An analysis [128] on van der waals heterostructures with multilayer TMDs sandwiched between two graphene electrodes identified WSe2 and MoSe2 as the ideal TMDs that can provide high TE conversion efficiency. Graphene/h-BN/graphene vertical heterostructures are synthesized [129] to study their TE performance. Measurements uncovered a significant Seebeck coefficient at the material interfaces, which makes such heterostrutures suitable candidates for TE applications. Twisted bilayer graphene vertical heterostructures are studied extensively due to their interesting properties like magic angles that can induce superconductivity, described in the previous section. With regards to their thermoelectric properties, measurements on twisted graphene bilayers [130] showed exceptional cross-plane thermopower, attributed to phonon drag from out-of-plane phonon modes (ZA/ZO’). Besides twisted graphene bilayers, MoS2 has been shown to also exhibit a phonon drag effect [131].

5 3D heterostructures and superlattices

With the discovery of NDR in narrow highly doped germanium p-n junctions by Esaki [132], a new era of electronics based on quantum tunneling devices was ushered in. These p-n junctions were called tunnel diodes as they demonstrated interband tunneling of carriers from the valence band on the p-side to conduction band on the n-side and vice versa. For this discovery, Leo Esaki shared the 1973 Nobel Prize in Physics with Brian Josephson. Around the same time, Tsu and Esaki [133] found that in the direction of a one-dimensional periodic potential, such as that found in a heterostructure (superlattice), electrons are localized to discrete energy states analogous to electrons in a two-dimensional electron gas. Further, due to the comparable dimensions of these potentials with electron wavelengths, the wave nature of electron leads to phenomena such as interference and tunneling. In such structures, when the energy of electrons coincides with one of the discrete energy states achieved by tuning the applied bias, electron can tunnel through the barriers and can have near-perfect transmission. With further increase in the applied bias, the amount of available electrons decreases, thus reducing the current flowing through the structure. Ultimately, the current increases again when high enough bias is applied to cause thermionic emission of electrons over the top of the potential barrier. This distinct feature is popularly known as negative differential resistance and has found numerous applications in devices such as high-frequency oscillators, frequency converters, and detectors and also exhibits great potential in high-speed logic devices and switches.

5.1 Dual barrier heterostructure: resonant tunneling diode

Shortly after the pioneering discoveries of Tsu and Esaki, it was shown [134] that the same effect can be achieved using a finite superlattice. With the advancements made in fabrication process, a finite superlattice with a double barrier is created [135] in GaAs by sandwiching it with Ga1-xAsx to form the potential barriers and demonstrated NDR phenomenon. A basic RTD is made up of a double-barrier heterostructure with nanoscale barriers made from a heavily doped semiconductors. Beside tunneling, there are other physical processes like interactions of carriers with lattice, impurities, surface roughness, other carriers and alloy disorder. A realistic physical model of RTDs including all these processes requires advanced quantum transport theory like density matrix, Wigner functions, and non-equilibrium Green’s functions (NEGF).

5.2 Progress in numerical simulation of 3D heterostructures

The semi-classical Boltzmann transport equation, which is used widely to simulate electron and thermal transport in electronic devices, is not adequate to capture quantum effects that are predominant in RTDs. Various quantum-transport frameworks such as the Wigner formalism and the related density matrix approach, as well as non-equilibrium Green’s functions (NEGF), have been employed to study transport in RTDs. NEGF was one of the first methods used to study the transport in RTDs. Lake and Datta [136] studied the effect of energy broadening and inelastic scattering and reported an enhanced valley current due to inelastic scattering with simultaneous enhancement of occupation of the resonant states. Time-dependent transport capabilities of NEGF were demonstrated by simulating RTD [137] in mesoscopic region. Nam Do et al. [138] used fully self-consistent non-equilibrium Green’s function approach to study the the impact of quantum-well width, the barrier thickness and the temperature, and showed their effect on PVR. PVR, as the name implies, is the ratio between the currents at the peak before the onset of the NDR and the valley after it, which helps to determine the feasibility of using the RTD for device applications.

Wigner formalism is one of the first methods [139] used to simulate the transport behavior of quantum-well RTD and their NDR. Biegel and Plummer [140] implemented a self-consistent Wigner-function-based quantum device simulation using the RTD as a test case. A comparison between different iterative methods was presented here to solve the Wigner function. It concluded that transient Gummel approach is reasonably accurate with low computational resource requirement for device simulations. A unified approach [141] was proposed to merge Wigner functions with semi-classical Boltzmann transport equation and treat the scattering term as a generating term. In that study, interference effects due to the Wigner potential were associated with particle generation having statistical weights.

5.3 Thermoelectric properties of superlattices

Quantum effects in heterostructures have a large impact on their thermoelectric properties, affecting both the Seebeck coefficient $S$, which is related to the average transport energy per carrier, and the conductivity $\sigma$. Together, the Seebeck coefficient and conductivity constitute the TE power factor $PF=S^{2}\sigma$, which is central to the TE conversion efficiency captured by the figure-of-merit $ZT=S^{2}\sigma/\kappa$. $\kappa$ is the sum of electronic and lattice contributions to the thermal conductivity. In their seminal work, Hicks and Dresselhaus[142] predicted a one order of magnitude increase in figure-of-merit (ZT) for quantum well superlattices made from Bi2Te3. High anisotropy along certain direction at low widths of quantum wells are shown to be the reason for this enhancement. This kicked off extensive research into understanding their thermoelectric properties and possible use in thermoelectric generators. Their predictions were later revised [143] when it was shown that the enhancement is controlled by period of the superlattice and not by the width of the well. Also, the overall $ZT$ is projected to be less than one when tunneling is introduced between quantum wells. This idea was experimentally validated using PbTe/Pb1-xEuxTe superlattices [144] with a five fold increase in $ZT$ compared to bulk PbTe. A more extensive analysis [145, 146] for the case which includes a non-zero barrier and effects of carrier tunneling, revealed that the in-plane thermoelectric performance is not superior to bulk systems. Later short-period Si/Ge superlattice structures were shown [147] to provide higher $ZT$ compared to thin-film SiGe and bulk Si-Ge alloys, which was attributed to the reduction in thermal conductivity. A further improvement [148] in $ZT$ is reported in Si/Ge superlattices by applying strain to further tune the conduction band structure. Experimental studies demonstrated the thin-film thermoelectric coolers using both single heterostructure and superlattice structures from SiGe/Si [149, 150, 151, 152, 153, 154].

Using an exact numerical solution, a complete treatment of powerfactor in PbTe quantum well and quantum wire superlattices showed [155] weaker dependencies on potential barrier. Extending the study to GaAs quantum wire superlattices [156], it was established that significant enhancements in power factor can be achieved from restricting phonon transport. PbTe/PbTeSe-based superlattices are experimentally shown [157, 158] to improve $ZT$ by 50% relative to bulk PbTe. This can be attributed to a reduction in thermal conductivity, an affect of the increase in scattering due to alloying. Later, majority of studies used this approach to improve $ZT$. Ge quantum dot in Si quantum dot superlattices were shown [159] to improve $ZT$ by reducing the thermal conductivity due to scattering of phonons. Theoretical calculations of TE properties in superlattice nanowires [160] based on lead salts (PbS, PbSe and PbTe) exhibited a significantly higher ZT values, especially for 5 nm diameter wires.

Even though significant improvement in $ZT$ is reported by restricting thermal conductivity, further improvement can be achieved by tuning the TE power factor. Many studies proposed mechanisms and ways to improve power factor in superlattices by tuning the band structure [161, 162, 163]. Energy filtering, a process of restricting movement of carriers with kinetic energy smaller than barrier height, has been consistently found to enhance the overall power factor. It raises the average energy of carriers and thereby, its Seebeck coefficient. Nanocomposites [164, 165, 166], superlattices [167] as well as single and multiple barrier structures have been used to introduce some form of energy filtering in a material. Studies [168, 169, 123, 170] tried to understand the effect of potential barrier smoothness and their structure on thermoelectric performance. Using a series of gates placed periodically in graphene, creating periodic potentials [123, 170] that resulted in a drastic increase in Seebeck coefficient. Thermionic emission over the barrier and tunneling through the barrier, control transport in these structures as shown in Figure 11. Semi-classical transport (diffusive transport) theory like Boltzmann transport equation is capable of recreating the thermionic emission over potential barriers. Quantum effects like tunneling are simulated using density matrix formulation, Wigner formalism and non-equilibrium Green’s functions. Here we aim to understand the strengths and limitations of these approaches in simulating the quantum effects along with the diffusive transport.

5.3.1 Wigner formalism

In the semi-classical Boltzmann transport equation (BTE), widely used in device simulation, electrons are treated as point particles [172, 173]. Despite its widespread use, this approach is unable to explain some device effects, such as underestimating the threshold voltage in ultra-thin body MOSFET [174, 175], and the carrier interactions with the rapid potential changes across heterojunctions. The impact of sharp spatially-varying potentials on transport in a superlattice can be simulated using the Wigner formalism. The Wigner formalism [176, 177] includes an additional quantum evolution term in addition to drift term that can capture the spatial variation in potential. The Wigner formalism has been the subject of an extensive recent review article [178]. BTE can be modified as follows to include Wigner formalism for rapid varying spatial potentials:

[TABLE]

where ${f}_{w}\left(r,k,t\right)$ is called Wigner distribution function, which is written as

[TABLE]

where the mixed state is represented by the density operator $\rho$. Here the coordinates $r$ and ${r}^{\prime}$ represent the center of mass and spread of the electron wave packet, respectively. The potential operator $Q{f}_{w}\left(r,k,t\right)$, often called the quantum evolution operator, is given as

[TABLE]

where the Wigner potential $V_{w}(r,k)$ is itself also obtained through a Wigner transform

[TABLE]

of the electrostatic potential $V$ across the material, expressed in the new Wigner coordinates $r$ and $r^{\prime}$. Eq. 9 can be simplified as

[TABLE]

where $\hat{V}(k)$ is spatial Fourier transform of $V$

[TABLE]

The potential operator $Q{f}_{w}\left(r,k,t\right)$ is usually decomposed into two components: a slowly-varying, often called classical, potential (such as the applied external bias) Vcl, and separate, rapidly-varying quantum-mechanical portion Vqm. Together, these two make up the potential according to

[TABLE]

By including the additional effect of rapidly varying potentials in the BTE, the resulting steady-state Wigner-Boltzmann transport equation (WBTE) can be written as (from Eq. 6)

[TABLE]

The collision term helps to further incorporate the semi-classical effects like energy and momentum relaxation from time-dependent perturbation theory. The Wigner formalism permits us to couple the interdependencies between quantum and semi-classical effects, as it uses a phase-space formulation. This allows us to study the carrier energy relaxation in the presence of both phonon scattering and periodic potential barriers. The use of conventional semi-classical boundary conditions for contacts in the Wigner formalism has been found to result in unphysical results, especially in the coherent regime. This is due to non-unique solutions [179, 180] and the nonlocal nature of Wigner approach. Dissipation and decoherence phenomena, even though they do lead to a unique solution, can also lead to unphysical results. Simplified local scattering models to define the dissipative transport can also cause negative probability-density [181, 182]. Alternate approaches using Lindblad-type scattering superoperators  [183, 184] and conditional wave functions [185] have been proposed to solve these problems within the Wigner approach. These limitations are avoided in periodic systems, where contacts are replaced with periodic boundary conditions [171].

In our recent study [171], we developed a numerical model that combines a full bandstructure with the WBTE, with the goal of studying the impact of quantum effects on thermoelectric transport in semiconductors. It was implemented by solving the full-band WBTE iteratively in a manner analogous to the widely-used Rode’s approach [186, 187] to study the influence of shape, $V_{0}$, and $L_{p}$ on thermoelectric performance. In Rode’s method, the accuracy of the solution to BTE is improved by including contributions from inelastic scattering processes, especially the in-scattering. This is achieved by splitting the BTE collision integral into two separate terms—the in-scattering and the out-scattering. Then the perturbation $g(k)$ to the equilibrium distribution function $f_{0}(k)$ (where $f_{w}(k)=f_{0}(k)+g(k)$) is solved iteratively until convergence. We refer the interested readers to Appendix A in [171] for a detailed derivation. The final iteration is expressed as

[TABLE]

The iteration is started from an equilibrium distribution, such that the initial in-scattering term is zero and the first update is straightforward to compute. The first iteration corresponds exactly to the relaxation time approximation (RTA). Next, Rode’s approach is applied to solving the WBTE by adding the Wigner potential contribution representing the additional forces arising from the potential variation in the structure. Then, the resultant perturbation to the distribution function in Eq. 14 can be modified as

[TABLE]

To study the effect of periodic potential structure (shape and size), a square barrier with smoothening parameter $\beta$ and a smooth cosine shaped potential barrier are used. A generalized potential that is spatially periodic can be represented as

[TABLE]

The quantum evolution (Eq. 8) (derived in Appendix B of [171]) of this generalized form for periodic potentials is

[TABLE]

The quantum evolution force for square barriers,

[TABLE]

of height $V_{0}$ with smoothening factor $\beta$ and width $2a$ or $L_{p}/2$ is

[TABLE]

Similarly for a cosine shaped potential (applied in $r$ direction) of form

[TABLE]

the quantum evolution is obtained as (Appendix B in [171]),

[TABLE]

[TABLE]

It should be noted that the quantum evolution force for a smooth potential is the first order approximation (m=1) to the quantum evolution force of a square barriers (Eq. 17).

The solution for the perturbation $g(k)$ to distribution function (Eq. 15) in WBTE can now be used to study the thermoelectric behavior of the system. First, to study the effect of potential barrier structure (apart from smoothening), the cosine shaped barriers are used. Introducing potential barriers in highly doped (ND = 1019 cm*-3*) silicon with Fermi level of 46 meV below conduction band edge $E_{c}$ ($E_{c}$ is set to zero), $g(k)$ is calculated. Cosine-shaped potential barriers with V0 = 1.6 $k_{B}$T $\cong$ 41 meV at T = 300 K ($k_{B}$ is the Boltzmann constant), show an decrease in Seebeck coefficient with increase in the period length (as seen in Figure 12(a)). Increase in period Lp from 3 nm to 9 nm, restricts the tunneling of electrons through the barrier causing a decrease in the Seebeck coefficient. This restriction in tunneling also reduces the electrical conductivity with increase in period length. At Lp = 9 nm, tunneling of carriers vanishes and the enhancement in Seebeck coefficient is a consequence of thermionic emission alone ($\cong$ V0/T $\cong$0.13 meV). Also, by nanostructuring (confining the structure along the directions apart from the transport direction) in the system, we observe a modest decrease in both the Seebeck coefficient and electrical conductivity due to the carrier scattering at the boundaries. Nonetheless, the reduction in thermal conductivity due to nanostructuring produces a significant boost to the overall TE performance.

When the period length is kept constant, increasing the potential barrier height (V0) results in a simultaneous increase in electrical conductivity and Seebeck coefficient, as shown in Figure 12(b). Even though the increase in Seebeck coefficient can be attributed to energy filtering, the simultaneous increase in electrical conductivity is a consequence of higher mobility electrons that drive most of the transport. The increase in mobility of electrons causes the conductivity increase with the potential barrier height. Again at wider period lengths (Lp = 9 nm), lack of tunneling and the presence of complete diffusive transport, transport occurs only through thermionic emission. It should be noted that irrespective of the potential barrier height, the uniform doping concentration in the system maintains the number of carriers available.

Apart from increasing the barrier height, another way of tuning the barrier seen by the carriers is by varying the position of Fermi level. The position of Fermi level EF is controlled by doping and it strongly influences the Seebeck coefficient, since the Seebeck is defined as the average energy of electrons that participate in the transport relative to the Fermi level. As shown in Figure 12(c), the system acts as bulk/no barrier structure when Lp = 9 nm, replicating the usual TE behavior (an increase in Seebeck coefficient and decrease in conductivity, with higher difference in V0 and Ef). Moving Fermi level towards the conduction band by doping the material gives an initial increase in Seebeck coefficient due to combined effect of filtering and an increase in carrier availability. A peak in Seebeck coefficient can be observed when V0-Ef $\cong$ 0.14 eV, which is consistent with previous studies[168]. It is followed by a reduction as the effective barrier height encountered by the carriers decreases and increase in the number of low-energy carriers that can overcome the barrier. An increase in doping results in higher electrical conductivity, which increases the overall power factor at higher doping conditions(smaller V0-Ef).

To study the effect of barrier smoothness, the parameter $\beta$ that controls the smoothness of a square potential is varied (Eq. 18). As shown in Figure 12(d), with increase in $\beta$, Seebeck coefficient increases due to the additional quantum reflections that are introduced in the quantum operator (as in Eq. 17). These additional reflections at the barrier helps to better filter the carriers compared to smooth barriers. The strength of these reflections increases with the weight for quantum evolution term, which itself is a function of $\beta$. Further, shape parameters like the duty cycle (ratio of potential barrier and the period) also effect the TE performance. Even for sharper barriers (higher $\beta$), Seebeck coefficient is maximum when the barriers are symmetric (50% duty cycle) as shown in the inset of Figure 12(d).

5.3.2 Non-equilibrium Green’s functions

The ability to simulate electron-electron correlations in the time domain is a central feature of the non-equilibrium Greens’s function (NEGF) formalism, which allows it to capture and simulate many-body quantum effects in nanoscale devices. A comprehensive review on nanoscale device modeling using NEGF was published by Datta  [188]. Here we provide a brief overview of the technique before delving into its recent applications in studying transport in RTDs and thermoelectric devices. The equations of motion in the NEGF formalism in the steady state are given by the Dyson and Keldysh relations [189, 190, 191] as

[TABLE]

[TABLE]

where $\hat{V_{c}}$ represents the coherently treated interactions that include the mean-field Coulomb potential $\hat{V_{c}}^{m-f}$ and the electronic disordered part $\hat{V_{e}}^{rand}$. The spectral self-energy used to treat the incoherent interactions is expressed as

[TABLE]

where the retarded self-energy used to include the electron-phonon interactions is defined as

[TABLE]

The expectation values of the observables and current operator are calculated using the relation between density matrix and Green’s functions.

The NEGF technique has found applications in understanding the impact of quantum effects on electronic and thermoelectric transport, as well as finding ways to further improve TE properties through confinement in heterostructures. The effect of confinement on the thermoelectric performance in Si/Ge/Si superlattices was simulated using NEGF. Strain-induced energy splitting in Si/Ge/Si superlattices is shown to improve power factor by four orders in magnitude [192]. However, the gains in TE performance expected in such structures was shown to be limited by the reduction in conductivity of superlattices with thin barriers [193]. Quantum wire superlattices with lateral confinement were also studied with NEGF, while including the scattering processes due to electron-phonon couplings, phonon anharmonicity, charged impurities, surface and interface roughness and alloy disorder [194]. Thermionic emission and tunneling of carriers in periodic superlattices were simulated in the NEGF formalism [168, 169], and their effect on TE properties has been found to improve the power factor, as shown in Figure 13. These studies also investigated the optimum size and shape of the potential barriers in Si-based superlattices in order to further improve TE performance. They found that hierarchically designed materials, containing heterostructures at varying scales from the atomic to the microscopic, can significantly impact transport and improve the TE power factor [195]. Double-barrier RTDs were simulated in the NEGF formalism, with self-consistent treatment of inelastic scattering, to show the enhanced valley current which is a consequence of enhanced occupation of resonant state  [136].

6 Conclusions

We have reviewed recent progress in understanding electronic and thermoelectric transport in heterostructures of 2D and 3D materials. Transport in both lateral and vertical heterostructures is strongly influenced by the band alignment between the two materials while the transmission of carriers between them depends on the simultaneous momentum and energy conservation in the transmission process. Consequently, homojunctions and, to a lesser degree, heterojunctions in 2D materials have conductances that strongly depend on the angle of rotational mismatch between the two adjacent materials in the heterostructure. Thermoelectric properties, including the Seebeck coefficients, are enhanced through heterostructuring in both the 2D and 3D case and offer promise for future on-chip waste-heat energy harvesting applications. Understanding quantum transport in 3D heterostructures and superlattices requires sophisticated modeling based on the Wigner or NEGF approaches. Open problems include simultaneously capturing quantum effects such as tunneling and semi-classical processes such as phonon and impurity scattering, as well as computational optimization of complex multi-layered heterostructures in order to enhance their properties. The latter involves simultaneous materials selection, doping, layer thickness, and, in the 2D case, angle of rotation; this multifaceted computational problem might benefit in the future from machine learning and artificial intelligence techniques in order to identify heterostructures with optimal electronic or thermoelectric properties.