Semi-analytical model of the contact resistance in two-dimensional semiconductors
Roberto Grassi Yanqing Wu Steven J. Koester Tony Low





rgrassi@umn.edu
Department of Electrical and Computer Engineering, University of Minnesota, 200 Union St. SE, Minneapolis, MN 55455, USA

Wuhan National High Magnetic Field Center and School of Optical and Electronic Information, Huazhong University of Science and Technology, Wuhan 430074, China


Department of Electrical and Computer Engineering, University of Minnesota, 200 Union St. SE, Minneapolis, MN 55455, USA



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.


Semi-analytical model of the contact resistance in two-dimensional semiconductors
    Tony Low
    December 30, 2023
=================================================================================



Introduction— 2D layered semiconducting materials, such as transition metal dichalcogenides (TMDs) and black phosphorus (BP), have many interesting electrical and optical properties <cit.>, 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 <cit.>.

Thermionic emission <cit.> is commonly assumed when extracting the SB height from temperature-dependent current measurements of field-effect transistors (FETs). In <cit.>, 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 <cit.>. Furthermore, experiments show that, as opposed to the insulating behavior of a SB contact, the two-terminal resistance <cit.> as well as contact resistance <cit.> can decrease with decreasing temperature at high gate voltage. The origin of this metallic behavior is debated and not yet clarified <cit.>.

Recently, a model for SB FETs has been proposed in <cit.> 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. <ref>a, transport is inherently 2D. To be precise, due to quantization, the SB height Φ_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 μ (see the schematic band profile in Fig. <ref>b 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. <ref>a). 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 <cit.>. A current can flow as indicated in Fig. <ref>a.

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 <cit.>. For the case of multilayer TMDs, one should employ a more fundamental tight-binding model <cit.>. 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. <ref>b, 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:

    F = V_G/a + (ϵ_s/ϵ_ox)t_ox  ,

where the workfunctions of the two metals are taken to be equal, t_ox is the back oxide thickness, and ϵ_s and ϵ_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 <cit.>

    E_i = - qF ( a - | ζ_i |/k_F) ,

where the energy reference is taken at the top of the barrier in Fig. <ref>b, q is the elementary electric charge, the wavevector k_F is defined as

    k_F = ( 2 m_z qF/ħ^2)^1/3 ,

and ζ_i are the zeros of Airy's function, i.e., Ai(ζ_i)=0, which can be approximated as <cit.>

    ζ_i ≈ - [ 3 π/8 (4 i-1) ]^2/3 .

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 <cit.>, 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 τ_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 <cit.>, which has been recently applied to describe tunneling in vertical heterostructures of 2D materials <cit.> and electron-hole bilayer tunnel FETs <cit.>. In the limit of large L_x, we get

    1/τ_i = {[ 1/hħ^2/2 m_zk_F^2/Ai'^2(ζ_i)4 √(- E_i (E_i - V_0))/-V_0 e^-2γ_0,                                                    V_0 < E_i < 0;                                                               0,                                                        otherwise ].

where h is Planck's constant, ħ=h/(2π), γ_0 is defined as

    γ_0    = 2/3ζ_0^3/2  , 
    ζ_0    = - k_F E_i/qF = k_F a - | ζ_i |   ,

and Ai'(ζ_i) is the derivative of Airy's function evaluated at ζ_i, which can be approximated as <cit.>

    |Ai'(ζ_i)| ≈1/√(π)[ 3 π/8 (4 i-1) ]^1/6 .

Note that τ_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 τ_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 ≪ 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

    G     = 2q^2/h∫_-∞^∞ d ε T(ε) (-∂ F_0/∂ε) , 
    
    F_0(ε)    = √(m_y k_B T/2 πħ^2)ℱ_-1/2(μ - ε/k_B T) ,

where ε is the total energy for electrons with k_y=0 and ℱ_-1/2 the Fermi-Dirac integral of order -12. The total transmission function T is defined as

    T(ε) = ∑_i T_i(ε)   ,

with the trasmission probability T_i of each subband given by

    T_i(ε) = {[                0,           ε < E_i; 1 - e^-2 L_x/λ_i,           ε > E_i ].   ,

where λ_i = |v_x| τ_i is the mean free path related to tunneling and |v_x| = √( 2 (ε-E_i)/m_x ) is the longitudinal carrier velocity. In the Supporting information, we provide a detailed derivation of the model.



Results— In Fig. <ref> we plot the tunneling rate 1/τ_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 (<ref>) is dominating. In this regime, an increase of F or a results in a decrease of the scattering rate 1/τ_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 (<ref>), 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 <cit.> and assuming the same non-self-consistent triangular potential profile as in Fig. <ref>b. Details on the GF calculation can be found in the Supporting information. Fig. <ref> shows the plot of the total transmission function T(ε) computed with the analytical expression in (<ref>)–(<ref>) 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 λ_i = |v_x| τ_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. <ref>, increasing the contact length tends to raise the transmission probability of each single subband to unity because the ratio L_x/λ_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. <ref>b and c, respectively. This is accompanied by a narrowing of the peaks, which is related to the decrease of 1/τ_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. <ref>). 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 μ, 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. <ref>a plots G, obtained by numerically computing the energy integral in (<ref>) 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. <ref>a shows that, at high electric field before the second subband starts to contribute significantly, the derivative ∂ G/∂ T is negative, i.e., contact resistance decreases with decreasing temperature. This has to do with the factor -∂ F_0/∂ε, which is plotted in Fig. <ref>b. It can be proved that its derivative with respect to temperature changes sign at the energy ε_0 ≈μ+0.857 k_B T, which has only a weak temperature dependence (ε_0=-0.243 and -0.228 eV at T=100 and 300 K, respectively) and appears as a crossover point in Fig. <ref>b. 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. <ref>b), more contribution to the integral in (<ref>) comes from the energy range where ∂(-∂ F_0/∂ε)/∂ T < 0 and eventually leads to ∂ G/∂ T<0.

The model in (<ref>) 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 μ, whichever is maximum. If L_x ≫λ_1 in this energy range, it follows that T_1(ε)≈ 1 (i.e., an almost transparent barrier) and (<ref>) simplifies to

    G ≈2q^2/h F_0(E_1)   .

If, in addition, Φ_1=E_1-μ≫ k_B T, then (<ref>) further reduces to

    G ∝  T^1/2exp( -Φ_1/k_B T) ,

which is the expression of the thermionic emission theory for a 2D system <cit.>. The prefactor is T^1/2 instead of T^3/2 because we are considering the low bias limit V_D ≪ k_B T/q. Expression (<ref>) implies that the SB height Φ_1 can be extracted from the slope of ln(G/T^1/2) vs 1/T. On the other hand, if L_x ≪λ_1 in most of the energy window for transport, one can derive (see Supporting information)

    G ≈2q^2/h√(m_x m_y)/ħL_x/τ_1 f_0(E_1)   ,

where f_0(E) = {exp[(E-μ)/(k_B T)]+1}^-1 is the Fermi-Dirac function. Fig. <ref>c compares the Arrhenius plot of G/T^1/2 computed with the rigorous model in (<ref>) and the approximated expressions in (<ref>) and (<ref>). E_1 is calculated according to (<ref>) in all three cases. Similar Arrhenius plots are commonly used to extract the SB height in experiments <cit.>. For the chosen set of parameter values, approximation (<ref>) is valid up to F≈ 0.08 V/nm. At higher electric field, the transmission function becomes increasingly peaked around the subband edge (see Fig. <ref>b) and (<ref>) becomes a better approximation. A positive slope, or metallic behavior, is predicted at high electric field similar to what has been reported in experiments <cit.>. In Fig. <ref>d, we plot the SB height obtained by fitting the data of the rigorous model in (<ref>) with the thermionic expression (<ref>), compared with the actual value of Φ_1=E_1-μ from (<ref>). 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 (<ref>) 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 τ_s (and a corresponding mean free path λ_s = |v_x| τ_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 <cit.>.

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.

34
natexlab#1#1bibnamefont#1#1bibfnamefont#1#1citenamefont#1#1url<#>1urlprefixURL 


[Novoselov et al.(2005)Novoselov, Jiang,
  Schedin, Booth, Khotkevich, Morozov, and Geim]novoselov2005two
authorK. Novoselov,
  authorD. Jiang,
  authorF. Schedin,
  authorT. Booth,
  authorV. Khotkevich,
  authorS. Morozov, and
  authorA. Geim,
  journalProceedings of the National Academy of Sciences of the
  United States of America volume102,
  pages10451 (year2005).

[Wang et al.(2012)Wang, Kalantar-Zadeh,
  Kis, Coleman, and Strano]wang2012electronics
authorQ. H. Wang,
  authorK. Kalantar-Zadeh,
  authorA. Kis,
  authorJ. N. Coleman,
  and authorM. S.
  Strano, journalNature nanotechnology
  volume7, pages699 (year2012).

[Xu et al.(2014)Xu, Yao, Xiao, and
  Heinz]xu2014spin
authorX. Xu,
  authorW. Yao,
  authorD. Xiao, and
  authorT. F. Heinz,
  journalNature Physics volume10,
  pages343 (year2014).

[Fiori et al.(2014)Fiori, Bonaccorso,
  Iannaccone, Palacios, Neumaier, Seabaugh, Banerjee, and
  Colombo]fiori2014electronics
authorG. Fiori,
  authorF. Bonaccorso,
  authorG. Iannaccone,
  authorT. Palacios,
  authorD. Neumaier,
  authorA. Seabaugh,
  authorS. K. Banerjee,
  and authorL. Colombo,
  journalNature nanotechnology volume9,
  pages768 (year2014).

[Low et al.(2016)Low, Chaves, Caldwell,
  Kumar, Fang, Avouris, Heinz, Guinea, Martin-Moreno, and
  Koppens]low2016polaritons
authorT. Low,
  authorA. Chaves,
  authorJ. D. Caldwell,
  authorA. Kumar,
  authorN. X. Fang,
  authorP. Avouris,
  authorT. F. Heinz,
  authorF. Guinea,
  authorL. Martin-Moreno,
  and authorF. Koppens,
  journalNature Materials  (year2016).

[Koppens et al.(2014)Koppens, Mueller,
  Avouris, Ferrari, Vitiello, and Polini]koppens2014photodetectors
authorF. Koppens,
  authorT. Mueller,
  authorP. Avouris,
  authorA. Ferrari,
  authorM. Vitiello, and
  authorM. Polini,
  journalNature nanotechnology volume9,
  pages780 (year2014).

[Sun et al.(2016)Sun, Martinez, and
  Wang]sun2016optical
authorZ. Sun,
  authorA. Martinez, and
  authorF. Wang,
  journalNature Photonics volume10,
  pages227 (year2016).

[Das and Appenzeller(2013)]das2013does
authorS. Das and
  authorJ. Appenzeller,
  journalNano letters volume13,
  pages3396 (year2013).

[Du et al.(2014)Du, Liu, Deng, and
  Ye]du2014device
authorY. Du,
  authorH. Liu,
  authorY. Deng, and
  authorP. D. Ye,
  journalACS nano volume8,
  pages10035 (year2014).

[Haratipour et al.(2015)Haratipour,
  Robbins, and Koester]haratipour2015black
authorN. Haratipour,
  authorM. C. Robbins,
  and authorS. J.
  Koester, journalElectron Device Letters, IEEE
  volume36, pages411 (year2015).

[Sze and Ng(2006)]sze2006physics
authorS. M. Sze and
  authorK. K. Ng,
  titlePhysics of semiconductor devices
  (publisherJohn wiley & sons, year2006).

[Allain et al.(2015)Allain, Kang,
  Banerjee, and Kis]allain2015electrical
authorA. Allain,
  authorJ. Kang,
  authorK. Banerjee, and
  authorA. Kis,
  journalNature Materials volume14,
  pages1195 (year2015).

[Yu et al.(2014)Yu, Lee, Ling, Santos,
  Shin, Lin, Dubey, Kaxiras, Kong, Wang et al.]yu2014graphene
authorL. Yu,
  authorY.-H. Lee,
  authorX. Ling,
  authorE. J. Santos,
  authorY. C. Shin,
  authorY. Lin,
  authorM. Dubey,
  authorE. Kaxiras,
  authorJ. Kong,
  authorH. Wang, et al.,
  journalNano letters volume14,
  pages3055 (year2014).

[Avsar et al.(2015)Avsar, Vera-Marun,
  Tan, Watanabe, Taniguchi, Castro Neto, and Özyilmaz]avsar2015air
authorA. Avsar,
  authorI. J. Vera-Marun,
  authorJ. Y. Tan,
  authorK. Watanabe,
  authorT. Taniguchi,
  authorA. H. Castro Neto,
  and
  authorB. Özyilmaz,
  journalACS nano volume9,
  pages4138 (year2015).

[Liu et al.(2015)Liu, Wu, Cheng, Yang,
  Zhu, He, Ding, Li, Guo, Weiss et al.]liu2015toward
authorY. Liu,
  authorH. Wu,
  authorH.-C. Cheng,
  authorS. Yang,
  authorE. Zhu,
  authorQ. He,
  authorM. Ding,
  authorD. Li,
  authorJ. Guo,
  authorN. O. Weiss,
  et al., journalNano letters
  volume15, pages3030 (year2015).

[Cui et al.(2015)Cui, Lee, Kim, Arefe,
  Huang, Lee, Chenet, Zhang, Wang, Ye et al.]cui2015multi
authorX. Cui,
  authorG.-H. Lee,
  authorY. D. Kim,
  authorG. Arefe,
  authorP. Y. Huang,
  authorC.-H. Lee,
  authorD. A. Chenet,
  authorX. Zhang,
  authorL. Wang,
  authorF. Ye, et al.,
  journalNature nanotechnology volume10,
  pages534 (year2015).

[Radisavljevic and
  Kis(2013)]radisavljevic2013mobility
authorB. Radisavljevic
  and authorA. Kis,
  journalNature materials volume12,
  pages815 (year2013).

[Penumatcha et al.(2015)Penumatcha,
  Salazar, and Appenzeller]penumatcha2015analysing
authorA. V. Penumatcha,
  authorR. B. Salazar,
  and
  authorJ. Appenzeller,
  journalNature communications volume6
  (year2015).

[Datta(1997)]datta1997electronic
authorS. Datta,
  titleElectronic transport in mesoscopic systems
  (publisherCambridge university press, year1997).

[Zhang et al.(2017)Zhang, Huang, Chaves,
  Song, Özçelik, Low, and Yan]zhang2017infrared
authorG. Zhang,
  authorS. Huang,
  authorA. Chaves,
  authorC. Song,
  authorV. O. Özçelik,
  authorT. Low, and
  authorH. Yan,
  journalNature Communications volume8,
  pages14071 (year2017).

[Kang et al.(2016)Kang, Zhang, and
  Wei]kang2016unified
authorJ. Kang,
  authorL. Zhang, and
  authorS.-H. Wei,
  journalThe journal of physical chemistry letters
  volume7, pages597 (year2016).

[Miller(2008)]miller2008quantum
authorD. A. Miller,
  titleQuantum mechanics for scientists and engineers
  (publisherCambridge University Press, year2008).

[Abramowitz and Stegun(1964)]abramowitz1964handbook
authorM. Abramowitz and
  authorI. A. Stegun,
  titleHandbook of mathematical functions: with formulas,
  graphs, and mathematical tables, number55
  (publisherCourier Corporation, year1964).

[Rudan(2015)]rudan2015physics
authorM. Rudan,
  titlePhysics of Semiconductor Devices
  (publisherSpringer, year2015).

[Bardeen(1961)]bardeen1961tunnelling
authorJ. Bardeen,
  journalPhysical Review Letters volume6,
  pages57 (year1961).

[Harrison(1961)]harrison1961tunneling
authorW. A. Harrison,
  journalPhysical Review volume123,
  pages85 (year1961).

[Duke(1969)]duke1969tunneling
authorC. B. Duke,
  titleTunneling in solids, vol. volume10
  (publisherAcademic Pr, year1969).

[Feenstra et al.(2012)Feenstra, Jena, and
  Gu]feenstra2012single
authorR. M. Feenstra,
  authorD. Jena, and
  authorG. Gu,
  journalJournal of Applied Physics volume111,
  pages043711 (year2012).

[Britnell et al.(2013)Britnell,
  Gorbachev, Geim, Ponomarenko, Mishchenko, Greenaway, Fromhold, Novoselov, and
  Eaves]britnell2013resonant
authorL. Britnell,
  authorR. Gorbachev,
  authorA. Geim,
  authorL. Ponomarenko,
  authorA. Mishchenko,
  authorM. Greenaway,
  authorT. Fromhold,
  authorK. Novoselov,
  and authorL. Eaves,
  journalNature communications volume4,
  pages1794 (year2013).

[Alper et al.(2013)Alper, Lattanzio,
  De Michielis, Palestri, Selmi, and Ionescu]alper2013quantum
authorC. Alper,
  authorL. Lattanzio,
  authorL. De Michielis,
  authorP. Palestri,
  authorL. Selmi, and
  authorA. M. Ionescu,
  journalElectron Devices, IEEE Transactions on
  volume60, pages2754 (year2013).

[Agarwal et al.(2014)Agarwal, Teherani,
  Hoyt, Antoniadis, and Yablonovitch]agarwal2014engineering
authorS. Agarwal,
  authorJ. T. Teherani,
  authorJ. L. Hoyt,
  authorD. A. Antoniadis,
  and
  authorE. Yablonovitch,
  journalElectron Devices, IEEE Transactions on
  volume61, pages1599 (year2014).

[Anwar et al.(1999)Anwar, Nabet, Culp,
  and Castro]anwar1999effects
authorA. Anwar,
  authorB. Nabet,
  authorJ. Culp, and
  authorF. Castro,
  journalJournal of applied physics volume85,
  pages2663 (year1999).

[Anugrah et al.(2015)Anugrah, Robbins,
  Crowell, and Koester]anugrah2015determination
authorY. Anugrah,
  authorM. C. Robbins,
  authorP. A. Crowell,
  and authorS. J.
  Koester, journalApplied Physics Letters
  volume106, pages103108
  (year2015).

[Stern(1972)]stern1972self
authorF. Stern,
  journalPhysical Review B volume5,
  pages4891 (year1972).







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















§ 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

    ℋ = -ħ^2/2∇·m̂^-1∇ + E_c(z) = 𝒯 + E_c(z),   m̂ = (
    [ m_x   0   0;   0 m_y   0;   0   0 m_z ]),

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. <ref>a):

    E_c(z) = {[     V_0,     z< 0; - q F z,  0 <z< a;       ∞,    z > a ].


Note that this is different from the Fowler-Nordheim field-emission problem <cit.> 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 <cit.> 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 <cit.>. 

More precisely, ℋ 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 ℋ_L that approximates well the true Hamiltonian ℋ on the metal side of the junction but whose eigenfunctions decay in the semiconductor. We take ℋ_L = 𝒯 + E_cL(z) with E_cL(z) a potential step (Fig. <ref>b):

    E_cL(z) = {[  V_0,  z< 0;    0, z > 0 ].

For t>t_P, one must choose a different unperturbed Hamiltonian ℋ_R which, conversely, approximates well ℋ on the semiconductor side of the junction but whose eigenfunctions decay in the metal. We take ℋ_R = 𝒯 + E_cR(z) with E_cR(z) a triangular well (Fig. <ref>c):

    E_cR(z) = {[ - q F z,    z < a;       ∞,    z > a ].

It is assumed that prior to the perturbation the electron wavefunction coincides with an eigenfunction ψ_L,α of ℋ_L. Since this is not an eigenstate of ℋ, 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 ψ_R,β of ℋ_R at t>t_P, is, to first order and per unit t_P,

    P_αβ = 2π/ħ| ⟨ψ_R,β | ℋ-ℋ_L | ψ_L,α⟩|^2 δ( E_L,α - E_R,β)   ,

where E_L,α and E_R,β are the eigenvalues corresponding to the initial state ψ_L,α and final state ψ_R,β, respectively, and δ is Dirac's delta function <cit.>. The conservation of energy is related to the perturbation being constant in time. In (<ref>), it is assumed that each set of eigenfunctions ψ_L,α and ψ_R,β 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 ψ_L,α/R,β(x=0,y,z)=ψ_L,α/R,β(x=L_x,y,z), ψ_L,α/R,β(x,y=0,z)=ψ_L,α/R,β(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 ψ_L,α(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 <cit.>, the matrix element in (<ref>) must be computed between eigenstates of different Hamiltonians. Also, the perturbation Hamiltonian must be evaluated with respect to the initial Hamiltonian. The validity of (<ref>) rests on the assumption that the two sets of eigenstates of ℋ_L and ℋ_L are “almost orthogonal” to each other, in particular that ⟨ψ_R,β | ψ_L,α⟩≪ 1 <cit.>. Similary, one has for the probability rate of the inverse transition

    P_βα = 2π/ħ| ⟨ψ_L,α | ℋ-ℋ_R | ψ_R,β⟩|^2 δ( E_L,α - E_R,β) = P_αβ  ,

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

From (<ref>)-(<ref>) we get

    ℋ = {[ ℋ_L,  z<0; ℋ_R,  z>0 ].

which means that ℋ satisfies the separability property of Bardeen's model Hamiltonian, i.e., that ℋ - ℋ_L  0 only in regions of space where ℋ - ℋ_R ≡ 0 <cit.> [Bardeen's theory is often introduced by writing the Hamiltonian in the form ℋ = ℋ_L + ℋ_R + ℋ_T, with ℋ_T being the “transfer” Hamiltonian, despite the fact that Bardeen himself did not make use of such decomposition in its original paper <cit.>. In our case, we get from (<ref>)

    ℋ_T = {[ -ℋ_R,   z<0; -ℋ_L,   z>0 ].

but this Hamiltonian does not correspond to either of the perturbation Hamiltonians that appear in (<ref>) or (<ref>). See also the discussion in <cit.>.]. Since ℋ - ℋ_L ≡ 0 for z<0, the matrix element in (<ref>) can be written as

    ⟨ψ_R,β | ℋ-ℋ_L | ψ_L,α⟩   = ∫_0^L_x dx ∫_0^L_y dy ∫_-∞^∞ dz  ψ_R,β^* ( ℋ-ℋ_L ) ψ_L,α
       = ∫_Ω_Rψ_R,β^* ( ℋ-ℋ_L ) ψ_L,α d^3 r   ,

where 𝐫=(x,y,z) and Ω_R = {𝐫;0<x<L_x,0<y<L_y,0<z<∞}. Noting that ℋ - ℋ_R ≡ 0 for z>0, we can get the symmetric expression

    ⟨ψ_R,β | ℋ-ℋ_L | ψ_L,α⟩   = ∫_Ω_R[ ψ_R,β^* ( ℋ-ℋ_L ) ψ_L,α -  ψ_L,α( ℋ-ℋ_R ) ψ_R,β^* ] d^3 r 
       = ∫_Ω_R[ ψ_R,β^* ( 𝒯-E_L,α) ψ_L,α -  ψ_L,α( 𝒯-E_R,β) ψ_R,β^* ] d^3 r   .

Due to the delta function in (<ref>), we are only interested in the case E_L,α = E_R,β, for which

    ⟨ψ_R,β | ℋ-ℋ_L | ψ_L,α⟩ = ∫_Ω_R[ ψ_R,β^* 𝒯ψ_L,α -  ψ_L,α𝒯ψ_R,β^* ] d^3 r   .

Applying Green's theorem, we finally get

    ⟨ψ_R,β | ℋ-ℋ_L | ψ_L,α⟩ = -iħ∫_Σ_R𝐧·𝐉_βα d^2 r   ,

where Σ_R is the surface of Ω_R, 𝐧 is the unit vector normal to Σ_R pointing in the outward direction, and

    𝐉_βα = -iħ/2m̂^-1[ ψ_R,β^* ∇ψ_L,α -  ψ_L,α∇ψ_R,β^* ]

is the matrix element of the probability current density operator between the states ψ_L,α and ψ_R,β.

Let us now start to pick up all the ingredients that we need to calculate (<ref>). The eigenfunctions and corresponding eigenvalues of ℋ_L are <cit.>

    ψ_L,α(𝐫)    ≡ψ_L(𝐤_L;𝐫) = 1/√(L_x L_y) e^i (k_xL x + k_yL y) b ×{[              0,          z<-L_z; cos(k_z z + φ),        -L_z<z<0;   cos(φ)e^-κ z,             z>0 ]. 
    
    E_L,α   ≡ E_L(𝐤_L) = ħ^2/2( k_xL^2/m_x + k_yL^2/m_y) + E_z   , 
    
    E_z    = V_0 + ħ^2 k_z^2/2 m_z  , 
    κ   = √(- 2 m_z E_z)/ħ  , 
    φ   = arctan(κ/k_z)   ,

where 𝐤_L = (k_xL,k_yL,k_z) and it is assumed that E_z < 0 [We limit the discussion to the case E_z < 0 because (<ref>) looses validity if ⟨ψ_R,β | ψ_L,α⟩ is not a small number.].
Because of the periodic boundary conditions, the transverse components of the wavevector are quantized as k_xL = 2 π l/L_x, k_yL = 2 π m/L_y (l,m integers). As for k_z, the allowed values are the roots of the transcendental equations k_z L_z - π (n-1/2) = φ in the interval 0 < k_z <√(-2 m_z V_0)/ħ. For large L_z, we have k_z ≈π n/L_z (n positive integer). The constant b can be obtained from the normalization condition

    1 = |b|^2 [ ∫_-L_z^0cos^2(k_z z + φ) dz + cos^2(φ) ∫_0^∞ e^-2 κ z dz ] = |b|^2/2( L_z+1/κ) ≈ |b|^2 L_z/2  ,

where only the leading term in L_z has been kept. Thus, up to an unimportant phase, b = √(2/L_z).

The solutions of the eigenvalue problem of ℋ_R are <cit.>

    ψ_R,β(𝐫)    ≡ψ_R,i(𝐤_∥ R; 𝐫) = 1/√(L_x L_y) e^i (k_xR x + k_yR y) c ×{[ Ai(ζ),    z<a;     0,    z>a ]. 
    
    E_R,β   ≡ E_R,i(𝐤_∥ R) = ħ^2/2( k_xR^2/m_x + k_yR^2/m_y) + E_i   , 
    
    E_i    = - qF ( a - | ζ_i |/k_F) , 
    ζ   = -k_F ( z+E_i/qF), 
    
    k_F    = ( 2 m_z qF/ħ^2)^1/3 ,

where 𝐤_∥ R = (k_xR,k_yR), Ai is Airy's function and ζ_i are its zeros, which can be approximated as (i positive integer) <cit.>

    ζ_i ≈ - [ 3 π/8 (4 i-1) ]^2/3 .

The constant c can be obtained from the normalization condition

    1 = |c|^2 ∫_-∞^aAi^2(ζ) dz = |c|^2/k_F∫_ζ_i^∞Ai^2(ζ) dζ  ,

where the last integral can be evaluated using integration by parts and the fact that Ai is a solution of Airy's equation Ai” = ζAi (the prime symbol indicates derivative with respect to ζ):

    ∫_ζ_i^∞Ai^2(ζ) dζ = - ∫_ζ_i^∞ 2 Ai(ζ) Ai'(ζ) ζ dζ = - ∫_ζ_i^∞ 2 Ai'(ζ) Ai”(ζ) dζ = Ai'^2(ζ_i)   .

Therefore, c = √(k_F)/|Ai'(ζ_i)|, in which we can use the approximated expression <cit.>

    |Ai'(ζ_i)| ≈1/√(π)[ 3 π/8 (4 i-1) ]^1/6 .


The surface Σ_R in (<ref>) is made up of six faces. By inserting (<ref>) and (<ref>) into (<ref>), 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 [For example, J_x,βα (x=L_x,y,z)- J_x,βα (x=0,y,z) ∝ e^i (k_xL -k_xR) L_x - 1=0 because of the periodic boundary conditions.]. The integral over the face at z=∞ is also zero because ψ_L,α is vanishingly small. We are only left with the integral over the face at z=0:

    ⟨ψ_R,β | ℋ-ℋ_L | ψ_L,α⟩   = ħ^2/2 m_z√(2/L_z)√(k_F)/|Ai'(ζ_i)|cos(φ) [ Ai(ζ) d/dz e^-κ z - e^-κ zd/dzAi(ζ) ]_z=0
         ×1/L_x∫_0^L_x e^i (k_xL -k_xR) x dx 1/L_y∫_0^L_y e^i (k_yL -k_yR) y dy 
       = - ħ^2/2 m_z√(2/L_z)√(k_F)/|Ai'(ζ_i)|cos(φ) [ κAi(ζ_0) - k_F Ai'(ζ_0) ] δ_k_xL,k_xRδ_k_yL,k_yR  ,

where

    ζ_0 ≡ζ(z=0) = - k_F E_i/qF = k_F a - |ζ_i|

and δ 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 ζ_0 = ( κ/ k_F )^2. Assuming ζ_0 ≫ 1 (which is consistent with ⟨ψ_R,β | ψ_L,α⟩≪ 1), we can substitute in (<ref>) the asymptotic expressions

    Ai(ζ_0)    ≈e^-γ_0/2 √(π)ζ_0^1/4  ,
    Ai'(ζ_0)    ≈ - ζ_0^1/4 e^-γ_0/2 √(π)  ,

where γ_0 = (2/3) ζ_0^3/2 <cit.>, to get

    ⟨ψ_R,β | ℋ-ℋ_L | ψ_L,α⟩   = - ħ^2/2 m_z√(2/L_z)√(k_F)/|Ai'(ζ_i)|cos(φ) [ κ + k_F √(ζ_0)] e^-γ_0/2 √(π)ζ_0^1/4δ_k_xL,k_xRδ_k_yL,k_yR
       = - ħ^2/2 m_z√(2/L_z)k_F/|Ai'(ζ_i)|2 k_z √(κ)/√(k_z^2+κ^2)e^-γ_0/2 √(π)δ_k_xL,k_xRδ_k_yL,k_yR  .

Finally, plugging (<ref>) into (<ref>), we obtain

    P_αβ = 1/h( ħ^2/2 m_z)^22 π/L_zk_F^2/Ai'^2(ζ_i)4 k_z^2 κ/k_z^2+κ^2 e^-2γ_0δ_k_xL,k_xRδ_k_yL,k_yRδ( E_L,α - E_R,β)   .




§ MODEL OF LONGITUDINAL DIFFUSION
 

Suppose that the states in the metal (L) are populated according to a Fermi-Dirac distribution with Fermi level μ_L:

    f_L(E_L,α) = 1/exp( E_L,α-μ_L/k_B T) + 1

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 <cit.>. 

Let f_i(x,𝐤_∥ R) be the distribution function in the four-dimensional phase space associated with the i-th subband [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

    v_x ∂ f_i/∂ x = C,    0<x<L_x

where v_x = ħ 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 [Other types of scattering, which could be responsible for a finite carrier mobility in the semiconductor, are here neglected.]:

    C    = ∑_α f_L(E_L,α) P_αβ ( 1-f_i) - f_i P_βα[ 1 - f_L(E_L,α) ] 
       = ∑_α P_αβ[ f_L(E_L,α) - f_i ]   .

Note that expression (<ref>) takes into account Pauli's exclusion principle. Using (<ref>), we get

    C = [ f_L(E_R,β) - f_i ] ∑_α P_αβ =  f_L(E_R,β) - f_i /τ_i  ,

where the relaxation time τ_i is defined as 

    1/τ_i = ∑_α P_αβ = ∑_k_z1/h( ħ^2/2 m_z)^22 π/L_zk_F^2/Ai'^2(ζ_i)4 k_z^2 κ/k_z^2+κ^2 e^-2γ_0δ( E_z - E_i )   .

Going to the limit of large L_z, we can replace

    ∑_k_z→∫L_z/π d k_z

so that

    1/τ_i = ∫_0^√(-2 m_z V_0)/ħ d k_z 1/h( ħ^2/2 m_z)^2 2 k_F^2/Ai'^2(ζ_i)4 k_z^2 κ/k_z^2+κ^2 e^-2γ_0δ( E_z - E_i )

and, with the change of variables k_z → E_z,

    1/τ_i   = ∫_V_0^0 d E_z 1/hħ^2/2 m_zk_F^2/Ai'^2(ζ_i)4 √(-E_z (E_z - V_0))/-V_0 e^-2γ_0δ( E_z - E_i ) 
       = {[ 1/hħ^2/2 m_zk_F^2/Ai'^2(ζ_i)4 √(-E_i (E_i - V_0))/-V_0 e^-2γ_0,                                                   V_0 < E_i < 0;                                                              0,                                                       otherwise ].

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 (<ref>) of ℋ_L by their WKB <cit.> approximation

    ψ_L,α^WKB(𝐫) = 1/√(L_x L_y) e^i (k_xL x + k_yL y)√(2/L_z)×{[                  0,              z<-L_z;   cos(k_z z + π/4),            -L_z<z<0; 1/2√(k_z/κ) e^-κ z,                 z>0 ].

the expression of τ_i simplifies to

    1/τ_i^WKB = {[ 1/hħ^2/2 m_zk_F^2/Ai'^2(ζ_i) e^-2γ_0,                         V_0 < E_i < 0;                                    0,                             otherwise ].

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^±(x,k_xR,k_yR) = f_i(x,± k_xR,k_yR) with k_xR>0. We impose the boundary conditions

    f_i^-(x = L_x,𝐤_∥ R)    = f_R(E_R,β)   , 
    
    f_i^+(x = 0,𝐤_∥ R)    = f_i^-(x = 0,𝐤_∥ R)   ,

where f_R is a Fermi-Dirac function similar to (<ref>) with μ_L replaced by μ_R (see Fig. <ref>). The latter condition makes sure that the longitudinal current vanishes at x=0. The Boltzmann equations (<ref>) for different wavevectors are independent of each other expect for ± k_xR. The solutions are

    f_i^±   = -[ f_L(E_R,β) - f_R(E_R,β) ] e^- L_x ± x/|v_x| τ_i + f_L(E_R,β)   .


The current (per unit width) from the semiconductor to the metal can be obtained by summing the net flux |v_x|(f_i^+-f_i^-) at x = L_x over all semiconductor states per unit area, multiplying by 2 for spin degeneracy and multiplying by the electronic charge q:

    I    = 2q/L_x L_y∑_k_xR>0∑_k_yR∑_i |v_x| (f_i^+-f_i^-)_x=L_x
       = 2 q/L_x L_y∑_k_xR>0∑_k_yR∑_i |v_x| ( 1 - e^-2 L_x/|v_x| τ_i) [ f_L(E_R,β) - f_R(E_R,β) ] .

Going to the limit of large L_x,L_y, we can replace

    ∑_k_xR→∫L_x/2π d k_xR  ,   ∑_k_yR→∫L_y/2π d k_yR

to get

    I = 2q/( 2 π)^2∫_0^∞ d k_xR∫_-∞^∞ d k_yR∑_i |v_x| ( 1 - e^-2 L_x/|v_x| τ_i) [ f_L(E_R,β) - f_R(E_R,β) ] .

Finally, with the change of variables k_xR→ε = ħ^2 k_xR^2 / (2 m_x) + E_i, we obtain the Landauer formula <cit.>

    I    = 2q/h∑_i∫_E_i^∞ d ε ( 1 - e^-2 L_x/|v_x| τ_i) 1/2 π∫_-∞^∞ d k_yR[ f_L( ħ^2 k_yR^2/2 m_y + ε) - f_R( ħ^2 k_yR^2/2 m_y + ε) ] 
       = 2q/h∑_i∫_E_i^∞ d ε ( 1 - e^-2 L_x/|v_x| τ_i) [ F_L(ε) - F_R(ε) ] 
       =  2q/h∫_-∞^∞ d ε T(ε) [ F_L(ε) - F_R(ε) ] ,

where the longitudinal velocity must be computed as |v_x| = √( 2 (ε-E_i) / m_x), the supply function F_L/R is defined as

    F_L/R(ε) = √(m_y k_B T/2 πħ^2)ℱ_-1/2( μ_L/R - ε/k_B T)

with ℱ_-1/2 the Fermi-Dirac integral of order -12 <cit.>, T is the transmission function

    T(ε) = ∑_i T_i(ε)   ,

with the transmission probability T_i given by

    T_i(ε) = {[                0,           ε < E_i; 1 - e^-2 L_x/λ_i,           ε > E_i ].

and λ_i = |v_x| τ_i. For ε > E_i, we can have the asymptotic behaviors 

    L_x ≫λ_i    :   T_i(ε) ≈ 1 
    
    L_x ≪λ_i    :   T_i(ε) ≈2 L_x/λ_i

Approximation (<ref>) holds, in particular, as ε→ E_i^+ (resonant tunneling). Note also that, when (<ref>) is satisfied, T_i decays as a function of energy as 1/√(ε).

The low bias conductance per unit width G can be evaluated from (<ref>). Let μ_L=μ and μ_R=μ-qV_D. We have

    G = .∂ I/∂ V_D|_V_D=0 = 2q^2/h∫_-∞^∞ d ε T(ε) (-∂ F_0/∂ε) ,

where F_0(ε) is given by (<ref>) with μ_L/R replaced by μ.

Assume for simplicity that only the first subband contributes to transport. In the two limiting cases when either (<ref>) or (<ref>) are satisfied over the whole energy range of interest for transport (from E_1 to few k_B T's above max{E_1,μ}), it is possible to derive analytical expressions for G. If L_x ≫λ_1 in this energy range, it follows immediately from (<ref>)

    G = 2q^2/h F_0(E_1)   .

To work out the expression of G in the other limiting case when L_x ≪λ_1 in most of the energy window for transport[The inequality does not hold for energies close to E_1 but their contribution becomes increasingly smaller as τ_1 increases.], it is convenient to go back to the double-integral formulation of the tunneling current in (<ref>) and do the change of variables k_xR→ E_R,β:

    I = 2q/h∫_E_1^∞ d E_R,β( 1/π∫_0^k_ymax d k_yR2 L_x/|v_x| τ_1) [ f_L(E_R,β) - f_R(E_R,β) ]   ,

where k_ymax = √(2 m_y (E_R,β-E_1))/ħ and |v_x| = √( 2 [E_R,β-ħ^2 k_yR^2/(2 m_y)-E_1] / m_x). The integral over the transverse wavevector can be easily computed with the change of variables k_yR→arcsin(k_yR/k_ymax) to give

    I    = 2q/h∫_E_1^∞ d E_R,βk_ymax L_x/√(2 (E_R,β-E_1)/m_x)τ_1[ f_L(E_R,β) - f_R(E_R,β) ] 
       = 2q/h√(m_x m_y)/ħL_x/τ_1∫_E_1^∞ d E_R,β[ f_L(E_R,β) - f_R(E_R,β) ]   .

Letting μ_L=μ and μ_R=μ-qV_D, we finally get

    G = .∂ I/∂ V_D|_V_D=0 = 2q^2/h√(m_x m_y)/ħL_x/τ_1∫_E_1^∞ d E_R,β(-∂ f_0/∂ E_R,β) = 2q^2/h√(m_x m_y)/ħL_x/τ_1 f_0(E_1)   ,

where f_0 is a Fermi-Dirac function similar to (<ref>) with μ_L replaced by μ.



§ GREEN FUNCTION ALGORITHM



We discretize the Hamiltonian in (<ref>) using finite differences on a two-dimensional rectangular grid (Fig. <ref>a). The same linear potential profile as in (<ref>) is assumed. The transmission function is computed as <cit.>

    T(ε) = tr[ Γ^L G^r Γ^R G^a ] ,

where the symbol tr indicates the trace, G^r is the retarded Green function, G^a=G^r†, Γ^L/R = i (Σ^r,L/R-Σ^a,L/R), with Σ^r,L/R the retarded self-energy representing the renormalization of the Hamiltonian of the channel region (black dots in Fig. <ref>a) due to the presence of the semi-infinite left/right lead, and Σ^a,L/R = ( Σ^r,L/R )^†. The channel region is partitioned into layers as shown in Fig. <ref>b. Using matrix block notation and noting that the only non-null block of Σ^r,L is Σ^r,L_0,0 and the only non-null block of Σ^r,R is Σ^r,R_N_x+1,N_x+1, (<ref>) can be rewritten as

    T(ε) = tr[ Γ^L_0,0 G^r_0,N_x+1Γ^R_N_x+1,N_x+1G^a_N_x+1,0] .

The self-energy of the left lead is computed analytically using the prescription given in <cit.>:

    Σ^r,L_0,0(i,i')    = ∑_m=1^N_xχ_m(i) σ_m χ_m(i')   , 
    χ_m(i)    = √(2/N_x+1)sin(k_x i)   , 
    σ_m    = t_z ×{[  λ-1+√(λ^2-2λ),           λ < 0;  λ-1-√(λ^2-2λ),           λ > 2; λ-1-i√(2λ-λ^2),       0 < λ < 2 ]. 
    λ   = ε - V_0 - 2 t_x (1 - cos k_x)/2 t_z  , 
    
    k_x    = π m/N_x+1  ,

where t_x = ħ^2/(2 m_x Δ_x^2) and similarly for t_z. The self-energy of the right lead is obtained numerically using a well-known iterative algorithm <cit.>. The matrix block G^r_0,N_x+1 is computed through a combination of the recursive and decimation algorithms <cit.>, modified so as to treat a non-tridiagonal-block Hamiltonian matrix. Let A = ε I - H_C - Σ^r,L - Σ^r,R, where H_C is the Hamiltonian matrix of the channel region alone, and define δ_1^(0) = A_0,0, δ_2^(0) = A_1,1, α^(0) = A_0,1, β^(0) = A_1,0. The algorithm consists in eliminating the layers from 1 to N_x with the formulas

    δ_1^(n)   = δ_1^(n-1) - α^(n-1)[ δ_2^(n-1)]^-1β^(n-1)  , 
    δ_2^(n)   = A_n+1,n+1 -   A_n+1,n[ δ_2^(n-1)]^-1 A_n,n+1  , 
    α^(n)   = -α^(n-1)[ δ_2^(n-1)]^-1 A_n,n+1 + A_0,n+1  , 
    β^(n)   = - A_n+1,n[ δ_2^(n-1)]^-1β^(n-1) + A_n+1,0

for n=1,…,N_x, where it is understood that A_0,N_x+1=A_N_x+1,0^† = 0. At the end, the required matrix block of the Green function can be obtained as

    G^r_0,N_x+1 = - [ δ_1^(N_x)]^-1α^(N_x){δ_2^(N_x) - β^(N_x)[ δ_1^(N_x)]^-1α^(N_x)}^-1 .


12
natexlab#1#1bibnamefont#1#1bibfnamefont#1#1citenamefont#1#1url<#>1urlprefixURL 


[Fowler and Nordheim(1928)]fowler1928electron
authorR. H. Fowler and
  authorL. Nordheim, in
  booktitleProceedings of the Royal Society of London A:
  Mathematical, Physical and Engineering Sciences (organizationThe
  Royal Society, year1928), vol. volume119, pp.
  pages173–181.

[Bardeen(1961)]bardeen1961tunnelling_S
authorJ. Bardeen,
  journalPhysical Review Letters volume6,
  pages57 (year1961).

[Harrison(1961)]harrison1961tunneling_S
authorW. A. Harrison,
  journalPhysical Review volume123,
  pages85 (year1961).

[Duke(1969)]duke1969tunneling_S
authorC. B. Duke,
  titleTunneling in solids, vol. volume10
  (publisherAcademic Pr, year1969).

[Oppenheimer(1928)]oppenheimer1928three
authorJ. R. Oppenheimer,
  journalPhysical review volume31,
  pages66 (year1928).

[Rudan(2015)]rudan2015physics_S
authorM. Rudan,
  titlePhysics of Semiconductor Devices
  (publisherSpringer, year2015).

[Miller(2008)]miller2008quantum_S
authorD. A. Miller,
  titleQuantum mechanics for scientists and engineers
  (publisherCambridge University Press, year2008).

[Abramowitz and Stegun(1964)]abramowitz1964handbook_S
authorM. Abramowitz and
  authorI. A. Stegun,
  titleHandbook of mathematical functions: with formulas,
  graphs, and mathematical tables, number55
  (publisherCourier Corporation, year1964).

[Messiah(1961)]messiah1961quantum
authorA. Messiah,
  titleQuantum mechanics, vol. 1
  (publisherNorth-Holland, Amsterdam, year1961).

[Datta(1997)]datta1997electronic_S
authorS. Datta,
  titleElectronic transport in mesoscopic systems
  (publisherCambridge university press, year1997).

[Sancho et al.(1985)Sancho, Sancho,
  Sancho, and Rubio]sancho1985highly
authorM. L. Sancho,
  authorJ. L. Sancho,
  authorJ. L. Sancho,
  and authorJ. Rubio,
  journalJournal of Physics F: Metal Physics
  volume15, pages851 (year1985).

[Low and Appenzeller(2009)]low2009electronic
authorT. Low and
  authorJ. Appenzeller,
  journalPhysical Review B volume80,
  pages155406 (year2009).



