Conductivity and size quantization effects in semiconductor \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\delta$$\end{document}δ-layer systems

We present an open-system quantum-mechanical 3D real-space study of the conduction band structure and conductive properties of two semiconductor systems, interesting for their beyond-Moore and quantum computing applications: phosphorus \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\delta$$\end{document}δ-layers and P \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\delta$$\end{document}δ-layer tunnel junctions in silicon. In order to evaluate size quantization effects on the conductivity, we consider two principal cases: nanoscale finite-width structures, used in transistors, and infinitely-wide structures, electrical properties of which are typically known experimentally. For devices widths \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$W<10$$\end{document}W<10 nm, quantization effects are strong and it is shown that the number of propagating modes determines not only the conductivity, but the distinctive spatial distribution of the current-carrying electron states. For \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$W>10$$\end{document}W>10 nm, the quantization effects practically vanish and the conductivity tends to the infinitely-wide device values. For tunnel junctions, two distinct conductivity regimes are predicted due to the strong conduction band quantization.


I. INTRODUCTION
Highly conductive δ -layer systems, i.e. thin, high-density layers of dopants in semiconductors are actively used as a platform for exploration of the future quantum and classical computing when patterned in plane with atomic precision [1][2][3].Such structures, with the dopant densities above to the solid solubility limit [4], have been shown to possess very high current densities [3,5] and thus have a strong potential for quantum computing applications [6,7] and advanced microelectronic devices [8,9].However, at the scale important for these applications [10], i.e. devices with sub-20 nm physical gate/channel lengths and/or sub-20 nm widths, that could compete with the future CMOS, the conductive properties of such systems are expected to exhibit a strong influence of size quantization effects.At the same time, experimental assessment of the conductivity of δ -layer systems is typically performed using Hall measurements on samples of macroscopic dimensions (> 1 µm) [5,[11][12][13].
The electronic structure and conductive properties of Si:P δ -layer systems have been a subject of previous studies based on either effective mass [14][15][16][17], tight-binding [18][19][20][21][22], den-FIG.1. Geometry of the Si:P δ -layer tunnel junction.Our device consists of a semi-infinite source and drain, in contact with the engineering channel.The channel is composed of a lightly doped Si body and Si cap and a very thin, highly P doped-layer with an intrinsic gap of length L gap .sity functional theory [23][24][25] formalisms or semiclassical Boltzmann theory [26].Recently it has been demonstrated in [15][16][17] that to accurately extract the conductive properties of highly-conductive, highly-confined systems, an open-system quantum-mechanical analysis is necessary.Such open-system treatment, that can be conducted for instance using the Non-Equilibrium Green's Function (NEGF) formalism [27,28], allows to compute the current and conductivity directly from the quantum-mechanical flux, thus avoiding semi-classical approximations, which are intrinsic to the traditional charge self-consistent closed-system or periodic boundary conditions band-structure calculation methods.It has been also demonstrated in [16] that an open-system treatment in highly P doped δ -layer in Si: i) permits to reveal the quantization in space and energy of the free electrons around the δ -layer; ii) permits to explain the existence of the shallow 3Γ sub-band, which has been observed experimentally [29], without any fitting parameters; iii) predicts significant quantum-mechanical dependence of the current on the δ -layer sheet thickness for a fixed dopant sheet density; and iv) provides a very good agreement with the experimental electrical measurements [5,[11][12][13]].
An accurate computational description of electron tunneling in semiconductor δ -layer tunnel junctions (such the one shown in Fig. 1) is additionally required because the tunneling rate at a δ -layer junction is affected not only by the gap length and the conductivity of the δ -layers, but also by quantization of the conduction electrons in energy and space [16].In this work we employ an efficient computational open-system quantum-mechanical treatment to explore the conductive band structure and the conductive properties of phosphorus δ -layer systems in silicon (Si:P δ -layer) for device widths, from nanoscale (< 20 nm) to macro-scale (> 1 µm) dimensions, and to analyze the influence of size quantization effects on the conductive properties for sub-12 nm device widths.All simulations are carried at the cryogenic temperature of 4K, in which we can neglect inelastic scattering events [5,30].The analysis of the influence of different kinds of non-idealities on the tunneling current is presented in [31].The open-system treatment is based on an application of Keldysh formalism [27], known as NEGF [28], and the effective mass theory.Generally for Si systems, the use of the effective mass approximation has been shown to be in a good agreement with tight-binding models for nanowire diameters down to 3 nm [32,33]; thus, the fidelity of this approximation start to decline for very narrow systems (< 3 nm).Here we have employed an efficient implementation of NEGF, refereed to as the Contact Block Reduction (CBR) method [34][35][36][37][38], which allows accurate computation of all quantum-mechanical quantities of interest (local density of states, transmission probability, current) and scales linearly O(N) with the system size.We find that at the scale most interesting for applications, i.e. for device widths W < 10 nm, quantization effects strongly affect both the conductivity and the spatial distribution of the current-carrying electron states, which, similarly to the charge distribution in the hydrogen atom, is determined by a "quantum number", i.e. the number of propagating modes.This strong spacial quantization of the current-carrying states can be utilized in novel electronic δ -layer devices.Conversely, for W > 10 nm, the quantization effects practically vanish and the conductivity tend to the values of infinitely-wide devices.Finally, two distinct conductivity regimes are predicted due to the strong conduction band quantization for tunnel junctions.

A. Effects of the device width on the conductive properties
The conductivity of infinitely-wide Si:P δ -layers has been first studied in [15][16][17] using the open-system quantummechanical approach for infinite-width systems (see Sect.III G).Here we apply the open-system treatment to more realistic/utilitarian devices of finite width and/or with an intrin-sic or lightly-doped tunnel gap as shown in Fig. 1.The device consists of a semi-infinite source and drain (represented by the NEGF open boundary conditions), in contact with a channel of length L, which is composed of a lightly doped Si body and Si cap and a very thin, highly P doped-layer with an intrinsic gap of length L gap , as shown in Fig. 1.In this work we consider tunnel junction devices with an intrinsic gap of length L gap , ranging from 0 nm to 12 nm, and width W , ranging from 2 nm to infinity, paying a particular attention to nano-scale dimensions (W, L gap < 12 nm) which are most interesting for quantum and beyond-CMOS computing applications.The length of the channel is assumed to be L = 30 nm + L gap to avoid the boundary effects, between the source and drain contacts, and the height of the device is chosen to be H = 12 nm.We also assume δ -layer thicknesses raging from monoatomic layer, t = 0.2 nm, to few atomic layers, t = 1.0 nm, a doping density of the δ -layer of N D = 1.0 × 10 14 cm −2 and an acceptor doping densities in the Si cap and Si body of N A = 5.0 × 10 17 cm −3 .Note that the doping density in the δ -layer is given in cm −2 (i.e.sheet doping density) to be consistent with experiments' nomenclature: D , where t is the δ -layer thickness, N (3D) D is the doping density in cm −3 and N (2D) D in cm −2 .Fig. 2 shows the computed conductivity in function of the device width W and for different gap lengths L gap .The dashed lines represent the conductivity values for infinitely wide devices, W → ∞.Details of the computational treatment for infinitely-wide systems are presented in Sect.III G.The simulations suggest that quantization due to the finite size of the device width starts to appear for widths below to 7-10 nm.This size quantization is reflected as non-monotonic increase of the conductivity.Interestingly, there is a peak in the conductivity at W = 5 nm and a dip at W = 3 nm.Conversely, as W increases, the conductivity tends to the values for infinitely wide devices that therefore can be seen as lower bound limits.Additionally, the size quantization effect on the conductivity of δ -layer tunnel structures is most notable for large tunnel gaps L gap > 7 nm.Indeed, in Fig. 2, one can see that the deviation from the conductivity values for infinitely wide devices is more significant for large tunnel gaps.As we will discussed later in Section II C, this is due to the consequence of the quantization of the low-energy conduction band in δ -layer systems and the wave-functions decoupling between the left and right δ -layers in large tunnel gaps.
The oscillations of the conductivity for narrow device widths W arise due to a small number of the propagating modes in a "waveguide", created by the finite-size width of the δ -layer.The corresponding dependence of the current on the device width is shown in Fig. 3 a for a gapless δ -layer structure L gap = 0.The existence of the conduction steps due to each new propagating mode is well known experimentally since 1980's [39].Here we show, however, that in highlyconfined, highly-conductive δ -layer systems, the quantum number m, representing the number of propagating modes, determines not just the total current, but also the spatial distribution of the corresponding current-carrying electrons.The total number of propagating modes m is determined by the number of peaks in the density of states (DOS) below the Fermi level FIG. 3. Propagation modes for Si: P δ -layer systems.a Current I vs device width W for δ -layer systems with L gap = 0 nm: the insets in blue color show the spatial distributions of current-carrying modes across a y-z plane, indicating the corresponding number of propagating modes m; Inset in green color shows the total electron density that includes all (not just current-carrying) occupied electron states for a device width of W = 12 nm (Note that n total ∼ 10 20 cm −3 n curr.−carr.∼ 10 16 cm −3 ).b Detail of the spatial distribution of current-carrying states n curr.−carr.(y, z) for a device width of W = 9 nm .For all calculations, N D = 1.0 × 10 14 cm −2 , N A = 5.0 × 10 17 cm −3 , t = 0.2 nm and an applied voltage of 1 mV.
or the number of steps in the electronic transmission function below the Fermi level.
The spatial distribution of the current-carrying electron states, n curr.−carr.(y, z), can be obtained by performing the energy integration of the local density of electron states (LDOS) weighted by the corresponding current spectrum i e (E) as: n curr.−carr.(y, z) = LDOS(y, z, E)i e (E)dE/ i e (E)dE.The current spectrum and the local density of states can be obtained by expression in Eqs. 1 and 7, respectively.The spatial current-carrying electrons for the different modes is shown in Fig. 3 a as insets in blue color, as well as the corresponding number of propagating modes.Additionally, the total electron density is also included in the figure as an inset in green color, demonstrating only weak spatial quantization along the ydirection.However, the specific portion of electrons with energies close to the Fermi level, i.e. the current-carrying states, do exhibit a strong spatial quantization.Indeed, for m = 1 the propagating mode reaches the maximum concentration at the center of the structure, the mode the corresponds to m = 2 is "excited" into the further penetration along the confinement direction (z-axis), leaving the center relatively depopulated (in terms of the current-carrying states), the mode m = 3 is again "pushed out" of the center along both z-and y-axis.One can further note that the modes m = 1, and m = 4, 6, etc. tend to form a regular "phase" distribution of the current-carrying states (i.e. the states distributed closer to the center of the δlayer along z-axis), while the modes m = 2, and m = 3, 5, etc. form "anti-phase" distributions (i.e. the states distributed further from the center of the δ -layer along z-axis) that have the maximum current being carried in the different regions of space, separated by a few nanometers.When W → ∞, the number of propagation modes in y-direction becomes infinite m → ∞ as expected.Fig. 3 b shows the density of currentcarrying electron states for the anti-phase case of δ -layer doping N D = 10 14 cm −2 , thickness t = 0.2 nm (to approximate a monoatomic δ -layer) and device width W = 9 nm, which corresponds to m = 5.As final remark, the number of propagating modes m in δ -layer structures is mainly determined by three factors: 1) the δ -layer doping level N D , 2) the the δ -layer doping thickness t and 3) the device width W .

B. Peculiarities of finite-width vs infinitely-wide devices
Next it is to get an insight into the transition from infinitewidth to finite-width devices, before the quantum effects arise on the conductivity properties, i.e. the transition marked with a black arrow in the right upper part of Fig. 3.In the following, therefore, we will compare the free electron energy distribution, the transmission function and the current spectrum for devices of 12 nm-width, in which the size quantum effects on the conductivity are minimum as was seen in Fig. 2, and the corresponding infinitely wide counterparts.All results of infinitely wide devices are scaled to an effective device width of 12 nm for comparison purposes.
Fig. 4 shows the free electron energy distribution of δ -layer structures of 12 nm-width (in continuous lines), together with the result of infinitly wide δ -layer structure with L gap = 0 nm (in dashed line), for different gap lengths L gap .The free electron energy distribution corresponds to the total density of states multiplied by the the Fermi-Dirac distribution func- tion, which determines the probability of a state to be occupied.First we notice that in contrast to the infinite-width δ -layer systems, the finite size along the y-direction manifests itself in splitting of the occupied 1Γ sub-band, that is located around -0.18 eV, and 2Γ sub-band, that is located around -0.033 eV, into additional modes.Interestingly, for tunnel junctions (L gap > 0), two of 2Γ modes are immediately dampened as manifested in the reduction of the corresponding peaks compared to the gapless case as shown in Fig. 4. In addition, as it is evident from Fig. 4, the free electron distribution in tunnel junctions is almost independent of the tunnel gap length L gap .
Fig. 5 shows the transmission function of infinite-width devices (in dashed lines) and 12 nm-width devices (in continuous lines) for different tunnel gap lengths L gap .The transmission function provides the sum of the probabilities for each mode of a carrier at certain given energy E to carry current from the source to drain.As can be seen, the transmission function is fairly similar for both systems, infinite and finite width, especially for wider tunnel gaps.The transmission function is reduced exponentially for the low-energy modes with the length of the tunnel gap.Additionally, the transmission function strongly dependent on the gap length: the energy window in which the carrier can be transmitted is quickly reduced with the increase of the tunnel gap length.
Fig. 6 shows the current spectrum, i e (E): in dashed lines, for infinite-width devices and, in continuous lines, for 12 nmwidth devices.It is evident that the current spectrum of infinite-width devices differs from the current spectrum of finite-width counterparts.The finite size along y-axis significantly affects the current spectrum limiting the currentcarrying states in all cases, including the gapless case, to only a narrow energy window around the Fermi level, with the size of the energy window being proportional to the applied voltage.The mechanism of such current spectra limitation for finite-width devices can be understood referring to the derivation in section III G: for finite-width devices, the sum in Eq. 41 contains only a limited number of quantized values of k y , therefore the contribution of the corresponding propagation modes with low energies E m is dampened, while only modes with the high energy (near Fermi level) can affect the current.The difference in the current spectra between finitely and infinitely-wide devices is diminished for very large tunnel junctions, as Fig. 6 indicates for L gap > 9 nm.Indeed, for sufficiently large tunnel gaps, the tunneling current spectrum is exponentially suppressed for all low-energy modes, since the effective tunnel barrier height is larger for them than for higher-energy modes.Consequently, only the modes in the immediate vicinity of the Fermi level contribute to the total current for both finite and infinite-width tunnel junctions.
Despite the seemingly significant differences between finite-width systems and the corresponding infinitely wide counterparts, evidenced in the free electron energy distribution (see Fig. 4) and in their respective current spectra (see Fig. 6), the computed total currents, obtained by integrating the current spectra, only differ by about 20%.We thus conclude that while the density of states and the current spectra significantly differ, the sheet conductance (i.e.conductivity) calculations are accurate enough for finite-size systems at least down to 12 nm-widths.At the same time, for tunnel junctions with sufficiently large gaps, both the total current and the current spectra are very similar between infinite-width and finite-width systems, since in both cases the current is carried only by the high-energy electrons in the vicinity of system's Fermi level.

C. Quantization effects in δ -layer tunnel junctions
The tunneling current vs the tunnel gap length L gap is shown in Fig. 7 for three different voltages V = 1 mV, 10 mV and 100 mV.In the range of gap lengths L gap = 0, ..., 7 nm, the I vs L gap trend is practically exponential for all three voltages, i.e. ln I ∼ ln I L gap =0 − L gap /B voltage , where I L gap =0 is the current when L gap = 0, L gap is the tunnel gap length and B voltage is a proportional constant related with the barrier height.Indeed, in this range of gap values we get B 1−10 mV ≈ 1.9 nm and B 100 mV = 2.9 nm, which can be translated into the effective barrier heights of an equivalent rectangular barrier (L gap ×h voltage ), where h voltage = d val h2 /(8B 2 voltage m t ), h is the reduced Planck's constant, d val is the valley degeneracy, and m t is the transverse effective electron mass.Substituting the corresponding B voltage values, we obtain h 1−10mV ≈ 80 meV and h 100mV ≈ 30 meV, which agrees well with the effective barrier height shown in Fig. 8 for a tunnel gap of L gap = 7 nm, i.e. the barrier height from the Fermi level to the maximum value of the effective 1D electrostatic potential, computed using our charge self-consistent scheme.However, a deviation from the exponential trend can be seen for large tunnel gaps, L gap > 7 nm, and it is the consequence of the quatization of the low-energy conduction band [16] and the resulting mismatch between occupied and unoccupied quasi-discrete electron states in the left and right δ -layers, respectively, at low bias.Fig. 8 shows the LDOS along x-direction for a tunnel junction of L gap = 7 nm, i.e. the available states which can be occupied by the free electrons in space-energy dimension; in the low-temperature regime, only the states below the Fermi level are occupied.As one might observe from Fig. 8, the lowenergy LDOS are quasi-quantized, highlighted with dashed lines in the figure, however, for high energies, the states are not quantized anymore, the LDOS are practically continuous in space-energy, as can be seen above the Fermi level for high energies.When the occupied quasi-discrete states from the left side overlap most with the unoccupied quasidiscrete states (corresponding to the ones above the Fermi level) from the right side, it results in a considerably increase of the tunneling current as shown in Fig. 7 for L gap = 10 nm; on the contrary, if the overlap is reduced, as it happens for L gap = 11 nm, there is a reduction of the current as the result of the mismatch.For low biases ( 10 mV), the mismatch can only exist for sufficiently large tunnel gaps, L gap > 7 nm, because the existing coupling of the left and right δ -layer wavefunctions for narrow tunnel gaps (L gap < 7 nm) effectively equalizes the electron states on both sides, thus reducing this mismatch.Conversely, the quantization effect of the conduction band on the conductivity vanishes for high applied voltages ( 100 mV) as shown in Fig. 7 for an bias of 100 mV.When a high bias is applied, e.g. to the right side (or drain) of the tunnel junction, it effectively makes the unoccupied highenergy continuous states of the right side available for the tunneling from the left side, thus negating the effects of the quantization of the occupied left side states on the current, as can be seen in Fig. 8 b.

D. Conclusions
In this work, we have used an efficient open-system quantum-mechanical treatment to explore the conductive band structure and the conductive properties in Si:P δ -layers for device widths ranging from nano-scale to macro-scale dimensions, and analyze the influence of size quantization effects on the conductive properties for sub-12 nm device widths.For device widths W < 10 nm, quantization effects strongly affect not only the conductivity, but also the spatial distribution of the current-carrying electron states.Conversely, for W > 10 nm, the quantization effects practically vanish and the conductivity tends to the values of infinitely-wide devices.Additionally, we have revealed and discussed the mechanism of two conductivity regimes in δ -layer tunnel junctions with L gap > 7 nm: a low-voltage regime, where conduction band quantization effects play a very significant role, and a highvoltage regime, where the quantum effects on the current are negligible.
Finally, we point out that the strong spacial quantization of the current-carrying states could be utilized in novel electronic δ -layer switches, where the number of propagating modes and their match/mis-match could be controlled by external electric fields, thus strongly affecting the current.In regular δlayer conductors the particular distribution of current-carrying states directly affects their penetration depth into Si body and cap, which typically has a large concentration of impurities (see e.g.[3]).Thus, the control of the number of propagating modes may give an additional degree of control over the rate of impurity scattering.

III. METHOD A. Open-system treatment
The open-system framework used in this work is based on an application of Keldysh formalism [27], known as Non-Equilibrium Green Function (NEGF) [28] (Sect.III B), and the effective mass theory (Sect.III E).The NEGF formalism, described in Sect.III B, defines the Green's function matrix through the inverse G(E) = [1 − H 0 − Σ(E)] −1 , which allows to compute the electron density and current through a charge self-consistent scheme as described in Sect.III F. However, the computation of an inverse matrix generally also scales as O(N 3 ), with N being the size of the device Hamiltonian matrix H 0 (Sect.III E).For instance, for a device of dimensions 12nm × 15nm × 50nm with a uniform grid of 0.2 nm, N > 10 6 and N 3 > 10 18 , making a direct computation of the inverse very challenging.Here we have employed an efficient implementation of NEGF, refereed to as the Contact Block Reduction (CBR) method, which allows accurate computation of all quantum-mechanical quantities of interest (the local density of states, transmission probability, current) and scales linearly O(N) with the system size.The CBR method in 3D real-space is summarized in Sects.III C and III D.

B. Non-Equilibrium Green's Function Formalism
The current for two-contact device (D) from source (s) to drain (d) can be computed within the Landauer-Buttiker formalism through the transmission function where i e (E) is the current spectrum, e is the electron charge, h is the Planck's constant, E is the energy, is the Fermi-Dirac distribution function within source (drain), to which a voltage V s(d) is applied, and T sd (E) is the electronic transmission from source to drain.Within the Green's function formalism [28], the transmission function is given by with where H 0 D + Σ s + Σ d is the non-Hermitian effective Hamiltonian (N D × N D )-matrix of the open device.Eq. 4 can be also expressed in terms of the Green's function of the decoupled device (i.e. the closed system) , where E + = I(E + iε), with ε → 0+, via the Dyson equation: The decoupled device Green function G D can be computed using its spectral representation: and E α and |ψ α are the eigenvalues and eigenvectors of where LDOS(r i , E) is the local density of states.Note that the dimension of G D and Γ λ is (N D × N D ).

C. Contact Block Reduction method
We next review the Contact Block Reduction (CBR) method presented in Mamaluy et al. [34][35][36], which allows a very efficient calculation of the density matrix, transmission function, etc. of an arbitrarily shaped, multiterminal twoor three-dimensoinal open device within the NEGF formalism.In the following we apply the CBR method to the twocontact device shown in Fig. 9.The device consists of two semi-infinite contacts, source (s) and drain (d), which are in contact to the engineering channel, named as device (D).
We first discretize the domain of the device in N D gridpoints, and subdivide them into N C = N C s + N C d boundary grid-points and into N D i = N D − N C interior grid-points.N C s(d) corresponds to the boundary grid-points between the device and source (drain).Furthermore, we assume that the realspace Hamiltonian matrix that corresponds to this discretization only couples sites within some finite range with one another, typically first nearest-neighbors.The total grid-points in the device domain is described by the following set With this discretization of the device domain, the selfenergy matrices Σ s(d) , which represent the coupling of the source (drain) to the device, is given by where W Ds(d) are the Hamiltonian coupling matrices between the device D and contact s(d) and G 0 s(d) are the retarded Green's function matrix of decoupled source (drain).The total self-energy matrix (N D × N D )-Σ can be expressed by the following block-diagonal matrix of the form where is of dimension of (N C × N C ). Analogously, the retarded Green's function matrix for the device can be also rewritten by the following submatrices where Using the new representation for the above matrices, the Dyson equation can be then rewritten as where the matrices have the following forms: and The submatrices Thus, the retarded Green's function of the open system is given by and the Green's function within the contact region by The electron transmission from source s to drain d can be rewritten as where ).
We note that the dimension of all submatrices involved in Eqs.16 and 17 is of ( Similarly, the electron density can also be written as where the density matrix due to source(drain), Ξ with and where 1 C is the identity matrix of dimension of (N C × N C ).Note that the dimension of Σ C and G 0 C is (N C × N C ) as well.

D. Incomplete set of CBR eigenstates
The major feature of the CBR method is the ability to use a greatly reduced, incomplete set of specially defined eigenstates to represent the true open-system solution.As was shown in Mamaluy et al. [34], it can be accomplished by imposing a special kind boundary condition to the decouple device.The idea is to be able to rewrite the Green's function matrix of the open system G D (E), Equation 4, as where H N D is independent of the energy, and Σ N s (E) and Σ N d (E) tend toward zero for values of E that lies not far from the band edge.This enables us to solve the Dyson equation with an incomplete basis.We start assuming that all nonzero coupling matrix elements of W s(d)D is equal to a real constant value, W s(d) .This leads to the following expression for the self-energy matrix Σ s(d) within each source(drain) contact, s(d), see e.g.Datta [28] where a s(d) = a is the constant parameter of the lattice, M T is the total number of propagating modes, and χ m s(d) are the modes of the Schrödinger equation . The wave vectors k m s(d) are functions of energy E as The exponential term in Equation 22 can be approximated in power of k m s(d) as which leads to Thus, the self-energy matrix in each contact, s and d, consists of a first term, which is an independent-energy diagonal matrix, and a second term, that vanishes not too far from the low energy band edge.The total self-energy matrix can be then rewritten as where and the retarded Green's function matrix as The self-energy Σ N (E) is small for E close to the band edge, and the Hamiltonian matrix of the device is defined by It can be shown that this Hamiltonian matrix corresponds to generalized Neumann boundary conditions for its eigenfunctions [36].Thus, the eigenfunctions of H N D are approximate solutions of the open-boundary problem in the low-energy limit.As a consequence it suffices to include an incomplete set in the spectral representation of H N D in the calculation of the open-system quantities of interest within some limited energy range.

E. Device Effective-Mass Hamiltonian
The effective mass Schrödinger equation for Γ-valley electrons in Si is given by where H 0 D is the effective-mass Hamiltonian operator and φ H (r i ) is the Hartree potential and φ XC (r i ) is the exchange-correlation potential.
We discretize the 3D domain in N x equidistant grid-points along the x-axis, in N y equidistant grid-points along the yaxis and in N z equidistant grid-points along the z-axis.For all directions, the separation between grid-points is chosen to be ∆x = ∆y = ∆z = 0.2 nm.A grid-point is defined by the following triple indices (i, j, k), with i = 1, ..., N x , j = 1, ..., N y and k = 1, ..., N z , and the global index n = k + N z (N x ( j − 1) + i − 1).Therefore, the effective mass Hamiltonian matrix for the decoupled device can be expressed with seven non-zero diagonals as follows where and The number of zero diagonals between the diagonals l and b is (N z − 2) and between the diagonals b and g is (N x − 1)N z − 1.The exchange-correlation potential φ XC (r i ) that accounts for electron-electron interaction is computed using the terms given in Perdew and Zunger [40].

F. Charge Self-Consistency
We solve charge self-consistently the open-system effective mass Schrödinger equation and the non-linear Poisson equation [37,38,41].We employ a single-band effective mass approximation with a valley degeneracy of d val = 6.For the charge self-consistent solution of the non-linear Poisson equation we use a combination of the predictor-corrector approach and Anderson mixing scheme [37,38].
Firstly, the Schrödinger equation is solved in a specially defined (see generalized Neumann BC in section III D) closedsystem basis taking into account the Hartree potential φ H (r i ) and the exchange and correlation potential φ XC (r i ).It is found that out of more than 1,000,000 eigenstates only about 700 (< 0.1%) of lowest-energy states is needed (see section III D), which is generally determined by the material properties (e.g.doping level), but not the device size.Then, the LDOS of the open system, ρ(r i , E), and the electron density, n(r i ), are computed using the CBR method.The potential and the carrier density are then used to calculate the residuum F of the Posisson equation where A is the matrix derived from the discretization of the Poisson equation and N D and N A are the total donor and acceptor doping densities arrays, respectively.If the residuum is larger than a predetermined threshold ε, the Hartree potential is updated using the predictor-corrector method, together with the Anderson mixing scheme.Using the updated Hartree potential and the corresponding carrier density, the exchangecorrelation is computed again for the next step, and an iteration of Schrodinger-Poisson is repeated until the convergence is reached with F [φ H (r i )] < ε = 10 −6 eV.A typical convergence rate for Si:P δ -layer structures of different thicknesses is shown in Fig. 10.It illustrates that the proposed charge self-consistent scheme has a robust convergence for different structures, providing the residuum error reduction of 5 orders of magnitude within 40 iterations.In all simulation in this work the standard values of electron effective masses m l = 0.98 × m e , m t = 0.19 × m e and the dielectric constant ε Si = 11.7 of Silicon were employed.
G. From finite-width to infinite-width systems In the scenario that the width of a conducting structure along the y-direction is assumed to be infinite (see Fig. 9) with a flat electrostatic potential, we can write the solutions of the Schrödinger equation as the product of plane-waves along yaxis and the solutions of the 2D Schrödinger equation in the X-Z plane of the device.
which is normalized to a length of L y , and The electron density can be expressed by summing over all (generally -open-system) states α as where k B = 8.617 × 10 −5 eV K −1 is the Boltzmann constant, T is the temperature of the system and f (E) is the Fermi-Dirac distribution function, which provides the probability that a state is occupied or unoccupied.The Fermi-Dirac distribution function is given by Using equations 37, 38 and 39 we get where we denote the term in parenthesis as the effective 2D Fermi-Dirac distribution function f 2D (E).Replacing the summation over k y with an integral as 1 2π , we can get the distribution function in an integral form as follows Finally, performing a variable change l = h2 k 2 y 2m t k B T we obtain Fig. 11a shows the Fermi-Dirac distribution function, Eq. 40, and the effective 2D Fermi-Dirac distribution function, Eq. 43, in equilibrium conditions and a temperature of 4.0 K.In both cases the occupied states exists only below the Fermi level, but the occupation rate for the effective 2D distributions increases as ∼ √ −E for E → −∞.We next examine the behaviour of the current integral given in the Eq. 1 at low temperatures.In general, the difference in the distribution functions f s (E) − f d (E) = f (E) − f (E + qV bias ) determines the energy range of non-zero contributions to the total current.The difference between Fermi-Dirac distributions is non-zero only in a small range of energies in the vicinity of the Fermi level as the result in red line in Fig. 11b indicates.Furthermore, the energy range for non-zero contribution is proportional to the applied voltage V bias = 0.1 V. On the contrary, the difference in the effective 2D distribution functions is non-zero for all energies below the Fermi level as the plot in blue color in Fig. 11b illustrates.This is a reflection of the asymptotic behaviour of f 2D (E) at E → −∞; for the difference in the distributions one gets f 2D s (E) − f 2D d (E) = f 2D (E) − f 2D (E + qV bias ) ∼ qV bias / √ −E.
We thus conclude that, within the Landauer-Buttiker/NEGF formalism, finite size structures in the transverse directions have a principally different current spectrum from the structures that are infinite along a transverse direction (e.g. when W → ∞ in δ -layer structures shown in Fig. 1).The employee owns all right, title and interest in and to the article and is solely responsible for its contents.The United States Government retains and the publisher, by accepting the article for publication, acknowledges that the United States Government retains a non-exclusive, paid-up, irrevocable, world-wide license to publish or reproduce the published form of this article or allow others to do so, for United States Government purposes.The DOE will provide public access to these results of federally sponsored research in accordance with the DOE Public Access Plan https://www.energy.gov/downloads/doe-public-access-plan.This version of the article has been accepted for publication, after peer review (when applicable) but is not the Version of Record and does not reflect post-acceptance improvements, or any corrections.The Version of Record is available online at: http://dx.doi.org/10.1038/s41598-022-20105-x

FIG. 2 .
FIG. 2. Conductivity for Si: P δ -layer devices.The conductivity is shown in function of the device width W for different tunnel gap lengths L gap .The corresponding conductivity values for infinitelywide devices are indicated in dashed lines.N D = 1.0 × 10 14 cm −2 , N A = 5.0 × 10 17 cm −3 , t = 1.0 nm and for an applied voltage of 1 mV.

FIG. 4 .
FIG. 4. Free electron energy distribution for Si:P δ -layer systems.It shows a comparison of the free electron energy distribution for infinite and 12 nm-width δ -layer devices.The free electron energy distribution for the infinite-width δ -layer is normalized to a width of the device of W = 12 nm for comparison purpose.N D = 1.0 × 10 14 cm −2 , N A = 5.0 × 10 17 cm −3 and t = 1.0 nm.

FIG. 5 .
FIG. 5. Electronic transmission for Si:P δ -layer systems.It shows a comparison of the transmission function (logarithmic scale) for an applied voltage of 1.0 mV for infinite and 12 nm-width δ -layer devices.The transmission function for the infinite-width δ -layer is normalized to a width of the device of W = 12 nm for comparison purpose.N D = 1.0 × 10 14 cm −2 , N A = 5.0 × 10 17 cm −3 and t = 1.0 nm.

FIG. 6 .
FIG. 6.Current spectrum for Si:P δ -layer systems.a Comparison of the current spectrum i e (logarithmic scale) for an applied voltage of 1.0 mV between infinite and 12 nm-width δ -layer devices.b Close-up of the current spectrum within the Fermi level.The currents for the infinite-width δ -layer systems are normalized to a width of the device of W = 12 nm for comparison purpose.N D = 1.0 × 10 14 cm −2 , N A = 5.0 × 10 17 cm −3 and t = 1.0 nm.

FIG. 7 .
FIG. 7. Characteristic tunneling current curves.The total current (logarithmic scale) I vs. tunnel gap length L gap for three different applied voltages of V = 1 mV, 10 mV and 100 mV is shown.N D = 1.0 × 10 14 cm −2 and N A = 5.0 × 10 17 cm −3 , W = 12 nm and t = 1 nm.Dashed lines represent least-square fits to the exponential trend.

FIG. 8 .
FIG. 8. Local Density of States for δ -layer tunnel junctions.The LDOS(E, x) for a tunnel junction of L gap = 7 nm is shown in a and b when a voltage of 1 mV and 100 mV is applied to the right side of the tunnel junction, respectively.Subfigure a indicates that the quantized states around the Fermi level might affect the conductivity for low applied biases due to a possible mismatch between the left and right quasi-discrete states, i.e., when the corresponding quasi-discrete peaks from the left and right sides are not aligned.Subfigure b shows that for high applied voltages the mismatch becomes impossible due to the availability of unoccupied continuum states on the right side for tunneling from the left side.In a and b, the corresponding effective 1D potentials are also shown, calculated by integrating over the (y,z)-plane the actual charge self-consistent 3D potentials weighted with the electron density.N D = 1.0 × 10 14 cm −2 , N A = 5.0 × 10 17 cm −3 , W = 12 nm and t = 1 nm.
) where Γ s(d) are the coupling (N D × N D )-matrices between the device and the source (drain), Σ s(d) are the self-energy (N D × N D )-matrices, which describe the effects of the source (drain) on the electronic structure of the device by providing the appropriate open-system boundary conditions, G D and G † D are the retarded and advanced Green's function (N D × N D )-matrices of the coupled device with the source and drain (open-system device), respectively.N D is the total number of grid points of the discretized device domain.The retarded Green's function matrix is defined by

FIG. 9 .
FIG. 9. Schematic computational model for the CBR method.The two-contacts model is composed of two semi-infinite contacts, source (s) and drain (d), and a channel or device (D).

FIG. 10 .
FIG. 10.Convergence evaluation of the method.The convergence rate of the charge self-consistent open-system Schrodinger-Poisson equations for different δ -layer thicknesses (t) is shown.

FIG. 11 .
FIG. 11.Fermi-Dirac distribution functions.a Comparison between the Fermi-Dirac distribution function, Eq. 40, and the effective 2D Fermi-Dirac distribution function, Eq. 43, for equilibrium conditions and T = 4.0 K. b Comparison between the Fermi-Dirac function difference between source and drain, f s (E) − f d (E) = f (E) − f (E + qV bias ), and the corresponding effective 2D Fermi-Dirac distribution function, f 2D s (E) − f 2D d (E) = f 2D (E) − f 2D (E + qV bias ), for an applied voltage of V bias = 0.1 V and T = 4.0 K.