Correlated impurity complex in the asymmetric tunneling contact: an ideal system to observe negative tunneling conductivity

We studied theoretically electron transport through the impurity complex localized between the tunneling contact leads by means of the generalized Keldysh diagram technique. The formation of multiple well pronounced regions with negative tunneling conductivity in the I-V characteristics was revealed. The appearance of negative tunneling conductivity is caused by the presence of both strong Coulomb correlations and the asymmetry of tunneling rates, which lead to the blockade of the electron transport through the system for a certain values of applied bias. The developed theory and obtained results may be useful for the application of impurity (dopant) atoms as a basic elements in modern nanoelectronic circuits.

We studied theoretically electron transport through the impurity complex localized between the tunneling contact leads by means of the generalized Keldysh diagram technique. the formation of multiple well pronounced regions with negative tunneling conductivity in the i-V characteristics was revealed. the appearance of negative tunneling conductivity is caused by the presence of both strong Coulomb correlations and the asymmetry of tunneling rates, which lead to the blockade of the electron transport through the system for a certain values of applied bias. the developed theory and obtained results may be useful for the application of impurity (dopant) atoms as a basic elements in modern nanoelectronic circuits.
Nowadays electron transport through impurity complexes attracts strong attention. It occurs because impurity (dopant) atoms are considered to be both promising candidates for the implementation in semiconductor nanoelectronic devices [1][2][3][4] and model systems for quantum transport phenomena investigation [5][6][7][8] . Application of individual atoms as building blocks of nanoelectronic devices is very perspective, as they have a stable well defined electronic structure. Individual atoms embedded in the semiconductor medium allow to fabricate unique single-atom single-electron tunneling devices as prototypes of quantum logic gates 9 , quantum bits [10][11][12] , charge pumps and turnstiles [13][14][15][16] , etc.
One of the most vital properties for electronic components functioning is the presence of negative tunneling conductivity 17 , which was observed for the first time in a highly doped tunneling diods 18 . Later negative tunneling conductivity has been revealed in a quite few low temperature experiments of electron transport through molecules 1,2,[19][20][21][22][23] , dopant atoms 24 and quantum dots 25 . For molecular systems the presence of negative tunneling conductivity was demonstrated both at low and room temperatures 1,2,19,21 . Usually, negative tunneling conductivity occurs at the back of sharp current peaks and is explained by the charge transfer between two energy levels which are in resonance at a certain value of applied bias and are out of resonance for other values of applied bias.
One of the first theoretical explanations of negative tunneling conductivity was proposed by Likharev with co-authors 8 . They attributed formation of negative tunneling conductivity with the enhancement of one of the two tunneling barriers of the transistor by the source-drain electric field. This mechanism made it possible to explain experimental results obtained in molecular systems 19,20 . Another possible mechanism for the appearance of negative tunneling conductivity in molecular electronic devices was proposed in 22 and originated from local orbital symmetry matching between an electrode and a molecule in a single molecular electronic device. However, mechanisms considered in 8,22 could not explain formation of negative tunneling conductivity in a small size strongly correlated systems 25 , where the presence of strong Coulomb interaction plays an important role and should be taken into account. Later it was demonstrated that in a low dimensional strongly correlated structures negative tunneling conductivity arises due to the presence of non-equilibrium Coulomb interaction [26][27][28] . Inter-particle correlations should be carefully taken into account as they could drastically affect the local charge distribution in the vicinity of impurity complexes in nanometer-size tunneling junctions. Moreover, interacting impurities is one of the most promising systems for the analysis of Coulomb correlations influence on the properties of the electron transport [28][29][30][31][32][33][34] .
Electronic circuits based on the individual impurities are recently under active investigation, therefore understanding the role of strong inter-particle correlations in such systems and analysis of their influence on the electron transport properties would be an important milestone on the way of single atom nanoelectronic devices creation. Here we theoretically analyze electron transport through the impurity cluster with strong Coulomb correlations between localized electrons by means of the generalized Keldysh diagram technique. We reveal that the presence of strong Coulomb correlations in the asymmetric tunneling contact (coupling strength between the impurity complex and metallic contacts strongly differ) leads to the formation of multiple well pronounced regions with negative tunneling conductivity in the I-V curves. This finding is very promising in the sense of dopant atoms application as basic elements in modern nanoelectronic circuits.

Model System and theoretical Approach
Model system. Further we will consider a model two-level impurity system localized between metallic tunneling contact leads. Two-level system is usually applied for the analysis of electron transport through the correlated impurities. Single electron energy levels ε i could be associated with two different impurities or they could both correspond to the single impurity. We would like to mention that energy spectrum of correlated impurity complex could be rather complicated depending on the strength of coupling between the impurities and the substrate, the value of inter-particle interaction, type of the substrate, etc. However, the situation when only two levels of the impurity complex contribute significantly to the electron transport could be realized by means of the applied bias and gate voltage tuning. Two-level system is relevant, when only two energy levels from the whole impurity complex spectrum are localized in the energy gap E F < ε i < E F −eV. Inequality means that only these two levels contribute to the tunneling current. The Hamiltonian of the two-level system could be written as Hamiltonian Ĥ imp describes impurity complex and includes on-site and inter-site Coulomb interaction in the Hubbard form: i i is localized electron occupation numbers operator and operator σ c i destroys electron with spin σ at the energy level ε i . σσ′ U ij is the Coulomb repulsion between localized electrons. Further we assume that on-site Coulomb repulsion σσ′ U ii exceeds inter-site Coulomb repulsion σσ′ U ij . Such an assumption corresponds to the situation, when impurity state localization radius is smaller than the distance between the impurities. Part Ĥ lead describes continuous spectrum states in the metallic tunneling contact leads  Tunneling transfer amplitudes t k(p)i between continuous spectrum states in the leads and the two-level system are considered to be independent on momentum and spin. Further we will analyze only elastic tunneling processes. As it was shown in 35 inelastic contribution to the tunneling current in the frame of adiabatic scheme at low temperature is much smaller, than the elastic one. Thus, it is sufficient to consider only elastic tunneling and take into account Coulomb correlations to obtain negative tunneling conductivity in asymmetric tunneling contact for impurity atoms with deep levels. It was demonstrated that inelastic processes lead to the additional peculiarities in the I-V characteristics, tunneling conductivity and current noise spectrum but it is not the scope of this paper 36 . In this work we analyze the effect of the two-level system intrinsic properties on the electron transport, while effects caused by the leads, such as the band width effect or the image charge effect, are not discussed.
We would like to mention, that among the most promising systems for the appearance of negative tunneling conductivity are impurity complexes or QDs embedded in a semiconductor matrix. The peculiarities of electron transport through such systems could be analyzed by means of the tunneling measurements. The most evident technique is the scanning tunneling microscopy/spectroscopy (STM/STS) 37,38 (see Fig. 1b for the scheme of the measurements geometry). The ability to image individual dopant atoms combined with scanning tunneling spectroscopy allows to directly study the transport mechanisms through the impurities. It is important, that impurity atoms could be localized directly on the semiconductor surface or several nanometers below the surface. In STM measurements semiconductor substrate and a tip of the scanning tunneling microscope form leads of the tunneling contact. In modern tunneling experiments typical current values can be of the order of 10 pA-10 nA. Tunneling rates could vary from 10 μeV to 100 meV depending both on the dopant atoms position and on the distance between the STM tip and the semiconductor surface. Another promising possibility to study properties of the electron transport through the impurity complexes deals with tunneling through a nanobridge where the Green's functions formalism. Several theoretical approaches are usually used for analysis of the electron transport through the atomic-scale devices. Among them are approaches based on the Hubbard operators 41,42 , pseudo-particle approach 43,44 , well known non-equilibrium Keldysh diagram technique or equations of motion formalism. The problem of Hubbard operators approach deals with the non-trivial commutation rules for Hubbard operators, so it is rather difficult to consider high order correlations only by means of this approach. That is why in 41,42 authors tried to combine non-equilibrium diagram technique with theoretical approach based on the Hubbard operators. Theoretical scheme based on the pseudo-particles seems to be more convenient in some cases as it allows to generalize Keldysh diagram technique with full account of constraint on the total number of pseudo-particles 43,44 . Here we tried to avoid all these difficulties introducing non-equilibrium Green functions in the operator form and performing averaging over localized electron states at the very last stage in the equations of motion. It allows to obtain the spectrum weight of Green functions with account of high order correlations, which are necessary for the proper treatment of strong Coulomb interaction. Such an approach gives us the possibility to analyze carefully electron transport through the two-level strongly correlated structure and to reveal multiple well pronounced areas with negative differential conductivity in the I-V curves.
Let us introduce operators for retarded ′ ( , ) Green functions as a products of fermion creation and annihilation operators 45 : , with indexes ν, μ = k, p, i, j. Indexes i, j describe impurity complex and indexes k, p correspond to the continuous spectrum states in the leads. Equations of motion for the impurity complex retarded Green functions operators where τ = t − t′ and square brackets […] denote the commutator, δ(τ) is the Dirac delta function and Hamiltonian Ĥ 0 is given by Eq. (1) without tunneling part. We consider the situation, when non-crossing approximation occurs, so the following ratios between the system parameters are valid Δε i is the difference between the energies of the single electron levels. In the stationary case total electron occupation numbers are time independent, so one can obtain the Fourier transformed retarded Green function operator ω σ G ( ) R 11 in the following form: 0 ( ) describe electron transitions between the impurity complex and the tunneling contact leads. ν 0L(R) is the continuous spectrum density of states in the leads.
could be obtained from Eq. (7) by the following indexes substitution 1 ↔ 2. Imaginary part of the retarded Green function operator given by Eq. (7) after averaging over localized electrons states directly determines the local density of states in one of the impurities depending on the charge and spin configuration of the whole system considering correlation effects. Poles of each Green function give the energy spectrum of the electrons localized in the vicinity of impurity. In the case of half-filling without taking into account second and high order correlation functions and Kondo effect all the states with the different number of electrons demonstrate the same spectral weight. Tunneling through the strongly correlated structure results in the difference between the spectral weights, as now second and high order correlation functions should be taken into account. For fermions the following relations take place , so all the terms in the perturbation series for the retarded Green function with the different numerators vanish in Eq. (7). Total retarded Green function is obtained as a sum of terms with uncorrelated denominators with the particular energies ε i , ε i + U i , ε i + U ij etc., and each term should be multiplied by a proper combination of the electron occupation number operators. Averaging of the lesser Green function operator allows to obtain the occupation numbers of the impurity complex energy levels. Moreover, non-equilibrium Green functions formalism makes it possible to take into account the non-equilibrium distribution of tunneling particles caused by the tunneling current flowing and the finite value of applied bias. The stationary equation of motion for the operator Ĝ ( ) ii ω < after decoupling localized and conduction electron states and averaging over the electron states in the leads has the form: Tunneling current operator can be expressed through the non-equilibrium lesser Keldysh Green function operator 45 : Considering Eq. (11) an expression for the tunneling current operator between the contact leads in the frequency representation can be obtained by means of the non-equilibrium diagram technique formalism: T ki ki ki ik During the averaging procedure over reservoir states we decouple electron operators in the leads from the localized electron operators. Thus all orders correlation functions for the localized electron operators in the impurity complex can be taken into account. Expression (12) for the reduced tunneling current operator (averaged over conduction electron states) could be re-written through the advanced, retarded and lesser impurity Green functions operators Considering Eqs. (8) and (9) and using averaging procedure over localized electron states one could re-write Eq. (13) as The retarded (advanced) Green function is obtained after averaging over localized electron states . Finally, averaged tunneling current through the impurity complex could be written as is the Fermi distribution function of electrons in the tunneling contact leads. Fermi levels in the L and R leads are shifted on the value of applied bias eV, so f L (ω) = f F (ω) and f R (ω) = f F (ω−eV). Expression (15) describes tunneling through the all electron states of the considered system (single-, two-three and four-electron states are available in the two-level system). Careful analysis of the tunneling processes through the multi-electron states beyond the mean-field approximation leads to the necessity of the second and high order correlation functions calculation between the electron occupation numbers. In the two-level system tunneling current depends on the localized electron correlation functions up to the third order. In the paramagnetic case when electron occupation numbers with the opposite spins have the same value ( = =  where matrix A depends on the functions N σ . The explicit form of the correlation functions is given in the Appendix section. Expression (15) has a simple form for tunneling through the single electron level The tunneling occupation numbers σ N X ( ) Expression (19) directly demonstrates that tunneling current flowing causes non-equilibrium distribution of the carriers in the intermediate structure. The first term in Eq. (18) describes tunneling through the single-electron states, while the second one corresponds to the tunneling current flowing through the two-electron states. For the applied bias being ε 1 < e V< ε 1 + U occupation numbers N Tσ (ε 1 + U) and are close to zero, so tunneling current is given only by the first term of Eq. (18). With the increasing of applied bias (ε 1 + U < eV) tunneling current is given by the electron transport through both single-and two-electron states. Considering the situation when impurity complex energy levels are well defined (the distance between energy levels strongly exceeds their widths >> ) one could get an expression for the current through the single -electron and two -electron states taking into account Coulomb correlations of localized electrons in all the orders [26][27][28] . When applied bias is smaller, than multi-electron energy levels eV < min(ε i + U ij ; ε i + U i ) only single electron states are available for tunneling. Tunneling current through the single electron states reads: Expression for the tunneling current through the two-electron states (ε 1 + U 12 < eV < min(ε i + 2U ij )) can be written as:

Results and Discussion
We will further analyze electron transport in the asymmetric tunneling contact. The energy spectrum of the two-level system including electron states contributing to the tunneling current depending on the ratio between the tunneling rates is shown in Fig. 1c,d. Further to avoid difficulties with the results presentation we will demonstrate only energy levels, which fall into the range E F < ε < E F + eV and, consequently, contribute to the tunneling current. We would like to stress, that Figs. 2 and 3 represent an increase in the voltage applied between the electrodes, not moving of the impurity complex energy spectrum. Let us start from the situation, when the single electron energy level ε 2 is asymmetrically coupled to the leads. In this case the following relation between the tunneling rates is realized: For the small values of applied bias 0 < eV < ε i (see Fig. 2a) occupation numbers ε N ( ) i T i are close to zero, consequently, following Eq. (20) tunneling current is negligibly small (I T → 0). With the increasing of applied bias the single -electron impurity level with the smallest energy starts to be localized between E F and E F + eV (see Fig. 2b). So, for the applied bias ε 1 < eV < ε 2 occupation numbers ε N ( ) www.nature.com/scientificreports www.nature.com/scientificreports/ Expression (23) reflects the growth of the tunneling current flowing through the impurity complex (see Fig. 4). Further increasing of applied bias (ε 2 < eV < ε 1 + U 12 ) makes the second single -electron energy level available for tunneling (see Fig. 2c). Considering Eq. (20) and relations between the tunneling rates (22) one could find that occupation numbers are ε  N ( ) 1/2 . So, the tunneling current through the single-electron states now reads: Analyzing expression (24), the formation of negative conductivity could be revealed, as the growth of applied bias corresponds to the decrease of the tunneling current (see Fig. 4). It happens due to the influence of both Coulomb correlations and the asymmetry between the tunneling rates. Ratio (22) means that single-electron state ε 2 fills much faster than the state ε 1 and blocks tunneling through the energy level ε 1 due to the Coulomb repulsion. Consequently, tunneling current value decreases.
For the higher values of applied bias (ε 1 + U 12 < eV < ε 2 + U 2 ) tunneling through the two-electron states becomes possible (see Fig. 2d). Now tunneling current is determined by the sum of Eqs. (20) and (21). Occupation numbers aspire to the values ε  N ( ) 1/2 and expression for the tunneling current through both single-and two-electron states (for our system parameters the lowest two-electron state has the energy ε 1 + U 12 ) takes the form (see Fig. 4): Figure 2. Scheme of the tunneling processes through the two -level system in the asymmetric contact with the following relation between the tunneling rates Γ L2 >> Γ L1 ≃ Γ R1 >> Γ R2 . The value of applied bias eV i increases from panel (b) to panel (f). Only energy levels contributing to the tunneling current are depicted for each value of applied bias. Full scheme of the two-level system energy levels contributing to the tunneling current for the different values of applied bias is shown in Fig. 1c.  Figure 3. Scheme of the tunneling processes through the two -level system in the asymmetric contact with the following relation between the tunneling rates Γ L1 >> Γ L2 ≃ Γ R2 >> Γ R1 . The value of applied bias eV i increases from panel (b) to panel (f). Only energy levels contributing to the tunneling current are depicted for each value of applied bias. Full scheme of the two-level system energy levels contributing to the tunneling current for the different values of applied bias is shown in Fig. 1d.  considers tunneling barrier transparency changing for the energy level ε 1 + U 12 as it is localized higher than the single-electron states. Further growth of applied bias (ε 2 + U 2 < eV < ε 1 + 2U 12 ) opens the possibility for the electrons with opposite spins to tunnel simultaneously through the two-electron state (see Fig. 2e). The tunneling rates asymmetry and the presence of Coulomb correlations cause blockade of the electrons, which tunnel through the single-electron state ε 1 . Occupation numbers aspire to the values ε  N ( ) 1/2 and tunneling current decreases. It is now determined only by the electrons tunneling through the two-electron state with the energy ε 2 + U 22 : Expression (26) corresponds to the formation of the second area with the negative conductivity, when the growth of applied bias results in the decrease of the tunneling current. The second minimum appears in the I-V curve (see Fig. 4). Further increasing of the applied bias (ε 1 + 2U 12 < eV < ε 1 + U 1 + 2U 12 ) opens the possibility for electrons to tunnel through the three-electron states (see Fig. 2f). Occupation numbers ε + N U ( 2 ) T 1 1 12 turn to 1/2, and expression for the current flowing through the system reads: takes into account tunneling barrier transparency changing for the three-electron states in comparison with the two-electron states. Further growth of applied bias leads to the formation of an additional step in the current voltage characteristic for eV > ε 1 + U 1 + 2U 12 , when tunneling through the four-electron states becomes possible (see Fig. 4).
We now consider the situation when strong asymmetry between the tunneling rates of the single electron energy level ε 1 takes place: The behavior of I-V characteristics is quite similar to the previously discussed situation. The difference exists in the tunneling current values and the ranges of applied bias where negative conductivity appears.
For the small values of applied bias 0 < eV < ε i (see Fig. 3a) tunneling current is negligibly small I T → 0. The first maximum in the I-V characteristic (see Fig. 5) corresponds to the applied bias range ε 1 < eV < ε 2 . The tunneling current amplitude is determined by the expression In this range of applied bias tunneling through the single-electron state with the lowest energy ε 1 occurs (see Fig. 3b). The increasing of applied bias (ε 2 < eV < ε 1 + U 12 ) opens the possibility for electrons to tunnel through both single-electron levels (see Fig. 3c). Consequently, current through the single-electron states reads: Expression (30) reveals the formation of negative tunneling conductivity (see Fig. 5). It happens due to the presence of both Coulomb correlations and the asymmetry between the tunneling rates. As a results, the blockade of electron transport through the single-electron level ε 2 takes place. For the higher values of applied bias (ε 2 + U 12 < eV < ε 1 + U 1 ) tunneling through the two-electron states becomes possible (see Fig. 3d). Tunneling www.nature.com/scientificreports www.nature.com/scientificreports/ current is then determined by the sum of Eqs. (20) and (21). It corresponds to the tunneling processes through the single-electron states and the two-electron states: 1 . Figure 5 demonstrates growth of the tunneling current amplitude and formation of the second maximum in the I-V curve. Further increasing of the applied bias (ε 1 + U 1 < eV < ε 2 + 2U 12 ) leads to the simultaneous tunneling of electrons with opposite spins through the two-electron states (see Fig. 3e). Transport through the single-electron states is now blocked and the second area with the negative conductivity in the I-V characteristic appears (see Fig. 5): With the increasing of applied bias tunneling through the three-(see Fig. 3f) and four-electron states becomes possible and two additional steps are formed in the I-V curves.
We would like to mention, that in the absence of Coulomb interaction between localized electrons only steps would be seen in the I-V characteristics. The presence of strong Coulomb interaction significantly modifies I-V curves. In asymmetric contact areas with negative tunneling conductivity could arise (peaks are well pronounced at particular values of applied bias). The width of the peaks depends on the values of the tunneling rates.
Tunneling rates values are determined by the the barriers widths and heights which are the functions of the particular geometry of the experiment and the symmetry of the localized states electronic orbitals. Moreover, tunneling barrier characteristics (width and height) are changed for the higher energy levels. As Coulomb interaction could be of order of several eV the barrier height for the high energy levels could be strongly changed.
Calculations demonstrated in Figs. 4 and 5 mostly correspond to the I-V curves available in the single-electron transistor geometry 4,39 . Typical current values obtained in the STM experiments are much smaller but all the features (peaks and steps) could be well resolved in both experimental schemes.

conclusions
To summarize, we studied the electron transport through the impurity cluster (modeled by the two-level system) localized between the tunneling contact leads taking into account strong Coulomb correlations and the asymmetry between the tunneling rates. By means of the generalized Keldysh diagram technique we derived general expressions for the tunneling current and obtained formulas, which define electron transport through the states with the different number of electrons. We demonstrated that the interplay between Coulomb correlations and the asymmetry between the tunneling rates result in the formation of multiple regions with negative tunneling conductivity in the I-V curves. We believe, that our results are very promising in the sense of single atoms transistors application in modern nanoelectronic circuits.

Appendix
Expression for the tunneling current (15) includes mean electron occupation numbers σ n i , pair and triple correlation functions for the localized electrons, which have to be determined for a proper treatment of correlation effects. Equations for the total electron occupation numbers σ n 1 and σ n 2 can be found from the conditions: Tunneling current σ I L2 can be easily obtained from expression for the σ I L1 by the following indexes changing 1 ↔ 2 and σ I Ri can be found by the indexes changing L ↔ R. Pair correlation functions can be found as: