Spectroscopy of multi-electrode tunnel barriers

TL;DR

This paper develops and experimentally validates a compact, physics-based model for multi-electrode tunnel barriers in nanoscale electronic devices, crucial for quantum device engineering and understanding electron transport phenomena.

Contribution

It introduces a new theoretical model for multi-electrode tunnel barriers, validated by transport spectroscopy experiments, capturing voltage-dependent barrier behaviors and effects on electron transport.

Findings

Barrier height exhibits linear and exponential voltage dependence.

Exponential dependence correlates with electron channel formation.

Model agrees qualitatively with semi-classical Poisson solutions.

Abstract

Despite their ubiquity in nanoscale electronic devices, the physics of tunnel barriers has not been developed to the extent necessary for the engineering of devices in the few-electron regime. This problem is of urgent interest, as this is the precise regime into which current, extreme-scale electronics fall. Here, we propose theoretically and validate experimentally a compact model for multi-electrode tunnel barriers, suitable for design-rules-based engineering of tunnel junctions in quantum devices. We perform transport spectroscopy at $T=4$ K, extracting effective barrier heights and widths for a wide range of biases, using an efficient Landauer-B\"uttiker tunneling model to perform the analysis. We find that the barrier height shows several regimes of voltage dependence, either linear or approximately exponential. The exponential dependence approximately correlates with the formation of an electron channel below an electrode. Effects on transport threshold, such as metal-insulator-transition and lateral confinement are non-negligible and included. We compare these results to semi-classical solutions of Poisson's equation and find them to agree qualitatively. Finally, we characterize the sensitivity of a tunnel barrier that is raised or lowered without an electrode being directly above the barrier region.

Figures

Tables

Table 1. Table 1: Extracted barrier parameters U 0 subscript 𝑈 0 U_{0} and w 𝑤 w for fits shown in Fig.  3 (a).

$V_L$ (V)	DOS of leads	$U_0$ (meV)	$w$ (nm)
0.90	0D	17$\pm$3	12$\pm$3
0.80	0D	18$\pm$4	11$\pm$3
0.70	0D	19$\pm$5	12$\pm$3
0.70	1D	27$\pm$8	12$\pm$3
0.70	2D	29$\pm$9	12$\pm$3
0.60	0D	20$\pm$5	12$\pm$3

Equations

$$E_f-E_C=\frac{\pi {\hslash }^{2}n}{m_t{g}_v},$$

$$n\left(V_{gate}\right)=\frac{C_{ox}}{e}(V_{gate}-V_{th})+n_{crit},$$

$$I={\int }_{{\mu }_0}^{{\mu }_S}\mathrm{d}E\frac{2e}{h}T(E)M(E).$$

$$V_{th}=V_0+V_{crit}.$$

$$V_0=V_{FB}+2|{\varphi }_p|+\frac{1}{C_{ox}}\sqrt{4{ϵ}_{s}e{N}_a|{\varphi }_p|},$$

$$E_0^L=m_L{V}_G,$$

$$V_{\infty }=\frac{E^0}{m_{\infty }},$$

$$\Delta V_L$$

$$\Delta V_L=V_{\infty }\left(\frac{m_{\infty }}{m_L}-1\right).$$

$$V_L=V_{\infty }\frac{m_{\infty }}{m_L}.$$

$$V_L=V_{\infty }\frac{m_{\infty }}{m_{\infty }+a{e}^{-bL}}.$$

$$I=\int dE\frac{2e}{h}\bar{T}(E)\left(f_s(E)-f_d(E)\right),$$

$$\bar{T}(E)$$

$$=T(E)M(E),$$

$${\psi }_1(x)$$

$${\psi }_2(x)$$

$${\psi }_3(x)$$

$${\psi }_1^{\prime }(x)$$

$${\psi }_2^{\prime }(x)$$

$${\psi }_3^{\prime }(x)$$

$$\text{mbox}Ai(z)$$

$$\text{mbox}Bi(z)$$

$${\text{Ai}}_{N=0}(a(0))$$

$${\text{Bi}}_{N=0}(a(0))$$

$${\psi }_2(x)$$

$$\equiv {\mathcal{A}}_2\mathcal{A}i\left(a(x)\right)+{\mathcal{B}}_2\mathcal{B}i\left(a(x)\right),$$

$$\underset{F\to 0}{\lim }\mathcal{A}i\left(a(x)\right)=e^{\frac{x}{\hslash }\sqrt{2m(U_0-E)}},$$

$$\underset{F\to 0}{\lim }\mathcal{B}i\left(a(x)\right)=e^{-\frac{x}{\hslash }\sqrt{2m(U_0-E)}},$$

$${\psi }_2(x)$$

$${\psi }_2^{\prime }(x)$$

$${\psi }_1(x)$$

$${\psi }_2(x)$$

$${\psi }_3(x)$$

$${\psi }_1^{\prime }(x)$$

$${\psi }_2^{\prime }(x)$$

$${\psi }_3^{\prime }(x)$$

$$1$$

$$i{k}_1$$

Full text

I Introduction

Electron tunneling is a phenomenon that is generally well-understood having been studied in numerous configurations over many decades Sze (1981); Davies (1997). In recent semiconductor devices, tunneling effects are becoming increasingly important as dimensions of devices are scaled to their ultimate limits. Many Beyond-Moore’s-Law device concepts, for example, invoke devices for which tunneling is a central element of the device (e.g., tunnel FETs) and the tunnel barrier is necessarily voltage tunable Ionescu and Riel (2011). In practice, numerous phenomenological tunnel models and interpretations of the models are invoked to design or extract barrier heights and widths (e.g., WKB Simmons (1963), Fowler-Nordheim Lenzlinger and Snow (1969); Huisman et al. (2009), and Non-equilibrium Green’s function Ilatikhameneh et al. (2015); Gao et al. (2014)). Despite the profound foundational understanding of tunneling, the details of the barrier and its geometry are known to produce non-trivial effects Friesen et al. (2013). This challenges both design and analysis of tunnel-based devices.

Voltage-tuned barriers represent a significant increase in complexity in contrast to barriers formed by static material barriers (e.g., heterostructures). In these circumstances the barriers deform according to the applied voltages and it is unclear how effective simple phenomenological models that assert, for example, two parameter descriptions like width and barrier height are. Furthermore, the functional dependence of these parameters on voltage is not well understood. Observations of linear dependence of barrier height on voltage have been conjectured for some geometriesMacLean et al. (2007), but it is not clear under what conditions this holds or what to expect in general.

In this work, we use tunneling spectroscopy at 4 K to characterize the voltage dependence of a barrier in a lateral, electrostatically gated nanostructure using a MOS gate stack. In particular, we methodically examine and parameterize the effect of multiple electrodes on how the barrier deforms in this geometry, which extends previous studies in this area. The lateral, electrostatically gated MOS nanostructure is a model system that applies abstractly to many other systems while also being directly informative to the ubiquitous MOS system. We show that a two-parameter model (width and barrier height) is sufficient to describe the observed tunneling over a wide range of voltages as long as the barrier height and width are dependent on electrode voltages. We find that the functional barrier height dependence has multiple regimes, ranging from linear to approximately exponential. The non-linear regime correlates with the formation of an electron channel in this MOS enhancement mode configuration. Changes in local bias from neighboring electrodes produce voltage shifts in the overall behavior, which otherwise can be described by similar parameterization. This suggests the possibility that tunnel barriers might be simulated with a relatively simple compact model, which would greatly simplify modeling of more complex laterally gated nanostructures (e.g., quantum dot networks Zajac et al. (2015)). Finally, we characterize the tunnel barrier control that can be achieved using a reservoir enhancement gate instead of a dedicated barrier gate (i.e., a more efficient single metal layer QD layout Rudolph et al. (2016)).

This paper is organized as follows: In Sec. II we describe the device fabrication and transport measurements of the tunnel barrier. In addition we discuss calculations of the electron density and Fermi energy in the leads and consider effects such as metal-insulator transition and lateral quantum confinement. Sec. III describes the 1D barrier model and we present two methods for computing transmission through the barrier. We then discuss the results of fitting the 1D barrier model to the data and compare these results to electrostatic simulations in Sec. IV.

II Experimental approach

The tunnel barriers are formed in a metal oxide semiconductor (MOS) field effect transistor (FET) inversion channel at a Si/SiO2 interface. The starting material is a p-type Si substrate produced by a float-zone process and with a doping concentration of $\approx 10^{14}$ cm*-3* boron. The structures are made by first forming heavily doped n+ ohmics using ion implantation of As. Next, a 35 nm Si$\text{O}_{2}$ gate dielectric is grown by thermal oxidation at 900*∘C. Degenerately doped n-type poly-Si is formed over the SiO2* to complete the MOS gate stack. The gate structures are patterned using electron beam lithography and a dry etch using a fabrication process described in more detail elsewhere Singh et al. (2016). A scanning electron microscope (SEM) image of the nanostructure gates and cross-sectional schematic view of the device are shown in Figs. 1(a-b). As seen in Fig. 1(a) the relevant gates are separated from each other by approximately 15 nm. The inversion channel is formed under the poly-Si enhancement gates and the barriers are formed at gaps between the poly-Si gates.

A simplified model of the conduction band at low temperature is shown in Fig. 1(c). A high-density electron channel forms when gates L, C, and R are biased above threshold. If we decrease the bias on one of the gates below threshold, electrons are locally depleted and the conduction band edge rises above the Fermi level forming a tunnel barrier with source and drain leads to the left and right of the barrier. The extent to which the Fermi level in the leads is above the conduction band edge is proportional to the electron density. Assuming a 2D density of states (DOS), the difference between the Fermi energy $E_{f}$ and conduction band edge $E_{C}$ is

[TABLE]

where $n$ is the 2D electron density, $g_{v}$ is the valley degeneracy, and $m_{t}$ is the electron transverse effective mass. Spin degeneracy, $g_{s}=2$, is already assumed. We use a valley degeneracy $g_{v}=2$ and an effective mass of $m_{t}=0.19m_{0}$, appropriate for 2D electron inversion layers in Si(100) Beenakker and van Houten (1991), where $m_{0}$ is the free-electron mass.

We measure the differential conductance $G=\mathrm{d}I/\mathrm{d}V_{SD}$ using a lock-in amplifier with an AC excitation voltage of 100 $\mu$V at a modulation frequency of 149 Hz in combination with a DC source-drain bias $V_{SD}$. We also monitor the DC source-drain current $I_{SD}$ using a digital multi-meter. All measurements are performed at 4 K by immersing the sample in liquid He. A plot of the differential conductance as a function of $V_{L}$ and $V_{SD}$ for $V_{R}=2.13$ V and $V_{C}=1.60$ V is shown in Fig. 1(d). Disorder such as charge defects in the gate oxide Nordberg et al. (2009), impurities Lansbergen et al. (2008), or strain from the gate stack Thorbeck and Zimmerman (2015) can cause sub-threshold resonances in these plots. The absence of resonances suggests a single tunnel barrier between the source and drain leads. We also note that the combination of the contact resistance and sheet resistance of the channel leads to approximately 10 k$\Omega$, which limits the maximum conductance.

To calculate the tunneling current, we first need to estimate the voltage dependence of the Fermi energy in the leads. We first estimate the electron density dependence on voltage in the field of the device. In the field we exclude effects due to the nanoscale gates, such as fringing fields and lateral quantum confinement (i.e., we assume a bulk poly-Si gate). We identify a field threshold of 0.48 V, Fig. 1(e). At cryogenic temperatures there is a metal-insulator transition for which there is a critical density of electrons necessary before appreciable conduction starts. We find a critical density of $n_{crit}=3.1\times 10^{11}$ cm*-2* measured in a Hall bar with nominally similar oxide/silicon interface as used in this tunnel barrier experiment. We estimate the critical density by extrapolating the low density mobility to zero mobility van Beveren et al. (2010); Borselli et al. (2011). Assuming that $n=n_{crit}$ at $V_{gate}=V_{th}$ where $V_{th}$ is the threshold voltage then the electron density as a function of gate voltage is

[TABLE]

Where $C_{ox}$ is the oxide capacitance per unit area and $e$ is the electron charge. The electron density and Fermi energy can be calculated in the field using Eqs. 2 and 1. The accounting of the metal-insulator transition accommodates fitting over the full range of voltages including the sub-threshold region using a simple 1D capacitance model to approximate the density in the leads. We note that this compact model approximation for the electron density is not strictly accurate in the sub-threshold region. Conductances at densities below $n_{crit}$ should be negligible. However, since the transmission is relatively weakly dependent on the lead density, we don’t expect significant error for extending the capacitance model for the densities below $n_{crit}$.

We also use the threshold voltage and critical density to extract an effective fixed charge density near the interface of $Q_{f}=3.9\times 10^{10}$ cm*-2*, assuming a standard MOS threshold calculation adjusting parameters for 4 K Muller and Kamins (2002). The details of this calculation are in Appendix A. The fixed charge would notably be below zero if the critical density was not considered. Hall measurements of mobility dependence on density can also be used to quantify a different but related charge density, the Coulomb scattering charge density. We find that the slope of the leading edge of the mobility dependence on density fits with a scattering charge density of $7.4\times 10^{10}$ cm*-2*, in near quantitative agreement with the extracted fixed charge density. The scattering charge density is extracted through interpolation between mobility vs. density curves calculated for the relevant range of scattering charge densities Tracy et al. (2009); Das Sarma and Hwang (2013); Stern (1980). The Hall measurements also show a peak mobility of 5950 cm2/(V$\cdot$s) and an estimated surface roughness and correlation length of 2.2 Å and 24 Å, respectively, from the falling edge dependence Mazzoni et al. (1999).

The electron density in the leads of the tunnel barrier will differ from that in the field due to two effects: fringing field and lateral quantum confinement. We calculated electrostatic modifications of threshold due to fringing fields numerically and were found to be a small, less than 10% effect. The effect of quantum confinement is more significant Lu et al. (2016). We calculated the confinement threshold shift for this geometry over various wire widths, see Fig. 1(g) and Appendix B. We measure a threshold of 0.60$\pm$0.05 V for the nanostructure (i.e., $V_{R}$ = 2.13 V, $V_{C}$ = 5.00 V, and $V_{L}$ is increased). This threshold is the linear extrapolation to zero current, Fig. 1(f). The threshold shift falls within a range that would be expected for the $\sim$70 nm width of the tunnel barrier lead. The electron density in the leads can, therefore, be estimated using the field threshold combined with an offset from quantum confinement.

III Theoretical model

We now examine how well a 1D barrier model and a voltage dependent parameterization of barrier height and width fits the observed tunneling. At zero source-drain bias, we assume a rectangular barrier with a width $w$, a left barrier height $U_{0}$, and a right barrier height $U_{1}$. The two barrier heights allow for asymmetric plateaus in the leads. When a source-drain bias $V_{SD}$ is applied, the chemical potentials or quasi-Fermi levels in the source $\mu_{S}$ and drain $\mu_{D}$ separate by an amount proportional to $V_{SD}$. The potential in the barrier region then varies linearly with a slope equal to $–F$, where $F$ is the electric field due to the source-drain bias, as shown in Fig. 1(c). Thus, for non-zero source-drain bias, the wave-functions in the barrier region are Airy functions. In the leads the wave-functions are propagating plane waves with an incident wave $\psi_{i}$, a reflected wave $\psi_{r}$, and a transmitted wave $\psi_{t}$.

We use the Landauer-Büttiker formalism to model the tunneling transport. Assuming low temperature, the Fermi functions in the source and drain leads can be approximated as step functions and the current through the barrier is given by Datta (1997)

[TABLE]

Here, $T(E)$ is the transmission coefficient of an electron with energy $E$, $M(E)$ is the total number of available transverse modes which depends on the DOS in the leads, and $\mu_{0}=\max\left(\mu_{D},0\right)$. Thus, the product $T(E)M(E)$ gives the total transmission summed over the transverse modes.

We examine both a numerical and approximate analytic solution to computing the transmission. Both are sufficiently efficient to find good fits for barrier height and width for each gate voltage. The numerical approach uses the analytical, piecewise solution to Schrödinger’s equation and then numerically solves the boundary matching problem as a system of equations. This approach gives us an exceptionally fast forward solve (contrasted against fully numerical approaches) to facilitate non-linear inference over a large experimental data set. For more details, see Appendix C.

Reasonably good quantitative agreement can also be provided through solving the transmission coefficient by approximating the linearly varying potential in the barrier region with two potential steps, see Fig. 1(c). One step has a fixed height of $U_{0}$ while the other has a variable step height that depends on $V_{SD}$. The transmission coefficient is then solved using transfer matrices and yields a simple, analytic formula for transmission that is a good approximation to the trapezoidal barrier problem over a wide parameter range. More details are provided in Appendix D. We find that both 1D models can match the measured current dependences on voltage, Fig. 2(a). We further find that the three-step model agrees well with the full numeric solution for lower $V_{SD}$ relative to the barrier $U_{0}$, Fig. 2(b-c). For the rest of this paper we will use the numerical solution to the trapezoidal model for the more accurate quantitative analysis.

We note that the three-step and trapezoidal barrier models are useful alternatives to traditional approximate models of tunneling phenomena, such as the WKB approximation and Fowler-Nordheim tunneling. Our approach avoids unnecessary (and often unjustified) assumptions by exactly solving Schrödinger’s equation and more accurately accounting for the lead geometry. Over the course of a numerical fit, the assumptions of these approximations can easily become violated, invalidating the parameter extraction. Both our numerical solution to the trapezoidal barrier and the three-step model are simple enough to allow rapid estimates of the barrier height and width from transport measurements using a modern computer.

IV Results and Discussion

Figure 3(a) shows line-cuts at $V_{L}=$ 0.90, 0.80, 0.70 and 0.60 V from the data in Fig. 1(d). We treat the barrier height $U_{0}$ and width $w$ as free parameters and fit the trapezoidal barrier model to the data by minimizing a chi-squared test statistic, resulting in maximum likelihood estimates (MLEs) of the parameters under the assumption of independent and identically distributed, Gaussian errors. These fits are shown in Fig. 3(a). The MLEs and expected errors for the parameters are shown in Table 1.

The expected errors were computed by constructing the profile likelihood function (separately, for each parameter), and numerically determining a confidence level of $95\%$. To do this, we make the standard assumption that the difference in the profile likelihood function from the MLE value should be $\chi^{2}(1)$ distributed. Unless otherwise noted, we assume throughout this paper that the experimental standard error in the observed currents is $\delta I=0.1\cdot\text{max}\left(I_{SD}\right)$, where $\text{max}\left(I_{SD}\right)$ is the maximum source-drain current for a given gate voltage. From Table 1, we see that the barrier width seems to be fixed at about 12 nm while the barrier height and its uncertainty increase as $V_{L}$ decreases.

We repeat the fitting procedure for different $V_{L}$ values. Fig. 3(b) shows a plot of the extracted barrier height as a function of $V_{L}$ for a fixed width of 12 nm and three values of $V_{C}$. In this case we assume $\delta I=0.3\cdot\text{max}\left(I_{SD}\right)$ because the errors were found to be smaller for a one-parameter fit. We find that at high $V_{L}$ the barrier height varies linearly with gate voltage. However for low $V_{L}$ the barrier height increases non-linearly.

Overall, we observe, that the dependence fits an exponential function of the form $f(x)=a\mathrm{e}^{-bx}+c$. In Fig. 3(c), we compare the exponent $b$ as a function of $V_{C}$, finding that it varies between 4.1 V*-1* to 3.7 V*-1* as $V_{C}$ changes from 1.60 V to 2.40 V. The average exponent is 3.9 V*-1*. This non-linear increase in barrier height approximately coincides with the regime for which charge density is changing exponentially as the electron channel is forming under the gate.

The width of the barrier is plotted as a function of $V_{L}$ for $V_{C}=1.60$ V in Fig. 3(d). In this case we show a two-parameter fit where the barrier height is not fixed. We see that the width appears to be a constant over the range of $V_{L}$ values considered. The average value is about 12 nm. We note that the extracted widths correspond well with our gate pattern.

In the context of compact models, the overall barrier height and width dependences on $V_{L}$ and $V_{C}$ have relatively simple forms that describe the behavior over very wide bias ranges.

In addition to profile likelihood error bars, we also compute full 2D confidence regions for the parameters $U_{0}$ and $w$. Similar to the 1D case, we assume that the difference between the MLE point and the likelihood function is distributed as $\chi^{2}(2)$. Fig. 3(e) shows the result for $V_{L}$ = 0.70 V, where the contours represent different confidence levels.

Besides performing the parameter extractions above, we consider cases where the DOS in the leads is 0D, 1D, or 2D to examine if there are any distinguishing signatures in the calculated dependences. The difference in DOS affects the total number of available modes $M(E)$ in Eq. 3. All three cases fit the data well when the barrier heights and widths are allowed to adjust to compensate the change in the DOS. We note that the estimated widths of the leads would be consistent with transversal quantization and the energy level splitting is too large to warrant 1D DOS for the transmission modes. We therefore use the 0D DOS for the extracted barrier heights and widths. More generally, it is unclear at this time what factors in the tunnel barrier geometry express clearer signatures due to the DOS dimensionality. This is a topic for future examination.

We lastly compare the results of the 1D barrier model to electrostatic simulations of the electron density and electric potential (i.e., the conduction band edge, $E_{C}$). The electrostatics of the device is modeled using COMSOL Multiphysics software. We use the Thomas-Fermi approximation to model the electric potential and the electron density at the oxide-semiconductor interface with a 2D density of states for the experimental gate voltages. The electron density is calibrated at a gate voltage of 0.48 V by applying an offset voltage to all electrodes while ramping gate L. The offset voltage is used to ensure that the electron density in the leads agrees with the definition of threshold given by Eq. 2.

Plots of the conduction band energy as a function of $V_{L}$ and position are shown in Fig. 4(a). As $V_{L}$ is decreased a single barrier forms between gates L and C. The barrier height is determined by calculating the saddle point in the potential, which is found by calculating the gradient of the potential. We note that the shape of the barrier will generally not be trapezoidal for a 1D cut of the conduction band energy along a curve through the saddle point. The barrier height as a function of $V_{L}$ for $V_{C}=1.60$ V, $V_{R}=2.13$ V, and an offset voltage of 0.50 V is shown in Fig. 4(b-c).

In the constriction where the barrier forms we expect quantum confinement to have a significant effect. These transverse confinement energies were computed for this device geometry and become appreciable at the lower densities where the channel is narrow, see Appendix B. The effect of lateral confinement is shown more clearly in the inset of Fig. 4(b). This produces a region of relatively linear dependence of the barrier height on voltage after the electron channel is completely formed but the gap is not completely saturated with electron density. A qualitatively similar region is seen in the extracted barrier heights in Fig. 3(b), although extended over a wider voltage range. We note that better quantitative agreement is beyond the scope of this paper and probably requires self-consistent Schrödinger-Poisson calculations.

Fig. 4(c) shows the barrier height as a function of $V_{L}$ for a variety of offset voltages (i.e., threshold offsets for the COMSOL calculation). Each curve shows the range from no accumulation (large barrier) to a non-existent barrier (barrier height is zero or lower).

V Conclusion

In conclusion, we have measured and modeled the transport spectroscopy of silicon MOS surface electrode-defined tunnel barriers at 4 K. Multiple electrodes are used to form the electrostatically defined barriers making for a large parameter space on which the tunnel barrier is dependent.

We examine a 1D barrier model using voltage dependent barrier height and width, which is found to fit the current-voltage dependence well over a large parameter range. The model is quasi-analytic providing both accurate 1D transmission values for the barrier model, while also being a fast extraction method to more readily enable its application to multi-electrode cases. Cryogenic and quantum confinement effects are included in the model to account for threshold and barrier height shifts.

The barrier height dependence on gate voltage was shown to be relatively linear in the high and low voltage regimes with an intermediate non-linear regime that is nearly exponential. The regimes and their voltage dependences can be fit in a compact way offering a potential approach to modeling more complex nanostructures that have multiple tunnel barriers and more electrodes.

We also compare the results of the 1D model to semi-classical solutions of Poisson’s equation using COMSOL and find that the simulations qualitatively agree well with the model predicted three regimes.

Transport spectroscopy, itself, is a relatively rapid way to characterize tunnel barriers relative to more time intensive and complex pulsing approaches used for these kinds of barriers Elzerman et al. (2004); Amasha et al. (2008); Thalakulam et al. (2011). The combination of this method with transport spectroscopy offers a relatively rapid way to build a compact model of a tunnel barrier for multiple gate electrodes and reasonably large bias ranges. This work should provide useful insight about the details of electrostatic barriers and how to characterize them.

Acknowledgements.

We gratefully recognize conversations with D.R. Ward and P.A. Sharma about this work and the manuscript. This work was performed, in part, at the Center for Integrated Nanotechnologies, a U.S. DOE, Office of Basic Energy Sciences user facility. The work was supported by the Sandia National Laboratories Directed Research and Development Program. Sandia National Laboratories is a multi-program laboratory operated by Sandia Corporation, a Lockheed-Martin Company, for the U. S. Department of Energy under Contract No. DE-AC04-94AL85000.

Appendix A Cryogenic threshold voltage calculation

We calculate the threshold voltage for a MOSFET at low temperatures by considering (1) the metal-insulator transition and (2) the consequences of low temperature on parameters used to calculate threshold, such as the metal-semiconductor work function, $\phi_{ms}$. First consider the metal-insulator transition.

One definition of threshold voltage for MOSFETs at high temperatures (such as room temperature) is that it is the gate voltage for which the mobile electron charge $Q_{n}=0$. As noted in the text, at cryogenic temperatures there is a metal-insulator transition where one must achieve a critical electron density $n_{crit}$ before appreciable conduction begins. Thus at low temperatures the threshold voltage is defined as the gate voltage for which $Q_{n}=-en_{crit}$. The expression for the threshold voltage is then given by

[TABLE]

In this expression $V_{0}$ is the standard MOSFET threshold voltage (discussed below) and $V_{crit}=en_{crit}/C_{ox}$, where $C_{ox}$ is the oxide capacitance per unit area. The critical density can be extracted from Hall measurements Tracy et al. (2009).

We next examine $V_{0}$. The usual formula for the threshold voltage, neglecting the body effect (that is, zero voltage difference between the source and bulk) and assuming a p-type substrate, is given by Muller and Kamins (2002)

[TABLE]

where $V_{FB}$ is the flat-band voltage, $\phi_{p}$ is the bulk potential (difference between the intrinsic and quasi-Fermi levels in the bulk of the semiconductor) and the third term is the voltage across the oxide due to the depletion layer charge. The flat-band voltage is given by $V_{FB}=\phi_{ms}-Q_{f}/C_{ox}$. The electrical permittivity of the semiconductor is $\epsilon_{s}$ and the acceptor concentration is $N_{a}$.

The parameters $\phi_{ms}$ and $\phi_{p}$ are most affected by changes in temperature. At cryogenic temperatures, the quasi-Fermi level in the semiconductor is pinned halfway between the acceptor energy levels and the top of the valence band due to freeze-out, as shown in Fig. 5. For boron, the acceptor energy levels are about 0.045 eV above $E_{V}$. It is well known that the band-gap energy $E_{g}$ changes with temperature; for Si, the band-gap energy is 1.17 eV at low temperatures. Thus we find $\phi_{p}=E_{g}/2-0.045/2=0.563$ V.

The parameter $\phi_{p}$ also affects $\phi_{ms}$ through the semiconductor work function $\phi_{s}$. Unlike $E_{g}$, the effect of low temperatures on the metal work function $\phi_{m}$ and the electron affinity $\chi$ is not clear. We use $\phi_{m}=4.05$ V for n*+* poly-silicon and assume $\chi=4.05$ V for Si. This gives $\phi_{s}=\chi+E_{g}/2+\phi_{p}=5.20$ V and a metal-semiconductor work function of $\phi_{ms}=\phi_{m}-\phi_{s}=-1.15$ V.

Appendix B Lateral quantum confinement

In a simple square potential well or particle in a box, the discrete energy levels increase as the width of the quantum well is reduced. Similarly, quantum confinement in a nanowire causes the lowest sub-band to increase in energy. In this appendix, we consider the effect of lateral quantum confinement on the threshold voltage. This effect causes a shift in the threshold voltage resulting in the threshold of a thin wire gate to be larger than the threshold of a much wider gate (field threshold). This theory was recently developed for nanowires in Si/SiGe heterostructures Lu et al. (2016); here, we apply it to the MOS geometry shown in Fig. 6(a).

To understand this effect, we consider a range of wire gate voltages and wire widths. For each wire width, we compute the electrostatic potential for each corresponding wire gate voltage by solving Poisson’s equation in 2D with the finite element method in COMSOL Multiphysics. We then calculate the ground state energy of the resulting confinement potential by solving Schrödinger’s equation. This gives us the energy of the lowest sub-band as a function of voltage. We find that the ground state energy is linear with gate voltage, that is,

[TABLE]

where $E_{0}^{L}$ is the ground state energy of a wire of width $L$, $m_{L}$ is the slope or lever arm, and $V_{G}$ is the wire gate voltage. As the wire width increases, the slopes $m_{L}$ saturate. We then define the lever arm for an infinitely wide wire as approximately equal to the lever arm for a 5000 nm wide wire, that is, $m_{\infty}\approx m_{5000}=$ -0.91 meV/mV for the MOS geometry shown in Fig. 6(a).

Assuming that conduction occurs when the ground state energy is aligned with some external reference energy level, $E^{0}$ (for example, the quasi-Fermi level), we define the field threshold voltage as,

[TABLE]

where $V_{\infty}$ is the threshold for an infinitely wide wire (field threshold).

If the threshold voltage of a wire of width $L$ is $V_{L}$ and we define the threshold shift as $\Delta V_{L}=V_{L}-V_{\infty}$, then we may write

[TABLE]

We can re-write the expression for the threshold shift explicitly in terms of the field threshold using Eq. 7 as,

[TABLE]

Moreover, using the definition $\Delta V_{L}=V_{L}-V_{\infty}$, we can write a direct expression for the wire threshold as,

[TABLE]

We find that the lever arm $m_{L}$ as a function of wire width $L$ is fit well by a decaying exponential function, that is, $m_{L}\approx m_{\infty}+ae^{-bL}$, where $a=0.32$ meV/mV and $b=0.003$ nm*-1*. The expression for the wire threshold Eq. 10 then becomes,

[TABLE]

Thus we can calculate the wire threshold given a field threshold voltage and wire width using this compact model. The wire threshold as a function of wire width for different field thresholds is shown in Fig. 6(b). We note this assumes that lateral quantum confinement is the dominant effect and the exact increase due to confinement is geometry dependent.

Appendix C Transport formalism for a general trapezoidal barrier

In this appendix, we present a detailed calculation of tunneling through a “1D” trapezoidal barrier. We will use the Landauer-Büttiker formalism, computing the transmission function via a semi-analytical approach.

The Landauer-Büttiker formalism is a general framework for computing tunneling quantum transport through a device. It can handle devices with many leads and high source-drain biases. The method breaks down when transport becomes non-ballistic, through “vertical” scattering events. This occurs when energy is not conserved throughout the scattering region; in this instance, more sophisticated techniques, such as non-equilibrium Green’s functions, must be employed.

The case we consider here is a simple, two-terminal device. Within the Landauer-Büttiker formalism, the total current through the scattering region is given by

[TABLE]

where $\bar{T}(E)$ is the total transmission function at a given energy (summed up over transverse modes), $f_{s}$ is the Fermi function of the source, and $f_{d}$ is the Fermi function of the drain. The difficult part of computing $I$ is to calculate the total transmission function $\bar{T}(E)$. To facilitate this, we decompose $\bar{T}$ as

[TABLE]

where $T_{m}(E)$ is the transmission of a given transverse mode, $M(E)$ is the total number of accessible modes at a given energy $E$. Hence, we have decomposed the problem into computing $T(E)$ and $M(E)$ indpendently. First, we examine the computation of $T(E)$.

C.1 Computing the transmission function

The transmission function gives the probability of an incident mode tunneling through the barrier. Schematically, we consider the case shown in Fig. 7. There, the transmission function is a function of the energy $E$ and is parameterized by the left-side barrier height $U_{0}$, the field $F$ between the source and drain, the barrier width $w$, and the right-side barrier height $U_{1}$.

C.1.1 Form of the Schrödinger equation

To compute the transmission function, we need to solve Schrödinger’s equation for an incoming plane wave at a specified energy. Within each region, we can write down the solution analytically:

[TABLE]

Here, we have $k_{1}=\frac{\sqrt{2mE}}{\hbar}$, $k_{3}=\frac{\sqrt{2m(E-(U0-U1-Fw))}}{\hbar}$, and $a(x)=\left[\frac{2mF}{\hbar^{2}}\right]^{1/3}\left[\frac{U_{0}-E}{F}-x\right]$. The transmission function is given by $T=\left|\frac{A_{3}}{A_{1}}\right|^{2}$. To solve for this quantity, we need to impose boundary conditions: $\psi_{1}(0)=\psi_{2}(0)$, $\psi_{1}^{\prime}(0)=\psi_{2}^{\prime}(0)$, $\psi_{2}(w)=\psi_{3}(w)$, $\psi_{2}^{\prime}(w)=\psi_{3}^{\prime}(w)$. To satisfy the boundary conditions, we will also need to know the formulas for the derivatives:

[TABLE]

C.1.2 Renormalizing the Airy functions

One issue that arises immediately is that Airy functions that form the solution in region 2 are not well-defined when $F\rightarrow 0$. This is a problem, since we want to recover the original square barrier problem in this case. Since the argument of the Airy functions is real, we can use the asymptotic form for $F\rightarrow 0$:

[TABLE]

These expressions are valid for real, large-magnitude $z$. We want to use this to work out expressions for the boundary conditions of $\psi_{2}(x)$. We will use these to renormalize the solution coefficients so that we don’t have undefined expressions. Numerically, we find that we need only go to $N=0$ to obtain good accuracy.

We write down the asympotic result at $x=0$:

[TABLE]

Using this, we rewrite the $\psi_{2}$ solution as:

[TABLE]

where the new functions $\mathcal{A}i$ and $\mathcal{B}i$ have the property that

[TABLE]

which are of course just the growing and decaying exponentials that span the solution for a flat barrier.

Hence, this renormalization allows us to make contact with the simple limiting case. We thus replace our region 2 wave function and derivative with following piecewise definitions:

[TABLE]

Here, we just relabeled the unknown coefficients $A_{2}$ and $B_{2}$ that we are solving for to include the renormalization. We identify a good value of a threshold field $F_{0}$ numerically.

C.1.3 Constructing the system of equations

Now that we have worked out the wave functions for all three regions, we are ready to compute the transmission function. However, notice that we have five unknown coefficients with only four constraints. This is because the full wave function needs to be normalized, which provides the fifth constraint. However, the transmission function only needs a ratio of coefficients, so we don’t need to worry about this additional constraint. Specifically, we divide through everywhere by $A_{1}$, giving:

[TABLE]

Here, the tilde variables indicate scaling by $A_{1}$, so $T=\left|\tilde{A_{3}}\right|^{2}k_{3}/k_{1}$, where the last factor comes from the difference in velocity between incoming and outgoing modes. We will now drop the tildes for convenience.

Considering first the case where the Airy functions are well-defined, we have the following system of equations to solve:

[TABLE]

If instead the Airy functions are not well defined and we have to use the asymptotic forms, we have:

[TABLE]

We cast these equations as a simple matrix problem solve for $A_{3}$.

C.2 Computing the total current

Now that we can reliably compute the transmission function, we need to use it to compute the total current within the Landauer-Büttiker formalism. To do this, we modify the schematic slightly to Fig. 8.

C.2.1 Density of transverse modes

The total transmission function $\bar{T}(E)$ is given by

[TABLE]

where $T_{m}(E)$ is the transmission of a given transverse mode, $M(E)$ is the total number of accessible modes at a given energy $E$, $g$ is the density of states, $L$ is the cross-sectional length scale of the scattering region, and $d$ is the dimension of the lead. The three cases of interest are for 0D (one mode), 1D (sheet contacts), and 2D (volume contacts) density of states, given by:

[TABLE]

where we have multiplied the usual formulas for the 1D and 2D DOS expressions by 2 for the extra valley degeneracy in silicon. Using these density of states formulas, we compute $M(E)$:

[TABLE]

C.2.2 Computing the current

Equipped with the number of modes, we can compute current:

[TABLE]

where

[TABLE]

The two cases here correspond to the case 1 and case 2 in Fig. 8. From the formulas above, we can write three different current models:

[TABLE]

where we assume that for $I_{0}$ the single mode is in the transport window.

Appendix D Three-step model

In the three-step barrier model we approximate a square barrier under bias using three step potentials as shown in Fig. 9. For $V_{SD}=0$, the barrier reduces to a square barrier of height $U_{0}$ and width $w$. The first step occurs at $x=-w/2$ and has a fixed height of $U_{0}$. The next step occurs at $x=0$ and has a variable step height that depends on the source-drain bias $V_{SD}$. The final step occurs at $x=w/2$ and also has a variable height that depends on $V_{SD}$. We use the transfer or $T$-matrix formalism to compute the transmission coefficient Davies (1997); Ando and Itoh (1987). For simplicity we only consider positive $V_{SD}$ but this approach can be easily applied to negative $V_{SD}$ as well.

First consider the case where $E>U_{0}$. There are four wavenumbers, each corresponding to a particular region of the barrier as indicated by the roman numerals in Fig. 9. The wavenumbers are given as follows:

[TABLE]

Here $m$ is the effective mass of the electron. Using the $T$-matrix formalism we find the transmission coefficient to be:

[TABLE]

where the coefficients $f_{1}$, $f_{2}$, $f_{3}$, and $f_{4}$ are functions of the wavenumbers given by

[TABLE]

The arguments of the trigonometric functions depend on the width of the barrier as $\theta=k_{2}w/2$ and $\phi=k_{3}w/2$. This is the general form of the transmission coefficient.

Now consider the large $V_{SD}$ case such that $V_{SD}>U_{0}-E$ and $U_{0}>E$. The form of the transmission coefficient and the wavenumbers remain the same except the wavenumber in region II is replaced in the following manner: $k_{2}\rightarrow i\kappa_{2}$. Then the trigonometric functions that have $\theta$ as the argument become hyperbolic functions, that is $\sin k_{2}w/2\rightarrow i\sinh\kappa_{2}w/2$ and $\cos k_{2}w/2\rightarrow\cosh\kappa_{2}w/2$. The wavenumber in region II is now given by $\kappa_{2}^{2}=2m\left(U_{0}-E\right)/\hbar^{2}$.

Then consider the small $V_{SD}$ case such that $V_{SD}<U_{0}-E$ and $U_{0}>E$, as shown in Fig. 9. Again the form of the transmission coefficient remains the same except now the wavenumbers in both region II and III are replaced: $k_{2}\rightarrow i\kappa_{2}$ and $k_{3}\rightarrow i\kappa_{3}$. The wavenumbers in region II and III are now given by $\kappa_{2}^{2}=2m\left(U_{0}-E\right)/\hbar^{2}$ and $\kappa_{3}^{2}=2m\left(U_{0}-E-eV_{SD}\right)/\hbar^{2}$ . Now all the trigonometric functions become hyperbolic functions as discussed above.

Finally Eq. 47 and the appropriate wavenumbers (based on the values of $V_{SD}$, $U_{0}$, and $E$) are used in Eq. 12 to compute the current through the device.