Josephson and Andreev transport through quantum dots

Josephson and Andreev transport through quantum dots


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.

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).

Figure 1: Schematic representation of a single spin-degenerate level QD connected to normal (left panel) and superconducting (right panel) leads with a large charging energy . In the normal case the local density of states (LDOS) in the dot exhibits the typical form corresponding to the Kondo regime with a narrow resonance at the Fermi energy and a a broad resonance (of width ) below it. In the superconducting case the Kondo resonance (assumed to be narrower than the superconducting gap) disappears due to the competition with the pairing correlations in the leads. Courtesy of C. Schönenberger.

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


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


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


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)


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


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


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


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


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


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


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


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

Figure 2: Schematic representation of the physical mechanism responsible for the formation of ABs in a generic nanostructure coupled to superconducting leads. Reprinted by permission from Macmillan Publishers Ltd: Nature Physics Pillet et al. (2010), copyright (2010).

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


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).

Figure 3: Andreev bound states for a non-interacting S-QD-S with different values: 1.0, 2.0 and 4.0. Left panel corresponds to the e-h symmetric case () and the right panel to a case with .

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).

Figure 4: Current spectral density for the non-interacting S-QD-S system with and for increasing values of 0.5 (red), 4 (green) and 16 (blue). A small broadening has been introduced to help visualizing the ABs contribution, which has been truncated for the sake of clarity. The inset shows the relative contribution of the continuous spectrum compared to the total value .

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


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


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


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


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.

Figure 5: Phase diagram in the plane for obtained using the HFA Rozhkov and Arovas (1999)). The phases are classified into , , and as explained in the text. Reprinted figure with permission from A.V. Rozhkov and D. Arovas, Physical Review Letters 82, 2788, 1999 Rozhkov and Arovas (1999). Copyright (1999) by the American Physical Society.

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.

Figure 6: Evolution of the ABSs and Josephson current in the “toy” model of Ref. Vecino et al. (2003) for and increasing parameter: 0.25 (upper panels), 0.75 (middle panels) and 1.5 (lower panels). Reprinted figure with permission from E. Vecino et al., Physical Review B 68, 035105, 2003 Vecino et al. (2003). Copyright (2003) by the American Physical Society.

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


In terms of these operators the terms and become


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


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


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).

Figure 7: Feynmann diagrams for the second-order self-energy in the Anderson model with superconducting electrodes. Reprinted figure with permission from E. Vecino et al., Physical Review B 68, 035105, 2003 Vecino et al. (2003). Copyright (2003) by the American Physical Society.

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).

Figure 8: Evolution of the DOS in the S-QD-S system in equilibrium within the second-order self-energy approximation for and electron-hole symmetric case with . The parameter takes the values 2.5 (a), 5 (b) and 10 (c). Reprinted figure with permission from E. Vecino et al., Physical Review B 68, 035105, 2003 Vecino et al. (2003). Copyright (2003) by the American Physical Society.

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


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).

Figure 9: Fermion (top) and boson (bottom) self-energy diagrams in the NCA approximation extended to the superconducting case. Reprinted figure with permission from G. Sellier et al., Physical Review B 72, 174502, 2005 Sellier et al. (2005). Copyright (2005) by the American Physical Society.

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.

Figure 10: Josephson current and subgap LDOS for the equilibrium S-QD-S model within the NCA for three values of and three different temperatures. Reprinted figure with permission from G. Sellier et al., Physical Review B 72, 174502, 2005 Sellier et al. (2005). Copyright (2005) by the American Physical Society.

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)


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


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


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.

Figure 11: Phase diagram of the S-QD-S system in the ZBW model for the leads obtained by exact diagonalization for , taken as unit of energy. The dashed lines correspond to the boundary between the 0 and regions within the HFA (lower line) and the finite-U SBMFA (upper line). The inset show the corresponding results for these boundaries within the full model. Reprinted figure with permission from F.S. Bergeret et al., Physical Review B 76, 174510, 2007 Bergeret et al. (2007). Copyright (2007) by the American Physical Society.

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.

Figure 12: Ground state energy for the and states of the S-QD-S system in the ZBW model for the leads. Full lines correspond to the exact results, dotted lines to the HFA and dashed lines to the second-order self-energy approximation. Reprinted figure with permission from E. Vecino et al., Physical Review B 68, 035105, 2003 Vecino et al. (2003). Copyright (2003) by the American Physical Society.

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)


where the initial Hamiltonian is given by


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


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

Figure 13: Phase diagram of the S-QD-S system in the plane for fixed and different values of obtained using the NRG method. Reprinted figure with permission from Y. Tanaka et al., New Journal of Physics 9, 115, 2007 Tanaka et al. (2007b). Copyright (2007) by IOP Publishing Ltd.

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

Figure 14: Current phase-relations for the S-QD-S system obtained using the fRG approach truncated at the HF level for different values of and . For comparison the results obtained using the NRG method are also plotted (indicated by the filled dots). Reprinted figure with permission from C. Karrash et al., Physical Review B 77, 024517, 2008 Karrasch et al. (2008). Copyright (2008) by the American Physical Society.

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,


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.

Figure 15: Experimental setup used in Ref. Cleuziou 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.

Figure 16: Color scale plots of the local density of states in a CNT QD coupled to SC leads as reported in Ref. Pillet et al. (2010) as a function of gate voltage and phase difference. The figure also shows the results obtained from a phenomenological model, as discussed in the text. Reprinted by permission from Macmillan Publishers Ltd: Nature Physics Pillet et al. (2010), copyright (2010).

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


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