# Josephson and Andreev transport through quantum dots

## Abstract

In this article we review the state of the art on the transport properties of quantum dot systems connected to superconducting and normal electrodes. The review is mainly focused on the theoretical achievements although a summary of the most relevant experimental results is also given. A large part of the discussion is devoted to the single level Anderson type models generalized to include superconductivity in the leads, which already contains most of the interesting physical phenomena. Particular attention is paid to the competition between pairing and Kondo correlations, the emergence of -junction behavior, the interplay of Andreev and resonant tunneling, and the important role of Andreev bound states which characterized the spectral properties of most of these systems. We give technical details on the several different analytical and numerical methods which have been developed for describing these properties. We further discuss the recent theoretical efforts devoted to extend this analysis to more complex situations like multidot, multilevel or multiterminal configurations in which novel phenomena is expected to emerge. These include control of the localized spin states by a Josephson current and also the possibility of creating entangled electron pairs by means of non-local Andreev processes.

###### Contents:

- I Introduction
- II Basic models and formalism
- III Equilibrium properties of quantum dots with superconducting leads
- IV Quantum dots with normal and superconducting leads
- V Voltage biased S-QD-S systems
- VI Beyond the single level model: multidot, multilevel and multiterminal systems
- VII Concluding remarks

## I Introduction

The field of electronic transport in nanoscale devices is experiencing a fast evolution driven both by advances in fabrication techniques and by the interest in potential applications like spintronics or quantum information processing. Within this context quantum dot (QD) systems are playing a central role. These devices have several different physical realizations including semiconducting heterostructures, small metallic particles, carbon nanotubes or other molecules connected to metallic electrodes. In spite of this variety a very attractive feature of these devices is that they can usually be described theoretically by simple ”universal-like” models characterized by a few parameters. In addition to their potential applications, these systems provide a unique test-bed for analyzing the interplay of electronic correlations and transport properties in nonequilibrium conditions.

Electron transport in semiconducting QDs
has been studied since the early 90’s when phenomena like Coulomb blockade (CB) was first observed Kastner (1993).
It soon became clear that QDs could allow to study the effect in transport properties of
basic electronic correlations phenomena like the
Kondo effect as suggested in early predictions Glazman and Raikh (1988); Lee and Ng (1988). These predictions were first
tested in metallic nanoscale junctions containing magnetic impurities Ralph and Burhman (1994).
However, a definitive breakthrough in the field came with the
observation of this effect in semiconducting QDs by Goldhaber et al. Goldhaber-Gordon
*et al.* (1998)
and Cronenwett et al. Cronenwett
*et al.* (1998).
A great advantage of these devices is to offer the possibility of controlling the
relevant parameters, thus allowing a more direct comparison with the theoretical predictions.
Since then the effect of Kondo correlations in electronic transport has been observed in
several physical realizations of QDs based on carbon nanotubes (CNTs) Nygard *et al.* (2000)
and big molecules like fullerenes Roch *et al.* (2009).

In parallel to these advances the study of superconducting (SC) transport in nanoscale devices has also experienced a great development. From a theoretical point of view, with the advent of mesoscopic physics, a more detailed understanding of superconducting transport was developed around the central concept of coherent Andreev reflection (AR) Andreev (1964); Beenakker (1997). This concept has allowed to unify the description of superconducting transport in different types of structures like normal metal-superconductor (N-S), S-N-S junctions and superconducting quantum point contacts (SQPC). Due to the multiple AR (MAR) mechanism the spectral density of systems like S-N-S or SQPCs is characterized by the presence of the so-called Andreev bound states (ABS) inside the superconducting gap. These states are sensitive to the superconducting phase difference and are thus current-carrying states which usually give the dominating contribution to the Josephson effect.

In recent years it has become feasible to produce hybrid systems combining different physical
realizations of QDs well contacted to superconducting electrodes (for a review see Franceschi
*et al.* (2010)).
Superconducting transport through QDs provides the interesting possibility to explore
the interplay of the AR mechanism and typical QD phenomena like CB and Kondo effect.
The central aim of this review article is to discuss the main advances which have
taken place on this issue during the last years.

A usual assumption in these studies is that a basic description of the main properties of
these hybrid systems can be provided by the Anderson model and its generalizations
to include SC leads, orbital degeneracy, etc.
The single level model applies when the dot
level spacing is larger than all other relevant energy scales.
In the normal state the model allows to describe in a unified
way CB and the Kondo effect both in and out of equilibrium conditions Kouwenhoven and Glazman (2001).
With two superconducting leads interesting new
physics already appear in the equilibrium case. Due to the Josephson effect, electron transport is
possible without an applied bias voltage and the model describes the
competition between Kondo effect and induced pairing within the dot.
Figure 1 illustrates this competition: depending on the ratio between
the Kondo temperature, , and the superconducting order parameter, , there is a
phase transition between a Kondo dominated spin-singlet ground state to
a degenerate magnetic ground state. This transition is accompanied by a
reversal of the sign of the Josephson current. Thus in the magnetic
case the S-QD-S system constitutes a realization of the so-called
-junction Glazman and Matveev (1989); Spivak and Kivelson (1964). A more detailed understanding of this
transition describing the appearance of intermediate phases with metastable states
was achieved more recently Rozhkov and Arovas (1999); Vecino *et al.* (2003).
The realization of a -junction in QD systems should distinguished from the
similar phenomena in SFS junctions, where F denotes a ferromagnetic material
Golubov *et al.* (2004).

Another basic situation which has been extensively explored (both theoretically and experimentally) is the N-QD-S case. This situation has been mainly analyzed in the linear transport regime in which it exhibits and interesting interplay between Kondo behavior and resonant Andreev reflection. In contrast to the S-QD-S case this system does not exhibit a quantum phase transition but there is instead a crossover from a Kondo dominated regime for large to a singlet superconducting regime in the opposite limit.

A third paradigmatic situation which has been studied is the voltage biased S-QD-S system. This situation constitutes a much more demanding task for the theory due to the need of describing properly the strong out-of-equilibrium distribution which is generated by the infinite series of multiple Andreev reflection (MAR) processes together with the effects of Coulomb interactions. The problem becomes simpler in two limiting cases: 1) when Coulomb interactions are small and treated in a mean field approximation thus allowing to analyze the interplay of MAR and resonant tunneling and 2) when the Coulomb energy is the larger energy scale in the problem and the contribution of MAR processes are largely suppressed.

More recent developments include the study of several QDs (connected either in series or in parallel) coupled to one or more SC electrodes. In these situations there is a competition between not only the Kondo and the SC correlations but also the possible magnetic coupling of the spins localized within the dots. This could open the possibility to control the spin state of the dots system by means of the Josephson current.

In addition there is a growing interest in
analyzing transport in these hybrid structures in a
multiterminal configuration. One of the aims of
these studies is the detection and
control of non-local Andreev processes, which offer the
possibility of producing entangled electron pairs Recher *et al.* (2001).

This review article is organized as follows: in Section II we introduce the basic theoretical models used to describe the hybrid QD systems, which are largely based on the Anderson model and its generalizations. In this section we also give a brief summary of the application of nonequilibrium Green functions techniques for the calculation of the electronic properties within these type of models. Section III is devoted to review the main results for the S-QD-S systems in equilibrium, i.e. in the dc Josephson regime. We first discuss the case of a non-interacting resonant level which is useful to illustrate the emergence of ABSs and its contribution to the Josephson current. In the subsequent subsections we give account of the different theoretical methods which have been used to analyze the effect of interactions in the dc Josephson regime. A main issue which is discussed in this section are the phase diagrams describing the transition to the -state as a function of the model parameters. We also give a brief account of the existing experimental results for S-QD-S devices in this regime. The case of a QD coupled to both a normal and a superconducting leads (N-QD-S) is addressed in Section IV. Most of the results obtained for this systems correspond to the linear regime with different levels of approximation to include the Coulomb interactions. We also briefly mention existing results for the non-linear regime and the few experiments which have been reported of this case up to date. Section V is devoted to the voltage biased S-QD-S system. We first give some technical details on the calculations for the non-interacting case in order to illustrate how to deal with the out of equilibrium MAR mechanism. We also comment in this section the few existing results including interactions in this regime and give an account of the related experiments. Finally, in Section VI we discuss several different situations which go beyond the single-level two-terminal case discussed in the previous sections. These include: multidot systems connected either in parallel or in series, the multilevel situation and setups in a multiterminal configuration. We conclude this article with a brief discussion of related issues not included in the present review and of topics which, in our view, deserve to be further analyzed in the near future.

## Ii Basic models and formalism

The minimal model for a QD coupled to metallic electrodes in the regime where is sufficiently large to restrict the analysis to a single spin-degenerate level is provided by the single level Anderson model Anderson (1961), with the Hamiltonian where corresponds to the uncoupled dot given by

(1) |

where creates and electron with spin on the dot level located at and is the local Coulomb interaction for two electrons with opposite spin within the dot (). On the other hand, describe the uncoupled left and right leads which can be either normal or superconducting. In this last more general case, they are usually represented by a BCS Hamiltonian of the type

(2) |

where creates an electron with spin at the single-particle energy level of the lead (usually referred to the lead chemical potential, i.e. ) and is the (complex) superconducting order parameter on lead . Finally, describes the coupling between the QD level to the leads and has the form

(3) |

In order to reduce the number of parameters it is usually assumed that the normal density of states of the leads is a constant in the range of energies around the Fermi level of the order of the superconducting gap and that the dependence of the hopping elements can be neglected within this range. The coupling to the leads is then characterized by a single parameter , which can be interpreted as the normal tunneling rate from the dot to the leads.

Within the above model would correspond to the normal state. For vanishing the model is in the so-called atomic limit which is characterized by sharp peaks in the spectral density at and . This limit corresponds to the Coulomb blockade regime in an actual QD where the conductance is strongly suppressed except at the charge degeneracy points. When the couplings to the leads increase ( become larger than temperature) virtual processes allow the charge and spin in the dot to fluctuate and a resonance at around the Fermi energy appears close to half-filling due to the Kondo effect. This simple model thus already captures the most relevant Physics of ultrasmall QDs with well separated energy levels, like the crossover from the Coulomb blockade to the Kondo regime as the temperature is lowered.

The simplicity of this model has allowed to obtain exact results in the equilibrium case by means of the Bethe ansatz Wiegmann and Tsvelick (1983). The most basic of these results is the expression for the Kondo temperature Hewson (1993)

(4) |

where .

This temperature characterizes the crossover from the so-called local moment regime for to the regime where Kondo correlations between the localized spin within the QD and the spin of the electrons in the leads sets in. Although this physics is basically well understood since the 70’s for the case of magnetic impurities in metals, its consequences for transport in artificial nanostructures has started to be developed much more recently specially driven by the advances in fabrication techniques. In this respect while the linear transport properties are well understood still open questions remain regarding the non-equilibrium regime.

In the superconducting case another energy scale, associated with the superconducting gap, appears bringing additional complexity to the problem, whose description is in fact the scope of this review. Even in the equilibrium situation the Anderson model with superconducting leads contains the non-trivial Physics associated to the Josephson effect. A relevant parameter is then provided by the superconducting phase difference .

In order to analyze the electronic and transport properties of a general superconducting system in the presence of interactions and in a non-equilibrium situation it is convenient to use Green function techniques. The Keldysh formalism provides the basic tools for this purpose.

Due to the presence of superconducting correlations it is convenient to introduce the Nambu spinor field operators , with where denotes the electrodes and the dot level respectively. The different terms in the model Hamiltonian of Eqs. (1), (2) and (3) can then be written as

(5) |

where , , , being the Pauli matrices defined in Nambu space.

Starting from these spinor field operators, generalized single-particle propagators can be defined along the Keldysh closed time loop as

(6) |

where denote the two branches in the Keldysh contour. These propagators allow to calculate in a straightforward way most of the relevant quantities like the mean charge and the superconducting order parameter within the dot as well as the mean current through it, which are given by

(7) | |||||

These expressions are formally exact but of little use unless the single-particle propagators are known. Fully analytical and exact results can only be obtained in the non-interacting case. In the presence of interactions numerical methods allow to obtain exact results in the equilibrium case. In a more general case one is bound to find reasonable approximations for these propagators valid for a restricted range of parameters. It is usually convenient to express these approximations in terms of a self-energy which is related to the propagators by the usual Dyson equation

(8) |

where denotes the unperturbed propagators corresponding to an appropriately defined non-interacting Hamiltonian (the symbol indicates matrix structure in Keldysh space) and where integration over internal times is implicitly assumed. Different approximations for the self-energy associated with electron-electron interactions are discussed in the forthcoming sections. The analysis of the problem is greatly simplified in the stationary case where Fourier methods can be applied both in the equilibrium and in the non-equilibrium situation. We shall start discussing in the next section the simplest possible case of an equilibrium situation. In this case further simplification arises from the possibility of expressing all Keldysh propagators in terms of the retarded-advanced propagators and the Fermi equilibrium distribution function, , where is the chemical potential and as

(9) |

Thus, for instance, the mean current can be written as

(10) |

## Iii Equilibrium properties of quantum dots with superconducting leads

For a nanoscale system coupled to superconducting electrodes a finite current can flow even in the absence of an applied bias voltage due to the Josephson effect. We shall illustrate this effect starting from the non-interacting situation within the single level Anderson model introduced above. In this case the advanced-retarded Green function of the coupled dot can be expressed as

(11) |

where and are the dimensionless BCS green functions of the uncoupled leads. For simplicity we focus below on the case where both leads are of the same material for which .

The spectral density associated with this model exhibits bound states within the superconducting gap (i.e ). Physically, the ABSs arise from virtual multiple Andreev reflection processes at the interface between the dot and each of the leads. In such processes and for energies inside the gap, the electrons (holes) incident towards the leads are reflected back as holes (electrons), as illustrated in Fig. 2. The condition for the appearance of ABSs is that the accumulated phase in the closed trajectory be a multiple of , which is equivalent to satisfying the equation

(12) | |||||

where we have used that the dimensionless BCS Green functions become real for energies inside the superconducting gap. It can be shown Rozhkov and Arovas (2000) that this equation has two real roots inside the gap , i.e. symmetrically located with respect to the Fermi energy.

An essential property of the ABSs is that they correspond to current-carrying states. In fact, due to its dependence on the superconducting phase difference they have associated a Josephson current . The total Josephson current is obtained by adding the contribution of all states with finite occupation. Thus, at zero temperature only the lower ABS contributes and there is an additional contribution from the continuous spectrum which we discuss below.

The ABS equation (12) becomes particularly simple for an electron-hole and left-right symmetric case (i.e. and ) when it can be reduced to

(13) |

which for has the simple solutions , where is a reduced gap parameter which for is given by . The ABSs for this case tend to the ones of a perfectly transmitting one channel superconducting contact , the main qualitative difference being the reduced amplitude of their dispersion detaching them from the gap edges at . This is illustrated in Fig. 3(a).

In the absence of electron-hole symmetry (i.e. ) a finite internal gap between the upper and lower ABSs appears as in the case of a non-perfect transmitting one channel contact. When the ABSs for this case are given by , where is the normal transmission at the Fermi energy. Outside this limiting case the ABSs exhibit both the internal gap and the detachment from the continuum states at , as it is illustrated in Fig. 3(b).

An interesting issue to comment is that the contribution from the states in the continuous spectrum becomes negligible when . In this limit the zero temperature current-phase relation (CPR) is simply given by

(14) |

which is the CPR of a one channel contact with transmission . For finite there is a contribution from the continuum states. The expression for the total Josephson current in the non-interacting case can be derived from Eq. (10), which yields the following compact form

(15) |

The current density in Eq. (15) contains both the contribution from the ABSs (region ) and from the continuous spectrum . The behavior of the current density as a function of is depicted in Fig. 4. As can be observed the contribution from the continuous spectrum has the opposite sign compared to the one arising from the ABS. The inset shows that this contribution becomes negligible in the limit and reaches a maximum for .

In the rest of this section we shall discuss the different theoretical approaches to include the effect of interactions in the dc Josephson effect through single level QD models. We also include a subsection on experimental results.

### iii.1 Cotunneling approach

From a theoretical point of view the simplest approach to account for the effect of interactions in the Josephson current is to perform a perturbative expansion to the lowest non-zero order in the tunnel Hamiltonian. This so-called cotunneling approach was first used by Glazman and Matveev Glazman and Matveev (1989), who predicted the onset of the -junction behavior by this method. More precisely they obtained for the limit

(16) |

where changes its value from 2 () to -1 (), thus describing the transition to the -phase, the function having the form

(17) |

This approximation is clearly not valid for describing the Kondo regime () which requires non-perturbative approaches like the ones discussed in following subsections.

### iii.2 Mean field and variational methods

Another simple approximate methods to deal with interactions are those of a mean field type like the Hartree-Fock approximation (HFA) or the slave-boson mean field (SBMF). In spite of their simplicity these approximations are able to capture important qualitative features due to interactions in certain limits.

We start by analyzing the HFA. In the context of magnetic impurities in
superconductors this method was first applied by Shiba Shiba (1973),
while for the analysis of the Josephson effect it was first considered in Ref. Rozhkov and Arovas (1999)
and further analyzed in Vecino *et al.* (2003). Within this approximation
electrons with a given spin “feel” a static potential due to the average
occupation of the opposite spin, which corresponds to a simple constant
self-energy and
. In principle within the same level
of approximation there appears a non-diagonal self-energy taking into
account the induced pairing within the dot due to proximity effect,
which can be written as .
The effect of this non-diagonal contribution,
which was not included in Ref. Rozhkov and Arovas (1999), was analyzed in Ref. Vecino *et al.* (2003).
The determination of the self-energy in the HFA requires a self-consistent
calculation by using Eqs. (II)
which cannot be performed analytically in general.

The most significant result within the HFA is the appearance of a broken symmetry state in which the dot acquires a finite magnetic moment (i.e. ) for certain ranges of parameters. In this respect one should be cautious in principle as the HFA is known to predict also broken symmetry states for the same model with normal leads Anderson (1961), which are known to be spurious. However, for the superconducting case ground states with a finite magnetization do exist for certain parameter range as commented in the introduction. As it is shown below the HFA gives a rather good estimate of the magnetic ground state energy in the regions where it is the most stable phase.

The general properties of the ground state within the HFA are most conveniently displayed by a phase-diagram like the one in Fig. 5. In this diagram the notation ”0”, ””, ”” and ”” corresponds to the different ground state symmetries. Thus, ”0” corresponds to the non-magnetic case for all values of (the absolute energy minimum being located at ), while the ”” denotes that the magnetic solution is the most stable for all values (with the absolute minimum at ). On the other hand, ”” and ”” refer to intermediate situations with mixed magnetic and non-magnetic solutions as a function of , the name indicating whether the absolute energy minimum corresponds to a non-magnetic or a magnetic solution. From Fig. 5 the broken symmetry ground states are predicted to appear around the line, which corresponds to the half-filled case for sufficiently large , i.e. . It is worth noticing that for normal leads this region corresponds to the deep Kondo regime, which anticipates an interesting interplay between both effects in the superconducting case beyond the HFA.

Further insight into the HFA solution can be provided by a ”toy” model
introduced in Ref. Vecino *et al.* (2003) (A similar model was analyzed in Ref. Benjamin *et al.* (2007)).
In this simplified model the finite
magnetization which appears in the HFA is simulated by means of an exchange
field parameter, , corresponding to the splitting of the diagonal dot levels for each
spin, i.e. .
The analysis of the Andreev states within this toy model is similar to the
one given at the beginning of this section for the non-interacting case, and becomes particularly simple in the limit
in which they adopt the analytical expression

where , and and is given by

This expression clearly shows that the effect of the exchange field is to break the spin degeneracy producing an splitting of the ABSs. Consequently for one could in principle observe up to four ABSs in the spectral density. The evolution of these states with increasing is shown in Fig. 6 together with the corresponding Josephson current. While for the splitting is small and states corresponding to different spin orientation do not cross, for increasing the upper and lower states closer to the Fermi energy begin to cross yielding a current-phase relation of or character. Eventually for sufficiently large these two states completely interchange position with a complete reversal of the sign of the current (-phase). It should be mentioned that although in this toy model the spin degeneracy is artificially broken, it nevertheless qualitatively simulates the behavior of the exact solution in the magnetic phase.

Another simple approach of a mean field type is provided by the slave boson mean field approximation (SBMFA). This method was introduced by Coleman for the normal Anderson model Coleman (1984); Hewson (1993). It is based in the introduction of auxiliary boson fields which act as projectors onto the empty impurity state. At the same time fermion creation and annihilation operators are introduced for describing the singly occupied states. In order to get rid of the doubly occupied states in the these operators should satisfy the completeness relation

(19) |

In terms of these operators the terms and become

(20) |

So far this transformation is exact in the limit. Specific diagrammatic methods to obtain the impurity self-energy in this slave boson formulation have been developed Bickers (1987). Within this formulation the simplest solution is provided by the mean field approximation in which the boson operator is treated as a c-number. In fact the dot Hamiltonian reduces in this case to

where is a Lagrange multiplier associated to the constraint (19). The problem becomes equivalent to a non-interacting impurity model with renormalized parameters and . Self-consistency is achieved by minimizing the system free energy.

Strictly speaking the mean field approximation is only valid in the limit
where is the level degeneracy of the Anderson model (for the single level case N=2).
However, the mean field approximation yields a reasonably good description of quantities like
the Kondo temperature in the normal case Hewson (1993).
When applied to the Anderson model with superconducting electrodes the SBMFA is only valid in the
regime as it is not able to describe the transition to the -phase when
. In the regime the ABSs as described by the SBMFA corresponds
to the non-interacting case with renormalized parameters and .
In this way the ABSs within the SBMFA in this regime would be given by with .
In principle, the self-consistent effective parameters in the superconducting state can differ from those in the
normal state. However, in the limit in which the approximation is supposed
to work this difference can be neglected. The SBMFA in the limit
has only been applied for the case of superconducting leads in a few references:
Avishai et al. Avishai *et al.* (2003) for analyzing the dc current with an applied bias voltage,
and in Ref. Yeyati *et al.* (2003) for studying the dynamics of Andreev states in the Kondo regime.
Both works correspond to the non-equilibrium situation which will be discussed in Section V.

For a proper description of the phase-diagram within a mean-field slave boson approach a finite- version of the method, like the one introduced by Kotliar and Ruckenstein Kotliar and Ruckenstein (1986), is necessary. Within this method the number of auxiliary boson fields is extended up to four, denoted by and , which project into the empty, singly occupied (with either spin orientation) and doubly occupied dot states respectively. These operators must satisfy the constraints

(22) |

For recovering the non-interacting limit it is necessary to introduce also an auxiliary operator , in terms of which the model Hamiltonian becomes

(23) | |||||

where and are the Lagrange multipliers associated with the constraints (22). Again, in a mean field approximation, the auxiliary fields are treated as c-numbers to be determined self-consistently.

The type of phase-diagram that is obtained within the finite-U SBMFA will be analyzed in subsection III.4, for the zero band-width limit which allows a comparison with exact diagonalizations. As it is shown in that subsection the finite-U SBMFA tends to underestimate the stability of the -phase in contrast with the HFA, which typically overestimates it.

Another relatively simple approach is provided by the use a variational wave-function. This approach was used by Rozhkov and Arovas Rozhkov and Arovas (2000) extending previous works Varma and Yafet (1976) in which variational wave-functions were proposed for analyzing the normal Kondo problem. In their work Rozhkov and Arovas propose different many-body variational states in the limit for the singlet and the doublet states, looking for the their relative stability. They find a transition between both ground states for and also predict the appearance of the intermediate phases and in addition to the pure and ones.

### iii.3 Diagrammatic approaches

#### Perturbation theory in the Coulomb interaction

A natural extension over the HFA is provided by applying diagrammatic perturbation theory in
the Coulomb parameter . Already at the level of second order one can obtain an approximation
for the self-energy which is able to capture part of the interplay between Kondo effect and
pairing interactions. This approximation has been applied both for the
S-QD-S case in equilibrium Matsumoto (2001); Vecino *et al.* (2003), as well as for the N-QD-S case Cuevas *et al.* (2001); Yamada *et al.* (2010).

The diagrams contributing to the second-order self-energy in the superconducting Anderson model are depicted in Fig. 7. The first diagrams (denoted as 11(a)) is equivalent to the one appearing in the normal case, describing interaction of an electron with an electron-hole pair with opposite spin. The other diagrams include anomalous superconducting propagators and are therefore characteristic of the superconducting state. The presence of these propagators gives several effects: the appearance of non-diagonal elements of the self-energy in Nambu space, and the presence of diagrams like 11(b) and 21(b) in Fig. 7 which corresponds to a double-exchange process. Finally, diagram 21(a) describes the interaction of a Cooper pair with fluctuations in the pairing amplitude within the dot.

Formally, these diagrams can be computed from the full Green functions of the non-interacting
case by means of the expressions Vecino *et al.* (2003)

The evaluation of these expressions for a general range of parameters requires a significant numerical effort. An efficient algorithm can be implemented to evaluate these expressions based on Fast Fourier transformations, as discussed in Schweitzer and Czycholl (1990).

In the limit Kondo correlations dominate over
pairing ones. The results of the second-order self-energy approach for the half-filled case
capture the main features of the onset of Kondo correlations in the spectral density when
. This is illustrated in Fig. 8 taken
from Ref. Vecino *et al.* (2003), which shows its evolution for increasing .
The spectral density is similar to the one found in the normal state except for the
superimposed features inside the superconducting gap.
The overall shape
evolves from the single Lorentzian broad resonance for to the three peaked
structure characteristic of the Kondo regime when . In this regime the
width of the central Kondo peak is set by the scale , which in the present
approximation is given by , where

(24) |

thus coinciding with the perturbative result in the normal state (see Ref. Yamada and Yoshida (1975)).
Although this perturbative approach fails to yield the exponential behavior of for
large , it provides a reliable description of the spectral density for moderate
values of this parameter Ferrer *et al.* (1986).

The second-order self-energy allows also to analyze the renormalization of the ABSs due to the presence of Coulomb interactions. For values of the renormalized ABSs maintain approximately the behavior of the non-interacting case but with a narrower dispersion set by , where .

On the other hand, when a transition to the -phase
is expected. Within the second-order self-energy approach the transition can be identified
by allowing for a breaking of the spin-symmetry in the initial non-interacting problem and
searching for self-consistency. Rather than imposing the consistency condition of the HFA, i.e.
in Ref. Vecino *et al.* (2003) it was
imposed that the effective dot level for each spin-orientation be determined by
the charge-consistency condition, i.e. , where
is the dot charge corresponding to the broken-symmetry non-interacting Hamiltonian.
Such a procedure was shown to eliminate the unstable behavior of perturbation theory
when developed from the HFA Yeyati *et al.* (1993).

#### NCA approximation

Within the diagrammatic approximations one can include the so-called non-crossing approximation
(NCA). In this case an infinite order resumation of the perturbation theory is performed starting
from the slave boson representation of the Anderson Hamiltonian
Bickers (1987). To the lowest order in (where is the spin degeneracy) the
family of diagrams in this resumation is represented in Fig. 9. The
dashed lines correspond to the fermion propagators and the wavy lines to the slave bosons.
In the normal case the NCA include only the first two diagrams in the Dyson equation for
the fermion and boson propagators. The extension to the superconducting case was proposed
in Ref. Clerk *et al.* (2000) and corresponds to including the anomalous propagators for describing
multiple Andreev reflection processes (last two diagrams in the fermion and boson self-energies
represented in Fig. 9).

In the normal case the NCA theory has been shown to yield reliable results for temperatures down below Bickers (1987); Hewson (1993) in spite of certain pathologies like its failure to fulfill the Friedel sum-rule. The self-consistent extension for the superconducting case by Clerk et. al. Clerk and Ambegaokar (2000) is also formally exact to order and it is thus expected to yield reasonable results even in the presence of MAR processes.

Further analysis of the NCA applied to the S-QD-S system was provided in Ref. Sellier *et al.* (2005).
Their results for the Josephson current and LDOS in the superconducting gap region are summarized in Fig. 10.
The fact that the calculations are performed for temperatures which are a quite large
fraction of yields very broad resonances for the subgap states. As can be observed in
Fig. 10, only one broad resonance can be clearly resolved within the gap.
These results are in contrast to what is obtained using exact numerical methods as
will be discussed in Sect. III.4.2.
In this approximation the transition to the -phase appears as a
smooth crossover which can be associated to the crossing of this resonance through the Fermi energy.

#### Real time diagrammatic approach

Another technique which has been applied to the study of quantum dots coupled
to superconducting leads is the real time diagrammatic approach first introduced
by König et al. König *et al.* (1996) for the normal Anderson model.
The main idea of this technique is to integrate out the fermionic degrees of freedom
of the electrodes leading to a reduced description of the density matrix
projected on the Hilbert space of the isolated dot states. In the superconducting case
this reduced density matrix also depends on the
number of Cooper pairs in the leads relative to some chosen reference.
The aim of the technique is to determine the
time evolution of this reduced density matrix in the Keldysh contour thus allowing
to consider both equilibrium and non-equilibrium situations.
In Ref. Governale *et al.* (2008) the method has been applied to the equilibrium S-QD-S
case obtaining results in agreement with those of Ref. Glazman and Matveev (1989) in the cotunneling
limit. The method has been mainly applied to analyze the properties of quantum dots
connected to both normal and superconducting leads in multiterminal configurations
out of equilibrium, an issue which will be commented in Sect. VI.

### iii.4 Diagonalization by numerical methods

Within this subsection we will review methods which attempt a direct diagonalization of the superconducting Anderson model, either by truncating the initial Hilbert space using physical arguments valid for certain parameter region or by using the Numerical Renormalization group (NRG) method.

#### Exact diagonalization for the large limit

An exact diagonalization of the model is possible in the limit .
In this limiting case the Hilbert space of the problem is automatically reduced to
states spanned by the different electronic configuration of the dot levels.
The effect of the superconducting leads appears
as a pairing term between the electrons within the dot. The effective Hamiltionian for the truncated
Hilbert space becomes Vecino *et al.* (2003); Tanaka
*et al.* (2007b); T. Meng and Simon (2009)

(25) | |||||

The eigenvalues of this reduced Hamiltonian can be determined straightforwardly by noting the decoupling of subspaces with even and odd number of electrons. The ground state for the even case (corresponding to total spin ) is a linear combination of the empty and doubly occupied dot state with an energy

(26) |

On the other hand, the odd number subspace simply corresponds to a single uncoupled spin with energy . The transition to the magnetic ground state thus occurs for . In the simpler electron-hole symmetric case this condition reduces to and thus the full state appears for . This simple model already gives a rough qualitative account of the quantum phase transition.

A further step in the idea truncating the Hilbert space is performed in the so called zero
band-width limit Affleck *et al.* (2000); Vecino *et al.* (2003); Bergeret *et al.* (2007). In this approximation the
superconducting leads are represented by a single
localized level (which formally corresponds to the limit of vanishing width of the leads
spectral density). This approximation is justified when the superconducting gap is large compared
to the other energy scales in the problem, and thus can be considered as a refinement with respect
to the previous approach. The corresponding Hamiltonian can be written as ,
with

(27) |

where denotes the left-right sites describing the leads in this approximation.

Although the total number of particles is not a good quantum number, their parity is conserved as in the previous case. This allows to reduce the initial 64 states in the Hilbert space to a subspace of 20 states for even parity with a total spin -component and 15 states for odd parity with . These values of are the ones corresponding to the ground state in each subspace with total spin and . In addition to providing a qualitative description of the phase diagram of the full model, this simplified calculation can furthermore be useful as a test for comparing different approximation methods.

The phase diagram obtained within this approximation was discussed in Refs. Vecino *et al.* (2003); Bergeret *et al.* (2007)
and is shown in Fig. 11 for , with .
As can be observed, the overall diagram is very similar to the one shown before for the
HFA (Fig. 5) exhibiting the four phases and in the
same sequence. For a more direct comparison it is necessary to perform the HFA of the ZBW model,
which as an exactly solvable model also provides a stringent test of the approximation. The
lower broken line in Fig. 11 indicates the boundary of the -phase
within the HFA. It can be noticed that the HFA overestimates the stability of this magnetic phase.
On the other hand, it is also possible to test the finite-U SBMF approximation in this ZBW model. The
corresponding boundary for the -phase is indicated by the upper broken line in Fig.
11. In opposition to the HFA this approximation overestimates the stability
of the phase, the exact boundary therefore lying in between the two different
mean-field approximations. It is interesting to point out that a similar difference between
both approximations is also found for the full model (see inset in Fig. 11).
One would then expect that the exact boundary for the full model lays in between these two.

The ZBW model can also be used to test approximate methods beyond the mean field ones. In Ref. Vecino *et al.* (2003)
this was done for the second-order self-energy approximation. As it is illustrated in Fig.
12 for the symmetric case
the second-order self-energy approach matches quite well the exact ground state energy
for values up to . It is interesting to note that for larger values
even when the HFA already yields a full -state, the second-order approximation
predicts a mixed ground state in agreement with the exact solution.

#### Numerical Renormalization Group (NRG)

The NRG method is based on the ideas of Wilson on logarithmic
discretization for magnetic impurity problems Wilson (1975) and was first applied to the
Anderson model by Krishna-murthy et al. Krishna-murthy
*et al.* (1980). The idea behind the method is
to discretize the energy levels in the leads on a logarithmic grid of energies
(with and ) with exponentially
high resolution on the low-energy excitations. This discretization allows then to map
the impurity model into a linear ”tight-binding” chain with hopping matrix elements
decaying as with increasing site index . The sequence of Hamiltonians
which is constructed by adding a new site in the chain is then diagonalized iteratively.
As the number of states grows exponentially an adequate truncation scheme is required.

The NRG scheme has been first generalized to the case of an Anderson impurity in a superconducting
host by Yoshioka and Ohashi Yoshioka and Ohashi (2000) and implemented by several authors to analyze
the S-QD-S model with a finite phase difference Choi *et al.* (2004); Oguri *et al.* (2004); Tanaka
*et al.* (2007b); Lim and Choi (2008).
For the left-right symmetric case (i.e. and )
the sequence of Hamiltonians can be written as Choi *et al.* (2004)

(28) | |||||

where the initial Hamiltonian is given by

(29) | |||||

The fermion operators correspond to an effective tight-binding chain resulting from the logarithmic discretization and the canonical transformation into the even-odd linear combination of original left-right states in the leads and

(30) | |||||

with and being an energy cut-off in the leads spectral density. The original Hamiltonian is recovered in the limit .

The NRG method was applied to analyzing the Josephson current
in a S-QD-S system in Ref. Choi *et al.* (2004), this work confirming the predicted quantum phase
transition at for the electron-hole symmetric case.
It should be mentioned that more recent calculations Karrasch *et al.* (2008) using NRG obtain
Josephson currents which are approximately a factor 2 larger than the ones of
Choi *et al.* (2004).
On the other hand, Oguri et al. Oguri *et al.* (2004) used NRG to analyze this model in the case of
. In this case the model can be exactly
mapped into a single channel model consisting on the right lead coupled to the Anderson
impurity with a local pairing ,
thus allowing a simpler implementation of the NRG algorithm. Further work for the case although with by Tanaka et al. Tanaka
*et al.* (2007b)
confirmed the presence of intermediate phases even in the left-right asymmetric
case. A characteristic phase diagram obtained for this case is shown in Fig.
13.

In addition to the ground state properties, NRG methods have been applied in an attempt to
clarify the structure of the subgap ABSs. In Ref. Lim and Choi (2008) the spectral density inside
the gap obtained from the NRG algorithm was analyzed, showing that a pair of ABSs
located symmetrically respect to the Fermi energy is present in the
case. This is in contrast to the NCA results discussed previously (shown in Fig.
10) where a single broad resonance appears. Similar conclusions are
obtained in Ref. Bauer *et al.* (2007) although for the single lead case and for finite .
A word of caution should be said regarding the analysis of the ABSs in this last
work in which the relation is assumed
in their Eq. (8) for the states inside the gap. This relation would not be strictly valid for
the doublet ground state when choosing a given spin orientation. In this case the quasi-particle
excitation energies would become spin dependent and the electron-hole symmetry would be broken.
This would allow in principle to have up to 4 ABSs inside the gap as predicted both by the Hartree-Fock
approximation and in the exact limit. Of course,
in the phase the spin is not frozen but is fluctuating. In this sense the above relation
between the self-energy components would be recovered when averaging over the and
states.
We believe in any case that a more detailed analysis of the ABSs using the NRG method is still lacking.

### iii.5 Functional renormalization group

The functional renormalization group (fRG) method is based on the application of an RG analysis to the diagrammatic expansions in terms of electronic Green functions. This is an approximate method whose accuracy depends on the initial diagrams used in the evaluation of the electron self-energies. The starting point is the introduction of an energy cut-off into the Matsubara non-interacting Green-functions

Using these propagators the particle vertex functions acquires a dependence.
The flow equations are determined differentiating these vertex functions with respect to
which are then solved iteratively for increasing . In Ref. Karrasch *et al.* (2008) the
method has been applied to the S-QD-S system employing a truncation scheme which keeps
only diagrams corresponding the the static Hartree-Fock approximation.
The -dependent Green function used in Ref. Karrasch *et al.* (2008) was of the form,

(31) |

where and , and being the dimensionless BCS Matsubara Green functions of the uncoupled leads. Within this approximation the flow equations lead to energy-independent self-energies, corresponding to an effective non-interacting model with renormalized parameters. It is important to notice that this approximation exhibits also the limitation already pointed out in the previous section as it imposes electron-hole symmetry which is not satisfied in the magnetic phase. Nevertheless the approximation allows to identify a transition to a phase with inversion of the Josephson current which is driven by an ”overscreening” of the induced pairing determined by . Fig. 14 shows the comparison of fRG results with those obtained with the NRG method. The agreement is rather good for large but it becomes poorer in the -phase. We believe that the agreement could be improved allowing for a broken symmetry state within the same fRG approach.

### iii.6 Quantum Monte-Carlo

The Josephson current in the S-QD-S system has also been analyzed using Quantum Monte Carlo (QMC) simulations
by Siano and Egger Siano and Egger (2004). The method used was the Hirsh-Fisher algorithm adapted to
this particular problem. They consider the deep Kondo regime and
and show that the results for the Josephson current
exhibit a universal dependence with provided that . They identify
the transition between the different phases at and for
, and respectively Siano and Egger (2005a). Being a finite temperature calculation
the resulting current-phase relations do not exhibit sharp discontinuities in the intermediate
phases. This smooth behavior was criticized in Ref. Choi
*et al.* (2005a) pointing out that the QMC
results did not match the NRG ones of Ref. Choi *et al.* (2004) at finite temperatures, which was attributed
by Siano and Egger in their reply Siano and Egger (2005b) to a possible limited accuracy
of the NRG calculation of Ref. Choi *et al.* (2004). More recent
NRG calculations of Ref. Karrasch *et al.* (2008) give a good agreement with QMC results at finite temperature
for values between and , whereas QMC underestimates the current in the range
. It is claimed in Ref. Karrasch *et al.* (2008) that the origin of the discrepancy lies in the
fact that the first excited state in this phase range is smaller or of the order of the
temperature values used in the calculations of Siano and Egger (2004).

The QMC method has more recently been applied to analyze the spectral properties
of this model in Ref. Luitz and Assaad (2010). The authors employ the so-called weak-coupling continuous-time
version of the method which is based on a perturbative expansion around the limit. They show
that the results for the spectral densities are in qualitative good agreement with the ones obtained
in the zero band-width approximation introduced in Ref. Vecino *et al.* (2003), which was discussed
before in this section.

### iii.7 Experimental results

Several physical realizations of the S-QD-S system have been obtained in the last few
years by means of contacting carbon nanotubes (CNT), C60 molecules or semiconducting nanowires
with superconducting electrodes (for a recent review see Ref. Franceschi
*et al.* (2010)). In view of
the existence of this review on the experiments in this section we give only a brief summary of the
main findings and its relation to the theoretical work.

CNTs have provided so far the most promising setups for a direct test of the theoretical predictions
concerning the Josephson effect through a QD. The first experiments detecting a supercurrent
through a CNT-QD strongly coupled to the leads (i.e. ) were performed
by Jarillo-Herrero et al Jarillo-Herrero
*et al.* (2006). These experiments were basically performed in the resonant-tunneling
regime with a single-level spacing . The results indicated a
strong correlation between the critical current and the normal conductance . However, the
product deviated from a constant value due to the effect of the electromagnetic environment
suppressing the critical current more strongly in off-resonance conditions due to phase fluctuations.

*et al.*(2006) for analyzing the transition in CNT QDs (upper panel). The middle panel shows the reversal of the oscillatory pattern of the current as a function of the magnetic flux through the SQUID across the transition. The lower panel shows the correlation of the -phase region with the Kondo ridges in the normal state (indicated by the dashed line). Reprinted by permission from Macmillan Publishers Ltd: Nature Physics Cleuziou

*et al.*(2006), copyright (2006).

In a subsequent work by the same group the first experimental evidence of -junction behavior
in a S-QD-S system was obtained using semiconducting InAs nanowires with Al leads in a SQUID
configuration van Dam *et al.* (2006). However, the observed features in this experiment could not be
explained completely on the basis of a single-level model but rather a multilevel description
was necessary. In particular the authors showed that in this case the -junction behavior
is not necessarily linked to the parity of the number of the electrons in the dot (this will be further
discussed in Sect. VI).
On the other hand, -junction behavior have also been demonstrated in Ref. Cleuziou *et al.* (2006)
using CNT-QDs in a SQUID geometry. The corresponding experimental setup
is depicted in Fig. 15. In this SQUID configuration the transition to the
-phase is directly demonstrated in the measured current-phase relation as a function
of one of the applied gate voltages (see Fig. 15). Remarkably, this experiment
allows to correlate the appearance of the -junction behavior with the presence of Kondo
correlations in the normal state. The dotted lines in the lower panel of Fig. 15 indicate the
Kondo ridge appearing in the normal state (which is recovered by applying a magnetic field).

Evidence of a transition has also been found in CNT-QD systems by analyzing the
current-voltage characteristic Jorgensen *et al.* (2007). In this work it was shown that the evolution of
the zero-bias conductance can be correlated to the behavior of the critical current
when transversing the boundary. As it corresponds to a non-equilibrium situation
this analysis will be discussed later in Sect. V. A similar observation holds for
other experimental works Vecino *et al.* (2004); Grove-Rasmussen
*et al.* (2007); Eichler *et al.* (2009)
which will be discussed in that section.

Finally, it is worth pointing out recent experimental developments which have allowed to directly
measure the Andreev bound states spectrum of a CNT-QD coupled to superconducting leads
in a SQUID configuration Pillet *et al.* (2010). In this experiment a weakly coupled lead was deposited at
the center of the CNT to allow for tunneling spectroscopy measurements. In this way
both the phase and gate voltage
variation of the Andreev bound states were determined. The results could be fitted satisfactorily
using a simplified phenomenological model corresponding to the superconducting Anderson model
within a mean field approximation, in which an exchange field is included to represent the magnetic
ground state. An example of the comparison between theory and experiment
is given in Fig. 16. This analysis showed that a double dot model was in
general necessary to fit the experimental results. Recent experimental work on graphene
QDs coupled to SC leads Dirks *et al.* (2011) provided also evidence on the crossing of ABs as a function of the gate
potential which is consistent with a magnetic ground state.

## Iv Quantum dots with normal and superconducting leads

A single-level quantum dot coupled to a normal and a superconducting lead provides a basic model system to study electron transport in the presence of Coulomb and pairing interactions. Compared to the S-QD-S situation, this case has a simpler response in non-equilibrium conditions due to the absence of the ac Josephson effect. For this reason this system has been widely analyzed theoretically.

As in any N-S junction the low bias transport properties are dominated by Andreev processes. This mechanism is in general highly modified by resonant tunneling through the localized levels in the dot. In addition, charging effects can strongly suppress the Andreev reflection in certain parameter ranges. Furthermore there is also an interesting interplay between Kondo and pairing correlations as in the case of the S-QD-S system discussed above.

As an illustration of the general formalism we derive here the linear transport properties of the non-interacting model. One can straightforwardly write the dot retarded Green function for this case from expression (11) by setting . Then, from the expression of the current in terms of Keldysh Green functions and using the corresponding Dyson equation one can write

(32) | |||||

where are the Keldysh Green functions of the uncoupled normal lead. By further using