#### Abstract

We study the *I*-*V* characteristic of mesoscopic systems or quantum dot (QD) attached to a pair of superconducting leads. Interaction effects in the QD are considered through the charging energy of the QD; that is, the treatment of current transport under a voltage bias is performed within a coupled Poisson nonequilibrium Green function (PNEGF) formalism. We derive the expression for the current in full generality but consider only the regime where transport occurs only via a single particle current. We show for this case and for various charging energies values and associated capacitances of the QD the effect on the *I*-*V* characteristic. Also the influence of the coupling constants on the *I*-*V* characteristic is investigated. Our approach puts forward a novel interpretation of experiments in the strong Coulomb regime.

#### 1. Introduction

The overall shape of the characteristic of a variety of systems (metals, semiconductors, and molecular conductors) in the nanometer scale sandwiched between metallic or superconductors leads has been recently a matter of study (see [1, 2] and references therein). In these systems, the energy level discreteness is quite important since level spacing is comparable with other energy scales [3, 4]. Indeed, the coupling with the bath modifies drastically the properties of an otherwise uncoupled nanometer system in a sharp contrast with similar nonequilibrium macroscopic systems [5–13]. They constitute hybrid systems. Theoretical studies [14–20] as well as experimental measurements have been done by many research groups [1–3] on such systems mostly at low enough temperature with negligible thermal and non-equilibrium fluctuations.

All the systems mentioned above underlay universal common features with the hybrid superconductor quantum dot devices we want to address in this work [2, 4, 21, 22]: (i) broadened energy levels of the quantum dot due to hybridization with the leads; (ii) spatial potential profile. (iii) a charging energy due to the potential profile. An insight behind these issues has been highlighted recently [23, 24] for molecular dots. The device we study in this work is shown in Figure 1. It constitutes a spin degenerated quantum dot level, which is coupled to a pair of biased superconductors contacts or leads (source and drain). When a source-drain voltage is applied, an electric current flows between the leads and across the quantum dot. The biasing defines a non-equilibrium steady state situation. Such situation is coming from the frustration to establish simultaneously an equilibrium configuration with both leads under a given bias. In addition, a gate voltage sets the quantum dot spectrum. However, the charge energy can modify it whenever the density of states is significant. In response to the applied voltages, an actual potential develops inside the dot; that is, an effective electrostatic profile potential inside the mesoscopic region exists in such a way, that it couples to both the electronic non-equilibrium state population and the non-equilibrium electric current. That approach, as introduced by Datta [4], links the electrostatic profile to the electronic population of the quantum dot [4, 25] via the non-equilibrium Keldysh formalism (NEGF) [26, 27]. The whole system is modeled by coupling capacitances which represents the drain, source, and gate contributions to the self-consistent electrostatic problem. Incoming electrons have to overcome an energy barrier (Coulomb blockade). On the other hand, gate or source-drain voltage can lower or increase this energy barrier. These source, drain, and gate electrodes capacitances (see Figure 2) constitute a simple capacitive model (in experiments [2, 28], these capacitances are measured) from which , the Laplacian part of the potential, can be obtained. In addition, the charge in dot can be expressed as the sum of the charges in the coupling capacitances. It yields the Poisson contribution to the total potential , as a function of the dot population. In other words, we solve the self-consistency (SC) of the total electrostatic potential together with the dot population. After that, the electric current is evaluated.

Previous to the self-consistent program, the non-equilibrium current through the dot and electronic occupation in the dot are worked out. We emphasize that the calculation is carried out in a general framework. However, we confine our attention to the single particle current contribution. We adapt the SC to two different approximation regimes. In Section 4, the equivalent capacitive circuit (Figure 2) is introduced; the spatial potential profile is calculated within the capacitive model. The SC scheme is applied to two cases [29, 30]. First, the so-called restricted case, where the gap is the bigger energy scale and the coupling QD-Leads is of the order of the charging energy (). In this case, quantitative results are expected to be accurate. We also make calculations for the so-called unrestricted case, where the charging energy is the dominant energy scale . In this case the results are quantitatively less accurate. The experiments of Ralph et al. [28] were done in this regime. Their characteristic shows that the spacing of the energy levels is subjected to strong fluctuations. According to our model, the fluctuations are due to complex multilevel charging effects. Our hybrid S/QD/S system has been studied in previous theoretical works [16–19]. However, to our knowledge, the coupled SC scheme which describes charging effects has not been considered so far. This is an important step; then, gauge invariant independence of the results as well independence of the zero reference voltage is fulfilled [31, 32]. Our model uses experimental values of the equivalent capacitances [4]. To this respect, pioneering work is done by Meir et al. [33, 34] for N/QD/N systems, considering the interatomic Coulomb term as a measure of the charging energy . Their purpose was to find the main object of the non-equilibrium formalism, namely, the QD Green-Keldysh function, in which the influence of the leads on the QD is taken into account. Due to the presence of the Coulomb term, its equation of motion generates a two-particle Green-Keldysh function. By ignoring correlations with the leads, the equation of motion for the QD Green function closes after truncation of higher order equations of motion. This solution (their Equation ) has two resonances, one at the energy level weighted by the probability that the other spin degenerate level (raised by ) is vacant and another one at the energy level raised by weighted by the probability that the level is occupied. It is correct for temperatures higher than the Kondo temperature and is exact in the noninteracting limit () and the isolated limit. Analogously, for S/QD/S hybrid systems Kang [16] has obtained an expression for the current through the QD (his Equation ), which is evaluated in the limit (his Equation ). The QD Green function from the very beginning does not contain off-diagonal terms that involve superconducting pairing, which excludes the possibility of Andreev reflection processes. The presence in the equation for the current (his Equation ) of terms proportional to affects the contribution to the current of the considered level. In order to complete the outlined program one has to calculate self consistently which is not carried out. Instead, Kang calculate the current (his Equation ) where the spectral function is calculated in the limit of zero coupling with the leads via a model taken from literature (his Reference [18]) and without taking into account the dependence of the contribution of one level to the current on the occupancy of the other. The point of view which neglects the unavoidable influence of the bath (the leads) on the small system (the QD) is accomplished by factorizing the density matrix () and integrating out the leads degrees of freedom which simplify the Liouville-von Neumann equation (Equation in [35]). This program is carried out by Kosov et al. for S/QD/S system [36]. In this way, a Markovian master equation is obtained and an expression for the current is calculated. In their Figure 2, they show the characteristic of a nondegenerated QD for a given set of parameters. In this case, the Cooper pair density in the QD is zero [37]. For the sake of comparison, we restrict our calculations to this case. A similar but not identical approach was done by Pfaller et al. [38]. Also, the approach of both Kosov et al. and Pfaller et al. misses the energy levels broadening as discussed in the introduction. This lack of broadening is a general deficit of quantum Markov approach [39]. In particular, Pfaller et al. [38] introduce a phenomenological broadening while our approach derives it from first principles. In fact, within the Keldysh formalism, this broadening appears naturally (see (58) below). Levy Yeyati et al. [17] writes an expression for the current (his Equation and Figure 2). They use that expression to explain the experimental results of Ralph et al. [28]. Their calculation was done in the limit. In addition, they include charging effects, although they do not say explicitly in which way these effects are included. In this respect, one has to realize that has important contributions to the QD mesoscopic charging effect. In , the leads and the QD maintain independent thermal equilibrium, that is, are uncoupled systems. When they become coupled, the Keldysh formalism yields the general behavior of the system. After a long enough time, this particular system reaches a steady state.

Our point of view is taken from the fact that the charging of the QD is the origin of the Coulomb repulsion between two electron is occupying a two-fold degenerate level. Therefore, we study the behavior of a noninteracting QD at where exact expressions are found. In this way, we obtain a formally similar expression (see (61) below) for the current as Equation in the work of Meir and Wingreen [33]. Later on, Coulomb repulsion is introduced via a self-consistent field (SCF) that depends dynamically on the applied bias () and, in consequence, on the actual number of electrons in the QD. This approach constitutes the coupled Poisson NEGF formalism that has been discussed in the context of molecular conductors by Datta et al. [4, 24]. We use a capacitive model in Section 4 to calculate and as discussed above, a numerical procedure is used to evaluate the current. Our approach has the known disadvantage of ignoring correlations in the QD (as pointed out in [39]). In that sense, there is a proposal by Datta (Equation in [4]) that improves the SCF method and permits more accurate quantitative results. In Section 4, we apply this improvement for the case when the Coulomb charging is greater than the value of the coupling constants. We discuss possible improvements of our approach in Section 6.

#### 2. Single Level QD-Model: Derivation of Nonequilibrium Currents

In macroscopic systems, the task of deriving transport equations or generalized Ginzburg-Landau equations relies on quasiclassical Green functions [7]. In addition, recently non-equilibrium transport in dirty Aluminium quasi-one-dimensional nanowires coupled with normal reservoirs [11] was studied experimentally and theoretically with quasi-classical Green functions [13]. As we want to include the possibility of particle interference effects, we do no resort to such objects. This point of view has been discussed in [14]. Instead, we use the equation of motion method (EOM) technique of Keldysh formalism for generating non-equilibrium states (see [8–10, 26]). We consider a spin degenerated single orbital as a quantum dot connected to superconductors leads. The Hamiltonian which describes this system is a generalized Anderson model [40]. It reads where , , and stand for the superconducting leads, the dot, and the tunneling term, respectively. Also , where and are the left and right lead Hamiltonians, respectively. They are given, within the BCS model [41], by with where is the conduction electron energy and is the superconductor gap of the lead . and are the Nambu spinors: Here denotes the creation (annihilation) operator for a conduction electron with wave vector and spin in the superconductor lead.

is the hamiltonian for the single-level quantum dot of energy : with

The model QD does not contain the Hubbard Coulomb repulsion interaction term. As explained in the introduction, Coulomb repulsion is modeled by means of the inclusion of capacitances, which are taken independent of the charge in the QD. The model also ignores possible superconducting correlations in the QD. For sufficiently small QDs, the discreteness of the single energy levels suppresses these correlations [37]. The position of the energy level will be treated first as fixed by the gate potential with respect to the left lead, while the effect of the applied voltage is taken into account by the coupled Poisson scheme. The tunneling hamiltonian is given by with connects the dot to the biased superconducting leads and it allows the electric charge flow. is the hybridization matrix element between a conduction electron in the superconductor lead and a localized electron on the dot with energy . and are the dot spinors: here, is the creation (annihilation) operator for an electron on the dot.

The flow of electric charge from the terminal is given by
where is the electron charge. is the thermodynamical average over the biased and leads at the temperature , taken at time , as indicated in the Keldysh contour in Appendix A:
and is the “number of particles” operator. Book-keeping calculation using (10) leads to
is the lesser Keldysh Green function,
and T_{K} is the time-ordering operator, the action of which is to rearrange product of operators, such that operator with later times, on the Keldysh contour are placed to the left of the product. Hereafter, for simplicity, we replace by an average at the Fermi surfaces () of the leads and . Using the scheme given in Appendix A for the rate of change of (13), we proceed to obtain the equation of motion:
which leads to
where
Note that is the QD single-particle Green function. Similarly, satisfies the equation of motion:
where
Here is the QD of two-particle Green’s function.

Equations (15) and (18) can be written in a compact form as follows (see Appendix B): We introduce the tilde Keldysh-Green functions: Consider the following matrix whose elements are the unperturbed Green-Keldysh functions, that is, defined for : where

According to Appendix A, their equations of motion are given by These equations can be written in matrix form as follows:

Equation (20) can be written as an integral along the Keldysh contour C_{K} (for an explanation see Appendix B)

From the last expression one can read for the following equation:

We now apply the procedure explained in Appendix C; in order to obtain the lesser component, we obtain Furthermore, the superscripts , , , and correspond to lesser, greater, retarded, and advanced Green’s functions, respectively.

Therefore, from (12), can be written as follows: with

When applying the Fourier transformations, (30) and (31) can be expressed as follows: with

In Appendices D to G, we evaluate the unperturbed Green’s functions , , , and in the wide band limit.

We summarize these results as follows:

All these expressions will used below.

#### 3. QD Green Function

We need to evaluate the most important objet for calculations, namely, QD Green’s functions given by (17) and (19), as well as their respective tilde functions:

Again using the scheme given in Appendix A, their equation of motion iswhich develops to This can be written as follows: When , one has with

The last two equations can be written as follows:

We write the last equation in its equivalent convolution integral along the Keldysh contour (see Appendix B):

An equivalent way to write the last equation (using (25)) as a convolution of and is with

We are interested in two regimes: a first regime in which and the Coulomb blockade effects are neglected because in this case the couplings to the leads are not extremely small and the dot capacitance is large enough, a second regime for where Coulomb blockade effects must be taken into account. For both regimes and from now on, we are interested in the case , where multiple Andreev reflection [42] processes are strongly suppressed. Therefore only the single particle current () has to be considered . From the above considerations, we have that the Keldysh Green function , which carries information of the quantum dot two-particle Green's function, can be neglected and all relevant information is contained in .

The Keldysh Green function becomes spin independent; . Element 11 of (43) is given by

Again, using the recipe given in Appendix C, we obtain for and the following: Taking the Fourier transform of (46) results in a set of algebraic equations: Dot Keldysh Green’s functions and are below straightforward evaluated. In this regime, quantities such as currents are independent of time. Therefore, we have Therefore (47) result in Solving (50), Moreover, we know that Resulting in Equation (49) is reduced to Here is the so-called quantum dot spectral function which is given in terms of the imaginary part () of the retarded Keldysh Green function : From (30), the single particle current () results in Substituting (51) and (55) in (57), with and . In our regime, ; therefore, in the above equations is zero. We use the expression for from Appendix D and obtain the single particle current : corresponds to the applied voltage between the superconductors electrodes with chemical potential . In the following, we fix the chemical potential and use as a measure of . In addition, the QD energy is measured with respect to . On the other hand, the limits of integration are given by the functions and . The extra factor arises from the dot Keldysh Green functions. and are given by

At steady state there is no net flow into or out of the mesoscopic channel or quantum dot which yields a stationary particle number in it. The population number , at the dot, is given by which becomes a weighted average over the and contacts: For the N/QD/N case, are just constants. This case was studied in the context of the generalized quantum master approach (section IV in [39]). That approach permits the inclusion of broadening in a natural way. They obtained Equations similar to (59)–(63).

#### 4. Coupled Poisson Nonequilibrium Green Function Scheme: The Capacitive Model

So far, we are not including the side effects of a potential profile inside the mesoscopic channel. On the one hand, its inclusion takes in order zero or Hartree approximation the electron-electron interaction in the QD. Its inclusion also guarantees current independence from the choice of zero potential [32]. Such potential is induced by the action of source, drain and gate applied voltages. In principle, we have to couple the number of population equations. Equation (62), with electric field . However, since the number of quantum levels in the channel is small and the particle number variation is negligible, the potential profile variation inside the channel is negligible. Then it is appropriate to visualize the channel as an equivalent circuit framework (Figure 2). In this framework, we associate capacitances , , and with the drain, source and gate, respectively. Whenever drain, source, and gate bias potentials , , and , respectively, are present, there is an electrostatic potential inside the QD, which induces an energy shift of the QD energy level , are channel electrostatic potential before we apply the source and drain biases, respectively.

The electronic populations before and after we apply the biases mentioned above are given by respectively. It leads us to where . Therefore, the energy shift is given by where In the expression for , represents a uniform shift for all levels, whereas the second term (the Poisson contribution denoted by in the introduction) represents a level of repulsion which is proportional to the averaged occupation of the QD level denoted by and proportional to the charging energy .

On the other hand, one has from (62) and (65) given by In the expression for (see (55)), the energy level shifts only () in the expression for the QD spectral function . Equations (66) and (68) are coupled nonlinear equations with unknowns and . We solve the coupled equations via an iteration procedure. First we guest a value for plug this value in , and then we calculate with the following equation: and so on until convergence is achieved. With the final value of obtained for a given bias voltage , is calculated via the equation:

In summary, the procedure for computing consists of the following steps. (i) Determine the spectral density. (ii) Specify , , , and coupling constants. (iii) Iteratively solve (69) and (66). (iv) Evaluate the current from (70) for the , , and . Once a converged has been found, the current is finally evaluated.

The way we consider electron-electron interactions imposes restrictions on the possible values of the charging energy . For the self-consistent scheme to be valid, we have to assume that . However, less precisely quantitative results, although qualitative correct results can be obtained if , when the so-called Coulomb Blockade energy dominates over the coupling constants. For this case, we use the improvement of the SCF method discussed in the introduction [4, 29, 30]. The self-consistent generalizes to where the up-spin level feels a potential due to the down-spin electrons and viceversa. Notice the different signs which reflects the Coulomb repulsion between otherwise degenerate levels: Here, and are the population of the spin-up and spin-down levels. Once the values of and are calculated, is calculated from As Datta has pointed out [4], the approach described above (called unrestricted SCF) can lead to a better quantitative agreement in comparison with a conceptually correct multi-level Master equation calculation.

#### 5. Numerical Results and Remarks on Experiments

##### 5.1. First Case:

In this regime, there are the multiple Andreev reflections [42] for voltages such that (MAR). Also, there is the possibility for quasiparticle co-tunneling current for energy levels far from . These cases will be considered in a future work that involves the whole expressions we have derived for the currents (see (30) and (31)) and eventually more accurate Green-Keldysh functions and the use of master equations [17]. This case was studied experimentally in [43]. For given values of capacitances and source voltage, we iterate (66) and (68) in order to find the potential . Then the single particle current is evaluated (see (70)). We put the charge before biasing , such that Coulomb repulsion with the QD-energy level is absent. Anyhow, in this regime, the effect of the second term in (74) is negligible. Consequently, the Laplace term essentially positions the QD degenerate energy level (with respect to ). In Figure 3, we show characteristic for gate voltage values and , whereas Figure 4 shows the occupation number . These curves are symmetric, due to the assumed equality of the coupling capacitances (). Otherwise, the shifts to right or to the left for or respectively. For this case, we show in Figure 7 the spectral density for . Notice that the position of the energy level is essentially , that is, just the sum of . Qualitatively, these results are similar to Levy Yeyati et al.’s [17]. Characteristic is the broadening of the BSC singularity. The effect of bigger values of is a more pronounced round off the BCS-type singularity. We discuss this issue below. For large enough bias, the current approaches the normal saturation value I_{Sat}.

##### 5.2. Second Case:

In this regime, the charging energy acts effectively in lifting the degeneracy of the otherwise single degenerate QD energy level. For this regime, we use the couple system defined by (71)–(75) and calculate the current according to (76). This is the unrestricted SCF method mentioned in the introduction. The transport begins through one level as long as there is in average less than one electron in it. For the given parameters, the onset of current is similar to the first case (no interaction with residual charge in the QD is considered). However, when the average occupation exceeds one, the other degenerate levels float according to the resulting values and . These values push down the position of this second level and push up the already occupied energy level. In Figures 5-6, we show and the number of electrons. In Figure 8, it is shown that the spectral density for . In this case, . The values obtained from the SCF calculation position the energy levels to and to (see (72)).

The experimental work of Ralph et al. [28] corresponds to this second case. In their Figure 3, they show the current for single state level peaks pA at a bias mV. According to Figures 5-6, for this sample at 30 mK, there should be another peak at mV with a current value pA. The charging energy for this sample is mV. However, in their samples with radios ~2.5 nm or greater, the level spacing of the energy levels is such, that ; therefore another current signal may occur before and a quantitative proper description would be a multilevel QD-model. One notice, however, the strong fluctuations in the spacing in Figure 2 [28], indicating complex charging for many levels of phenomena. Notice that the theoretical explanation of Levy Yeyati et al. [17, 18] does not contain our prediction.

##### 5.3. Influence of Coupling Constants

The influence of the coupling constant is to broaden the otherwise sharp energy QD level. However, the broadening is not equally strong and depends on the relative values of . Ralph et al. [28] determined for the sample in his Figure 3 . Figure 10 shows the characteristic for the restricted case when both coupling constants are equal. For larger values of the coupling constants, the broadening is stronger. If they are dissimilar in value, the broadening is stronger when . This effect is shown in Figures 9, 10, 11, 12, 13, and 14. This effect is due to the stronger involving of the BSC-DOS singularity of the left lead in the integral expression for the current (see (59)) and the particular choice of the zero bias voltage (see Figure 2).