Single Charge Current in a Normal Mesoscopic Region Attached to Superconductor Leads via a Coupled Poisson Nonequilibrium Green Function Formalism

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 U 0 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.


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][6][7][8][9][10][11][12][13]. They constitute hybrid systems. Theoretical studies [14][15][16][17][18][19][20] as well as experimental measurements have been done by many research groups [1][2][3] on such systems mostly at low enough temperature with negligible thermal and nonequilibrium 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 0 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 [ 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 sourcedrain 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 L , 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 = L + 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 socalled 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. Theircharacteristic 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][17][18][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 2 / . 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 (8)) 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 ( = 0) 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 (8)), which is evaluated in the → ∞ limit (his Equation (13)). The QD Green function from the very beginning does not contain offdiagonal terms that involve superconducting pairing, which excludes the possibility of Andreev reflection processes. The presence in the equation for the current (his Equation (14)) of terms proportional to (1−⟨ − ⟩) 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 (8)) 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 ( ( ) = QD ⊗ Baths ) and integrating out the leads degrees of freedom which simplify the Liouville-von Neumann equation (Equation (3.140) 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 Scientific World Journal 3 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 (2) 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 (12) 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 ( QD + SCF ) 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 SCF 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 (3.4.9) 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.

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-onedimensional nanowires coupled with normal reservoirs [11] was studied experimentally and theoretically with quasiclassical 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][9][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 , QD , 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 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.
QD is the hamiltonian for the single-level quantum dot of energy 0 : 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 4 The Scientific World Journal 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 0 . † 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 0 → −∞, as indicated in the Keldysh contour in Appendix A: and = † ⃗ ⃗ is the "number of particles" operator. Book-keeping calculation using (10) leads to 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 ( ⃗ ≡ √⟨| ⃗ | 2 ⟩ FS ) 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 G ( , ) is the QD single-particle Green function. Similarly, F ⃗ ( , ) satisfies the equation of motion: where Here G ( , ) 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 2 × 2 matrix whose elements are the unperturbed Green-Keldysh functions, that is, defined for = 0: where The Scientific World Journal 5 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 F ⃗ ( , ) the following equation: We now apply the procedure explained in Appendix C; in order to obtain the F ⃗ ( , ) lesser component, we obtain Furthermore, the superscripts (<), (>), (r), and (a) correspond to lesser, greater, retarded, and advanced Green's functions, respectively. Therefore, from (12), ( ) can be written as follows: When applying the Fourier transformations, (30) and (31) can be expressed as follows: 6 The Scientific World Journal In Appendices D to G, we evaluate the unperturbed ⃗ ( , ), and f < ⃗ ( , ) in the wide band limit. We summarize these results as follows: (34) All these expressions will used below.

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 is which develops to This can be written as follows: 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 G ( , ) is ) , ) .
We are interested in two regimes: a first regime in which 0 ∼ Γ < Δ 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 0 ∼ Δ > Γ 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 (SP) has to be considered SP . From the above considerations, we have that the Keldysh Green function G ( ), which carries information of the quantum dot two-particle Green's function, can be neglected and all relevant information is contained in G ( ).
The Keldysh Green function becomes spin independent; G ( ) ≡ G( ). Element 11 of (43) is given by Again, using the recipe given in Appendix C, we obtain for G < ( , ) and G (a) ( , ) the following: The Scientific World Journal Taking the Fourier transform of (46) results in a set of algebraic equations: 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 (I) of the retarded Keldysh Green function G (r) ( ): From (30), the single particle current ( SP ) results in Substituting (51) and (55) in (57), with Γ ( ) = −IΣ (r) ( ) = Γ ( , Δ ) and Γ( ) = ∑ Γ ( ).
In our regime, > Δ; therefore, RΣ (r) ( ) in the above equations is zero. We use the expression for Σ (r) ( ) from The Scientific World Journal 9 Appendix D and obtain the single particle current SP ≡ ( ,SP − ,SP )/2: (59) − = − corresponds to the applied voltage between the superconductors electrodes with chemical potential . In the following, we fix the chemical potential = 0 and use as a measure of . In addition, the QD energy 0 is measured with respect to . On the other hand, the limits of integration are given by the functions Γ ( ) and Γ ( + ). The extra 2 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).

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 QD inside the QD, which induces an energy shift of the QD energy level = − ( QD − 0 ), 0 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 , L 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 0 and proportional to the charging energy 0 = 2 / .
On the other hand, one has Δ from (62) and (65) given by In the expression for G < ( , ) (see (55)), the energy level shifts only ( 0 ⇒ ( 0 + )) 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:

10
The Scientific World Journal and so on until convergence is achieved. With the final value of obtained for a given bias voltage , SP 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 0 . For the self-consistent scheme to be valid, we have to assume that Δ ≫ Γ , ≃ 0 . However, less precisely quantitative results, although qualitative correct results can be obtained if Δ ≃ 0 ≫ Γ , , 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 downspin 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 spindown levels. Once the values of ↑ and ↓ are calculated, SP 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 multilevel Master equation calculation.

Numerical Results and Remarks on Experiments
In this regime, there are the multiple Andreev reflections [42] for voltages such that < Δ (MAR). Also, there is the possibility for quasiparticle cotunneling 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 SP ( , ) is evaluated (see (70)). We put the charge before biasing 0 = 0, 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 L essentially positions the QD degenerate energy level (with respect to = 0). In Figure 3, we show -characteristic for gate voltage values = 0 and = Δ, whereas Figure 4 shows the occupation number Δ . These curves are symmetric, due to the assumed equality of the coupling capacitances ( / = 0.5). Otherwise, the -shifts to right or to the left for / = 0.5 > 0.5 or < 0.5 respectively. For this case, we show in Figure 7 the spectral density ( ) for = −6Δ. Notice that the position of the energy level is essentially −4.5Δ, that is, just the sum of 0 + L . 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 .

Second
Case: Δ ≃ 0 ≫ Γ , . 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 = −8Δ. In this case, 0 + L = −5.5Δ. The values obtained from the SCF calculation position the energy levels to −4.846Δ and to −6.396Δ (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 ∼ 5 pA at a bias ∼ 2.4 mV. According to Figures 5-6, for this sample at 30 mK, there should be another peak at ∼ 4.4 mV with a current value SP ≃ 10 pA. The charging energy for this sample is 0 ≃ 2.0 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.

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 Γ /Γ ≫ 1. 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).

Summary and Perspectives
We have studied the single particle current through a quantum dot coupled with two superconductor leads via a coupled Poisson Nonequilibrium Green function (PNEGF) formalism. In a systematic and self-contained way, we derived the expressions for the current in full generality. In this work we focused only on the weak coupling regime where single particle current is the dominant one. The QD is a single degenerate energy level system modeled via a capacitive circuit. The influence of the potential on the QD and on the - The Scientific World Journal characteristic is calculated for relevance values of the coupling and capacitances and the implication of experiments is discussed. This was done in the weak coupling regime and for Δ ≫ Γ , ≃ 0 . A second case when Δ ≃ 0 ≫ Γ , also in the weak coupling regime was analyzed. Admittedly, our model of a hybrid system S/QD/S possesses potentially physical extensions. One important missed point is dephasing. This physical effect due to scattering of transport electrons can be incorporated in the self energy phenomenologically [44,45], or in a stochastic fashion [46]. Another point is to consider a QD with many energy levels and within the self-consistent scheme, to consider the strong and intermediate regimes and many body correlations due to different kinds of electronelectron interaction. Here we have to notice that it is not just to scale the level spacing by the charging energy [47]. It is a genuine many body problem. But the most important missed point was correlations. As pointed out by Datta [4,chapter III], there has been much effort in order to find a suitable SCF that considers correlations. For example, to modify (66) to consider occupancies probabilities. As discussed in the introduction, Kang [16] and Meir and Wingreen [33] find a solution for the QD Green function that contains this type of correlation. In other words, one could go to scheme where a more accurate Green function for the QD is used together with a multielectron picture and corresponding master equation. These would mean to use the Anderson model with a Coulomb interaction that is obtained from a SCF. We want to check if this point of view is correct. Work in this direction is in progress.

A. Equation of Motion for the Keldysh-Green Function
In general, any Green-Keldysh function of two operators ( ) and ( ) function is given by where the operator T K acts on the Keldysh contour shown in Figure 15. A Heaviside function on the Keldysh contour is given by whereas be derivative of the Heaviside function on the Keldysh contour is given by In general a function F( , ) defined on the Keldysh contour is given by  Figure 15: The contour C K = C K − ∪C K + runs on the real axis, but for clarity its two branches C K − and C K + are shown slightly away from the real axis. The contour C K runs from 0 and returns to 0 . and one obtains i G ( , ) If there is not an applied potential, that is, if = = 0, (see Figure 2), the whole system is at equilibrium, and one can use the usual commutator Green functions or equivalently the Matsubara or Temperature Green function to quantify correlations functions. In that case, (A.1) depends only on the time difference ( − ). However, for times > 0 , once the potential difference has been applied, the simple dependence on the time differences no longer holds which signalizes a nonequilibrium situation. In that case, the Keldysh method applies.

B. Equivalent Integral Equation on the Keldysh Countour on Figure 15
We encounter two cases (see (26) and (42)) where the general strategy to find an integral expression for the Keldysh-Green functions is the following: (1) one first considers the equation of motion for for the Keldysh-Green function of each system (two leads and the QD) (see (25) and (38)). (2) The resulting equation of motion in each case has a delta function as inhomogeneity.
(3) The systems are connected in at 0 , we have as a result two coupled equations of motion (see (20) and (41)).

C. Calculation of Convolutions
In this appendix, we explain how we evaluate the some important convolutions used to calculate lesser Keldysh-Green functions [48,49]. In general, a function given as a convolution on the Keldysh contour poses the definition given in Appendix A. One encounters situations where a Keldysh-Green function is given by a convolution of two other functions: However, for evaluation of quantities like, for example, the current, one needs One sees that the relative position of and divides the contour in three regions of integration: (1) < C K , (2) < C K < , and (3) > C K . This traduces into the following integrals:

D. Calculation of Unperturbed Green Functions and Retarded Self-Energy
Below, we proceed to evaluate the unperturbed Green functions g ⃗ and f ⃗ . To achieve this, it is necessary to introduce the chemical potential shift in each superconductor, because [ , ] = 0. Consequently, the first unperturbed retarded Green function g (r) ⃗ ( , ) is given by where the fermion operators † ⃗ , ⃗ create and annihilate the "Bogoliubov quasi-particles" and The WBL means that the width of the electronic energy bands of the leads is the largest energy. The density of states in the contacts varies on a scale of Fermi energy. These scales are of order 1-10 eV (∼10 4 -10 5 K) which are much larger than the energies involved in the quantum dot ∼ meV ∼ 10 K. Furthermore, Γ ( − ) varies slowly with − and the prefactor D ( − ) varies in the range of wide band and changes − on the average of | ⃗ | 2 occur on the order of meV. Therefore, we ignore the dependence of − in Γ ( − ). The WBL establishes that an electron in the dot decays in an continuum of states of the leads and is sufficient condition for the existence of a stationary state, as has been shown rigorously in [50].