Semi-analytical model of the contact resistance in two-dimensional semiconductors

TL;DR

This paper presents a semi-analytical model for contact resistance in 2D semiconductors that accounts for 2D transport and Schottky barrier modulation, explaining experimental metallic behavior in contact resistance.

Contribution

The work introduces a semi-analytical model based on Bardeen's approach that captures key effects neglected in previous models, improving understanding of contact resistance in 2D semiconductors.

Findings

Model reproduces metallic behavior in contact resistance at high gate voltage.

Accounts for temperature dependence and Schottky barrier modulation.

Explains experimental observations of contact resistance in 2D devices.

Abstract

Contact resistance is a severe performance bottleneck for electronic devices based on two-dimensional layered (2D) semiconductors, whose contacts are Schottky rather than Ohmic. Although there is general consensus that the injection mechanism changes from thermionic to tunneling with gate biasing, existing models tend to oversimplify the transport problem, by neglecting the 2D transport nature and the modulation of the Schottky barrier height, the latter being of particular importance in back-gated devices. In this work, we develop a semi-analytical model based on Bardeen's transfer Hamiltonian approach to describe both effects. Remarkably, our model is able to reproduce several experimental observations of a metallic behavior in the contact resistance, i.e., a decreasing resistance with decreasing temperature, occurring at high gate voltage.

Figures

Equations

$$F=\frac{V_G}{a+({ϵ}_s/{ϵ}_{\mathrm{ox}})t_{\mathrm{ox}}},$$

$$E_i=-qF\left(a-\frac{|{\zeta }_i|}{k_F}\right),$$

$$k_F={\left(\frac{2m_{z}qF}{{\hslash }^2}\right)}^{1/3},$$

$${\zeta }_i\approx -{\left[\frac{3\pi }{8}(4i-1)\right]}^{2/3}.$$

$$\frac{1}{\tau_{i}}=\left\{\begin{array}[]{ll}\frac{1}{h}\frac{\hbar^{2}}{2m_{z}}\frac{k_{F}^{2}}{\mathrm{Ai}^{\prime 2}(\zeta_{i})}\frac{4\sqrt{-E_{i}(E_{i}-V_{0})}}{-V_{0}}e^{-2\gamma_{0}},&\,V_{0}<E_{i}<0\\ 0,&\,\text{otherwise}\end{array}\right.$$

$${\gamma }_0$$

$${\zeta }_0$$

$$|{\mathrm{Ai}}^{\prime }({\zeta }_i)|\approx \frac{1}{\sqrt{\pi }}{\left[\frac{3\pi }{8}(4i-1)\right]}^{1/6}.$$

$$G$$

$$F_0(\epsilon )$$

$$\overline{T}(\epsilon )=\sum _i{T}_i(\epsilon ),$$

$$T_{i}(\varepsilon)=\left\{\begin{array}[]{ll}0,&\quad\varepsilon<E_{i}\\ 1-e^{-\frac{2L_{x}}{\lambda_{i}}},&\quad\varepsilon>E_{i}\end{array}\right.\,,$$

$$G\approx \frac{2q^2}{h}{F}_0(E_1).$$

$$G\propto T^{1/2}\exp \left(-\frac{{\Phi }_1}{k_{B}T}\right),$$

$$G\approx \frac{2q^2}{h}\frac{\sqrt{m_x{m}_y}}{\hslash }\frac{L_x}{{\tau }_1}{f}_0(E_1),$$

$$\mathcal{H}=-\frac{\hbar^{2}}{2}\nabla\cdot\hat{m}^{-1}\nabla+E_{c}(z)=\mathcal{T}+E_{c}(z),\quad\hat{m}=\left(\begin{array}[]{ccc}m_{x}&0&0\\ 0&m_{y}&0\\ 0&0&m_{z}\end{array}\right),$$

$$E_{c}(z)=\left\{\begin{array}[]{ll}V_{0},&\quad z<0\\ -qFz,&\quad 0<z<a\\ \infty,&\quad z>a\end{array}\right.$$

$$E_{cL}(z)=\left\{\begin{array}[]{ll}V_{0},&\quad z<0\\ 0,&\quad z>0\end{array}\right.$$

$$E_{cR}(z)=\left\{\begin{array}[]{ll}-qFz,&\quad z<a\\ \infty,&\quad z>a\end{array}\right.$$

$$P_{\alpha \beta }=\frac{2\pi }{\hslash }{\left|⟨{\psi }_{R,\beta }|\mathcal{H}-{\mathcal{H}}_L|{\psi }_{L,\alpha }⟩\right|}^2\delta \left(E_{L,\alpha }-E_{R,\beta }\right),$$

$$P_{\beta \alpha }=\frac{2\pi }{\hslash }{\left|⟨{\psi }_{L,\alpha }|\mathcal{H}-{\mathcal{H}}_R|{\psi }_{R,\beta }⟩\right|}^2\delta \left(E_{L,\alpha }-E_{R,\beta }\right)=P_{\alpha \beta },$$

$$\mathcal{H}=\left\{\begin{array}[]{ll}\mathcal{H}_{L},&\quad z<0\\ \mathcal{H}_{R},&\quad z>0\end{array}\right.$$

$$\mathcal{H}_{T}=\left\{\begin{array}[]{ll}-\mathcal{H}_{R},&\quad z<0\\ -\mathcal{H}_{L},&\quad z>0\end{array}\right.$$

$$⟨{\psi }_{R,\beta }|\mathcal{H}-{\mathcal{H}}_L|{\psi }_{L,\alpha }⟩$$

$$={\int }_{{\Omega }_R}{\psi }_{R,\beta }^{\ast }\left(\mathcal{H}-{\mathcal{H}}_L\right){\psi }_{L,\alpha }{d}^{3}r,$$

$$⟨{\psi }_{R,\beta }|\mathcal{H}-{\mathcal{H}}_L|{\psi }_{L,\alpha }⟩$$

$$={\int }_{{\Omega }_R}\left[{\psi }_{R,\beta }^{\ast }\left(\mathcal{T}-E_{L,\alpha }\right){\psi }_{L,\alpha }-{\psi }_{L,\alpha }\left(\mathcal{T}-E_{R,\beta }\right){\psi }_{R,\beta }^{\ast }\right]d^{3}r.$$

$$⟨{\psi }_{R,\beta }|\mathcal{H}-{\mathcal{H}}_L|{\psi }_{L,\alpha }⟩={\int }_{{\Omega }_R}\left[{\psi }_{R,\beta }^{\ast }\mathcal{T}{\psi }_{L,\alpha }-{\psi }_{L,\alpha }\mathcal{T}{\psi }_{R,\beta }^{\ast }\right]d^{3}r.$$

$$⟨{\psi }_{R,\beta }|\mathcal{H}-{\mathcal{H}}_L|{\psi }_{L,\alpha }⟩=-\mathrm{i}\hslash {\int }_{{\Sigma }_R}\mathbf{n}\cdot {\mathbf{J}}_{\beta \alpha }{d}^{2}r,$$

$${\mathbf{J}}_{\beta \alpha }=-\frac{\mathrm{i}\hslash }{2}{\hat{m}}^{-1}\left[{\psi }_{R,\beta }^{\ast }\nabla {\psi }_{L,\alpha }-{\psi }_{L,\alpha }\nabla {\psi }_{R,\beta }^{\ast }\right]$$

$${\psi }_{L,\alpha }(\mathbf{r})$$

$$E_{L,\alpha }$$

$$E_z$$

Full text

Introduction— 2D layered semiconducting materials, such as transition metal dichalcogenides (TMDs) and black phosphorus (BP), have many interesting electrical and optical properties Novoselov et al. (2005); Wang et al. (2012); Xu et al. (2014); Fiori et al. (2014); Low et al. (2016); Koppens et al. (2014); Sun et al. (2016), but tend to form Schottky barriers (SB) at the interfaces with metal contacts, resulting in a large contact resistance that severely degrades the device performance Das and Appenzeller (2013); Du et al. (2014); Haratipour et al. (2015).

Thermionic emission Sze and Ng (2006) is commonly assumed when extracting the SB height from temperature-dependent current measurements of field-effect transistors (FETs). In Allain et al. (2015), it was pointed out that this procedure is correct only at the gate voltage corresponding to the flat-band condition. Considering an n-type device, for example, above the flat-band voltage, the conduction band edge in the channel is higher than at the interface with the contact, hence electrons traversing the channel see a larger barrier than the SB height. Below the flat-band voltage, tunneling starts to contribute and the thermionic emission theory loses validity. As a result, this can lead to unphysical negative SB heights Yu et al. (2014); Avsar et al. (2015). Furthermore, experiments show that, as opposed to the insulating behavior of a SB contact, the two-terminal resistance Liu et al. (2015) as well as contact resistance Cui et al. (2015) can decrease with decreasing temperature at high gate voltage. The origin of this metallic behavior is debated and not yet clarified Radisavljevic and Kis (2013); Allain et al. (2015).

Recently, a model for SB FETs has been proposed in Penumatcha et al. (2015) and applied to extract the SB height and bandgap of BP devices. This model assumes one-dimensional transport and a bias-independent SB height. However, in a typical geometry with a top contact to a multilayer 2D semiconductor as in the sketch of Fig. 1a, transport is inherently 2D. To be precise, due to quantization, the SB height $\Phi_{1}$ to a 2D semiconductor should be defined as the difference between the edge $E_{1}$ of the first energy subband in the semiconductor and the Fermi level of the metal $\mu$ (see the schematic band profile in Fig. 1b for an n-type device). Transport occurring at energies above (below) $E_{1}$ is generally referred to as “thermionic” (“tunneling”). However, in the presence of a back gate, the subband edge and thus the SB height are expected to be modulated by the vertical electric field. When $E_{1}$ is lower than the bulk band edge at the interface with the metal, even the electrons traversing the junction at energies above $E_{1}$ see a tunneling barrier in the vertical (i.e., $z$) direction and a new transport regime arises.

In this paper, we present a semi-analytical model of the contact resistance to multilayer 2D semiconductors in this “vertical tunneling” regime. The model is based on a triangular barrier approximation of the vertical potential profile in the semiconductor underneath the contact. 2D transport is separated into a sequence of two 1D mechanisms: (i) quantum tunneling through the SB at the metal-to-semiconductor interface, followed by (ii) semiclassical “diffusive” transport across the semiconductor (the source of scattering being the in- and out-tunneling across the SB). The model is benchmarked against numerical solutions of the 2D quantum transport problem and employed to study the dependence of the contact resistance on vertical electric field and temperature. We show that, when the SB height is sufficiently lowered by the vertical electric field, contact resistance shows a metallic behavior with temperature, as observed in experiments. The model predicts a smooth transition from a thermionic-like regime at low electric field, where the tunneling barrier is almost transparent, to a true vertical tunneling regime at high electric field. In the former case, the extraction method of the SB height based on thermionic emission theory can still be applied.

Model— We consider a single planar junction of length $L_{x}$ between a metal and a multi-layer 2D semiconductor with thickness $a$ (Fig. 1a). A vertical electric field is created inside the semiconductor by the presence of a back gate. Let $x$ and $z$ be the longitudinal and vertical directions, respectively. The device is uniform in the $y$ direction. We focus on the portion of the semiconductor covered by the metal (contact region) assuming that the uncovered part (channel) simply acts as a “reflectionless” contact or semi-infinite lead Datta (1997). A current can flow as indicated in Fig. 1a.

We neglect hole transport and discuss only injection of electrons from the metal to the conduction band of the semiconductor. A simple single-valley effective mass Hamiltonian, with the same values of the effective masses $m_{x,y,z}$ in the metal and in the semiconductor, is adopted. Such effective mass model has been shown to provide an accurate description of the out-of-plane quantization in multilayer BP Zhang et al. (2017). For the case of multilayer TMDs, one should employ a more fundamental tight-binding model Kang et al. (2016). Whereas the form of the equations will be different, the general trends are not expected to change significantly. The bulk band edge profile is assumed to be uniform along the $x$ direction and is approximated with a triangular barrier along the $z$ direction, as shown in Fig. 1b, where the value of the band edge in the metal $V_{0}$ is chosen low enough so that it provides significant density of states at the Fermi level. Within this non-self-consistent approximation, which is valid at low carrier concentration, the magnitude of the vertical electric field $F$ in the semiconductor is simply proportional to the voltage $V_{G}$ applied between the back gate and the top metal:

[TABLE]

where the workfunctions of the two metals are taken to be equal, $t_{\mathrm{ox}}$ is the back oxide thickness, and $\epsilon_{s}$ and $\epsilon_{\mathrm{ox}}$ are the dielectric constants of the semiconductor and oxide, respectively.

Due to vertical confinement, the energy spectrum in the semiconductor splits into a set of discrete 2D subbands. Within a triangular well approximation, the subband edges $E_{i}$ ($i$ positive integer) can be computed as Miller (2008)

[TABLE]

where the energy reference is taken at the top of the barrier in Fig. 1b, $q$ is the elementary electric charge, the wavevector $k_{F}$ is defined as

[TABLE]

and $\zeta_{i}$ are the zeros of Airy’s function, i.e., $\mathrm{Ai}(\zeta_{i})=0$, which can be approximated as Abramowitz and Stegun (1964)

[TABLE]

We limit the discussion to the case $E_{i}<0$. Indeed, the subband description looses validity above the barrier.

We assume that transport within the semiconductor can be described by a set of decoupled 1D Boltzmann’s transport equations Rudan (2015), one for each subband, where the tunneling from the metal to the semiconductor and vice versa is included as a scattering mechanism. The corresponding relaxation time $\tau_{i}$, or inverse of the probability rate that an electron originally in the $k$-space state $(k_{x},k_{y})$ of the $i$-th subband tunnels into the metal, is computed according to Bardeen’s transfer Hamiltonian theory Bardeen (1961); Harrison (1961); Duke (1969), which has been recently applied to describe tunneling in vertical heterostructures of 2D materials Feenstra et al. (2012); Britnell et al. (2013) and electron-hole bilayer tunnel FETs Alper et al. (2013); Agarwal et al. (2014). In the limit of large $L_{x}$, we get

[TABLE]

where $h$ is Planck’s constant, $\hbar=h/(2\pi)$, $\gamma_{0}$ is defined as

[TABLE]

and $\mathrm{Ai}^{\prime}(\zeta_{i})$ is the derivative of Airy’s function evaluated at $\zeta_{i}$, which can be approximated as Abramowitz and Stegun (1964)

[TABLE]

Note that $\tau_{i}$ is independent of both $k_{x}$ and $k_{y}$. The tunneling current is computed from the $x$-dependent distribution function of each subband, which is obtained by solving Boltzmann’s transport equation with $\tau_{i}$ as the scattering relaxation time and with appropriate boundary conditions. In particular, we assume that the electrons are backscattered at the left end of the contact region. Differentiating the tunneling current with respect to the applied bias $V_{D}$ and taking the limit $V_{D}\ll k_{B}T/q$ ($k_{B}$ is Boltzmann’s constant and $T$ the temperature) gives us the low bias conductance $G$ or inverse of contact resistance (per unit width). We obtain the semi-analytical expression

[TABLE]

where $\varepsilon$ is the total energy for electrons with $k_{y}=0$ and $\mathcal{F}_{-1/2}$ the Fermi-Dirac integral of order $-\tfrac{1}{2}$. The total transmission function $\overline{T}$ is defined as

[TABLE]

with the trasmission probability $T_{i}$ of each subband given by

[TABLE]

where $\lambda_{i}=|v_{x}|\tau_{i}$ is the mean free path related to tunneling and $|v_{x}|=\sqrt{2(\varepsilon-E_{i})/m_{x}}$ is the longitudinal carrier velocity. In the Supporting information, we provide a detailed derivation of the model.

Results— In Fig. 2 we plot the tunneling rate $1/\tau_{i}$ of the first two subbands as a function of electric field and semiconductor thickness. The predicted tunneling rate goes to zero at small $F$ or $a$ because the Bardeen model does not account for the above-the-barrier regime at $E_{i}>0$. At high electric field or large semiconductor thickness, the exponential term in (5) is dominating. In this regime, an increase of $F$ or $a$ results in a decrease of the scattering rate $1/\tau_{i}$. This can be understood by noting that, since both $k_{x}$ and $k_{y}$ are conserved in the tunneling process, an electron can tunnel from the metal to the semiconductor only if its vertical energy is equal to $E_{i}$. However, according to (2), $E_{i}$ shifts to lower energies with increasing $F$ or $a$. Because of that shift, the tunneling distance, which is equal to $|E_{i}|/(qF)$ at the vertical energy $E_{i}$, becomes longer as $F$ or $a$ increase.

In order to benchmark the proposed model, we solve numerically the 2D Schrödinger equation with open boundary conditions using the Green function (GF) method Datta (1997) and assuming the same non-self-consistent triangular potential profile as in Fig. 1b. Details on the GF calculation can be found in the Supporting information. Fig. 3 shows the plot of the total transmission function $\overline{T}(\varepsilon)$ computed with the analytical expression in (11)–(12) and with GF for different sets of parameter values. The two models are in good general agreement. The transmission function increases by one at each energy corresponding to a subband edge $E_{i}$, indicating a resonant tunneling regime, and shows a decaying behavior between two successive subband edges. Indeed, different energies correspond to different $k_{x}$ states. Since the length of the contact $L_{x}$ is finite and the transfer length or average distance traveled by an electron in the semiconductor before tunneling into the metal is equal to the mean free path $\lambda_{i}=|v_{x}|\tau_{i}$, the probability of escaping into the metal is larger for the states closer to the subband edge which have smaller velocity. As shown in Fig. 3, increasing the contact length tends to raise the transmission probability of each single subband to unity because the ratio $L_{x}/\lambda_{i}$ between the contact length and the average distance before tunneling increases, which means that the electrons have more chances to enter the contact. The shift of $E_{i}$ to lower energy as $a$ or $F$ increase is clearly seen from the shift of the transmission peaks in Fig. 3b and c, respectively. This is accompanied by a narrowing of the peaks, which is related to the decrease of $1/\tau_{i}$ discussed above.

It should be noted that, with reference to the generic subband of index $i$, our model predicts no vertical tunneling contribution at energies below $E_{i}$ (see Eq. 12). This process could be possible in the case of a realistic band bending between the channel and the contact region. However, since tunneling is decreasing exponentially with the tunneling distance, the effect would be concentrated at the contact edge. Therefore, if $E_{i}$ is sufficiently close to the metal Fermi level $\mu$, the contribution from energies below $E_{i}$ (lateral tunneling) is negligible compared to energies above $E_{i}$ (vertical tunneling) because of the large surface-to-edge ratio of the contact.

Fig. 4a plots $G$, obtained by numerically computing the energy integral in (9) and resolved for the first two subbands, as a function of $F$ and temperature $T$. It can be seen that $G$ is a non-monotonic function of the electric field. This is consistent with our previous observations: the transmission probability reduces with increasing $F$ because the tunneling distance increases. Fig. 4a shows that, at high electric field before the second subband starts to contribute significantly, the derivative $\partial G/\partial T$ is negative, i.e., contact resistance decreases with decreasing temperature. This has to do with the factor $-\partial F_{0}/\partial\varepsilon$, which is plotted in Fig. 4b. It can be proved that its derivative with respect to temperature changes sign at the energy $\varepsilon_{0}\approx\mu+0.857k_{B}T$, which has only a weak temperature dependence ($\varepsilon_{0}=-0.243$ and $-0.228$ eV at $T=100$ and $300$ K, respectively) and appears as a crossover point in Fig. 4b. As the transmission function shifts to lower energy with increasing $F$ (compare the plots at $F=0.08$ and $0.1$ V/nm in Fig. 4b), more contribution to the integral in (9) comes from the energy range where $\partial(-\partial F_{0}/\partial\varepsilon)/\partial T<0$ and eventually leads to $\partial G/\partial T<0$.

The model in (9) allows for an analytical solution in two limiting cases. In order to simplify the discussion, we assume that only the first subband contributes to transport. Note that the energy range relevant for transport goes from $E_{1}$ to few $k_{B}T$’s above $E_{1}$ or $\mu$, whichever is maximum. If $L_{x}\gg\lambda_{1}$ in this energy range, it follows that $T_{1}(\varepsilon)\approx 1$ (i.e., an almost transparent barrier) and (9) simplifies to

[TABLE]

If, in addition, $\Phi_{1}=E_{1}-\mu\gg k_{B}T$, then (13) further reduces to

[TABLE]

which is the expression of the thermionic emission theory for a 2D system Anwar et al. (1999). The prefactor is $T^{1/2}$ instead of $T^{3/2}$ because we are considering the low bias limit $V_{D}\ll k_{B}T/q$. Expression (14) implies that the SB height $\Phi_{1}$ can be extracted from the slope of $\ln(G/T^{1/2})$ vs $1/T$. On the other hand, if $L_{x}\ll\lambda_{1}$ in most of the energy window for transport, one can derive (see Supporting information)

[TABLE]

where $f_{0}(E)=\{\exp[(E-\mu)/(k_{B}T)]+1\}^{-1}$ is the Fermi-Dirac function. Fig. 4c compares the Arrhenius plot of $G/T^{1/2}$ computed with the rigorous model in (9) and the approximated expressions in (13) and (15). $E_{1}$ is calculated according to (2) in all three cases. Similar Arrhenius plots are commonly used to extract the SB height in experiments Yu et al. (2014); Anugrah et al. (2015); Avsar et al. (2015); Cui et al. (2015). For the chosen set of parameter values, approximation (13) is valid up to $F\approx 0.08$ V/nm. At higher electric field, the transmission function becomes increasingly peaked around the subband edge (see Fig. 4b) and (15) becomes a better approximation. A positive slope, or metallic behavior, is predicted at high electric field similar to what has been reported in experiments Yu et al. (2014); Avsar et al. (2015). In Fig. 4d, we plot the SB height obtained by fitting the data of the rigorous model in (9) with the thermionic expression (14), compared with the actual value of $\Phi_{1}=E_{1}-\mu$ from (2). It is seen that the extraction method based on the thermionic emission theory can provide good results at low electric field values, where the tunneling barrier is almost transparent. In the high-field regime, a fitting based on (15) would provide a more physical result.

We conclude by noting that the model presented in this work can be easily extended to account for a finite carrier mobility in the semiconductor by introducing an additional relaxation time $\tau_{s}$ (and a corresponding mean free path $\lambda_{s}=|v_{x}|\tau_{s}$) related to elastic scattering. The main effect of scattering would be, for each subband, a shorter transfer length and a transmission probability that does saturate to unity in the limit of a long contact length. In the regime when only one subband is populated, the model could also be extended to include self-consistent electrostatics using the variational approach in Stern (1972).

Conclusions— In summary, we have demonstrated that the metallic behavior of the contact resistance observed in recent experiments can be explained by taking into account the modulation of the vertical tunneling due to the SB lowering with increasing electric field in back-gated devices. To the best of our knowledge, this transport regime has not been discussed before. The model also suggests a non-monotonic behavior of the contact resistance with respect to vertical electric field and semiconductor thickness. Our semi-analytical model provides a reasonable description of contact resistance in 2D semiconductors and could be useful for contact engineering in future 2D electronics.

I Model of vertical tunneling

We compute the tunneling rate in the limit of an infinite contact length and assume that the electric potential does not depend on the longitudinal position $x$, which implies translational invariance along $x$. Let $z=0$ be the vertical position of the metal-to-semiconductor interface. We consider a simple effective mass Hamiltonian

[TABLE]

where the values of the effective masses $m_{x,y,z}$ are taken to be the same in the metal and in the semiconductor and the conduction band edge profile $E_{c}(z)$ is modeled as a triangular barrier ($F>0$ is the magnitude of the vertical electric field, see Fig. S1a):

[TABLE]

Note that this is different from the Fowler-Nordheim field-emission problem Fowler and Nordheim (1928) because of the presence of the hard wall at $z=a$. Suppose that an electron is launched from $z<0$ towards the interface. The electron wavefunction will be totally reflected at $z=a$, resulting in a reflection coefficient, measured as the ratio between the probability currents of reflected and incident waves, identically equal to one at all energies. Does it mean that the tunneling probability is zero? The Bardeen Transfer Hamiltonian method Bardeen (1961); Harrison (1961); Duke (1969) provides a way to overcome this difficulty: the tunneling process across the barrier is thought of as a scattering event between states localized on different sides of the junction and the corresponding transition probability is computed through Oppenheimer’s version of time-dependent perturbation theory Oppenheimer (1928).

More precisely, $\mathcal{H}$ is taken as the perturbed Hamiltonian acting in the time interval $0<t<t_{P}$. For $t<0$, the unperturbed Hamiltonian must be identified with an Hamiltonian $\mathcal{H}_{L}$ that approximates well the true Hamiltonian $\mathcal{H}$ on the metal side of the junction but whose eigenfunctions decay in the semiconductor. We take $\mathcal{H}_{L}=\mathcal{T}+E_{cL}(z)$ with $E_{cL}(z)$ a potential step (Fig. S1b):

[TABLE]

For $t>t_{P}$, one must choose a different unperturbed Hamiltonian $\mathcal{H}_{R}$ which, conversely, approximates well $\mathcal{H}$ on the semiconductor side of the junction but whose eigenfunctions decay in the metal. We take $\mathcal{H}_{R}=\mathcal{T}+E_{cR}(z)$ with $E_{cR}(z)$ a triangular well (Fig. S1c):

[TABLE]

It is assumed that prior to the perturbation the electron wavefunction coincides with an eigenfunction $\psi_{L,\alpha}$ of $\mathcal{H}_{L}$. Since this is not an eigenstate of $\mathcal{H}$, the electron wavefunction will evolve during the time interval $0<t<t_{P}$ according to the time-dependent Schrödinger equation. If $t_{P}$ is sufficiently large, the probability that the electron is subsequently found in the eigenstate $\psi_{R,\beta}$ of $\mathcal{H}_{R}$ at $t>t_{P}$, is, to first order and per unit $t_{P}$,

[TABLE]

where $E_{L,\alpha}$ and $E_{R,\beta}$ are the eigenvalues corresponding to the initial state $\psi_{L,\alpha}$ and final state $\psi_{R,\beta}$, respectively, and $\delta$ is Dirac’s delta function Oppenheimer (1928). The conservation of energy is related to the perturbation being constant in time. In (S5), it is assumed that each set of eigenfunctions $\psi_{L,\alpha}$ and $\psi_{R,\beta}$ is discrete and orthonormal. We consider a rectangular domain with finite sides $L_{x}$ and $L_{y}$ in the plane parallel to the junction and prescribe the periodic boundary conditions $\psi_{L,\alpha/R,\beta}(x=0,y,z)=\psi_{L,\alpha/R,\beta}(x=L_{x},y,z)$, $\psi_{L,\alpha/R,\beta}(x,y=0,z)=\psi_{L,\alpha/R,\beta}(x,y=L_{y},z)$. In addition, we consider a finite length $L_{z}$ of the metal region in the $z$ direction with the hard-wall boundary condition $\psi_{L,\alpha}(x,y,z=-L_{z})=0$. This way, both energy spectra are discrete and the corresponding eigenfunctions normalizable. Later, we will take the limit as $L_{x},L_{y},L_{z}$ go to infinite in order to recover the continuous case. It should be noted that, contrary to the standard time-dependent perturbation theory Rudan (2015), the matrix element in (S5) must be computed between eigenstates of different Hamiltonians. Also, the perturbation Hamiltonian must be evaluated with respect to the initial Hamiltonian. The validity of (S5) rests on the assumption that the two sets of eigenstates of $\mathcal{H}_{L}$ and $\mathcal{H}_{L}$ are “almost orthogonal” to each other, in particular that $\langle\psi_{R,\beta}|\psi_{L,\alpha}\rangle\ll 1$ Oppenheimer (1928). Similary, one has for the probability rate of the inverse transition

[TABLE]

where the last equality follows from the delta function and the Hermiticity of the various Hamiltonians.

From (S2)-(S4) we get

[TABLE]

which means that $\mathcal{H}$ satisfies the separability property of Bardeen’s model Hamiltonian, i.e., that $\mathcal{H}-\mathcal{H}_{L}\neq 0$ only in regions of space where $\mathcal{H}-\mathcal{H}_{R}\equiv 0$ Bardeen (1961) 111Bardeen’s theory is often introduced by writing the Hamiltonian in the form $\mathcal{H}=\mathcal{H}_{L}+\mathcal{H}_{R}+\mathcal{H}_{T}$, with $\mathcal{H}_{T}$ being the “transfer” Hamiltonian, despite the fact that Bardeen himself did not make use of such decomposition in its original paper Bardeen (1961). In our case, we get from (S7)

\mathcal{H}_{T}=\left\{\begin{array}[]{ll}-\mathcal{H}_{R},&\quad z<0\\ -\mathcal{H}_{L},&\quad z>0\end{array}\right.

(S8)

but this Hamiltonian does not correspond to either of the perturbation Hamiltonians that appear in (S5) or (S6). See also the discussion in (Duke, 1969).. Since $\mathcal{H}-\mathcal{H}_{L}\equiv 0$ for $z<0$, the matrix element in (S5) can be written as

[TABLE]

where $\mathbf{r}=(x,y,z)$ and $\Omega_{R}=\{\mathbf{r};0<x<L_{x},0<y<L_{y},0<z<\infty\}$. Noting that $\mathcal{H}-\mathcal{H}_{R}\equiv 0$ for $z>0$, we can get the symmetric expression

[TABLE]

Due to the delta function in (S5), we are only interested in the case $E_{L,\alpha}=E_{R,\beta}$, for which

[TABLE]

Applying Green’s theorem, we finally get

[TABLE]

where $\Sigma_{R}$ is the surface of $\Omega_{R}$, $\mathbf{n}$ is the unit vector normal to $\Sigma_{R}$ pointing in the outward direction, and

[TABLE]

is the matrix element of the probability current density operator between the states $\psi_{L,\alpha}$ and $\psi_{R,\beta}$.

Let us now start to pick up all the ingredients that we need to calculate (S12). The eigenfunctions and corresponding eigenvalues of $\mathcal{H}_{L}$ are Rudan (2015)

[TABLE]

where $\mathbf{k}_{L}=(k_{xL},k_{yL},k_{z})$ and it is assumed that $E_{z}<0$ 222We limit the discussion to the case $E_{z}<0$ because (S5) looses validity if $\langle\psi_{R,\beta}|\psi_{L,\alpha}\rangle$ is not a small number.. Because of the periodic boundary conditions, the transverse components of the wavevector are quantized as $k_{xL}=2\pi l/L_{x}$, $k_{yL}=2\pi m/L_{y}$ ($l,m$ integers). As for $k_{z}$, the allowed values are the roots of the transcendental equations $k_{z}L_{z}-\pi(n-1/2)=\varphi$ in the interval $0<k_{z}<\sqrt{-2m_{z}V_{0}}/\hbar$. For large $L_{z}$, we have $k_{z}\approx\pi n/L_{z}$ ($n$ positive integer). The constant $b$ can be obtained from the normalization condition

[TABLE]

where only the leading term in $L_{z}$ has been kept. Thus, up to an unimportant phase, $b=\sqrt{2/L_{z}}$.

The solutions of the eigenvalue problem of $\mathcal{H}_{R}$ are Miller (2008)

[TABLE]

where $\mathbf{k}_{\parallel R}=(k_{xR},k_{yR})$, $\mathrm{Ai}$ is Airy’s function and $\zeta_{i}$ are its zeros, which can be approximated as ($i$ positive integer) Abramowitz and Stegun (1964)

[TABLE]

The constant $c$ can be obtained from the normalization condition

[TABLE]

where the last integral can be evaluated using integration by parts and the fact that $\mathrm{Ai}$ is a solution of Airy’s equation $\mathrm{Ai}^{\prime\prime}=\zeta\mathrm{Ai}$ (the prime symbol indicates derivative with respect to $\zeta$):

[TABLE]

Therefore, $c=\sqrt{k_{F}}/|\mathrm{Ai}^{\prime}(\zeta_{i})|$, in which we can use the approximated expression Abramowitz and Stegun (1964)

[TABLE]

The surface $\Sigma_{R}$ in (S12) is made up of six faces. By inserting (S17) and (S25) into (S12), it can be shown that the integrals over the two faces at $x=0$ and $x=L_{x}$, as well as the integrals over the two faces at $y=0$ and $y=L_{y}$, cancel out each other exactly 333For example, $J_{x,\beta\alpha}(x=L_{x},y,z)-J_{x,\beta\alpha}(x=0,y,z)\propto e^{\mathrm{i}(k_{xL}-k_{xR})L_{x}}-1=0$ because of the periodic boundary conditions.. The integral over the face at $z=\infty$ is also zero because $\psi_{L,\alpha}$ is vanishingly small. We are only left with the integral over the face at $z=0$:

[TABLE]

where

[TABLE]

and $\delta$ is Kronecker’s delta function. Conservation of transverse momentum is a consequence of the translational symmetry along $x$ and $y$. Combined with energy conservation, it implies that $E_{z}=E_{i}$ and thus $\zeta_{0}=(\kappa/k_{F})^{2}$. Assuming $\zeta_{0}\gg 1$ (which is consistent with $\langle\psi_{R,\beta}|\psi_{L,\alpha}\rangle\ll 1$), we can substitute in (S34) the asymptotic expressions

[TABLE]

where $\gamma_{0}=(2/3)\zeta_{0}^{3/2}$ Abramowitz and Stegun (1964), to get

[TABLE]

Finally, plugging (S38) into (S5), we obtain

[TABLE]

II Model of longitudinal diffusion

Suppose that the states in the metal ($L$) are populated according to a Fermi-Dirac distribution with Fermi level $\mu_{L}$:

[TABLE]

with $k_{B}$ Boltzmann’s constant and $T$ the temperature. As for the semiconductor, we cannot assume that the states are in equilibrium because a current has to flow in the $x$ direction as shown in Fig. 1 of the main text. In order to compute the population of such states, we assume semiclassical diffusive transport and make use of Boltzmann’s transport equation Rudan (2015).

Let $f_{i}(x,\mathbf{k}_{\parallel R})$ be the distribution function in the four-dimensional phase space associated with the $i$-th subband 444$f_{i}$ is independent of $y$ because of the translational symmetry along $y$.. Under the assumption that the electric potential is uniform along $x$, Boltzmann’s equation reads

[TABLE]

where $v_{x}=\hbar k_{xR}/m_{x}$ is the longitudinal carrier velocity and the transitions from the metal to the semiconductor and vice versa due to vertical tunneling are included through a collision term $C$ 555Other types of scattering, which could be responsible for a finite carrier mobility in the semiconductor, are here neglected.:

[TABLE]

Note that expression (II) takes into account Pauli’s exclusion principle. Using (S39), we get

[TABLE]

where the relaxation time $\tau_{i}$ is defined as

[TABLE]

Going to the limit of large $L_{z}$, we can replace

[TABLE]

so that

[TABLE]

and, with the change of variables $k_{z}\rightarrow E_{z}$,

[TABLE]

Besides the subband index $i$, the relaxation time depends on the parameters $m_{z}$, $F$, $a$, and $V_{0}$. It can be shown that, replacing the eigenfunctions (S17) of $\mathcal{H}_{L}$ by their WKB Messiah (1961) approximation

[TABLE]

the expression of $\tau_{i}$ simplifies to

[TABLE]

where $V_{0}$ appears only as an energy cut-off. This last formulation, which does not depend on the precise bandstructure of the metal, could be useful for treating injection from the metal to the valence band of the semiconductor.

Let $f_{i}^{+}$ and $f_{i}^{-}$ denote the distribution functions of right-going and left-going states, respectively, i.e., $f_{i}^{\pm}(x,k_{xR},k_{yR})=f_{i}(x,\pm k_{xR},k_{yR})$ with $k_{xR}>0$. We impose the boundary conditions

[TABLE]

where $f_{R}$ is a Fermi-Dirac function similar to (S40) with $\mu_{L}$ replaced by $\mu_{R}$ (see Fig. S2). The latter condition makes sure that the longitudinal current vanishes at $x=0$. The Boltzmann equations (S41) for different wavevectors are independent of each other expect for $\pm k_{xR}$. The solutions are

[TABLE]

The current (per unit width) from the semiconductor to the metal can be obtained by summing the net flux $|v_{x}|\left(f_{i}^{+}-f_{i}^{-}\right)$ at $x=L_{x}$ over all semiconductor states per unit area, multiplying by $2$ for spin degeneracy and multiplying by the electronic charge $q$:

[TABLE]

Going to the limit of large $L_{x},L_{y}$, we can replace

[TABLE]

to get

[TABLE]

Finally, with the change of variables $k_{xR}\rightarrow\varepsilon=\hbar^{2}k_{xR}^{2}/(2m_{x})+E_{i}$, we obtain the Landauer formula Datta (1997)

[TABLE]

where the longitudinal velocity must be computed as $|v_{x}|=\sqrt{2(\varepsilon-E_{i})/m_{x}}$, the supply function $F_{L/R}$ is defined as

[TABLE]

with $\mathcal{F}_{-1/2}$ the Fermi-Dirac integral of order $-\tfrac{1}{2}$ Rudan (2015), $\overline{T}$ is the transmission function

[TABLE]

with the transmission probability $T_{i}$ given by

[TABLE]

and $\lambda_{i}=|v_{x}|\tau_{i}$. For $\varepsilon>E_{i}$, we can have the asymptotic behaviors

[TABLE]

Approximation (S62) holds, in particular, as $\varepsilon\rightarrow E_{i}^{+}$ (resonant tunneling). Note also that, when (S63) is satisfied, $T_{i}$ decays as a function of energy as $1/\sqrt{\varepsilon}$.

The low bias conductance per unit width $G$ can be evaluated from (II). Let $\mu_{L}=\mu$ and $\mu_{R}=\mu-qV_{D}$. We have

[TABLE]

where $F_{0}(\varepsilon)$ is given by (S59) with $\mu_{L/R}$ replaced by $\mu$.

Assume for simplicity that only the first subband contributes to transport. In the two limiting cases when either (S62) or (S63) are satisfied over the whole energy range of interest for transport (from $E_{1}$ to few $k_{B}T$’s above $\max\{E_{1},\mu\}$), it is possible to derive analytical expressions for $G$. If $L_{x}\gg\lambda_{1}$ in this energy range, it follows immediately from (S64)

[TABLE]

To work out the expression of $G$ in the other limiting case when $L_{x}\ll\lambda_{1}$ in most of the energy window for transport666The inequality does not hold for energies close to $E_{1}$ but their contribution becomes increasingly smaller as $\tau_{1}$ increases., it is convenient to go back to the double-integral formulation of the tunneling current in (S57) and do the change of variables $k_{xR}\rightarrow E_{R,\beta}$:

[TABLE]

where $k_{y\mathrm{max}}=\sqrt{2m_{y}(E_{R,\beta}-E_{1})}/\hbar$ and $|v_{x}|=\sqrt{2[E_{R,\beta}-\hbar^{2}k_{yR}^{2}/(2m_{y})-E_{1}]/m_{x}}$. The integral over the transverse wavevector can be easily computed with the change of variables $k_{yR}\rightarrow\arcsin(k_{yR}/k_{y\mathrm{max}})$ to give

[TABLE]

Letting $\mu_{L}=\mu$ and $\mu_{R}=\mu-qV_{D}$, we finally get

[TABLE]

where $f_{0}$ is a Fermi-Dirac function similar to (S40) with $\mu_{L}$ replaced by $\mu$.

III Green function algorithm

We discretize the Hamiltonian in (S1) using finite differences on a two-dimensional rectangular grid (Fig. S3a). The same linear potential profile as in (S2) is assumed. The transmission function is computed as Datta (1997)

[TABLE]

where the symbol $\mathrm{tr}$ indicates the trace, $G^{r}$ is the retarded Green function, $G^{a}=G^{r\dagger}$, $\Gamma^{L/R}=\mathrm{i}(\Sigma^{r,L/R}-{\Sigma^{a,L/R}})$, with $\Sigma^{r,L/R}$ the retarded self-energy representing the renormalization of the Hamiltonian of the channel region (black dots in Fig. S3a) due to the presence of the semi-infinite left/right lead, and $\Sigma^{a,L/R}=(\Sigma^{r,L/R})^{\dagger}$. The channel region is partitioned into layers as shown in Fig. S3b. Using matrix block notation and noting that the only non-null block of $\Sigma^{r,L}$ is $\Sigma^{r,L}_{0,0}$ and the only non-null block of $\Sigma^{r,R}$ is $\Sigma^{r,R}_{N_{x}+1,N_{x}+1}$, (S69) can be rewritten as

[TABLE]

The self-energy of the left lead is computed analytically using the prescription given in Datta (1997):

[TABLE]

where $t_{x}=\hbar^{2}/(2m_{x}\Delta_{x}^{2})$ and similarly for $t_{z}$. The self-energy of the right lead is obtained numerically using a well-known iterative algorithm Sancho et al. (1985). The matrix block $G^{r}_{0,N_{x}+1}$ is computed through a combination of the recursive and decimation algorithms Low and Appenzeller (2009), modified so as to treat a non-tridiagonal-block Hamiltonian matrix. Let $A=\varepsilon I-H_{C}-\Sigma^{r,L}-\Sigma^{r,R}$, where $H_{C}$ is the Hamiltonian matrix of the channel region alone, and define $\delta_{1}^{(0)}=A_{0,0}$, $\delta_{2}^{(0)}=A_{1,1}$, $\alpha^{(0)}=A_{0,1}$, $\beta^{(0)}=A_{1,0}$. The algorithm consists in eliminating the layers from $1$ to $N_{x}$ with the formulas

[TABLE]

for $n=1,\ldots,N_{x}$, where it is understood that $A_{0,N_{x}+1}=A_{N_{x}+1,0}^{\dagger}=0$. At the end, the required matrix block of the Green function can be obtained as

[TABLE]