Uncovering anisotropic effects of electric high-moment dipoles on the tunneling current in \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

The precise positioning of dopants in semiconductors using scanning tunneling microscopes has led to the development of planar dopant-based devices, also known as \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-based devices, facilitating the exploration of new concepts in classical and quantum computing. Recently, it has been shown that two distinct conductivity regimes (low- and high-bias regimes) exist in \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 due to the presence of quasi-discrete and continuous states in the conduction band of \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. Furthermore, discrete charged impurities in the tunnel junction region significantly influence the tunneling rates in \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. Here we demonstrate that electrical dipoles, i.e. zero-charge defects, present in the tunnel junction region can also significantly alter the tunneling rate, depending, however, on the specific conductivity regime, and orientation and moment of the dipole. In the low-bias regime, with high-resistance tunneling mode, dipoles of nearly all orientations and moments can alter the current, indicating the extreme sensitivity of the tunneling current to the slightest imperfection in the tunnel gap. In the high-bias regime, with low-resistivity, only dipoles with high moments and oriented in the directions perpendicular to the electron tunneling direction can significantly affect the current, thus making this conductivity regime significantly less prone to the influence of dipole defects with low-moments or oriented in the direction parallel to the tunneling.


Introduction
Atomic precision advanced manufacturing (APAM) has enabled the creation of 2D doped regions, also known as δ -layers, in semiconductors with single-atom precision [1][2][3][4][5] and high conductivity [6][7][8][9][10][11] .The APAM is a process to incorporate dopants, such as P or B, at the atomic scale onto Si surface using chemistry surface 10,12 .In a simplified way, this process consists in several steps as follows for phosphorus-doped planar structures embedded in silicon (Si: P δ -layer systems) 13 : We start with a Si surface, normally (100), fully passivated with H; with the tip of a Scanning Tunneling Microscope (STM), it gives the capability to remove H atom by H atom in the exact locations where we want to incorporate the dopants; then, the surface is exposed with a precursor gas, such as phosphene (PH 3 ) 12 , followed by a annealing process to incorporate the dopants into the surface; finally, an epitaxial Si is overgrown through a series of annealing processes to protect the planar structure.
APAM has various applications, including the exploration of novel electronic devices such as nanoscale diodes or transistors for classical computing and sensing systems 10,[14][15][16] .But, most importantly, this technology has been used to explore dopantbased qubits in silicon, with recent advancements in understanding the exchange-based 2-qubit operations 17 , the limits to qubit fidelity from environmental noise 18 , the advantages of leveraging the number of dopants as a degree of freedom 19,20 , and the exploration of many body 21 and topological 22 effects in dopant chains.One of the building block of these devices is the δ -layer tunnel junction (see Fig. 1 a), which consists of two highly doped thin layers separated by an intrinsic gap and embedded in a semiconductor material.These devices require precise control of the tunneling rate for their functioning, as they are very sensitive to tunneling rates.Additionally, it is known that imperfections near the tunnel junction strongly alter the tunneling rate.
In our previous work 23,24 , we reported the existence of two different Ohmic regimes in the conductivity of δ -layer tunnel junctions, namely low-and high-bias regimes, connected by a transition regime.For convenience, we include the characteristic IV curve for δ -layer tunnel junctions in Fig. 2 reported in Ref. 24.One can clearly observe these regimes from the result: the first conductivity regime, characterized by a resistance of approximately 5 − 6 MΩ, is observed for an applied bias in the range of 0 − 50 mV between the source and drain; the second one, characterized by a resistance of 0.2 − 0.3 MΩ, occurs at voltages above 80 mV; and, these two Ohmic regimes are separated by a transition regime, approximately between 50 − 80 mV, where the resistance decreases as the bias is increased.We can understand this phenomenon because of the existence of quasi-discrete states and continuous states in δ -layer systems [23][24][25] .The threshold voltages of these regimes are determined by the doping density N D and the δ -layer thickness t.Further in Ref. 24, we investigated how the tunneling rate in δ -layer tunnel junctions can be affected by the presence of charged impurities, specifically focusing on single n-type and p-type impurities in the middle of the tunnel junction.The results, for convenience, are shown in Figs. 3 and 4 (see continuous and dashed black lines) for distinct tunnel junction lengths.It was found that the tunneling rate is strongly affected by the presence of these charges: the tunneling rate significantly increases by the presence of n-type impurities, while it decreases for p-type impurities.Additionally, the effect of n-type impurities was found to be stronger than for p-type impurities, specially for large tunnel gaps, in which the effect can be up to one order of magnitude higher.Therefore, these previous works also motivates the present analysis, in which we investigate the influence of electrical dipoles, i.e. zero-charge defects, on the tunneling rate in δ -layer tunnel junction systems.The formation of a dipole occurs when two impurities, one n-type and the other p-type, come into proximity either during the device fabrication process or throughout the device's operational lifespan, as a result of dopant diffusion; Additionally, another possibility is when the formation of a charge-neutral point defect with a non-zero dipole moment.In this work, we focus on the former scenario, and similar conclusions can be drawn for the latter case as well.
In this work, we have employed an efficient implementation of the Non-Equilibrium Green's Function (NEGF) method, referred to as the Contact Block Reduction (CBR) method [26][27][28][29][30] , to investigate how electric dipoles might alter the tunneling rate in phosphorus δ -layer tunnel junctions in silicon.We have found quantum-mechanical effects, which can not be described by classical methods: electrical dipoles located near the gap in δ -layer tunnel junctions can significantly alter the current, and its orientation indeed matters.The effect of dipoles oriented in the directions perpendicular to the electron tunneling direction is significantly stronger than for dipoles oriented in the direction parallel to the tunneling direction, revealing an anisotropic effect.This intriguing effect is only present for sufficiently 'large' (a few nanometers or more) dipole moments.For smaller dipole moments, the anisotropic effect vanishes.Similarly, at low bias levels, the anisotropic effect is less pronounced.

Simulation Setup
To explore the impact of electric dipoles in the intrinsic gap of a P δ -layer tunnel junction embedded in silicon, we adopt the structure shown in Fig. 1 a.This type of structure is normally referred to as Si: P δ -layer tunnel junction.In the open-system NEGF framework, the computational device consists of a semi-infinite source and drain, in contact with the channel of length L, which is composed of a lightly doped Si body and Si cap and two very thin, highly P doped layers (referred to as left and right δ -layers) separated by an intrinsic gap of length L gap .The channel length L is chosen to be 30 nm + L gap to avoid boundary effects between the source and drain contacts, the device height H is 8 nm, and the device width W is chosen to be 15 nm, with an effective width of 13 nm for the δ -layers, to minimize size quantization effects on the conductive properties 23 .In our analyses, we have considered a thickness of δ -layer of t = 1 nm, a sheet doping density of N D = 10 14 cm −2 (N for the δ -layers, and a doping density of N A = 5 × 10 17 cm −3 in the Si cap and Si body.These doping densities and dimensions are of the order of published experimental work 10,12,[31][32][33] .Furthermore, all simulations are carried out at the cryogenic temperature of 4.2 K, for which we can neglect inelastic scatterings 6,34 . In this work, we only focus on the effect of dipoles located in the middle of the tunnel gap, with three different orientations, along x-, y-and z-directions, as shown in Fig. 1 b, c, and d, respectively, and two distinct electric dipole moments, corresponding to two different scenarios, i.e. with low moment and high moment.The moment of a dipole p is defined as p = q × l dipole , where q is the electric charge and l dipole is the distance between both charges.For this work, we use a dipole length l dipole of 0.8 and 4 nm, corresponding to the low moment and the high moment, respectively.When dipoles are oriented along x-direction, we refer to them as dipoles oriented in the propagation direction (i.e.electron tunneling direction), whereas along y-and z-directions we refer to them as dipoles oriented in the transverse directions (i.e.perpendicular to the electron tunneling direction).In these simulations, dipoles are modeled as two point charges of opposing electrical sign, in which each charge is modeled by approximating a point charge with a density of (positive or negative) 4.6 × 10 21 cm −3 homogeneously distributed in a total volume of (0.6 nm) 3 .The final spatial distribution of the total charge will be dictated by the self-consistent solution of the open-system Schrödinger and Poisson equations.Additionally, while our analysis in this work is restricted to the center-gap location, exploring the influence of other locations could be intriguing.This is because free electrons in δ -layer systems form distinct conducting layers perpendicular to the confinement direction, thus signaling a highly non-homogeneous electron density distribution 25,35 .

Results and Discussion
We begin by examining the current ratio, defined as the ratio of the current of the system with the dipole to the current of the system without the dipole, and the current change ratio, defined as the ratio of the current change magnitude between the system with and without the dipole to the current of the system without the dipole.Figs. 3 and 4 show these analyses for an applied bias of 1 and 100 mV, respectively, representing the low-and high-bias regimes.Both figures include the analyses for both dipole moments, l dipole = 0.8 nm for low-moment and l dipole = 4 nm for high-moment, and all considered orientations, along the propagation direction and along the transverse propagation directions.From these results, we can notably discern very interesting results: i) The trend of the current ratio with the tunnel gap length in the low-bias regime shows a considerable oscillation, and this oscillation vanishes in the high-bias regime; ii) High-moment dipoles exhibit anisotropic behaviour since dipoles oriented in the transverse propagation directions behave differently than dipoles along the propagation direction; in contrast, low-moment dipoles does not show that strong anisotropic behaviour since the impact of the dipole orientation is minimum; additionally, high-moment dipoles oriented along the propagation direction behave very similar to low-moment

4/12
dipoles; iv) In the high-bias regime, high-moment dipoles oriented in the transverse propagation directions significantly affect the current, whereas the contribution is minimum for all low-moment dipoles and high-moment dipoles oriented along the propagation direction; v) In the low-bias regime, specially for large tunnel gap lengths, dipoles of all orientations and moments noticeably affect the tunneling rate in an indistinguishable manner; for narrow tunnel gap length, the effect of all dipoles on tunneling rate is of the same order of magnitude and less pronounced.In the following, we discuss these intriguing results in more details.
To have a deeper understanding of our analyses, we need to examine the local density of states (LDOS) of the free electrons, which represents the conduction band in real space for the free electrons.Figs. 5 and 6 show the LDOS(x, E) along the x-direction for a δ -layer tunnel junction of L gap = 10 nm with the presence of a dipole of length l dipole = 0.8, 4 nm, respectively, in the middle gap: a and b for an applied voltage of 1 and 100 mV to the drain contact while the source is grounded, respectively.In the figures, the states between x = 0 − 15 nm and between x = 25 − 40 nm correspond to the left and right δ -layers, respectively; the states between x = 15 − 25 nm correspond to the intrinsic gap region.At 0 K, the states below Fermi level are occupied, and above the Fermi level are unoccupied.As the temperature increases, some states above the Fermi level start to be occupied in detriment of states below Fermi level.As one can observe from the LDOS(x, E), within the δ -layer regions, the low-energy LDOS(x, E) are strongly quantized, i.e. for energies below the Fermi level and approximately up to 50 meV above the Fermi level, highlighted with dashed lines in Fig. 5; for higher energies, approximately above 50 meV, the LDOS(x, E) are practically continuous in the energy-space dimension, thus signaling that these states are not quantized.Both regions are accordingly marked in the figures.The quantization of the conduction band arises from the confinement of dopants in one direction, resulting in very sharp doping profile and leading to the size quantization of the δ -layer 23,24 .The presence of discrete states have been observed experimentally in several high resolution Angle-resolved Photoemission Spectroscopy (ARPES) measurements for δ -layers in silicon 36,37 .Our previous simulations revealed that the number of quantized conduction bands and their corresponding energy splittings are strongly influenced by both the δ -layer thickness t and the doping density N D 25 .When a positive voltage is applied to the drain contact while the source is grounded, the Fermi level corresponding to the drain contact is shifted down, resulting in lowering the energies of all states in the right side as well.Consequently, new unoccupied states, either in the right δ -layer or those created by the dipole, become available to be occupied by tunneling electrons coming from the left δ -layer.Thus, we can distinguish two tunneling processes: the first one corresponds to the direct tunneling, i.e. electrons tunneling from occupied states in the left δ -layer to available states in the right δ -layer; and, the second one is the defect-mediated tunneling, i.e. electrons first tunnel from the left δ -layer to the available state in the defect (dipole), and from there to the right δ -layer.Since we are only considering elastic scattering, it is required for both tunneling processes that the energy and momentum are conserved.
The most noticeable observation from Figs. 3 and 4 is that the current ratio considerably oscillates with the tunnel gap length in the low-bias regime (1 mV), while it becomes very smooth in the high-bias regime (100 mV).This is as a consequence of the strong quantization of the low-energy conduction band in δ -layer systems 24 .When a low voltage is applied, e.g. 1 mV, see panel a in Figs. 5 and 6, only the unoccupied quantized states in the right δ -layer (those just above the Fermi level), and all quasi-bounded states created by the dipole will play an important role in the tunneling process.If the occupied quasi-discrete states in the left δ -layer or the quasi-bounded states in the dipole align with the unoccupied quasi-discrete states in the right δ -layer, it will result in a considerable increase of the tunneling current.Conversely, if the overlap is minimum, the tunneling current will be reduced.For low biases, this alternating mismatch can only exist for sufficiently large tunnel gaps, because the coupling of the left and right δ -layer wave-functions for narrow tunnel gaps equalizes the electron states on both sides, increasing the overlap and thus eliminating the mismatch.On the other hand, when an high bias is applied, e.g. for 100 meV, see panel b in Figs. 5 and 6, it makes the continuous unoccupied high-energy states in the right side available for tunneling from the left side, thus diminishing the influence of the conduction band quantization on the current and easing the tunneling process.Thus, this difference in the tunneling mechanism explains the oscillating behaviour of the tunneling ratio in the low-bias regime, specially for large tunnel gaps, and the disappearance of these oscillations in the high-bias regime.Likewise, it also explain the existence of the two conductivity regimes: the first one, between 0 and 50 mV, with high tunneling resistance, where only the quasi-discrete states play a role in the tunneling process; and, the second one, with low tunneling resistance, for voltages above 80 mV, where quasi-discrete and continuous states contribute to the tunneling process.
In the high-bias regime, the effect of low-moment dipoles on the tunneling current is minimum.This finding can be deduced from Fig. 4, in which the tunneling ratio for all orientations of low-moment dipole (l dipole =0.8 nm) is approximately unity (cf. the results in color dashed lines in a) and the change of the tunneling current is between 1-30% (cf. the results in color dashed lines in b) for all studied tunnel junction lengths.Specifically, the tunneling change increases almost exponentially with the tunnel gap length, from a change of 1-3% for a tunnel gap length of L gap = 3 nm up to around 20-30% for L gap = 12 nm.In addition, our results also indicate negligible impact of the dipole orientation on the tunneling since the three studied orientations show similar change in the tunneling rate.These results also indicate that when both impurities are very close to each other, their effects on the tunneling balance out, resulting in minimal changes to the tunneling current, as well as explaining the 5/12   diminished effect on the tunneling for low-moment dipoles.Furthermore, it is evident from the results that the change in the tunneling current exhibits exponential growth as the tunnel junction length increases.This is primarily due to the dominant effect of the n-type impurity, which exhibits exponential growth (cf.continuous line in b), while p-type impurity exhibit lower order of growth with the tunnel junction length (cf.dashed line in b).
A different story occurs for high-moment dipoles in the high-bias regime.From Fig. 4, we can observe that the orientation of high-moment dipoles in the high-bias regime plays an important role: our results reveals an intriguing anisotropic effect on the tunneling current due to their orientation.When the dipole is oriented along the transverse propagation directions (i.e.y-and z-directions), the tunneling current is strongly affected by the dipole; on contrary, when the dipole is oriented along the propagation direction (x-direction), the tunneling current is only weakly affected.Additionally, it behaves very similar to low-moment dipoles (see red continuous line).In the same manner as low-moment dipoles and for the same reason, the change of the tunneling current increases exponentially with the tunnel gap length.However, for dipoles oriented in the transverse propagation directions, the tunneling change goes from 5-6% for L gap = 3 nm up to 2000% for L gap = 12 nm, while for dipoles oriented in the propagation direction, the change goes from 3% for L gap = 3 nm up to 100% for L gap = 12 nm.Furthermore, we can also observe that high-moment dipoles oriented in the transverse propagation directions almost reproduces the dotted line in Fig. 4 a, which means that the total effect for these dipoles on the current can be approximated as the arithmetic average of independent single n-type and p-type impurities.To understand why the effect of dipoles oriented in the transverse propagation directions is overall stronger than those oriented along the propagation direction, we might need to examine the effective 1D electrostatic potentials in equilibrium, shown in Fig. 7.We can see that the barrier height in the electrostatic potential for dipoles in the transverse propagation directions (cf.red lines in panels b and c) is reduced overall, leading to an increase in the tunneling current.Conversely, in the case of dipoles oriented in the propagation direction (cf.red line in panel a), we can see that the barrier height is increased near the p-type impurity and decreased near the n-type impurity.In other words, an increase in the barrier height results in a reduction of the current, while a decrease in the barrier height leads to an increase in the current.Thus, the total tunneling current is a combination of these two opposing effects, resulting in only a minimal change in the tunneling current in this particular case.Comparing the effective 1D electrostatic potentials for low-and high-moment dipoles, we can also note that the peak and dip in the electrostatic potential are less pronounced for low-moment dipoles, revealing again the cancellation of the effects when both impurities are in close proximity to each other.
Interestingly, in the low-bias regime, specially for larger tunnel junction lengths, dipoles of all orientations and moments affect noticeably the tunneling current in an almost indistinguishable manner.This is because only the quantized conduction band and the quasi-bounded states of the electric dipole are involved on the tunneling process in the low-bias regime (see panel a in Figs. 5 and 6).Unlike in high-bias regime, the tunneling change is exponential with the tunnel junction length for both dipole moments in a very similar manner.The tunneling change goes from 2-4% for L gap = 3 nm up to 620-1800% for L gap = 11 nm in low-moment dipoles, and from 2-5% for L gap = 3 nm up to 620-2200% for L gap = 11 nm in high-moment dipoles.Additionally, it is evident from the results that the anisotropic effect of the orientation of high-moment dipoles is much less significant compared to the high-bias regime.
Finally, we investigate the dwell time of the electrons in the tunnel junction, which is the average time that the electrons spend in the barrier, regardless of whether the electrons are reflected or transmitted [38][39][40] .It is defined mathematically as t dwell time = Ω n(r)∂ Ω I tun.curr./q , where I tun.curr. is the tunneling current, q is the elementary charge, n(r) is the electron density and Ω refers to the tunnel junction domain.From this definition, an alternative physical meaning of the dwell time can be readily deduced, which corresponds to the ratio of the total (average) number of electrons in the junction to the flux (number of electrons per second) going through the system in the steady-state.This ratio is the time necessary to fill the tunnel junction domain Ω with free electrons to achieve the charge neutrality from the fully depleted state.In other words, it corresponds to the time necessary to switch a device (e.g.FET-type device) from fully "off-" to "on-" state and/or vice versa.Thus, the dwell time can provide insights into the maximum operating frequency of δ -layer tunnel junction devices, which is important for applications.Fig. 8 shows the dwell time for the free electrons in the tunnel junction in the presence of single charged impurities and single dipole impurities, as well as without defects, in a and b for an applied voltage of 1 and 100 mV, respectively.We can observe several interesting results.Firstly, the dwell time of electrons in the barrier almost decreases inversely proportionally with the applied voltage.This can be explained by the weak increase in the number of electrons within the tunnel junction region with voltage, alongside a corresponding strong (proportional) increase in the current.Secondly, for the low-bias regime, the dwell time grows exponentially with the tunnel junction length, and it is almost identical for all cases.However, for the high-bias regime, we observe a strong dependence of the dwell time on the type of the impurity in the tunnel junction.Indeed, Fig. 8 b shows that a single p-type impurity significantly increases the dwell time, while a n-type impurity significantly decreases the dwell time.The corresponding microscopic interpretation is that a single additional acceptor atom in the gap binds a free electron (with a certain life-time), and thus effectively increases the average time that free electrons spend inside the barrier, whereas a donor has the opposing effect.Similarly, dipole impurities can also significantly affect dwell time of electrons in the barrier.From Fig. 8 b, it follows that, in the presence of low-moment dipoles and high-moment dipoles oriented along the propagation direction, the dwell times of the electrons are very similar to the dwell time for the ideal device.However, for high-moment dipoles oriented along the transverse propagation directions, the dwell times are noticeably reduced, indicating that electrons spend less time in the barrier.

Conclusions
In this work we analyze the influence of electric dipole impurities on the tunneling current in δ -layer systems.We employ an efficient implementation of the Non-Equilibrium Green's Function (NEGF) framework, referred to as the Contact Block Reduction (CBR) method [26][27][28][29][30] to carry out the simulations.Our analysis reveal several interesting results: i) The trend of the current ratio with the tunnel gap length in the low-bias regime shows a considerable oscillation; this oscillation vanishes in the high-bias regime.
ii) High-moment dipoles exhibit anisotropic behaviour since dipoles oriented in the transverse propagation directions behave differently than dipoles in the propagation direction; in contrast, low-moment dipoles does not show that strong anisotropic behaviour since the impact of the dipole orientation is minimum; Additionally, high-moment dipoles oriented in the propagation direction behave very similar to low-moment dipoles.iv) In the high-bias regime, high-moment dipoles oriented in the transverse propagation directions significantly affect current, whereas the contribution is minimum for all low-moment dipoles and high-moment dipoles oriented in the propagation direction.
v) In the low-bias regime, specially for large tunnel gap lengths, dipoles of all orientations and moments noticeably affect the tunneling rate in an indistinguishable manner; for narrow tunnel gap lengths, the effect of all dipoles on tunneling rate is of the same order of magnitude and less pronounced.
vi) Finally, we also perform the analysis of the dwell times for δ -layer tunnel junctions in the presence of single charged impurities and electrical dipole impurities in the tunnel gap, as well as without defects.Our results indicates the general suitability of δ -layer tunnel junction devices for TeraHertz applications in the "high-conductivity" regime (V ≥ 80 mV ).

Method
The simulations in this work are conducted using the open-system charge self-consistent Non-Equilibrium Green's Function (NEGF) Keldysh formalism 41 , together with the Contact Block Reduction (CBR) method [26][27][28][29][30]42 and the effective mass theory. The BR method allows a very efficient calculation of the density matrix, transmission function, etc. of an arbitrarily shaped, multiterminal two-or three-dimensional open device within the NEGF formalism and scales linearly O(N) with the system size N.The numerical details of the efficient open-system charge self-consistent treatment in 3D real-space are given in Refs.23-25.Within this framework, we solve self-consistently the open-system effective mass Schrödinger equation and the non-linear Poisson equation 26,28,30 .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 29,30 .First, the Schrödinger equation is solved in a specially defined closed-system basis taking into account the Hartree potential φ H (r r r i ) and the exchange and correlation potential φ XC (r r r i ) 43 .Second, the LDOS of the open system, ρ(r r r i , E), and the electron density, n(r r r i ), are computed using the CBR method for each iteration.Then the electrostatic potential and the carrier density are used to calculate the residuum F of the Poisson equation where A A A is the matrix derived from the discretization of the Poisson equation, and N N N D and N N 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 exchange-correlation is computed again for the next step, and an iteration of Schrödinger-Poisson is repeated until the convergence is reached with F F F[φ φ φ H (r r r i )] < ε = 10 −6 eV.In our simulations we have utilized a 3D real-space model, with a discretization size of 0.2 nm along all directions, thus with about 10 6 real-space grid points, and around 3,000 energy points were used.The CBR algorithm automatically ascertains that out of more than 1,000,000 eigenstates only about 700 (< 0.1%) of lowest-energy states is needed, which is generally determined by the material properties (e.g.doping level), but not the device size.We have also employed the standard values of electron effective masses, m l = 0.98 × m e , m t = 0.19 × m e , the dielectric constant of Silicon, and ε Si = 11.7.Further details of the methodology can be found in our previous publications [23][24][25]42 .

Validation
To validate our computational framework for δ -layer tunnel junctions, we computed the tunneling resistance for different tunnel gaps, L gap , and compared the calculations against recently measured data from M. Donnelly et al. 31 .Fig. 9 shows our computed tunneling resistance from our previous work 23 for an effective width of the δ -layer of 7 nm.The figure also includes the resistance measurements and tight-binding calculations for tunnel junctions of δ -layer width of 7 nm from M. Donnelly et al. 31 .One can observe that our predicted tunneling resistances (full black circles) are very close to both the experimental data (empty red circles) and the parameter-fitting tight-binding simulations (blue crosses).Fig. 9 demonstrates the true predictive power of our open-system effective mass framework for highly-conductive highly-confined systems.We also emphasize that in our calculations there is no fitting parameters, we use only the standard values of electron effective masses and dielectric constant for silicon.Finally, the slight differences between our computed tunneling resistances and experimental measurements can also be accounted for by the following reasons: i) Certain variations in the width, thickness, and doping density of the δ -layer (note that a doping density of N D = 10 14 was assumed in Ref. 23, while a higher doping density, N D = 2 × 10 14 , was employed in Ref. 31); ii) The possible presence of impurities and/or defects near the tunnel gap.

Figure 1 .
Figure 1.Computational model used in our simulations.A Si: P δ -layer tunnel junction device is shown in a, which consists of a semi-infinite source and drain, in contact with the channel.The channel is composed of a lightly doped Si body and Si cap and two very thin, highly P doped layers with an intrinsic gap of length L gap .The representation of an electric dipole oriented along x-direction (propagation direction) is shown in b, y-direction (transverse direction) in c, and z-direction (transverse direction) in d.The negative and positive charged impurities are represented as a blue and red spheres, respectively.

Figure 3 .
Figure 3. Tunneling ratio and tunneling rate change for an applied bias of 1 mV.Current ratio, I dipole /I ideal , between the δ -layer tunnel junction with and without the electric dipole in the intrinsic gap, oriented along the x-, y-and z-directions, for an applied bias of 1 meV in a, and the corresponding tunneling current change in b. t = 1 nm, N D = 10 14 cm −2 and N A = 5 × cm −3 .

Figure 4 .
Figure 4. Tunneling ratio and tunneling rate change for an applied bias of 100 mV.Current ratio, I dipole /I ideal , between the δ -layer tunnel junction with and without the electric dipole in the intrinsic gap, oriented along the x-, y-and z-directions, for an applied bias of 100 meV in a, and the corresponding tunneling current change in b. t = 1 nm, N D = 10 14 cm −2 and N A = 5 × cm −3 .

Figure 5 .
Figure 5. Local Density of States for a δ -layer tunnel junction with a low-moment dipole.It shows the LDOS(E, x) with a dipole of length 0.8 nm located in the middle of the tunnel junction and oriented along x-direction: a and b for an applied voltage of 1 and 100 mV to the drain contact while the source contact is grounded, respectively.The Fermi levels indicated in the figures correspond to the Fermi levels of the source and drain contacts.The corresponding effective 1D potentials for the ideal device and the device with the dipole are shown as blue dashed and red lines, respectively.L gap = 10 nm, N D = 10 14 cm −2 , N A = 5 × 10 17 cm −3 , and t = 1 nm.

Figure 6 .
Figure 6.Local Density of States for a δ -layer tunnel junction with a high-moment dipole.It shows the LDOS(E, x) with a dipole of length 4 nm located in the middle of the tunnel junction and oriented along x-direction: a and b for an applied voltage of 1 and 100 mV to the drain contact while the source contact is grounded, respectively.The Fermi levels indicated in the figures correspond to the Fermi levels of the source and drain contacts.The corresponding effective 1D potentials for the ideal device and the device with the dipole are shown as blue dashed and red lines, respectively.L gap = 10 nm, N D = 10 14 cm −2 , N A = 5 × 10 17 cm −3 , and t = 1 nm.

Figure 7 .
Figure 7. Electrostatic potential in equilibrium.It shows the effective 1D electrostatic potential in equilibrium along the x-direction for a dipole oriented along the x-direction in a, along the y-direction in b, and along the z-direction in c.The results are presented for both dipole lengths, l dipole = 0.8, 4 nm.The figures also include the electrostatic potential for the ideal δ -layer tunnel junction device represented by the black line.t = 1 nm, L gap = 10 nm, N D = 10 14 cm −2 and N A = 5 × 10 17 cm −3 .

Figure 8 .
Figure 8. Dwell time.It shows the dwell time of electrons in the tunnel junction in the presence of a dipole in a and b for an applied bias of 1 and 100 mV, respectively.The results are presented for both dipole lengths, l dipole = 0.8, 4 nm, and all considered orientations.The figures also include the dwell time for the device without defects, as well as with n-type and p-type impurities.t = 1 nm, L gap = 10 nm, N D = 10 14 cm −2 and N A = 5 × 10 17 cm −3 .

9 / 12 Figure 9 .
Figure 9. Tunneling resistance for Si: P δ -layer tunnel junctions.Comparison of the resistances computed in our previous work 23 , and measured and computed by M. Donnelly et al. 31 for diverse tunnel gap lengths.In our study, the dimensions of the δ -layer are 7 nm in width, 1 nm in thickness, and a doping density of N D = 10 14 cm 2 .In M. Donnelly et al.'s work, the dimensions of the δ -layer are approximately 7 nm in width, 0.25 nm in thickness, and a doping density of N D = 2 × 10 14 cm 2 .