Many-body perturbation theory for the superconducting quantum dot: Fundamental role of the magnetic field
MMany-body perturbation theory for the superconducting quantum dot: Fundamentalrole of the magnetic field
V´aclav Janiˇs ∗ and Jiawei Yan Institute of Physics, Academy of Sciences of the Czech Republic,Na Slovance 2, CZ-18221 Praha 8, Czech Republic (Dated: February 26, 2021)We develop the general many-body perturbation theory for a superconducting quantum dot rep-resented by a single-impurity Anderson model attached to superconducting leads. We build ourapproach on a thermodynamically consistent mean-field approximation with a two-particle self-consistency of the parquet type. The two-particle self-consistency leading to a screening of the bareinteraction proves substantial for suppressing the spurious transitions of the Hartree-Fock solution.We demonstrate that the magnetic field plays the fundamental role in the extension of the pertur-bation theory beyond the weakly correlated 0-phase. It controls the critical behavior of the 0 − π quantum transition, lifts the degeneracy in the π -phase, where the limits to zero temperature andzero magnetic field do not commute. The response to the magnetic field is quite different in 0- and π -phases. While the magnetic susceptibility vanishes in the 0-phase it becomes of the Curie typeand diverges in the π -phase at zero temperature. PACS numbers: 72.10.Fk73.63.Rt,74.40.Kb,74.81.-g
I. INTRODUCTION
Nanostructures with well separated localized energylevels are objects that can be isolated in regions of a fewnanometers or microns. They can be experimentally re-alized by magnetic impurities on metallic surfaces [1–6],semiconducting quantum dots [7] nanowires [8, 9], car-bon nanotubes [10–28] or single C molecules [29]. Theyare ideal systems for studying elementary quantum me-chanical phenomena according to the substrates on whichthey are grown or in which they are embedded due toa detailed control of the relevant microscopic parame-ters. When the impurity atoms with unpaired correlatedelectrons are placed in metals one observes the Kondoeffect [4–6]. The correlated quantum nanostructres at-tached to superconductors represent tunable microscopicJosephson junctions [10, 12, 30]. The simultaneous pres-ence of strong electron correlations on semiconductingimpurities and proximity of superconductors allow us toobserve and analyze the interplay between the Kondo ef-fect and the formation of the Cooper pairs carrying theJosephson current through the semiconducting nanode-vices [8, 15–17, 21, 31–45].Strong Coulomb repulsion on quantum dots attachedto superconducting leads may cause a local quantum crit-ical point at which the lowest many-body eigenstates ofthe system cross and a spin-singlet ground state withthe positive supercurrent (0-phase) goes over to a spin-doublet state with a small negative supercurrent ( π -phase) [8, 15, 16, 25, 31–33, 37, 38, 41, 44, 46–49]. Thistransition is associated with crossing of the Andreevbound states (ABS) at the Fermi energy as has also beenobserved experimentally [9, 23, 26]. ∗ [email protected] A number of theoretical techniques have been usedto address the 0 − π transition and related propertiesof superconducting quantum dots. A very good quan-titative agreement with the experiments [26, 28, 44]can be obtained in a wide range of parameters usingheavy numerics such as numerical renormalization group(NRG) [26, 38, 39, 45, 46, 50–53] and quantum MonteCarlo (QMC) [28, 37, 44, 54, 55]. However, both NRGand QMC demand extensive time and computational re-sources and they do not disclose the microscopic origin ofthis quantum critical behavior. They are also unable todistinguish the physically different properties of the in-gap states in the 0- and π -phases. Alternatively, analyticapproaches have been used mostly based on perturba-tion expansions, either in the strength of the Coulombrepulsion [56–61] or around the atomic limit [62–65].The perturbation expansion in a small parameter can-not describe any collective behavior. A self-consistentsummation of infinite series must be included to inter-polate between weak and strong couplings needed todescribe the 0 − π transition. Summations via self-consistences are used both in the expansion in theCoulomb repulsion and around the atomic limit [66]. Per-turbation expansions around the atomic limit miss thestrong-coupling Kondo effect for narrow superconductinggaps. The expansion in the Coulomb repulsion is well de-fined only at zero temperature and in the weakly-coupledspin-symmetric state of the 0-phase. The standard wayto include critical behavior and phase transitions is to usea mean-field approximation with spin-polarized states asa starting point for the perturbation expansion [57]. Al-though the mean-field, Hartree-Fock approximation maygive reasonably good quantitative predictions for weakand moderate coupling [67] it is conceptually unaccept-able, since the real 0 − π transition is a consequence of aspurious critical transition to the magnetic state [66].There is a way to improve upon the improper start of a r X i v : . [ c ond - m a t . s t r- e l ] F e b the perturbation expansion in the interaction strength.One has to replace the weak-coupling mean-field approx-imation with an advanced one that is able to interpo-late consistently between the weak and strong couplings.It must be a self-consistent theory suppressing the spu-rious transition to the magnetic state and reproducingthe Kondo strong-coupling regime in the impurity mod-els. One of the authors proposed a mean-field theorywith a two-particle self-consistency that is free of any un-physical behavior and qualitatively correctly reproducesthe Kondo limit of the single-impurity Anderson model(SIAM) [68–72]. This mean-field approximation is a ther-modynamically consistent extension of the weak-couplingtheory to the whole range of the input parameters.It is the aim of this paper to apply the mean-fieldapproximation from references [68–72] on the Andersonimpurity attached to superconducting leads. The super-conducting leads induce a gap on the impurity with nostates at the Fermi energy. Instead, discrete in-gap statesemerge the position of which depends on the interactionstrength and the phase difference between the attachedsuperconducting electrodes. The theory developed forthe SIAM with non-zero density of states at the Fermienergy will be appropriately modified to offer a reliabledescription of the models with a gap. The extension ofthe mean-field approach to the singlet phase of the su-perconducting quantum dot seems straightforward, sinceit is the ground state in weak coupling. An extension tothe doublet phase with a degenerate ground state and noweak-coupling regime appears to be more elaborate.The many-body perturbation theory for low-energyexcitations can be used only with a unique, non-degenerate many-body ground state. It means that a de-generacy of the ground state in the doublet phase must belifted before we can apply the many-body Green functiontechnique and the diagrammatic expansion. The doubletground state is degenerate with respect to the spin reflec-tion. We must then use a small magnetic field on the im-purity to lift the degeneracy. We hence need to formulatethe mean-field approximation for the dot in an externalmagnetic field. The properties of the superconductingquantum dot with a Zeeman field were recently studiedexperimentally [73–76] and also theoretically [65, 77–83].The role of the Zeeman field in the perturbation theoryof the superconducting quantum dot is crucial. It al-lows us to circumvent the quantum critical point of the0 − π transition and to extend the many-body approachfrom weak to strong coupling regimes at zero tempera-ture. The limits to zero filed and zero temperature donot commute in the π -phase. Moreover, the response tothe magnetic field is crucial for distinguishing betweenthe 0-phase with bound singlet Cooper pairs and the π -phase with in-gap fermionic excitations carrying a localmagnetic moment. This will be demonstrated on the be-havior of the magnetic susceptibility. This feature hasnot yet been disclosed because a full consistent many-body theory of the superconducting quantum dot withthe Zeeman field and at arbitrary temperature is still missing.The layout of the paper is the following. We intro-duce the model and the Nambu formalism of the corre-lated impurity attached to superconducting leads in Sec.II. We introduce the basic ingredients of standard many-body perturbation theory in Sec. III. The core of ourthermodynamically consistent mean-field approximationwith a two-particle self-consistency is presented in Sec.IV. We apply our mean-field approximation to study thebehavior of the in-gap states and the 0 − π transition inSec. V. Explicit calculations are performed in the asymp-totic atomic limit of the infinite superconducting gap inSec. VI. Numerical results are presented in Sec. VII andSec. VIII brings concluding remarks. Less important andelucidating technical details are presented in AppendicesA-C. II. MODEL HAMILTONIAN AND ANDREEVBOUND STATES
Standardly, a single impurity is used to simulate thenanowire with separated energy levels connecting super-conducting leads in the experimental setup. The Hamil-tonian of the system consisting of a single impurity at-tached to BCS superconducting leads is H = H dot + (cid:88) s = R,L ( H slead + H sc ) , (1)where the impurity Hamiltonian is a single-level atomwith the level energy ± (cid:15) for single electron (hole) andCoulomb repulsion U in the Zeeman magnetic field h H dot = (cid:88) σ = ± ( (cid:15) − σh ) d † σ d σ + U d †↑ d ↑ d †↓ d ↓ , (2)where σ = ± H slead = (cid:88) k σ (cid:15) ( k ) c † s k σ c s k σ − ∆ s (cid:88) k ( e i Φ s c † s k ↑ c † s − k ↓ + H.c.) (3)with s = L, R denoting left, right lead. Finally, the hy-bridization term for the contacts reads H sc = − t s (cid:88) k σ ( c † s k σ d σ + H.c.) . (4)We use the Nambu spinor formalism to describethe Cooper pairs and the anomalous functions relatedwith the superconducting order parameters and break-ing charge conservation. The Nambu spinors in the su-perconducting leads are (cid:98) ϕ s k σ = (cid:32) c s k σ c † s ¯ k ¯ σ (cid:33) , (cid:98) ϕ † s k σ = (cid:16) c † s k σ c s ¯ k ¯ σ (cid:17) , (5)where we introduced ¯ k = − k and ¯ σ = − σ .Due to the hybridization the Cooper pairs can pene-trate onto the impurity giving rise to anomalous impurityGreen functions. Hence, we introduce the Nambu spinorsalso for the impurity (local) operators (cid:98) φ σ = (cid:18) d σ d † ¯ σ (cid:19) , (cid:98) φ † σ = (cid:0) d † σ d ¯ σ (cid:1) . (6)The individual degrees of freedom of the leads areunimportant for the impurity quantities and we integratethem out leaving only the impurity variables dynamical.The fundamental function after projecting the lead de-grees of freedom is the one-electron impurity Green func-tion measuring (imaginary) time fluctuations that in theNambu formalism is a 2 × (cid:98) G σ ( τ − τ (cid:48) )= − (cid:32) (cid:104) T (cid:2) d σ ( τ ) d † σ ( τ (cid:48) ) (cid:3) (cid:105) , (cid:104) T [ d σ ( τ ) d ¯ σ ( τ (cid:48) )] (cid:105)(cid:104) T (cid:104) d † ¯ σ ( τ ) d † σ ( τ (cid:48) ) (cid:105) (cid:105) , (cid:104) T (cid:104) d † ¯ σ ( τ ) d ¯ σ ( τ (cid:48) ) (cid:105) (cid:105) (cid:33) = (cid:18) G σ ( τ − τ (cid:48) ) , G σ ( τ − τ (cid:48) )¯ G σ ( τ − τ (cid:48) ) , ¯ G σ ( τ − τ (cid:48) ) (cid:19) (7)correlating appearance of electrons and holes with spe-cific spin on the impurity. We introduced normal parti-cle and hole G σ , ¯ G σ propagators that conserve spin andanomalous G σ and ¯ G σ Green functions that create andannihilate singlet Copper pairs in the spin-polarized so-lution.The electron and hole functions are connected bysymmetry relations ¯ G σ ( τ ) = − G ¯ σ ( − τ ) = − G ¯ σ ( τ ) ∗ ,¯ G σ ( τ ) = G ¯ σ ( − τ ) ∗ , and G σ ( τ ) = −G ¯ σ ( − τ ).The problem can be exactly solved for an impuritywithout the onsite interaction, U = 0. In this case the in-verse unperturbed propagator for the spin-polarized sit-uation can be represented in the Nambu formalism as amatrix. We use identical left and right hybridizations tosuperconductors, t L = t R = t without loss of general-ity. The asymmetric situation can be transformed to asymmetric one [84]. Due to energy conservation it is con-venient to use Fourier transform from (imaginary) timeto frequency (energy) where the Green function can an-alytically be continued to complex values. The matrix ofthe inverse Green function for a complex energy z reads (cid:98) G − σ ( z ) = (cid:18) z [1 + s ( z )] + σh − (cid:15) , ∆ cos(Φ / s ( z )∆ cos(Φ / s ( z ) , z [1 + σ ( z )] + σh + (cid:15) (cid:19) , (8)where s ( z ) = i Γ ζ sgn( (cid:61) z ) . (9)is the “hybridization self-energy” σ ( z ), that is, a dynam-ical renormalization of the impurity energy level due tothe hybridization to the superconducting leads. We ap-proximated the Green function in the leads by its value at the Fermi energy and denote Γ = 2 πt ρ being the ef-fective hybridization strength. We denoted Φ = Φ L − Φ R the difference between the phases of the attached super-conducting leads and ρ the density of states of the leadelectrons at the Fermi energy. To represent explicitlythe hybridization self-energy we introduced a new com-plex number ζ = ξ + iη derived from the complex energy z = x + iy by a quadratic equation ζ = z − ∆ . Therebythe following convention for the complex square root hasbeen used ξη = xy, sgn( ξ ) = sgn( x ) , sgn( η ) = sign( y ) . (10)The renormalized energy ζ along the real axis z = x ± i − ∆ , ∆) and imaginarywithin it ζ = sgn( x ) (cid:112) x − ∆ for | x | > ∆ ,ζ = ± i (cid:112) ∆ − x for | x | < ∆ . (11)Accordingly the hybridization self-energy is purely imag-inary outside the gap and real within it s ( x ± i
0) = ± i Γ sgn( x ) √ x − ∆ for | x | > ∆ ,s ( x ± i
0) = Γ √ ∆ − x for | x | < ∆ . (12)With the above definitions the unperturbed ( U = 0)impurity Green function is (cid:98) G (0) σ ( z ) = 1 D σ ( z ) × (cid:18) z [1 + s ( z )] + σh + (cid:15) , − c Φ ∆ s ( z ) − c Φ ∆ s ( z ) , z [1 + s ( z )] + σh − (cid:15) (cid:19) . (13)where we denoted c Φ = cos(Φ /
2) and introduced D σ ( z ) = [ z (1 + s ( z )) + σh ] − (cid:15) − c ∆ s ( z ) the determinant of the matrix of the inverse unperturbedimpurity Green function. It is decisive for the determi-nation of the gap states. This determinant is real withinthe gap and can goes through zero determining the gapstates that are simultaneously the Andreev states. Theyare four of them ± ω σ in the external magnetic field. Wedenote the two independent ω σ (1 + s σ ) = − σh ± (cid:113) (cid:15) + c ∆ s σ . (14)where we used Eq. (13) and denoted s σ = s ( ω σ ). III. PERTURBATION EXPANSION:DIAGRAMMATIC REPRESENTATION
The best way to represent the many-body perturba-tion expansion is to use a graphical, diagrammatic rep-resentation that can be introduced also in the Nambuformalism. We start with the diagrammatic represen-tation of the Nambu spinor of the impurity propagatorto which we assign solid lines decorated with arrows asfollows (cid:18) G σ ( iω n ) , G σ ( iω n )¯ G σ ( iω n ) , ¯ G σ ( iω n ) (cid:19) = (cid:18) G σ ( iω n ) , G σ ( iω n ) G ∗ ¯ σ ( − iω n ) , − G ¯ σ ( − iω n ) (cid:19) = (15)We used the symmetry relations of the unperturbedGreen functions that remain generally valid and read in(complex) energy representation¯ G σ ( iω n ) = − G ¯ σ ( − iω n ) = − G ∗ ¯ σ ( iω n ) , (16a)¯ G σ ( iω n ) = G ∗ ¯ σ ( − iω n ) = −G ∗ ¯ σ ( iω n ) . (16b)We keep the time (charge) propagation (from left toright) in the diagrammatic representation and attach thespin up/down to the upper/lower line. Anomalous prop-agators do not conserve charge by annihilating two elec-trons with opposite spins (arrows against each other)or create a Cooper pair (arrows from each other). Wecan construct standard Feynman many-body diagramsfor processes induced by the Coulomb interaction of theelectrons on the impurity between two superconductingleads. The Coulomb interaction will be represented via awavy line. Since the interaction is static, the interactionwavy line is always vertical. Before we start to analyzethe diagrammatic contributions from the perturbationexpansion we resume basic exact relations.The impact of the Coulomb repulsion on the one-electron Green function is included in a matrix self-energy ˆΣ( z ) so that the full inverse propagator in thespin-polarized situation reads (cid:98) G − ( iω n ) = (cid:98) G − ( iω n ) − (cid:98) Σ( iω n ). Its explicit component representation is (cid:98) G σ ( iω n )= 1 D σ ( iω n ) (cid:18) − X ¯ σ ( − iω n ) , − c Φ Y ( iω n ) − c Φ Y ∗ ( − iω n ) , X σ ( iω n ) (cid:19) (17)with X σ ( iω n ) = iω n [1 + s ( iω n )] + σ [ h − ∆Σ( iω n )] − (cid:15) − Σ( iω n ) , (18a) Y ( iω n ) = s ( iω n )∆ − S ( iω n ) . (18b)We denoted Σ( iω n ), ∆Σ( iω n ) the even and odd partsof the normal self-energy with respect to the magneticfield and S ( iω n ) the anomalous superconducting part of the interaction-induced self-energy. The even, spin-symmetric self energy, Σ( iω n ) and the anomalous one, S ( iω n ), will be determined form the dynamical spin-symmetric Schwinger-Dyson equation. The odd self-energy ∆Σ( iω n ) generalizes the classical order parame-ters and will be related with the two-particle irreduciblevertex via a linearized Ward identity [72].The spin-dependent determinant of the inverse of thematrix propagator in this notation is D σ ( iω n )= − X σ ( iω n ) X ¯ σ ( − iω n ) − c Y ( iω n ) Y ∗ ( − iω n ) , (19)with the electron-hole symmetry D σ ( iω n ) = D − σ ( − iω n ).The normal spin-dependent impurity propagators are G σ ( iω n ) = − X ¯ σ ( − iω n ) D σ ( iω n ) , (20a)¯ G σ ( iω n ) = X σ ( iω n ) D σ ( iω n ) , (20b)while the anomalous propagators are G σ ( iω n ) = − c Φ Y ( iω n ) D σ ( iω n ) , (21a)¯ G σ ( iω n ) = − c Φ Y ∗ ( − iω n ) D σ ( iω n ) . (21b)The existence and positions of the Andreev states areagain determined from zeros of determinant D σ ( iω n ).They depend on the behavior of the normal and anoma-lous self-energies for which we introduce a diagrammaticexpansion. We first formulate the perturbation expan-sion in the thermodynamic language using the Matsub-ara representation. Only after having constructed con-tributions to the perturbation expansion and within theselected approximations we perform analytic continua-tion to real frequencies so that to control the behavior ofthe Andreev bound states (ABS). IV. PERTURBATION EXPANSION: REDUCEDPARQUET EQUATIONS
The basic element of the many-body perturbation ex-pansion is the one-particle propagator. Knowing it wedetermine all the physical quantities. The Dyson equa-tion introduces the self-energy containing the whole im-pact of the particle interactions on the one-particle prop-agator. That is why most of the theoretical approachesfocus on the self-energy. It is, however, not the best wayto control the critical and crossover behavior from weakto strong coupling regimes. Although it is more elaborateand complex in its analytic structure, perturbation the-ory applied directly to two-particle functions has gainedon popularity in recent years. The idea to extend the per-turbation theory and its renormalizations to two-particlefunctions is rather old [85, 86]. Presently, this general ap-proach is used within the so-called parquet equations thatadd a two-particle self-consistency [87, 88]. Generally,the full unrestricted approximations at the two-particlelevel can be solved only numerically and in the Matsu-ubara formalism at non-zero temperatures. One has toresort to simplifications if the critical behavior should becontrolled analytically. We developed the so-called re-duced parquet equations to reach this objective [68–72].The fundamental idea of this two-particle approach is totreat approximate two-particle vertex functions and one-particle self-energy separately and match them at the endso that to keep the theory thermodynamically consistentand conserving.
A. Two-particle vertex: Effective interaction
The fundamental element in the two-particle pertur-bation theory is the two-particle vertex Γ. It has gen-erally three dynamical variables, two fermionic iω n , iω n (cid:48) ,one bosonic iν m , and two spin indices σ, σ (cid:48) . An irre-ducible vertex Λ plays the role of the two-particle self-energy. The two-particle irreducibility is not uniquelydefined and hence there is not a unique way to selectthe irreducible vertex [87]. The most important one is,however, that from the two-particle scattering channelleading to a singularity and a critical behavior in inter-mediate coupling. It is the spin-singlet electron-hole scat-tering channel. The full vertex then can be decomposedinto its irreducible Λ and reducible, K , partsΓ ↑↓ ( iω n , iω n (cid:48) ; iν m ) = Λ ↑↓ ( iω n , iω n (cid:48) ; iν m )+ K ↑↓ ( iω n , iω n (cid:48) ; iν m ) , (22)where ω n and ω n (cid:48) are energies of the incoming and outgo-ing electron, respectively, and ν m is the energy differencebetween the electron and the hole that is conserved inthe multiple singlet electron-hole scatterings.Generally, the reducible vertex in one scattering chan-nel becomes irreducible in the other scattering chan-nels. The parquet equations self-consistently intertwinethem to determine both irreducible and reducible partsof the full vertex. Our approximation resorts to a two-channel version of the parquet equations with only sin-glet electron-hole and electron-electron multiple scatter-ings. The sum of the series of the repeated scatteringsof particle pairs are mathematically represented by theBethe-Salpeter equations. The Bethe-Salpeter equationin the electron-hole channel determines the reducible ver-tex K as a functional of the irreducible one Λ and is dia-grammatically represented in Fig. 1. The irreducible ver-tex in our approximation is determined from a reducedBethe-Salpeter equation that is diagrammatically repre-sented in Fig. 2. The reduction of the full Bethe-Salpeterequation in the electron-electron channel consists in sup-pressing convolution of two diverging reducible vertices K GG K so that not to destroy the possible quantum crit-icality in the strong-coupling regime of the full vertex Γ. The suppressed term is expected to be compensatedby higher-order terms not included in the two-channelapproximation [72].The mean-field approximation enters these reducedparquet equations by replacing the irreducible vertex bya frequency and spin-independent constant Λ that thenplays a role of an effective interaction. The reduciblevertex determined by the equation of Fig 1 is K σ ( iν m ) = − Λ φ σ ( iν m )1 + Λ φ σ ( iν m ) , (23)where the fermionic frequencies are ω n = (2 n + 1) πT and ω n (cid:48) = (2 n (cid:48) + 1) πT , and the bosonic is ν m = 2 mπT . Wedenoted the full electron-hole bubble φ σ ( iν m ) = 1 β (cid:88) ω n [ G ¯ σ ( iω n + iν m ) G σ ( iω n )+ G ¯ σ ( iω n + iν m ) G σ ( iω n )] . (24)The reduced parquet equations are justified in thecritical region of the magnetic transition, that is, in thespin symmetric case where G ↑ = G ↓ . The mean-fieldapproximation must be, however, defined in the wholerepresentation space, including the spin-polarized state.Since we introduced only a spin-independent renormal-ization of the bare interaction strength, we replace thespin-dependent bubble with its symmetric form, that is, φ σ ( iν m ) → φ ( iν m ) = ( φ ↑ ( iν m )+ φ ↓ ( iν m )) / (cid:34) β (cid:88) ν m K ( − iν m ) G ↑ ( iω n + m ) × G ↓ ( iω n (cid:48) − m ) (cid:21) Λ =
U , (25)which cannot, however, be satisfied for all fermionic fre-quencies. An approximate treatment of this equation isnecessary to close the mean-field scheme.The dominant contribution in metallic systems to ver-tex Λ comes from the lowest Matsubara frequencies closeto the Fermi energy, that is | n | ≈ | n (cid:48) | ≈
0. We canthen take the lowest values near the Fermi energy atlow-temperatures as we did in the SIAM [70–72]. TheFermi energy of the superconducting quantum dot lies inthe gap and there is no contribution from small fermionicfrequencies to screening of the interaction. The fluctu-ations in the fermionic Matsubara frequencies may shiftthe value of the critical interaction but do not affect theuniversal critical behavior. We can use averaging over thefermionic Matsubara frequencies to obtain a mean-field(static) renormalization of the bare interaction strengthat any temperature within the same universality class[68]. The averaging is not uniquely defined and the opti-mal one, producing the most accurate result, depends onthe studied problem. We found that the most suitable K = − Λ + Λ + K FIG. 1. Diagrammatic representation of the Bethe-Salpeter equation for the reducible vertex in the electron-hole channel. Theelectron-hole propagator contains simultaneous normal and anomalous propagators. The lines from the central part in thebrackets are attached to the left and right vertices (the three parts separated by brackets are mathematically multiplied) toform two connected diagrams.
Λ = − K Λ FIG. 2. The reduced Bethe-Salpeter equation as explainedin the text for the irreducible vertex from the electron-holescattering channel. The electron-electron propagator does notcontain an anomalous part due to conservation laws. averaging scheme here is to multiply Eq. (25) by a prod-uct G ↑ ( − iω n (cid:48) ) exp( − iω n (cid:48) + ) G ↓ ( − iω n ) exp( − iω n + ) andsum over the fermionic frequencies. The resulting equa-tion for the effective interaction Λ then isΛ = U n ↑ n ↓ n ↑ n ↓ + Λ X (26)where n σ is the density electrons with spin σ and X = − β (cid:88) ν m ψ ( iν m ) ψ ( − iν m ) φ ( − iν m )1 + Λ φ ( − iν m ) . (27)We introduced the electron-electron bubble ψ ( iν m ) = 1 β (cid:88) ω n G ↓ ( iω m + n ) G ↑ ( − iω n )= 1 β (cid:88) ω n G ↑ ( iω m + n ) G ↓ ( − iω n ) , (28)which is spin independent.Equation (26) determines the effective interaction forthe known densities n σ and the screening integral X . Theexplicit solution for Λ can be obtained by a substitutionwith an auxiliary variable w Λ = w − n − m w X , (29a)where we used n σ = ( n + σm ) / n and m being thetotal charge and spin density, respectively. The cube ofthe new variable w satisfies a quadratic equation with asingle positive root w = U (cid:0) n − m (cid:1) X (cid:34) (cid:114) n − m U X (cid:35) . (29b) The known value of w determined the effective interac-tion Λ from Eq. (29a). The consistency condition forpositivity of the effective interaction is12 w X ≥ n − m . (30)Equation (29) does not, however, determine the effec-tive interaction explicitly since integral X depends onthe solution. The final solution can be reached only viaiterations. B. Thermodynamic propagators: Staticself-energies
One cannot close the equation for the two-particle ver-tex Λ without connecting it with the one-particle den-sities. It means that we must determine how the self-energy in the one-particle propagators determining thecharge and spin densities is related with the two-particlevertex from Eq. (26). We introduce two self-energies ac-cording to their symmetry with respect to the spin reflec-tion to keep the theory conserving and thermodynami-cally consistent. We split the self-energy to two. Onewith odd and the other with even symmetry with re-spect to the symmetry-breaking field of the critical pointof the two-particle vertex. They will be related to thetwo-particle vertex differently in approximate schemes.The odd self-energy stands for the order parameteremerging below the critical point with a diverging vertex.The system with the repulsive interaction is driven in in-termediate coupling towards a magnetic order. The oddself-energy must then enter the Ward identity in order tokeep thermodynamic consistency between the criticalityin the two-particle vertex and the order parameter in thesymmetry-broken phase. We argued in previous publica-tions [70–72] that it is sufficient to obey the Ward identityonly in the leading linear order in the symmetry-breaking(magnetic) field to describe the critical behavior qualita-tively correctly. The odd self-energy determined from thestatic irreducible vertex Λ satisfying the linearized Wardidentity is ∆Σ = − Λ2 m . (31)There is no critical behavior in the charge sector andthe even self-energy does not affect the critical behaviornear the transition to the magnetic state. It need notbe related to the two-particle irreducible vertex via theWard identity. It is responsible for the charge dynamicsand should obey the Schwinger-Dyson equation of mo-tion. Its mean-field (static) version is just the Hartree-Fock spin symmetric approximation. We then haveΣ ( ω ) = U n . (32a)Analogously the anomalous self energy, that has no oddpart, is S ( ω ) = U ν , (32b)is proportional to the density of the Cooper pairs on theimpurity.The components determining the one-electron Greenfunction of the superconducting quantum dot in themean-field approximation are X σ ( ω ) = ω [1 + s ( ω )]+ σ (cid:18) h + Λ2 m (cid:19) − (cid:18) (cid:15) + U n (cid:19) , (33a) Y ( ω ) = ∆ s ( ω ) − U ν . (33b)The charge and spin densities are n = 1 β (cid:88) ω n e iω n + [ G ↑ ( iω n ) + G ↓ ( iω n )] , (34a) m = 1 β (cid:88) ω n e iω n + [ G ↑ ( iω n ) − G ↓ ( iω n )] , (34b)and the density of the Cooper pairs on the dot is νc Φ = 12 β (cid:88) ω n e iω n + [ G ↑ ( iω n ) + G ↓ ( iω n )] . (34c)The equations for the effective interaction Λ, the den-sity of Cooper pairs ν , the charge density n , and mag-netization m close our mean-field approximation. It isfree of the unphysical and spurious finite-temperaturetransition to the magnetic state due to the two-particleself-consistency. It qualitatively correctly describes thebehavior of the quantum dot in weak as well as in strongcoupling, including the Kondo regime for the dot at-tached to metallic leads. It can be applied at all tem-peratures and also in an arbitrary magnetic field. Thismean-field approximation serves as the proper startingpoint for the perturbation expansion to include dynami-cal corrections. The mean-field one-particle Green func-tions replace the bare propagators in the perturbationexpansion around the mean-field solution. We call themthermodynamic propagators. C. Spectral representation
The whole mean-field approximation can be fullysolved in the Matsubara formalism. What cannot be de-termined from the Matsubara frequencies are the spectral properties of the one and two-particle Green functions.To determine also the spectral properties one has to per-form analytic continuation to the real frequencies. Oneneeds to rewrite the sum over Matsubara frequencies tointegrals with Fermi and Bose distribution functions.The one-electron Green functions have a gap aroundthe Fermi energy. Since the hybridization self-energy s ( z )has a square-root singularity at the gap/band edges, thegap is fixed in the one-electron Green function and doesnot depend on the interaction strength. The poles andthe band edges of the higher-order Green functions do,however, depend on the interaction strength. We hencemust be careful when treating the two-particle functionsin the spectral representation.The sum over the fermionic Matsubara frequenciesfor the one-particle function can then be rewritten in thespectral representation1 β (cid:88) n F ( iω n ) e iω n + → (cid:88) i f ( ω i ) Res[ F, ω i ] − (cid:34)(cid:90) − ∆ −∞ + (cid:90) ∞ ∆ (cid:35) dωπ f ( ω ) (cid:61) F ( ω + i . (35)Functions with bosonic symmetry have no gap in theirspectra at non-zero temperatures with discrete Matsub-ara frequencies.The Andreev bound states are determined from zerosof the denominator D σ from Eq. (19). The frequenciesof the poles of the one-electron Green function in themean-field approximation are ω σ (1 + s σ ) = − σ (cid:18) h + Λ2 m (cid:19) + (cid:115)(cid:18) (cid:15) + U n (cid:19) + c ( s σ ∆ − U ν ) . (36)The other two frequencies of the gap states are symmet-rically situated on the other side of the Fermi energy.The spectral representation for the densities are n = n g + n b = (cid:88) α,σ f ( ασω σ ) Res[ G σ , ασω σ ] − (cid:88) σ (cid:34)(cid:90) − ∆ −∞ + (cid:90) ∞ ∆ (cid:35) dωπ f ( ω ) (cid:61) G σ ( ω + ) , (37a) m = m g + m b = (cid:88) α,σ σf ( ασω σ ) Res[ G σ , ασω σ ] − (cid:88) σ σ (cid:34)(cid:90) − ∆ −∞ + (cid:90) ∞ ∆ (cid:35) dωπ f ( ω ) (cid:61) G σ ( ω + ) , (37b) c Φ ν = c Φ ( ν g + ν b ) = 12 (cid:88) α,σ f ( ασω σ ) Res[ G σ , ασω σ ] − (cid:88) σ (cid:34)(cid:90) − ∆ −∞ + (cid:90) ∞ ∆ (cid:35) dωπ f ( ω ) (cid:61)G σ ( ω + ) . (37c)We abbreviated the notation of the frequency with aninfinitesimal imaginary part ω + i + = ω + . We split thecontributions to the densities to those from the in-gapstates, subscript g and from the band states, subscript b .Notice that the density of the Cooper pairs from the gapand band states is now spin dependent.The residues of the one-electron Green function areRes [ G σ , σσ (cid:48) ω σ (cid:48) ] = 1 K σ (cid:48) (cid:20) X σ (cid:48) + σσ (cid:48) (cid:18) (cid:15) + U n (cid:19)(cid:21) , (38a)Res [ G σ , σσ (cid:48) ω σ (cid:48) ] = − σ (cid:48) c Φ K σ (cid:48) [ s σ (cid:48) ∆ − U ν ] , (38b)with X σ = (cid:115)(cid:18) (cid:15) + U n (cid:19) + c ( s σ ∆ − U ν ) , (39a) K σ = 2 X σ (cid:20) s σ ∆ − ω σ (cid:21) − c ( s σ ∆ − U ν ) ω σ s σ ∆∆ − ω σ . (39b)The analytic representation of the two-particle Greenand vertex functions is more complex. The integrand ofthe screening integral has no gap at non-zero tempera-tures with a simple analytic representation of the sumover bosonic Matsubara frequencies X = − P (cid:90) ∞−∞ dxπ b ( x ) (cid:61) (cid:20) ψ ( x + ) ψ ( − x + ) φ ( − x + )1 + Λ φ ( − x + ) (cid:21) . (40)The explicit analytic representations separating the gapand band contributions of the electron-hole φ σ ( ω + ) andelectron-electron ψ ( ω + ) are presented in Appendices Aand B. D. Full Green function and the spectral self-energy
The spectral representation is necessary not only todetermine the positions of the in-gap states. It is gen-erally needed to disclose the whole spectral structure ofthe interacting system when we go beyond the mean-fieldapproximation in the perturbation expansion. The firststep beyond the static theory are dynamical correctionsto the static self-energy. The even self-energy is deter-mined from the dynamical Schwinger-Dyson equation ofmotion. Its form with the static irreducible vertex Λ isfor the normal partΣ Sp ( ω + ) = − U (cid:90) ∞−∞ dxπ (cid:26) f ( x ) (cid:61) ¯ G Sp ( x + )1 + Λ φ ( x − ω + ) − b ( x ) ¯ G Sp ( ω + + x ) (cid:61) (cid:20)
11 + Λ φ ( x + ) (cid:21)(cid:27) (41a) and analogously for the anomalous self-energy c Φ S Sp ( ω + ) = − U (cid:90) ∞−∞ dxπ (cid:26) f ( x ) (cid:61) ¯ G Sp ( x + )1 + Λ φ ( x − ω + ) − b ( x ) ¯ G Sp ( ω + + x ) (cid:61) (cid:20)
11 + Λ φ ( x + ) (cid:21)(cid:27) . (41b)The integrand in the Schwinger-Dyson equation con-tains two parts, the two-particle and one-particle ones.The former part, consisting of the electron-hole bub-ble φ ( ω + ) and vertex Λ, controls the thermodynamicresponse and the critical behavior. It hence must bethe same as used to determine the two-particle irre-ducible vertex Λ and the odd self-energy ∆Σ. Theone-particle propagators G Sp ( ω + ) and G Sp ( ω + ) in theSchwinger-Dyson equation carry information about thespectral properties. Its odd self-energy ∆Σ must beidentical with that used to determine the two-particlevertex. Its noncritical even self-energy Σ Sp ( ω + ) canbe selected self-consistently containing the spectral self-energy, a solution of the Schwinger-Dyson equation.Since the Schwinger-Dyson equation determines only thespin-symmetric self-energy we used the spin-averagedpropagators ¯ G Sp ( x + ) = (cid:16) G Sp ↑ ( x + ) + G Sp ↓ ( x + ) (cid:17) / G Sp ( x + ) = (cid:16) G Sp ↑ ( x + ) + G Sp ↓ ( x + ) (cid:17) / G Spσ and G Spσ used in theSchwinger-Dyson equation then are G Spσ ( ω + ) = ω + σ ( h − ∆Σ) + (cid:15) + Σ Sp ( − ω + ) D Spσ ( ω + ) , (42a) G Spσ ( ω + ) = − c Φ s ( ω + )∆ − S Sp ( ω + ) D Spσ ( ω + ) (42b)with the denominator D Spσ ( ω + ) = (cid:2) ω + + σ ( h − ∆Σ) − (cid:15) − Σ Sp ( ω + ) (cid:3) × (cid:2) ω + + σ ( h − ∆Σ) + (cid:15) + Σ Sp ( − ω + ) (cid:3) − c (cid:2) s ( ω + )∆ − S Sp ( ω + ) (cid:3) , (43)where ∆Σ = − Λ m T / m T is the magnetization cal-culated with the thermodynamic propagator determiningthe effective interaction Λ.The normal dynamical self-energy Σ Sp ( ω + ) and theanomalous one S Sp ( ω + ) from the Schwinger-Dyson equa-tion (41) and the odd one ∆Σ from the Ward iden-tity, Eq. (31) are the physical self-energies. It meansthat a mean-field approximation at the two-particlelevel generates nontrivial dynamical contributions to theone-particle self-energy in an analogous manner as therandom-phase approximation generates a dynamical self-energy for the Hartree-Fock mean-field thermodynam-ics. We will analyze the dynamical corrections from theSchwinger-Dyson equation in a separate paper. V. GAP STATES AND − π TRANSITION
The spectral representation is needed for the deter-mination of the positions of the in-gap states and findingthe point of their crossing signaling the 0 − π transition atzero temperature. We need to keep the applied magneticfield positive in order to be able to continue the solutionfrom the weak-coupling 0-phase to the strong-coupling π -phase. We resort to the static mean-field approximationto determine the 0 − π transition.We split the contributions from the band and gapstates and introduce the following abbreviations (cid:15) U = (cid:15) + U n , (44a)Γ σ = s σ ∆ − U ν . (44b)The one-electron parameters are n − n b = n g = 1 K ↑ K ↓ (cid:88) σ K σ [ X ¯ σ − (cid:15) U ∆ f ¯ σ ] , (45) m − m b = m g = 1 K ↑ K ↓ (cid:88) σ K σ [ X ¯ σ ∆ f ¯ σ − (cid:15) U ] , (46) ν − ν b = ν g = 12 K ↑ K ↓ (cid:88) σ K σ Γ ¯ σ ∆ f ¯ σ , (47)where the subscripts b, g refer to the band and gap con-tributions, respectively. We denoted ∆ f σ = f ( − ω σ ) − f ( ω σ ). We recall that the poles of the mean-field prop-agators G σ ( ω ) and G σ ( ω ) are ω σ and − ω ¯ σ . The equa-tions for the in-gap-state frequencies are determined inEq. (36).The 0 − π transition in the external magnetic field ina spin-polarized state happens at ω ↑ = 0, that is h + Λ2 m = (cid:115)(cid:18) (cid:15) + U n (cid:19) + c ( s σ ∆ − U ν ) . (48)This equation tells us that the effective interaction Λ af-fects the transition only in the spin-polarized solutionwith h >
0. The transition in our mean-field approxima-tion with no spectral self-energy at h = 0 coincides withthe Hartree-Fock result. It is, however, important to re-alize that unlike the Hartree-Fock solution the mean-fieldapproximation with an effective interaction Λ is free ofthe spurious transition to the magnetic state at non-zerotemperatures and is thermodynamically consistent in thewhole range of the input parameters. Notice, however,that the effective interaction does affect the position ofthe 0 − π transition in the spin-symmetric state if weemploy the spectral self-energy from Eq. (41). A. Spin-symmetric state
We first approach the 0 − π transition from the weak-coupling regime in the spin-symmetric state. We then have X = (cid:113) (cid:15) U + c ( s ∆ − U ν ) with s = Γ / (∆ − ω )and ω (1 + s ) = X . Further on, K = 2 κ X = 2 (cid:26) s ∆(1 + s ) (∆ − ω ) × (cid:2) ∆ + c U ν + s ∆ (cid:0) − c (cid:1)(cid:3)(cid:9) X . (49)The equation for the positive frequency of the gap stateis (cid:20) (1 + s ) κ ω + U (cid:18) βω (cid:19)(cid:21) = (cid:18) κ (cid:15) b + U (cid:19) + c κ Γ b , (50)where we used an identity ∆ f ( ω ) = tanh ( βω/ (cid:15) b = (cid:15) + U n b / b = s ( ω )∆ − U ν b . Thereis aways a solution for ω > U at non-zero temperature. There is hence no crossing of the in-gap states at non-zero temperature in the spin-symmetricstate as already observed in Ref. [59].The in-gap-state frequency reaches the Fermi energy,that is ω = 0, at the 0 − π transition only at zero tem-perature. The spin-symmetric state can reach the criticalinteraction strength U c of the 0 − π transition only for βω = ∞ defining a quantum critical point. The equa-tion for the critical interaction reads U c (cid:20)(cid:18) (cid:19) (cid:15) b + U c (cid:21) + c (cid:18) (cid:19) Γ b . (51)The equilibrium spin-symmetric solution must be sta-ble with respect to the perturbations caused by a smallmagnetic field. Its local stability is determined from thestatic magnetic susceptibility. It is critically dependenton the effective interaction of the mean-field approxima-tion. The mean-field static susceptibility has the Stonerform χ = − φ (0)1 + Λ φ (0) . (52)The denominator on the right-hand side of Eq. (52)is positive at any temperature and non-diverging atnon-zero temperatures due to the appropriately chosenscreening of the interaction strength in the self-consistentequation (26). The susceptibility can diverge only at zerotemperature in the π -phase as we demonstrate later. B. Magnetic state
We introduce an effective magnetic field containingthe entire effect of the applied magnetic field in the mean-field approximation to simplify the notation h Λ = h (cid:18) m h (cid:19) . (53)0The crossing of the in-gap states takes place when ω ↑ = 0 for which ω ↓ (1 + s ↓ ) = 2 h Λ . Consequently, K ↑ =4 h Λ (1 + Γ / ∆), X ↑ = (cid:15) U + c (Γ − U ν ) and X ↓ = (cid:15) U + c Γ∆ (cid:113) ∆ − ω ↓ − U ν , (54a) K ↓ = 4 h Λ (cid:16) ∆ − ω ↓ (cid:17) / × ∆ − c Γ∆ (cid:113) ∆ − ω ↓ − U ν (54b)We further have ∆ f ↑ = 0 at the crossing point atnon-zero temperatures and hence (cid:15) U = U [ K ↓ ∆ + 2 X ↓ (Γ + ∆)] + 4 X ↓ (Γ + ∆) (cid:15) b (Γ + ∆) (2 K ↓ + U ∆ f ↓ ) , (55a)Γ − U ν = 2 K ↓ Γ b + U ∆ f ↓ (Γ − s ↓ ∆)2 K ↓ + U ∆ f ↓ . (55b)The equation for frequency ω ↓ at the crossing is∆ (cid:18) K ↓ + U tanh (cid:18) βω ↓ (cid:19)(cid:19) = (cid:20) U K ↓ + 4 X ↓ (cid:18) (cid:19) (cid:18) (cid:15) b + U (cid:19)(cid:21) + c (cid:20) K ↓ Γ b + U tanh (cid:18) βω ↓ (cid:19) (Γ − s ↓ ∆) (cid:21) . (56)The crossing leads to the 0- π transition only at zero tem-perature and zero magnetic field and it is a quantum crit-ical point with a diverging magnetic susceptibility whenapproached from the spin-symmetric state.The solution of Eq. (56) shows a universal behaviorfor non-zero magnetic field. We can divide all energyvariables T, U, Λ , ∆ , (cid:15), K, X, Γ, by 2 h > ω ↓ = 1 +¯Λ m will then become universal, independent of the actualvalue of the applied magnetic field. VI. ASYMPTOTIC ATOMIC LIMIT
The important test of reliability of the approxima-tions is a comparison with the existing exact solu-tions in specific limiting situations. The present two-particle approximation with the effective interaction Λ from Eq. (26) and the spectral self-energy from Eq. (41)was shown to reproduce qualitatively correctly the Kondoregime of the exact solution of the SIAM for ∆ = 0.The opposite asymptotic atomic limit ∆ → ∞ can alsobe exactly solved [33, 49, 50, 57, 66]. Here we comparethe predictions of the presented mean-field approxima-tion with the exact results of the atomic limit with nohybridization to the band electrons. The exact results ofthe atomic limit are summarized in Appendix C.The normal and anomalous Green functions in theatomic limit are G σ ( ω ) = 12 X (cid:20) X + (cid:15) U ω − ω σ + X − (cid:15) U ω + ω ¯ σ (cid:21) (57a) G σ ( ω ) = − c Φ Γ2 X (cid:20) ω − ω σ − ω + ω ¯ σ (cid:21) . (57b)where ¯ σ = − σ and ω σ = ¯ σh Λ + X , (58a) X = (cid:113) (cid:15) U + c (Γ − U ν ) . (58b)The explicit value of the electron-hole bubble is φ ( ω + ) = φ ↑↓ ∆ ωω − ∆ ω , (59)where we denoted φ ↑↓ = (∆ f ↓ − ∆ f ↑ ) / − φ (0)∆ ω = m and ∆ ω = ω ↓ − ω ↑ . Further on, φ ( ω + )1 + Λ φ ( ω + ) = φ ↑↓ ∆ ωω − ω φ , (60)where the poles ± ω φ of this function are ω φ = ∆ ω [∆ ω − Λ φ ↑↓ ] = ∆ ω [1 + Λ φ (0)] . (61)The particle-particle bubble in the atomic limit is ψ ( ω )= (∆ f ↑ + ∆ f ↓ )8 X (cid:34) ( X − (cid:15) U ) ω + 2 X − ( X + (cid:15) U ) ω − X (cid:35) . (62)We further use the following identities11 + Λ φ (0) = 1 + Λ m h , (63a)∆ ω = 2 h (cid:18) m h (cid:19) (63b) ω φ = 2 h (cid:114) m h (63c)to represent the screening integral1 X = (∆ f ↓ + ∆ f ↑ ) m X (cid:16) X − ω φ (cid:17) (cid:114) m h (cid:40) ( X + (cid:15) U ) + ( X − (cid:15) U ) X (cid:20) X coth (cid:18) βω φ (cid:19) − ω φ coth( βX ) (cid:21) + (cid:0) X − (cid:15) U (cid:1) (cid:16) X − ω φ (cid:17) (cid:0) X + ω φ (cid:1) coth (cid:18) βω φ (cid:19) − X ω φ coth( βX ) − βω φ (cid:16) X − ω φ (cid:17) ( βX ) . (64)We used ∆ b ( ω ) = b ( ω ) − b ( − ω ) = coth( βω/ (cid:15) U = 0, reduces to X = (∆ f ↓ + ∆ f ↑ ) m X (cid:114) m h (cid:34) X coth (cid:32) βh (cid:114) m h (cid:33) − h (cid:114) m h cosh ( βX ) + βX ( βX ) (cid:35) (65)and in the spin-symmetric state, ω ↑ (cid:37) ω ↓ = X to X = T χ X tanh (cid:18) βX (cid:19) (cid:104) ( X + (cid:15) U ) + ( X − (cid:15) U ) + 2 (cid:0) X − (cid:15) U (cid:1) (cid:105) , (66)where χ = lim h → mh = 2 f ( X )(1 − f ( X )) T − f ( X )(1 − f ( X )) (67)is the magnetic susceptibility.The spin symmetric solution is locally stable at allnon-zero temperatures. If we introduce a generalizedKondo scale a = 1 + Λ φ (0) and assume that the solutionapproaches the critical point a (cid:38) T f (1 − f ) , (68a) a = T C U X n f (1 − f ) (68b)and C = 1 X tanh (cid:18) βX (cid:19) (cid:104) ( X + (cid:15) U ) + ( X − (cid:15) U ) + 2 (cid:0) X − (cid:15) U (cid:1) (cid:105) . (68c)The magnetic susceptibility can diverge only at zerotemperature. The critical region can, however, bereached only if the product f (1 − f ) ≥ T /U with thedecreasing temperature so that U ≥ Λ > U c (cid:18) (cid:15) + U (cid:19) + c Γ , (69) which is the exact result for the 0 − π transition in theatomic limit [49, 57].Resolving the particle and Cooper-pair densities inthe spin-polarized state we obtain n = 2 2 X − (cid:15) (∆ f ↓ + ∆ f ↑ )4 X + U (∆ f ↓ + ∆ f ↑ ) , (70a) ν = Γ (∆ f ↓ + ∆ f ↑ )4 X + U (∆ f ↓ + ∆ f ↑ ) . (70b) (cid:15) U = 4 X δ X + U (∆ f ↑ + ∆ f ↓ )= δ (cid:20) − U U c (∆ f ↑ + ∆ f ↓ ) (cid:21) , (71a)Γ − U ν = 4 X Γ4 X + U (∆ f ↑ + ∆ f ↓ ) , (71b)where we denoted δ = (cid:15) + U/ U c = 2 (cid:112) δ + c Γ .The equation for X needed to obtain n and ν is (cid:20) X + U f ↑ + ∆ f ↓ ) (cid:21) = 4 (cid:0) δ + c Γ (cid:1) = U c . (72)Since X ≥ X = U c − U f ↑ + ∆ f ↓ ) ≥ . (73)Using this solution we obtain n = 1 − δU c (∆ f ↑ + ∆ f ↓ ) , (74a) ν = Γ (∆ f ↑ + ∆ f ↓ )2 U c , (74b)and m = 12 (∆ f ↓ − ∆ f ↑ ) . (75)Since ω ↑ = − ( h +Λ m/ X and ω ↓ = ( h +Λ m/ X , we have always three independent variables to de-termine self-consistently, Λ , X , and m . The three cou-pled equations determining these variables are Eq. (29),2Eq. (73), and Eq. (75). The charge density and thedensity of the Cooper pairs are then calculated fromEqs. (74).The equation determining the crossing of the in-gapstates in the applied magnetic field at non-zero temper-ature is ω ↑ = 0 and ∆ f ↑ = 0, hence U c − U (cid:18) βω ↓ (cid:19) = 2 h (cid:18) m h (cid:19) = ω ↓ . (76)There is always a crossing of the in-gap states for arbi-trary interaction U at non-zero temperature at an appro-priate magnetic field.The zero-temperature solution behaves differently.An infinitesimally small magnetic perturbation generatesa fully polarized magnetic state in the π -phase, U > U c .The π -phase is bounded by ω ↑ <
0, where m = 1 fromEq. (75), n = 1, and ν = 0 from Eqs. (74) for β = ∞ .The effective interaction Λ = U at zero temperature andthe screening integral from Eq. (26) is proportional to n ↑ n ↓ = ( n − m ) / π -phase. It then means that ω ↓ = U c + U while ω ↑ = U c − U with h (cid:38)
0. TheHartree-Fock solution at zero temperature is exact alsoin the π -phase at zero temperature.There is a fundamental difference between the 0-phaseand the π -phase and we can distinguish the two phasesby the low-temperature asymptotics of the magnetic sus-ceptibility. We have 2 X = U c − U in the leadinglow-temperature asymptotics in the 0-phase, that is, for U < U c . Consequently, the magnetic susceptibility fromEq. (67) at low temperatures is χ . = 8 T e − β ( U c − U ) , (77)which reflects the Meissner effect due to the presence ofthe singlet Bound states (ABS).The situation in the π -phase, U > U c , is quite differ-ent. Equation (73) leads at low temperatures to βX = arctan (cid:18) U c U (cid:19) (cid:20) − T UU − U c (cid:21) . (78)The magnetic susceptibility is χ = U − U c U a (79)with the effective interaction and the Kondo scale fromEqs. (68) areΛ = 2 U TU − U c , (80) a = U U c U − U c ) arctan − (cid:18) U c U (cid:19) T . (81)The magnetic susceptibility follows the Curie law dueto the presence of the local magnetic moment of the fermionic excitations on the dot. The non-universal Curieconstant C = (cid:0) U − U c (cid:1) U U c arctan (cid:18) U c U (cid:19) (82)is in the static mean-field approximation overestimatedand grows with the increasing interaction strength. Itlies above the exact value and will be corrected by thedynamical spectral self-energy.The observed behavior of the in the strong-couplinglimit for U > U c discloses another feature of the solu-tion in the external magnetic filed. We proved that thelimit to zero magnetic field and to zero temperature donot commute. If we keep the magnetic field non-zeroand limit the temperature to zero the solution behavesanalytically and continuously reaches the fully saturatedstate at zero temperature. On the other hand, if we keepthe temperature non-zero and switch off the magneticfield we stay in the spin-symmetric states down to zerotemperature where the magnetic susceptibility divergesand an infinitesimally small magnetic field lifts the de-generacy to a magnetically saturated state. VII. NUMERICAL RESULTS
We apply the mean-field approximation in the atomiclimit to show its similarities and stress the substantialdifferences to the Hartree-Fock solution and to test thereliability of our mean-field approximation in differentregimes. The major asset of the mean-field approxima-tions is that they can be used in the whole range of themodel parameters. They are qualitatively reliable if theydo not lead to an unphysical and spurious behavior. Weknow that the reduced parquet equations with a two-particle self-consistency reproduce qualitatively correctlythe limit of the zero gap. It is instructive to apply itin the opposite limit of the infinite superconducting gapwhich is the least fitting situation for the application ofthe many-body construction. Many of the qualitativefeatures of the solution in the atomic limit are genericand mimic the behavior of the finite-gap model exceptfor the Kondo limit of the vanishing gap.
A. Spin-symmetric solution
Our static mean-field solution in the spin-symmetricstate, that is in the absence of the magnetic field, coin-cides with the Hartree-Fock approximation. This mayseem a limiting factor, but it holds only at the one-particle level when the dynamical corrections in the spec-tral self-energy from the Schwinger-Dyson equation areneglected. The first, and the most important differencebetween our mean-field theory and the Hartree-Fock ap-proximation in the spin-symmetric state is the stabilitywith respect to small fluctuations of the magnetic field,which is reflected in the magnetic susceptibility.3 T / U = ExactHFRPE T / U = 4 ExactHFRPE
FIG. 3. Magnetic susceptibility as a function of tempera-ture in the spin-symmetric state at half-filling in the 0- phase( U = Γ) and the π -phase ( U = 4Γ) for the phase differ-ence Φ = 0. Hartree-Fock (HF), reduced parquet equations(RPE) and exact (EXACT) solutions are compared. The un-physical instability with the diverging susceptibility makes theHartree-Fock mean-field solution in strong coupling unreliableat low temperature@‘s. T / U = RPEHF T / U = 4 RPEHF
FIG. 4. Different temperature behavior of the vertex renor-malization in the spin-symmetric state at half filling in the0-phase and the π -phase for the phase difference Φ = 0. We plotted the magnetic susceptibility of the spin-symmetric state at half filling as a function of tempera-ture in Fig. 3. We compared the two mean-field approx-imations, our, based on the reduced parquet equations(RPE), and the Hartree-Fock one (HF), with the exactsolution in the atomic limit. There is no big difference inthe 0-phase where all solutions asymptotically approachzero at zero temperature. Quite a different behavior is,however, observed in the π -phase. Both the exact andRPE solutions lead to a divergent susceptibility at zerotemperature, while the HF solution predicts an unphys-ical critical point with diverging susceptibility at a tem-perature of order of the hybridization strength Γ that istaken as the energy unit. The magnetic susceptibility isa physical, measurable quantity being able to distinguishthe character of the in-gap states. The in-gap states inthe 0-phase are bound pairs, ABS singlets being insensi-tive to small magnetic perturbations. The in-gap statesin the π -phase carry a local magnetic moment and reactstrongly on magnetic perturbations. The whole π -phasedisplays a Curie susceptibility diverging at zero temper-ature.The reason why the RPE suppress the HF instabil-ity is the two-particle self-consistency renormalizing thebare interaction strength U to a screened effective one Λ. T / U = ExactHF & RPE T / U = 4 ExactHF & RPE
FIG. 5. Density of the Cooper pairs ν in the spin-symmetricstate at half filling as a function of temperature in the 0-phaseand the π -phase for the phase difference Φ = 0. T / U = U = 4 U = 8 FIG. 6. Temperature dependence of the positive energy ofthe in-gap state in the spin-symmetric solution for weak andstrong interactions for the phase difference Φ = 0.
We plotted its temperature dependence at half filling inFig. 4. The renormalization gets stronger with the de-creasing temperature but starts abating in the zero phaseand dies out at zero temperature. The effective interac-tion approaches zero, maximizing the renormalization ofthe interaction, in the π -phase consistent the divergenceof the magnetic susceptibility of the exact solution.The thermodynamic mean-field solution with a staticself-energy produces good results for quantities with oddsymmetry and sensitive to the symmetry-breaking field.It is quantitatively less accurate in determining the spin-symmetric one-particle quantities with even symmetrywith respect to spin flips. We plotted the temperaturedependence of the Cooper-pair density ν at half filling inthe RPE/HF approximation together with the exact re-sult in Fig. 5. We can see how the static spin-symmetricvalue deviates from the exact value at low temperaturesof the π -phase. Unlike the HF mean-field, the RPE offera direct improvement by including the dynamical cor-rections from the Schwinger-Dyson equation (41). Ituses the renormalized interaction and this interactionis strongly renormalized at low temperatures of the π -phase. The Cooper-pair density will then better followthe dependence of the effective interaction as observedin the 0-phase in Fig. 4. We discuss the impact of thedynamical corrections to the static self-energy in detailelsewhere.4 T =0 ++ T =0.05 ++ T =0.1 ++ T = ++ FIG. 7. In-gap-state energies as a function of the phase differ-ence Φ between the superconducting leads in a weak magneticfield h = 0 .
01Γ for different temperatures at half filling and U = Γ. The critical angle of the crossing increases with tem-perature. The spin-symmetric solution of the many-bodyGreen-function approach cannot be extended to the π -phase at zero temperature since one has to cross thequantum critical point. One can nevertheless circumventthe critical point in that one extends the spin-symmetricsolution to non-zero temperatures. There is no criti-cal point at non-zero temperature and the solution canbe extended continuously from weak to strong coupling[59]. No crossing happens and the energy of the in-gapstate remains positive and approaches the quantum criti-cal point at zero temperature, as demonstrated in Fig. 6.The spin-symmetric solution in the π -phase becomes un-stable there and decays into the degenerate spin doubletwith a saturated local magnetic moment. B. Zeeman field
The magnetic field acting on the spin of the electrons(Zeeman field) plays an essential role in the applicationof the many-body Green functions in the superconduct-ing quantum dot. It is needed to approach the zero-temperature solution in the π -phase and to see the cross-ing of the in-gap states at non-zero temperatures. Thedoublet ground state, π -phase, is degenerate and the Zee-man field is the means to lift the degeneracy. Here weanalyze the properties of the low-temperature solutionwith an applied magnetic field.We plotted the dependence of the in-gap-state ener-gies on the phase difference between the attached super-conducting leads in Fig. 7 and on the interaction strengthin Fig. 8 for a very small magnetic field h = 0 .
01Γ at dif-ferent temperatures, The value of the magnetic field atthe crossing increases with temperature. The curves ofthe in-gap-state energies are continuous due the presenceof the symmetry-breaking magnetic field. We also plot- T = 0 ++ T = 0.1 ++ U / T = 0.5 ++ U / T = ++ FIG. 8. In-gap-state energies as a function of the Coulombrepulsion U in a weak magnetic field h = 0 .
01Γ for differenttemperatures at half filling and phase difference Φ = 0. Thecritical angle of the crossing increases with temperature.
15 10 5 0 5 10 15( + U /2)/1050510 ++ FIG. 9. In-gap-state energies as a function of the impurityenergy level (cid:15) for a low temperature, T = 0 . h = 0 . U = 8Γ, and the phase differenceΦ = π/ ted the dependence of the in-gap-state energies on theimpurity energy level (cid:15) + U/ U = 8Γand phase difference Φ = π/
2. We used a small magneticfield h = 0 .
2Γ to demonstrate the expected behavior inthe π -phase, cf. Fig. 9. Note that the RPE solution re-produces the exact positions of the in-gap states in thelimit T → h → m shows just the opposite dependence on temperature thanthe Cooper-pair density ν in the weak Zeeman field ( h =0 . π -phase, see Fig. 10.The magnetization vanishes and ν saturates in the 0-phase in both mean-field approximations as well as in theexact one. In the π -phase the temperature asymptoticsis inverted in both quantities in all solutions. The HFmean-field better simulates the exact dependence of theCooper-pair density ν while the RPE mean-field thenbetter fits the exact magnetization curve.The renormalized interaction Λ in the Zeeman field5 m U = ExactHFRPE U =4 ExactHFRPE T / U = ExactHFRPE T / U =4 ExactHFRPE
FIG. 10. Temperature dependence of magnetization m andthe Cooper-pair density ν at half filling, h = 0 . U = Γ) and the π -phase ( U = 4Γ) for themean-field solutions RPE and HF compared with the exactbehavior. We see that the exact behavior of the magnetiza-tion, with odd symmetry with respect to spin flip, is betterreproduced by the RPE while the Cooper-pair density, witheven symmetry, is better reproduced by the HF approxima-tion. T / U = RPEHF T / U = 4 RPEHF
FIG. 11. The temperature dependence of the renormalized in-teraction strength Λ at half filling in a Zeeman field h = 0 . has a different low-temperature asymptotics in the π -phase than in the spin-symmetric case, see Fig. 11. Oncethe magnetic field is kept non-zero down to zero temper-ature the effective interaction approaches the bare valueand the exact solution is reproduced.Although the HF solution quantitatively better ap-proximates the exact behavior of the thermodynamicquantities with even symmetry with respect to spin flips,it fails to reproduce the exact limit of the vanishing Zee-man field since it predicts a non-zero magnetization atzero field below its critical transition to the magneticstate as documented in Fig. 12. It is our mean-field thatqualitatively correctly reproduces all the limits of bothone and two-particle thermodynamic quantities. h / m U = ExactHFRPE h / U = 4 ExactHFRPE
FIG. 12. Magnetic-field dependence of magnetization m inweak ( U = Γ) and strong ( U = 4Γ) couplings for Φ = 0 and T = 0 .
5Γ and half filling. The strong-coupling HF solution isbelow its magnetic critical point is completely off in the limitto zero magnetic field unlike the RPE solution.
VIII. CONCLUSIONS
The quantum dot attached to superconducting leadsposes a challenging problem to the perturbation theorywith many-body Green functions. First, electron corre-lations on the dot lead to a line of first-order transitionsfrom the spin-singlet to the spin-doublet state that endsup at a quantum critical point at zero temperature wherethe in-gap states reach the Fermi energy. Moreover, thedoublet state is degenerate and it cannot be continuouslyapproached from the weak-coupling spin-singlet state.The basic assumptions of the applicability of the many-body perturbation theory is a non-degenerate groundstate and the existence of an analyticity region from theweak-coupling limit within which it can be applied. Itcannot cope with two independent many-body equilib-rium states with a first order transition between them. Itmay lead to a new equilibrium state, phase, only if thereis a divergence accompanied by a continuous symmetrybreaking. It means that the many-body perturbation ex-pansion cannot reliably be applied to the superconduct-ing quantum dot at low temperatures around the quan-tum critical point unless a self-consistency is introduced.A self-consistency must be introduced in the many-bodyperturbation expansion in order to deal with the quan-tum critical behavior. A static mean-field approximationis the simplest way to achieve this goal.It has been known for long that the weak-couplingHartree-Fock self-consistency is not appropriate to dealwith the quantum critical point of the superconductingquantum dot since it fails at non-zero temperatures whereit leads to a spurious transition to a magnetically orderedstate without the external magnetic field. We added atwo-particle self-consistency to the HF solution in thatwe replaced the bare interaction with a renormalized,screened one. We thereby suppressed the HF spurioustransition to the magnetic state and produced a fullythermodynamically consistent mean-field approximationapplicable in the whole range of the input parameters.We demonstrated that it is able to deal qualitatively cor-rectly with the quantum critical behavior of the 0 − π − π transi-tion and in distinguishing the spin-singlet from the spin-doublet. The magnetic field not only lifts the degener-acy of the π -phase it allows us to determine the differentcharacter of the in-gap states in the two phases. The0 − π transition is signaled by a crossing of the energiesof the in-gap states. The in-gap states in the spin-singletphase are the genuine Andreev bound states of two elec-trons with opposite spins that are insensitive to smallmagnetic perturbations. The low-lying excitations in thespin-doublet phase are fermions, carry a local magneticmoment, and are sensitive to the magnetic filed. Themagnetic susceptibility vanishes in the 0-phase and di-verges in the π -phase at zero temperature. The equi-librium state at non-zero temperatures turns magneticonly when the Zeeman field is applied. Consequently, theweak-coupling spin-symmetric solution can continuouslybe extended to strong coupling at non-zero temperatureswithout crossing any critical point. The limit to zero fieldleads to a magnetic state only at zero temperature abovethe critical interaction strength of the 0 − π transition.The limits to zero magnetic field and to zero temperaturedo not commute and the results depend on their orderin which they are performed. The Zeeman field playsthe role of a symmetry-breaking field we know from con-tinuous phase transitions and the π -phase beyond thequantum critical point mimics the ordered phase in thelattice models. It means that the 0 − π transition happensonly at zero temperature and zero magnetic field and itis a true local quantum phase transition. The crossingof the in-gap states in the magnetic field or at non-zerotemperatures is noncritical with no phase transition.The mean-field theory presented in this paper is thefirst fully consistent analytic approximation that can de-scribe not only the critical behavior of the 0 − π transitionbut it can qualitatively correctly and reliably reproducethe behavior of the quantum dot attached to both super- conducting and normal leads in the whole range of themodel parameters. It was derived within the perturba-tion expansion for two-particle vertices to control theircritical behavior. It offers a starting point for adding dy-namical corrections in a systematic way. The first step,without changing the static irreducible vertex, is to usethe Schwinger-Dyson equation to determine the spectralproperties of the model. This opens a new way to in-clude dynamical fluctuation in the thermodynamic andspectral properties with the controlled renormalizationsof the one- and two-particle functions. ACKNOWLEDGMENT
Research on this problem was supported in part byGrants 19-13525S of the Czech Science Foundation. VJthanks the INTER-COST LTC19045 Program of theCzech Ministry of Education, Youth and Sports for finan-cial support. We thank Tom´aˇs Novotn´y for illuminatingdiscussions.
Appendix A: Spectral representation - Electron-holebubble
The thermodynamic quantities including the effec-tive interaction can be calculated entirely in the Mat-subara formalism without the necessity to continue ana-lytically to real frequencies. If we want, however, to usethe Schwinger-Dyson equation and determine the spec-tral properties we need a spectral representation of thetwo-particle bubbles, at least the electron-hole one.We decompose the imaginary part of the electron-holebubble φ ( ω + ) to a sum of three contributions (cid:61) φ ( ω + ) = (cid:61) φ bb ( ω + ) + (cid:61) φ bg ( ω + ) + (cid:61) φ gg ( ω + ), according to whetherthe arguments of the Green functions of the integrandlie both within the band, one within the band and onewithin the gap, and both within the gap, respectively.We have for ω > (cid:61) φ bb ( ω + ) = − (cid:88) σ (cid:34)(cid:90) − ∆ − ω −∞ + (cid:90) − ∆min(∆ − ω, − ∆) + (cid:90) ∞ ∆ (cid:35) dx π [ f ( x ) − f ( x + ω )] [ (cid:61) G σ ( x + + ω ) (cid:61) G − σ ( x + )+ (cid:61)G σ ( x + + ω ) (cid:61)G − σ ( x + )] , (A1a) (cid:61) φ bg ( ω + ) = − (cid:88) σ (cid:34)(cid:90) min(∆ − ω, − ∆) − ∆ − ω + (cid:90) ∆max(∆ − ω, − ∆) (cid:35) dx π [ f ( x ) − f ( x + ω )] [ (cid:61) G σ ( x + + ω ) (cid:61) G − σ ( x + )+ (cid:61)G σ ( x + + ω ) (cid:61)G − σ ( x + )] , (A1b) (cid:61) φ gg ( ω + ) = − (cid:88) σ (cid:90) max(∆ − ω, − ∆) − ∆ dx π [ f ( x ) − f ( x + ω )] [ (cid:61) G σ ( x + + ω ) (cid:61) G − σ ( x + ) + (cid:61)G σ ( x + + ω ) (cid:61)G − σ ( x + )] , (A1c)7and for ω < (cid:61) φ bb ( ω + ) = − (cid:88) σ (cid:34)(cid:90) − ∆ −∞ + (cid:90) max(∆ , − ∆ − ω )∆ + (cid:90) ∞ ∆ − ω (cid:35) dx π [ f ( x ) − f ( x + ω )] [ (cid:61) G σ ( x + + ω ) (cid:61) G − σ ( x + )+ (cid:61)G σ ( x + + ω ) (cid:61)G − σ ( x + )] , (A2a) (cid:61) φ bg ( ω + ) = − (cid:88) σ (cid:34)(cid:90) min(∆ , − ∆ − ω ) − ∆ + (cid:90) ∆ − ω max(∆ , − ∆ − ω ) (cid:35) dx π [ f ( x ) − f ( x + ω )] [ (cid:61) G σ ( x + + ω ) (cid:61) G − σ ( x + )+ (cid:61)G σ ( x + + ω ) (cid:61)G − σ ( x + )] , (A2b) (cid:61) φ gg ( ω + ) = − (cid:88) σ (cid:90) ∆min(∆ , − ∆ − ω ) dx π [ f ( x ) − f ( x + ω )] [ (cid:61) G σ ( x + + ω ) (cid:61) G − σ ( x + ) + (cid:61)G σ ( x + + ω ) (cid:61)G − σ ( x + )] . (A2c)The subscript at ω + = ω + i + denotes the way the realaxis is reached from the complex plane.The real part of the bubble is then determined fromthe Kramers-Kronig relation (cid:60) φ ( ω ) = P (cid:90) ∞−∞ dxπ (cid:61) φ ( x + ) x − ω + φ ( ∞ ) (A3) Appendix B: Spectral representation -Electron-electron bubble
The electron-electron bubble has a simpler spectralrepresentation. It is not needed for the spectral self- energy, but its spectral representation is useful the de-termination of the effective interaction Λ at low temper-atures with a high precision. It has no contribution fromanomalous Green functions. Its imaginary part can berepresented as (cid:61) ψ ( ω + ) = − (cid:90) ∞−∞ dxπ [ f x ) − f ( x − ω )] × (cid:61) G ↑ ( ω + − x ) (cid:61) G ↓ ( x + ) (B1)Taking into account the induced gap on the dot we canrepresent the three contribution from the band and gapstates. (cid:61) ψ bb ( ω + ) = − (cid:34)(cid:90) min( − ∆ ,ω − ∆) −∞ + (cid:90) ∞ max(∆ ,ω +∆) (cid:35) dxπ [ f ( x ) − f ( x − ω )] (cid:61) G ↑ ( ω + − x ) (cid:61) G ↓ ( x + ) , (B2a) (cid:61) ψ bg ( ω + ) = − (cid:34)(cid:90) max( − ∆ ,ω − ∆)min( − ∆ ,ω − ∆) + (cid:90) max(∆ ,ω +∆)min(∆ ,ω +∆) (cid:35) dxπ [ f ( x ) − f ( x − ω )] (cid:61) G ↑ ( ω + − x ) (cid:61) G ↓ ( x + ) , (B2b) (cid:61) ψ gg ( ω + ) = − (cid:90) min(∆ ,ω +∆)max( − ∆ ,ω − ∆) dxπ [ f ( x ) − f ( x − ω )] (cid:61) G ↑ ( ω + − x ) (cid:61) G ↓ ( x + ) . (B2c) Appendix C: Atomic limit - Exact solution
We summarize the basic results of the exact solutionof the atomic limit with infinite superconducting gap.The atomic Hamiltonian is a matrix H = c Φ (cid:15) d + h (cid:15) d − h c Φ (cid:15) d + U . (C1) One can simply diagonalize the Hamiltonian matrix,Eq.(C1), by observing that the central 2 × E − d = (cid:15) d − h with theeigenstates | , (cid:105) ; (ii) E + d = (cid:15) d + h with the eigenstates | , (cid:105) ; (iii) E − s = (cid:104) (cid:15) d + U − (cid:112) (2 (cid:15) d + U ) + 4Γ c (cid:105) with the eigenstates √ Γ c +( E − s ) (Γ c Φ | , (cid:105) + E − s | , (cid:105) );and (iv) E + s = (cid:104) (cid:15) d + U + (cid:112) (2 (cid:15) d + U ) + 4Γ c (cid:105) with8the eigenstates √ Γ c +( E + s ) (Γ c Φ | , (cid:105) + E + s | , (cid:105) ).The general thermodynamic quantity is Q = 1 Z (cid:88) i (cid:104) E i | ˆ Q | E i (cid:105) e − β ( E i − µN i ) , (C2)where E i and | E i (cid:105) are the eigenvalues and correspond-ing eigenstates of the Hamiltonian, β = 1 /k B T , and thepartition function is Z = e − βE − s + e − βE + s + e − β ( (cid:15) d + h ) + e − β ( (cid:15) d − h ) . (C3)The thermodynamic properties can, alternatively, bederived from the derivatives of the grand potential F = − β ln Z . The charge and spin densities are n = 1 Z (cid:34)(cid:88) σ e − β ( (cid:15) d − σh ) + 2( E − s ) Γ c φ + ( E − s ) e − βE − s + 2( E + s ) Γ c φ + ( E + s ) e − βE + s (cid:35) , (C4a) m = 1 Z (cid:34)(cid:88) σ σe − β ( (cid:15) d − σh ) (cid:35) . (C4b)The density of the Cooper pairs is ν = 1 Z (cid:20) Γ E − s Γ c + ( E − s ) e − βE − s + Γ E + s Γ c + ( E + s ) e − βE + s (cid:21) . (C5) [1] D. C. Ralph and R. A. Buhrman, Physical Review Letters , 3401 (1994).[2] V. Madhavan, W. Chen, T. Jamneala, M. Crommie, andN. Wingreen, Science , 567 (1998).[3] J. Li, W.-D. Schneider, R. Berndt, and B. Delley, Phys-ical Review Letters , 2893 (1998).[4] D. Goldhaber-Gordon, H. Shtrikman, D. Mahalu,D. Abusch-Magder, U. Meirav, and M. A. Kastner, Na-ture , 156 (1998).[5] D. Goldhaber-Gordon, J. G¨ores, M. A. Kastner,H. Shtrikman, D. Mahalu, and U. Meirav, Physical Re-view Letters , 5225 (1998).[6] S. Cronenwett, T. Oosterkamp, and L. Kouwenhoven,Science , 540 (1998).[7] G. Katsaros, P. Spathis, M. Stoffel, F. Fournel,M. Mongillo, V. Bouchiat, F. Lefloch, A. Rastelli,O. Schmidt, and S. De Franceschi, Nature Nano , 458(2010).[8] J. A. van Dam, Y. V. Nazarov, E. P. A. M. Bakkers,S. De Franceschi, and L. P. Kouwenhoven, Nature ,667 (2006).[9] W. Chang, V. E. Manucharyan, T. S. Jespersen,J. Nyg˚ard, and C. M. Marcus, Physical Review Letters , 217005 (2013).[10] A. Y. Kasumov, R. Deblock, M. Kociak, B. Reulet,H. Bouchiat, I. I. Khodos, Y. B. Gorbatov, V. T. Volkov,C. Journet, and M. Burghard, Science , 1508 (1999).[11] A. F. Morpurgo, J. Kong, C. M. Marcus, and H. Dai,Science , 263 (1999).[12] A. Kasumov, M. Kociak, M. Ferrier, R. Deblock,S. Gu´eron, B. Reulet, I. Khodos, O. St´ephan, andH. Bouchiat, Phys. Rev. B , 214521 (2003).[13] P. Jarillo-Herrero, J. A. van Dam, and L. P. Kouwen-hoven, Nature , 953 (2006).[14] H. I. Jørgensen, K. Grove-Rasmussen, T. Novotn´y,K. Flensberg, and P. E. Lindelof, Phys. Rev. Lett. ,207003 (2006). [15] J. P. Cleuziou, W. Wernsdorfer, V. Bouchiat, T. Ondar-cuhu, and M. Monthioux, Nat. Nano , 53 (2006).[16] H. I. Jørgensen, T. Novotn´y, K. Grove-Rasmussen,K. Flensberg, and P. E. Lindelof, Nano Lett. , 2441(2007).[17] K. Grove-Rasmussen, H. I. Jørgensen, and P. E. Lindelof,New Journal of Physics , 124 (2007).[18] E. Pallecchi, M. Gaass, D. A. Ryndyk, and C. Strunk,Applied Physics Lett. , 072501 (2008).[19] Y. Zhang, G. Liu, and C. Lau, Nano Research , 145(2008-08-01).[20] H. I. Jørgensen, K. Grove-Rasmussen, K. Flensberg, andP. E. Lindelof, Phys. Rev. B , 155441 (2009).[21] A. Eichler, R. Deblock, M. Weiss, C. Karrasch, V. Meden,C. Sch¨onenberger, and H. Bouchiat, Phys. Rev. B ,161407 (2009).[22] G. Liu, Y. Zhang, and C. N. Lau, Phys. Rev. Lett. ,016803 (2009).[23] J.-D. Pillet, C. H. L. Quay, P. Morfin, C. Bena, A. L.Yeyati, and P. Joyez, Nat. Phys. , 965 (2010).[24] E. J. H. Lee, X. Jiang, R. Aguado, G. Katsaros, C. M.Lieber, and S. De Franceschi, Phys. Rev. Lett. ,186802 (2012).[25] R. Maurand, T. Meng, E. Bonet, S. Florens, L. Marty,and W. Wernsdorfer, Phys. Rev. X , 011009 (2012).[26] J. D. Pillet, P. Joyez, R. ˇZitko, and M. F. Goffman,Phys. Rev. B , 045101 (2013).[27] A. Kumar, M. Gaim, D. Steininger, A. L. Yeyati,A. Mart´ın-Rodero, A. K. H¨uttel, and C. Strunk, Phys.Rev. B , 075428 (2014).[28] R. Delagrange, D. J. Luitz, R. Weil, A. Kasumov,V. Meden, H. Bouchiat, and R. Deblock, Phys. Rev.B , 241401(R) (2015).[29] C. B. Winkelmann, N. Roch, W. Wernsdorfer, V. Bouch-iat, and F. Balestro, Nature Phys. , 876 (2009).[30] S. De Franceschi, L. Kouwenhoven, C. Sch¨onenberger,and W. Wernsdorfer, Nat. Nano , 703 (2010). [31] T. Matsuura, Prog. Theor. Phys. , 1823 (1977).[32] L. I. Glazman and K. A. Matveev, JETP Lett. , 659(1989).[33] A. V. Rozhkov and D. P. Arovas, Physical Review B ,6687 (2000).[34] M. R. Buitelaar, T. Nussbaumer, and C. Sch¨onenberger,Phys. Rev. Lett. , 256801 (2002).[35] T. Aono, A. Golub, and Y. Avishai, MolecularNanowires and Other Quantum Objects , 203 (2004).[36] M. R. Gr¨aber, T. Nussbaumer, W. Belzig, andC. Sch¨onenberger, Nanotechnology , S479 (2004).[37] F. Siano and R. Egger, Phys. Rev. Lett. , 047002(2004).[38] M.-S. Choi, M. Lee, K. Kang, and W. Belzig, Phys. Rev.B , 020502 (2004).[39] Y. Tanaka, A. Oguri, and A. C. Hewson, New. J Phys. , 115 (2007).[40] J. S. Lim and M.-S. Choi, J. Phys.: Condens. Matter ,415225 (2008).[41] C. Karrasch, A. Oguri, and V. Meden, Phys. Rev. B ,024517 (2008).[42] Y. Yamada, Y. Tanaka, and N. Kawakami, Journal ofthe Physical Society of Japan , 043705 (2010).[43] Y. Yamada, Y. Tanaka, and N. Kawakami, Physical Re-view B , 075484 (2011).[44] D. J. Luitz, F. F. Assaad, T. Novotn´y, C. Karrasch, andV. Meden, Phys. Rev. Lett. , 227001 (2012).[45] A. Oguri, Y. Tanaka, and J. Bauer, Phys. Rev. B ,075432 (2013).[46] T. Yoshioka and Y. Ohashi, J. Phys. Soc. Jpn. , 1812(2000).[47] G. Sellier, T. Kopp, J. Kroha, and Y. S. Barash, Phys.Rev. B , 174502 (2005).[48] T. Novotn´y, A. Rossini, and K. Flensberg, Phys. Rev.B , 224502 (2005).[49] T. Meng, S. Florens, and P. Simon, Physical Review B , 224521 (2009).[50] J. Bauer, A. Oguri, and A. C. Hewson, J. Phys.: Con-dens. Matter , 486211 (2007).[51] T. Hecht, A. Weichselbaum, J. von Delft, and R. Bulla,J. Phys.: Condens. Matter , 275213 (2008).[52] A. Oguri, Y. Tanaka, and A. C. Hewson, J. Phys. Soc.Jpn. , 2494 (2004).[53] A. Mart´ın-Rodero and A. L. Yeyati, J. Phys.: Condens.Matter , 385303 (2012).[54] D. J. Luitz and F. F. Assaad, Phys. Rev. B , 024509(2010).[55] V. Pokorn´y and M. ˇZonda, Physica B: Condensed Matter , 488 (2018).[56] A. T. Alastalo, R. J. Joynt, and M. M. Salomaa, Journalof Physics: Condensed Matter , L63 (1998).[57] E. Vecino, A. Martin-Rodero, and A. L. Yeyati, PhysicalReview B , 035105 (2003).[58] M. ˇZonda, V. Pokorn´y, V. Janiˇs, and T. Novotn´y, Sci-entific Reports , 8821 (2015).[59] V. Janiˇs, V. Pokorn´y, and M. ˇZonda, European PhysicalJournal B , 197 (2016).[60] M. ˇZonda, V. Pokorn´y, V. Janiˇs, and T. Novotn´y, Phys-ical Review B , 024523 (2016).[61] T. Doma´nski, M. ˇZonda, V. Pokorn´y, G. G´orski,V. Janiˇs, and T. Novotn´y, Physical Review B , 045104(2017).[62] M. Governale, M. G. Pala, and J. K¨onig, Physical Re- view B , 134513 (2008).[63] T. Meng, S. Florens, and P. Simon, Phys. Rev. B ,224521 (2009).[64] S. Droste, S. Andergassen, and J. Splettstoesser, Journalof Physics: Condensed Matter , 415301 (2012).[65] N. Wentzell, S. Florens, T. Meng, V. Meden, and S. An-dergassen, Physical Review B , 085151 (2016).[66] V. Meden, Journal of Physics: Condensed Matter ,163001 (2019).[67] A. Martin-Rodero and A. L. Yeyati, Journal of Physics-Condensed Matter , 385303 (2012).[68] V. Janiˇs and P. Augustinsk´y, Physical Review B ,165108 (2007).[69] V. Janiˇs and P. Augustinsk´y, Physical Review B ,085106 (2008).[70] V. Janiˇs, A. Kauch, and V. Pokorn´y, Physical Review B , 045108 (2017).[71] V. Janiˇs, V. Pokorn´y, and A. Kauch, Physical Review B , 165113 (2017).[72] V. Janiˇs, P. Zalom, V. Pokorn´y, and A. Kl´ıˇc, PhysicalReview B , 195114 (2019).[73] L. Cornils, A. Kamlapure, L. Zhou, S. Pradhan, A. Kha-jetoorians, J. Fransson, J. Wiebe, and R. Wiesendanger,Physical Review Letters , 197002 (2017).[74] T. Dvir, M. Aprili, C. H. Quay, and H. Steinberg, Phys-ical Review Letters , 217003 (2019).[75] A. G. Corral, D. M. T. van Zanten, K. J. Franke, H. Cour-tois, S. Florens, and C. B. Winkelmann, Physical ReviewResearch , 012065 (2020).[76] A. M. Whiticar, A. Fornieri, A. Banerjee, A. C. C.Drachmann, S. Gronin, G. C. Gardner, T. Lindemann,M. J. Manfra, and C. M. Marcus, Eprint Archive ,arXiv:2101.09706 (2021).[77] Y. Tanaka, S. Kashiwaya, and H. Takayanagi, JapaneseJournal of Applied Physics , 4566 (1995).[78] Y. Yamada, Y. Tanaka, and N. Kawakami, Physica E-Low-Dimensional Systems and Nanostructures , 265(2007).[79] L. Li, B.-B. Zheng, W.-Q. Chen, H. Chen, H.-G. Luo,and F.-C. Zhang, Physical Review B , 245135 (2014).[80] G. Kirˇsanskas, M. Goldstein, K. Flensberg, L. I. Glaz-man, and J. Paaske, Physical Review B , 235422(2015).[81] R. ˇZitko, J. S. Lim, R. L´opez, and R. Aguado, PhysicalReview B , 045441 (2015).[82] W. V. van Gerven Oei, D. Tanaskovi´c, and R. ˇZitko,Physical Review B , 085115 (2017).[83] J. S. Lim and R. L´opez, Physical Review B , 245427(2020).[84] A. Kadlecov´a, M. ˇZonda, and T. Novotn´y, Physical Re-view B , 195114 (2017).[85] C. D. Dominicis and P. C. Martin, Journal of Mathemat-ical Physics , 14 (1964).[86] C. D. Dominicis and P. C. Martin, Journal of Mathemat-ical Physics , 31 (1964).[87] N. E. Bickers and S. R. White, Physical Review B ,8044 (1991).[88] G. Rohringer, H. Hafermann, A. Toschi, A. Katanin,A. Antipov, M. Katsnelson, A. Lichtenstein, A. Rubtsov,and K. Held, Reviews of Modern Physics90