Spectral and transport properties of a half-filled Anderson impurity coupled to phase-biased superconducting and metallic leads
SSpectral and transport properties of a half-filled Anderson impurity coupled tophase-biased superconducting and metallic leads
Peter Zalom, ∗ Vladislav Pokorn´y, † and Tom´aˇs Novotn´y ‡ Institute of Physics, Czech Academy of Sciences,Na Slovance 2, CZ-18221 Praha 8, Czech Republic Department of Condensed Matter Physics, Faculty of Mathematics and Physics,Charles University, Ke Karlovu 5, CZ-12116 Praha 2, Czech Republic (Dated: June 12, 2020)We derive and apply a general scheme for mapping a setup consisting of a half-filled single levelquantum dot coupled to one normal metallic and two superconducting phase-biased leads onto anordinary half-filled single impurity Anderson model with single modified tunneling density of states.The theory allows for the otherwise unfeasible application of the standard numerical renormalizationgroup and enables to obtain phase-dependent local spectral properties as well as phase-dependentinduced pairing and Josephson current. The resulting transport properties match well with the nu-merically exact continuous-time hybridization-expansion quantum Monte Carlo. For weakly couplednormal electrode, the spectral properties can be interpreted in terms of normal-electrode-broadenedAndreev bound states with phase-dependent position analogous to the superconducting Andersonmodel, which coexist in the π -like phase with a Kondo peak whose phase-dependent Kondo temper-ature follows the form conjectured previously in Doma´nski et al., Phys. Rev. B , 045104 (2017). I. INTRODUCTION
Gradual advance in experimental techniques over thepast decades allowed to study electronic transport in in-creasingly sophisticated nanoscale systems with various,competing correlations. A prototype experiment typi-cally includes a strongly interacting mesoscopic systemattached to a reservoir with well defined properties. Intheory, the mesoscopic system is frequently described interms of one or multiple quantum dots (QDs) in theCoulomb blockade regime while the reservoir consists ofnormal metallic and/or superconducting leads. Exper-imental realizations of such QDs include, for example,carbon nanotubes [1–9] or semiconductor nanowires [10–13].The case of normal metallic electrodes attached to oneQD can be modeled microscopically by the single impu-rity Anderson model (SIAM) which is one of the mostunderstood models in the many-body physics [14]. Here,free conduction electrons of the reservoir can completelyor partially screen the magnetic doublet of QD depend-ing on the parameters under the study. When the screen-ing is effective, an emergent Kondo singlet becomes theground state of the system [14]. The Kondo singlet is acoherent many-body state composed of a screening cloudand QD states. It is specific due to its logarithmic energyscaling which can be fully understood only by applyingrenormalization group (RG) techniques [15–17] or at leasteffective renormalization schemes [18, 19].Once superconducting correlations are considered inthe reservoir, the electron transport is altered by the An-dreev scattering on the interface between the leads and ∗ [email protected] † [email protected] ‡ [email protected]ff.cuni.cz QD [20–24], but it is still well described theoretically interms of the Anderson impurity model with supercon-ducting leads [25]. For purely superconducting reservoirs,both the theoretical [26, 27] and experimental [28] under-standing is fairly complete. In particular, a sufficientlylarge gap depopulates electrons around the Fermi energyto such an extent that the screening cloud around theimpurity becomes disrupted. The ground state of thesystem changes then from a singlet (effective screeningat small sized gaps) to a doublet, which is an exampleof a impurity quantum phase transition (QPT). This so-called 0- π transition is accompanied by the reversal of thesupercurrent, which is positive in the singlet and negativein the doublet phase [5, 10, 29]. At the transition, onealso observes the crossing of the Andreev bound states(ABSs) at the Fermi energy.Hybrid systems incorporating simultaneously normalas well as superconducting reservoirs, lead to even moreintricate interplay of quantum correlation effects, wherethe understanding is limited both theoretically and ex-perimentally [26]. In a simplest realization, one metal-lic and one superconducting lead have been studied ex-perimentally in N-QD-S heterostructures [30, 31]. Fromthe theoretical perspective, such problem is of only two-channel nature and thus well tractable by the numericalrenormalization group (NRG) which offers unbiased in-sights and thus complements the purely numerical quan-tum Monte Carlo (QMC) simulations with reliable spec-tral properties [32]. However, having just one supercon-ducting lead does not allow superconducting phase dif-ference across the QD. Consequently, such systems lackany supercurrent flow and no interplay of Kondo andJosephson effects takes place. Thus, the more interestingscenario includes one normal and two phase-biased su-perconducting leads. The resulting three-terminal struc-ture is, however, beyond the reach of standard NRGtechniques since the corresponding discretized reservoir a r X i v : . [ c ond - m a t . m e s - h a ll ] J un model corresponds to three spin dependent and mutu-ally interconnected hopping chains. The standard NRGscheme has thus been so far employed only to two channelproblems [33, 34].In the standard, computationally intractable approachto NRG, one would first discretize the bath of the threeterminals. The resulting semi-infinite hopping chainwould be then transformed via the Bogolyubov-Valatintransformation at each chain site. Subsequently, in thesecond step, particle-hole transformations on odd andeven sites separately need to be carried out [33, 34]. Theresulting Wilson chain would, however, consist of threemutually interconnected spin-polarized chains which isbeyond the present computational power [32]. To cir-cumvent the problem, one may apply a general proce-dure of Ref. [35] or introduce suitable unitary transfor-mations to diagonalize the problem [36, 37]. The secondapproach has already been successfully applied to thehybrid normal-superconductor reservoir problem in thelimit of infinite superconducting gap. Although, the au-thors of Refs. [36, 37] make explicit reference to Wilsonchains corresponding to discretized versions of the modelunder the study, as shown in the present paper, it is pos-sible to limit the transformations in the case of half-fillingto just the local electrons of QD and map the infinite-gapmodel onto an ordinary asymmetric SIAM directly.Generalizing the approach of Refs. [36, 37], we are ableto treat the present general three-terminal problem atthe half-filling with finite superconducting gap and mapit onto a single impurity Anderson model of fermionswith tunneling density of states (DOS) in the reservoirthat corresponds to the standard one-channel-lead casetractable by NRG in the scheme of Ref. [38]. Since it isbelieved that such an approach is not feasible [39], wepresent the details of the transformation in Sec. II wherealso a detailed microscopic formulation of the problem isformulated. In Sec. III, we proceed to the ∆ → ∞ casetreated previously in Refs. [36, 37] and show that theirapproach is in complete equivalence with the one pre-sented here once half-filled case is considered. However,as opposed to Refs. [36, 37], no reference to NRG dis-cretization scheme is required at the microscopic level.Finally in Sec. IV, the general three-terminal problemwith finite superconducting gap is solved at the half-filling. To this end, the mapping of the finite-gap problemonto a single-channel SIAM with altered tunneling DOSis performed. Subsequently, standard NRG approach ofRef. [38] is employed utilizing the NRG Ljubljana code[40]. Using the corresponding backwards transforma-tions, all spectral and transport properties of the orig-inal three-terminal setup are then determined. The mostimportant conclusions are finally summarized in Sec. V.Technical calculations regarding the transformation ofthe interaction term in Sec. II and the effect of the fi-nite band width are discussed in the Appendices A andB, respectively. The comparison with QMC is shown inthe Appendix C. II. MAPPING ONTO SIAM-LIKE MODELSA. Microscopic formulation
The hybrid three-terminal setup consists of a meso-scopic system modeled as a usual Anderson magneticimpurity connected to one normal metallic and two su-perconducting electrodes. The superconducting elec-trodes follow the Bardeen-Cooper-Schriefer (BCS) the-ory of conventional s -wave superconductivity. One of theBCS leads is then referred to as the left ( L ) and the otheras the right ( R ) lead, see Fig. 1. They are considered tohave arbitrary BCS phase parameters ϕ L and ϕ R but thesame gap parameter ∆ L = ∆ R ≡ ∆. The total Hamilto-nian of the system is then the sum of the dot Hamiltonian H d , the Hamiltonian of the normal lead H N , two BCSHamiltonians for superconducting leads H L and H R , andthree tunneling Hamiltonians H T,α with α ∈ { N, L, R } which connect each lead separately to the dot. The con-stituent Hamiltonians read as H d = (cid:88) σ ε d d † σ d σ + U d †↑ d ↑ d †↓ d ↓ , (1) H α = (cid:88) k σ ε k α c † α k σ c α k σ − ∆ α (cid:88) k (cid:16) e iϕ α c † α k ↑ c † α − k ↓ + H.c. (cid:17) , (2) H T,α = (cid:88) α k σ (cid:16) V ∗ α k c † α k σ d σ + V α k d † σ c α k σ (cid:17) , (3)where c † α k σ creates an electron of spin σ and quasi-momentum k in the lead α while c α k σ annihilates it. Inanalogy, d † σ creates a dot electron of spin σ while d σ an-nihilates it. The QD is characterized by the Coulombrepulsion U and the level energy ε d which in the mostgeneral case is arbitrary but we will later concentrateonly at ε d = − U/
2. The QD hybridizes with the leadsvia V α k and the gap parameter vanishes in the normallead, thus ∆ N = 0.To allow for an unbiased treatment of the Kondo andsuperconducting correlations on the same footing, RGtechniques should be applied. However, in the basis ofthe d and c electrons, superconducting leads mix togetherthe spin up and spin down electrons of the bath whichrequires to solve a general problem of three mutuallyinteracting spin-polarized-bath channels. Applying theoriginal numerical RG scheme [15] is therefore extremelydemanding from the point of view of the computationaltime [32]. In general, multichannel problems pose a gen-uine challenge when applying NRG and have been tackledonly quite recently within the interleaved formulation ofthe generalized Wilson chain [41].Nevertheless, in the present case Hamiltonians of allthree leads are quadratic because lead electrons are con-sidered to be mutually non-interacting. Consequently,the lead electrons can be integrated out to obtain a gen-uine one channel impurity problem. To this end, it is ad-vantageous to reformulate Eqs. (1)-(3) using the Nambuformalism. B. Hamiltonian in the Nambu basis
Nambu formalism represents a convenient way of re-arranging Hamiltonians involving BCS superconductiv-ity in a way where the lead electrons are all treated onthe same footing. To this end, let us combine the spinup and spin down component of the corresponding fieldsdescribing the c electrons into spinors C † α k = (cid:16) c † α k ↑ , c α − k ↓ (cid:17) (4)with α ∈ { N, L, R } , while the spinor D of the dot elec-trons d is constructed in complete analogy as D † = (cid:16) d †↑ , d ↓ (cid:17) . (5)The Hamiltonians (2) and (3) become then H α = (cid:88) k C † α k E α k C α k , (6) H T,α = (cid:88) k D † V α k C α k + C † α k V α k D (7)with E α k = ∆ α C α σ x − ∆ α S α σ y + (cid:15) k σ z , (8) V α k = V α k σ z , (9)where V α k has been considered to be a real function only, σ i , i ∈ { x, y, x } , are Pauli matrices while C α ≡ cos ϕ α , S α ≡ sin ϕ α . The blackboard bold typeface is from nowon used to distinguish matrices from scalars.Before applying the Nambu formalism to the Hamilto-nian (1), let us first separate it into a quadratic part H d, = (cid:88) σ ε d d † σ d σ + U D † ( σ x + σ z ) D = D † E d D (10)with E d = U σ x + (cid:18) U (cid:15) d (cid:19) σ z (11)and a mixed quadratic and quartic interaction part H U = U d †↑ d ↑ d †↓ d ↓ − U D † ( σ x + σ z ) D. (12)The advantage of this non-standard partitioning will bediscussed in Sec. II C.Taking together, in the Nambu formalism the non-interacting quadratic part of the present problemspanned by the D and C spinors reads as H = H d, + (cid:88) α ( H α + H T,α ) (13)while the modified interaction part H U is given byEq. (12). C. Bogolyubov-Valatin transformations in thespace of local electrons
Let us now study the effect of unitary transformations T and ˜ T α onto the spinors D and C α respectively. Weintroduce the spinors W and ˜ C via T D ≡ W, D † T † ≡ W † , (14)˜ T α C α k ≡ ˜ C α k , C † α k ˜ T † α ≡ ˜ C † α k . (15)At this point, we consider no other constraints onthe transformations T , ˜ T α except of unitarity so thatthe many-body energy spectra of the problem remainthe same in both spinor bases D and W . On theother hand, such transformations may crucially affect theform of the one-particle operators in the correspondingnon-interacting Hamiltonians thus allowing for compu-tationally more suitable non-interacting Green functionsand/or self-energy contributions from the integrable de-grees of freedom for the problem under the study.The effect of the transformations T and ˜ T α on thequadratic part H d, reads as H d, = W † (cid:0) TE d T † (cid:1) W , (16)while the tunneling Hamiltonians change as H T,α = (cid:88) k W † (cid:16) TV α k ˜ T † α (cid:17) ˜ C α k + (cid:88) k ˜ C † α k (cid:16) ˜ T α V α k T † (cid:17) W. (17)The kinetic Hamiltonians are affected only by thetransformations ˜ T α as they involve no operators of thelocal electrons: H α = (cid:88) k ˜ C † α k (cid:16) ˜ T α E α k ˜ T † α (cid:17) ˜ C α k . (18)Since the non-interacting ( U = 0) Hamiltonian isquadratic we can easily obtain the retarded Green func-tion G W ( ω + ) which corresponds to the W spinors. Em-ploying the equation of motion technique for the Greenfunctions in an exact analogy to Ref. [42], we introducean infinite-dimensional vectorΨ † = ( W † , ˜ C † N k , ˜ C † L k , ˜ C † R k ) , (19)where Ψ is its Hermitian conjugate and the spinors ˜ C † α k are understood to be repeated in Ψ † for all possible quasi-momenta of lead electrons. This allows us to rearrangethe non-interacting Hamiltonian as H = Ψ † E W Ψ , (20)with E W = TE d T † TV N k ˜ T † N TV L k ˜ T † L TV R k ˜ T † R ˜ T N V L k T † ˜ T N E N k ˜ T † N T L V L k T † T L E L k ˜ T † L T R V R k T † T R E R k ˜ T † R , (21) FIG. 1. ( a ) Scheme of Y -shaped three-terminal junction where a QD (black disc) couples via the hybridization strength Γ N toone normal electrode (red pointed teardrop) and left (L) and the right (R) superconducting electrode (blue pointed teardrop)which are phase-biased by ϕ = ϕ L − ϕ R . The BCS leads hybridize with quantum dot with strengths Γ L and Γ R , respectively,and are considered here as having the same BCS gap parameter ∆ L = ∆ R ≡ ∆. ( b ) Equivalent scheme consisting of one-terminal reservoir containing Bogolyubov-like quasiparticles (purple pointed teardrop) which are hybridized with the QD viaa structured hybridization function Γ W ( ω ). In Sec. II C we show that by employing Γ W ( ω ) according to Eq. (30) the schemes( a ) and ( b ) may be in terms of spectral properties mapped onto each other at the half-filling. To obtain transport propertiescorresponding to the Y shape geometry of panel ( a ) transformations according to Sec. IV C are required. where the upper index W was introduced to clearly dis-tinguish the underlying spinor basis W for the formula-tion of the infinite-dimensional matrix E W .In general, the non-interacting problem can be solvedby finding the retarded Green function G ( ω + ) where ω + ≡ ω + iη while ω is a real frequency and η is aninfinitesimally small positive number. To this end, stan-dard equation of motion technique formulated in the ma-trix form requires one to solve the resolvent equation G ( ω + ) = ( ω + − H ) − with being the unit matrix.However, for the present problem we only need to ob-tain the solution for the local Green function of the dotelectrons which corresponds the left upper 2 × G W ( ω + ) in the spinor basis W directly as G W ( ω + ) = (cid:0) ω + − TE d T † − TΣ D T † (cid:1) − (22)with Σ D ( ω + ) = (cid:88) α k V α k ˜ T † α (cid:16) ω + − ˜ T α E α k ˜ T † α (cid:17) − ˜ T α V α k = (cid:88) α k V α k (cid:0) ω + − E α k (cid:1) − V α k , (23)which thus represents the self-energy contribution fromthe leads expressed with respect to the spinor basis D (as denoted by the upper index D ). Moreover, one mayalso obtain the self-energy contribution Σ W as Σ W ( ω + ) = TΣ D ( ω + ) T † , (24)which not only defines Σ W ( ω + ), but also gives us thetransformation rule to easily interchange the spinor bases D and W when required. By exploiting the unitarity of T , we may also extract an analogous transformation rulefor the non-interacting retarded Green functions G W ( ω + ) = T (cid:0) ω + − E d − Σ D (cid:1) − T † = TG D ( ω + ) T † (25) which yields the transformation rule as well as the def-inition of the non-interacting ( U = 0) retarded Greenfunction with respect to the spinor basis D . Clearly, inboth bases the effect of the leads is fully integrated outand only enters the G W ( ω + ) and G D ( ω + ) via the corre-sponding self-energy contributions Σ W ( ω + ) and Σ D ( ω + ),respectively. Green function as well as the self-energycontributions in different bases relate to each other viathe local dot transformation T since the transformations˜ T α are canceled out in Eqs. (24) and (25).It is now our aim to construct a suitable transformation T , so that the self-energy contribution Σ W is diagonal.To this end, we first perform all summations in Eq. (23),which is a quite straightforward with details given in theAppendix B. The resulting expression for Σ D ( ω + ) hasthe following matrix structure: Σ D ( ω + ) = Σ DN ( ω ) + Σ DA ( ω ) σ x , (26)where Σ DN ( ω ) and Σ DA ( ω ) are functions of frequency witha for now unimportant form which is given in the Ap-pendix B. We insert now Σ D ( ω + ) back into Eq. (25) toobtain the matrix structure of (cid:0) G D (cid:1) − . We first concen-trate exclusively on the half-filled case where (cid:0) G D (cid:1) − hasa much simpler structure. Afterwards, the more general ε d (cid:54) = − U/ (cid:0) G D (cid:1) − is onlyproportional to the unit matrix which remains unal-tered under unitary transformation. On the other hand,the off-diagonal parts of Σ D and E d are both proportionalto σ x and may thus be simultaneously diagonalized byenforcing the condition T σ x T † = ± σ z onto the transfor-mation T . The resulting non-interacting Green function G W is then also diagonal as intended. Each sign option ofthe condition T σ x T † = ± σ z , is solved by two linearly in-dependent transformations which we denote T ± and T ± with the given superscript indicating the sign of σ z in thecondition. Explicitly, we obtain T ± = √ ( σ x ± σ z ) , T ± = √ ( ± iσ y ) . (27)Note, that all transformations fulfill T † = T − = T andare of Bogolyubov-Valatin type but unlike in BCS super-conductors they do not involve momentum or frequencydependence.Now it also becomes obvious, why the half-filled casewas selected first. Generalizing the previous approach to ε d (cid:54) = − U/ (cid:0) G D (cid:1) − con-siderably by adding an extra diagonal term proportionalto ( ε d + U/ σ z . Diagonalizing (cid:0) G D (cid:1) − by T would nowrequire to simultaneously fulfill not only T σ x T † = ± σ z but also T σ z T † = ± σ z or T σ z T † = ± . Neither ofthe two combined conditions is however solvable. Con-sequently, outside of the half-filled case there exists nounitary transformation to diagonalize (cid:0) G D (cid:1) − . For ex-ample, applying T in the form of Eq. (27) away from thehalf-filling rotates the superconducting terms onto the di-agonal elements of (cid:0) G W (cid:1) − which is however traded offfor rotating the originally diagonal terms proportionalto ( ε d + U/ σ z into off-diagonal terms proportional to( ε d + U/ σ x . From now on, we therefore concentrate exclusively on the half-filled case .Although, for the non-interacting ( U = 0) half-filledcase the corresponding Green function can be diagonal-ized by transformations (27), in the end, in the full in-teracting case the action of transformations (27) on theinteraction part H U needs to be considered. As shownin the Appendix A, using the interaction term mixed ofquadratic and quartic terms according to Eq. (12) allowsto obtain the Hubbard interaction term when the trans-formation T − is used. Explicitly: H U = U w †↑ w ↑ w †↓ w ↓ , (28)The remaining options for T produce all an extraquadratic terms in H U which is proportional to σ x , thusspoiling the diagonal form of G W ( ω + ). We thus selectthe transformation T − in what follows and denote it as T ≡ T − . Nevertheless, we stress that the remainingchoices would also be possible when different partition-ings of the full Hamiltonian are considered. Moreover,the transformation T +2 fully corresponds to the transfor-mation used in the standard NRG treatments of mag-netic impurities coupled to superconducting reservoirs(see Refs. [33, 34] for more detail) while T +1 relates tothe transformations applied in the Refs. [36, 37] to solvethe ∆ → ∞ limit of the present model. However, in all ofthe aforementioned references, even at the half-filling thelead transformations corresponding to ˜ T α are explicitlyperformed within the applied NRG algorithms.Since the resulting G W ( ω + ) is proportional to σ z , onemay actually drop the Nambu formalism and employ di-rectly the w †↑ and w †↓ fields which constitute the W spinor. In such a case, we obtain: G W ↑ ( ω + ) = G W ↓ ( ω + ) = 1 ω + U/ − Σ W ( ω ) , (29)and the independence of Green functions of the spinindex of the field w becomes explicit. Assuming themomentum-independent hybridizations V α k ≡ V α , thecalculation of Σ W ( ω ) can be performed for arbitraryband-width as shown in the Appendix B. In the limitof infinitely wide band, it takes the following form Σ W ( ω + ) = − i Γ N − Γ S (cid:20) i sgn( ω ) √ ω − ∆ Θ( ω − ∆ )+ Θ(∆ − ω ) √ ∆ − ω (cid:21) (cid:16) ω + ∆ cos ϕ (cid:17) , (30)where Θ is the Heaviside step function and Γ N ≡ πρ V N and Γ S ≡ πρ V S with ρ being the constant DOS of theband and V S = V L = V R . As shown in Ref. [43], the V L (cid:54) = V R case can be directly obtained form the symmet-ric scenario of V L = V R introduced explicitly above. Theimaginary part of Eq. (30) is traditionally referred to asthe hybridization function, see also Sec. IV A, while thereal part is connected to the imaginary one via Kramers-Kronig relations.Taking together, the unitary transformation T al-lows to map the original non-diagonal model expressedvia D spinor onto a model described in terms ofBogolyubov-type quasiparticles w σ coupled to a singlenormal lead which has an altered tunneling DOS due tothe frequency-dependent self-energy Σ W ( ω ). Moreover,the Bogolyubov quasiparticles w σ interact locally via theordinary Hubbard interaction term. This allows us to re-define the three-terminal setup as one-channel-lead prob-lem similar to the ordinary SIAM and apply NRG in astraightforward way as described in Ref. [38]. This way,all spectral properties in the spinor basis W can be ob-tained.Although the normal spectral function A WσN ( ω ) and thetrivially vanishing anomalous spectral function A WσA ( ω )in the basis of the w fields can be determined, in theend, we need to transform back to the original basis ofthe d electrons. Therefore, one needs to relate the Greenfunctions and their corresponding spectral functions be-tween both bases. Since the unitary transformation T mixes the original d σ and d † σ fields only in a linear way,the corresponding normal Green function G DσN and thecorresponding anomalous Green function G DσA in the ba-sis of the d electrons are given by G DσN ( ω + ) = 12 (cid:2) G WσN ( ω + ) − G WσN ( − ω + ) (cid:3) , (31) G DσA ( ω + ) = 12 (cid:2) G WσN ( ω + ) + G WσN ( − ω + ) (cid:3) , (32)where − ω + ≡ − ω − iη . Recall that G Wσ ( ω + ) = G W ∗ σ ( ω − )with ω − = ω − iη and the imaginary part of the Greenfunction equals the spectral function up to a multiplica-tive prefactor − /π . Therefore, we can directly constructthe normal spectral function A DσN ( ω ) and the anomalousspectral function A DσA ( ω ) in the d basis in the d basisusing A DσN ( ω ) = 12 (cid:2) A WσN ( ω ) + A WσN ( − ω ) (cid:3) , (33) A DσA ( ω ) = 12 (cid:2) A WσN ( ω ) − A WσN ( − ω ) (cid:3) . (34)We will often refer to the backwards transformationsof the normal spectral and anomalous function as sym-metrization and antisymmetrization, respectively. III. ∆ → ∞ CASE
The ∆ → ∞ case of the present model is alreadyquite well understood with a general phase-dependentNRG solution given in Refs. [32, 36, 37]. Interestingly, inRefs. [36, 37] the authors have mapped the problem ontothe standard asymmetric SIAM with constant tunnelingDOS by employing similar Bogolyubov-Valatin transfor-mations as discussed above. However, instead applyingthe transformations on the level of the Hamiltonian, theytransform the Wilson chain corresponding to the dis-cretized version of the ∆ → ∞ model. Such a procedurehowever requires not only the T transformation of the dotelectrons but also an analog of the ˜ T α transformation ofthe leads which is applied to each site of the Wilson chainrepresenting the lead electron contributions. Unnoticedby the authors of [36, 37], at the half-filling only thetransformation T is, however, explicitly required whichwe will now show. To this end, we briefly review themost crucial points made in Refs. [36, 37]. The ∆ → ∞ Hamiltonian simplifies to H ∆ →∞ = H d, + H U + H N + H T,N − ∆ d ( ϕ ) (cid:16) d †↑ d †↓ + d ↓ d ↑ (cid:17) , (35)where ∆ d ( ϕ ) ≡ Γ S cos( ϕ/
2) and H N , H d, , H T,N fol-low our previous notations and the authors have appliedtheir approach to arbitrary values of ε d . Obviously, theeffective Hamiltonian is just of the one-channel-lead na-ture. However, the BCS effects are manifested by non-zero off-diagonal terms which cannot be treated directlyby the original NRG algorithm [15]. In general, whenoff-diagonal terms emerge, as for example in the super-conducting Anderson model (SCIAM), the constructionof the Wilson chain requires to implement Bogolyubov-Valatin transformation to each site of the Wilson chainwhich is then followed by particle-hole transformationsto each of the site except of the initial one. Such a se-ries of transformations ensures then the diagonality ofthe emerging Wilson chain [33]. Due to the presence ofthe off-diagonal terms, a similar approach is then alsoapplied for the ∆ → ∞ case [36, 37]. In detail, the log-arithmically discretized version of ∆ → ∞ model givesrise to a semi-infinite Wilson hopping chain with the firstnode populated by the local d σ electrons ( σ ∈ {↑ , ↓} ) of the QD and remaining sites labeled by i ∈ N and popu-lated by fermions c iσ which represent the bath degrees offreedom. In Refs. [36, 37], each site of the Wilson chainis then rotated using unitary transformation O Θ given by O Θ = (cid:32) cos Θ2 − sin Θ2 sin Θ2 cos Θ2 (cid:33) , (36)with cos Θ2 = (cid:118)(cid:117)(cid:117)(cid:116) (cid:32) ε d (cid:112) ˜ ε d + ∆ d (cid:33) , (37)sin Θ2 = (cid:118)(cid:117)(cid:117)(cid:116) (cid:32) − ˜ ε d (cid:112) ˜ ε d + ∆ d (cid:33) , (38)which acts on the Nambu spinors of the given site of theWilson chain, δ = (cid:112) ˜ ε d + ∆ d and ˜ ε d = ε d + U/
2. Thefactor δ is later shown to correspond to the particle-holeasymmetry factor of the ordinary SIAM with constanttunneling DOS. For example, the D spinor of the firstsite of the Wilson chain transforms as O Θ D = O Θ (cid:32) d ↑ d †↓ (cid:33) = (cid:32) w ↑ w †↓ (cid:33) ≡ W, (39)where W follows the notation of the previous section.Afterwards, all sites corresponding to the lead electrons(indexed by i ) are subjected to additional particle-holetransformations which are required to transform away theoff-diagonal BCS terms in the Wilson hoping chain. Onlythen a semi-infinite chain feasible for standard NRG al-gorithms is constructed. The resulting coefficients of thesemi-infinite Wilson chain are then noticed to be identicalto those of the ordinary asymmetric SIAM and becauseof the Hausholder transformation [36, 37] equivalence ofthe ∆ → ∞ model to the asymmetric SIAM is finallynoticed even at the microscopic level [36, 37]. We stressthat the findings in Refs. [36, 37] hold at arbitrary filling.However, once only the half-filling is considered, theunitary transformation O Θ may be restricted to act onlyin the space of local electrons while the transformationson the additional sites turn out to be merely a technicaltool required if NRG-based discretizations are introducedbefore lead electrons are transformed accordingly. Theproof is completely analogous to that of the finite-gaphalf-filled problem discussed above. We thus skip thedetails and discuss only the crucial steps.We first notice, that the half-filled case in Refs. [36, 37],i.e. ˜ ε d = ε d + U/ π/ O Θ transformation becomes the T +2 transformationfound previously for the general finite-gap model. Fur-thermore, the three-terminal structure of the leads en-ters the QD equation of motion only via the correspond-ing self-energy contribution, which for ∆ → ∞ is purelyimaginary and proportional to the 2 × d basis of the local electrons. Thus, its form isinvariant under any unitary transformation. Explicitly:Σ DN ( ω ) = Σ WN ( ω ) = − i Γ N . (40)On the other hand, the off-diagonal parts of the dotHamiltonian turn out to have a more intricate form. Per-forming then the same partitioning of the Hamiltonian asin Eqs. (10) and (12), we may apply the transformation O π/ and transform the non-interacting part of the dotHamiltonian into a diagonal form H d, = W † (cid:18) δ − U − δ + U (cid:19) W. (41)Applying then the transformation O π/ to the interac-tion term gives us subsequently the Hubbard repulsionterm typical for SIAM. Consequently, at the half-fillingtransformation (39) maps the ∆ → ∞ model onto a gen-eral asymmetric SIAM with the particle-hole asymmetryrelated to the phase difference ϕ which remained unno-ticed by the authors of Refs. [36, 37]. Applying the NRGalgorithm according to Ref. [38] to the diagonal formula-tion of the ∆ → ∞ gives Wilson chain coefficients whichare the same as those obtained by Refs. [36, 37]. Theresulting spectral and transport properties are thereforethe same and we review these only briefly below.Before proceeding further, let us first clarify the role ofthe QD filling in the transformation procedures. Clearly,the transformations T derived in this paper have beendesigned to work exclusively at the half-filling and fail todiagonalize the finite ∆ as well as the ∆ → ∞ case awayfrom the half-filling. On the other hand, the transfor-mation O Θ given by Eq. (39) and accompanied by cor-responding particle-hole transformations of the Wilsonchain diagonalizes only the ∆ → ∞ model but at an ar-bitrary filling. Nevertheless, at the half-filling O Θ is ac-tually exactly the T +2 transformation which under properchoice of partitioning of the Hamiltonian allows to mapeven the finite-gap problem at the half-filling onto a sim-ple SIAM. However, even for the ∆ → ∞ model the O Θ nor any other unitary transformation of the dot electronsis sufficient to solely diagonalize the problem away fromthe half-filling. This regarding, it is important to noticethat the off-diagonal self-energy contribution of the hy-brid reservoir of normal and superconducting electronsis proportional to σ x in the basis of the d fields. Suchcontributions transform under O Θ as O † Θ (cid:18) (cid:19) O Θ = (cid:18) sin Θ cos Θcos Θ − sin Θ (cid:19) (42)and become diagonal only at the half-filling in completeanalogy to the finite-gap model discussed above.Since both transformations T and O Θ accomplish thesame for the ∆ → ∞ model at the half filling, we willbriefly discuss the results and interpretations obtainedvia the transformation T while all details on transportproperties can be found in Refs. [36, 37]. The transfor-mation T maps the H ∆ →∞ Hamiltonian onto that of an ordinary particle-hole asymmetric SIAM with the asym-metry parameter δ = Γ S cos( ϕ/ A WN ( ω ) of the asymmetric SIAM with the asymmetry pa-rameter δ matching the parameters of the ∆ → ∞ case.Then, symmetrization procedure (33) is applied to obtainthe normal spectral function A DN ( ω ) of the three-terminal∆ → ∞ model in the original basis of the d fields. Forthe anomalous functions (subscript A ) one proceeds ana-logically by antisymmetrization (34).In Fig. 2, we present the spectral functions of theasymmetric SIAM for selected values of δ (see the toprow of panels) which have been obtained using the opensource NRG Ljubljana code [40] in the one-channel modewith intertwined z -discretization [45] with z = n/ n ∈ { , . . . } . The middle row of panels shows the nor-mal spectral functions of the ∆ → ∞ model obtained bythe symmetrization of those from the top row. In the bot-tom row, the anomalous spectral functions of the ∆ → ∞ model are presented. Panels in the same columns corre-spond to each other via the relation δ = Γ S cos( ϕ/
2) andare ordered from left to right with the increasing phasedifference ϕ and thus decreasing asymmetry parameter δ of the underlying SIAM. The particle-hole symmetriccase of the effective SIAM is then realized at ϕ = π .The mapping T makes thus obvious that the competitionbetween the Kondo and superconducting correlations, asexplicitly seen in the basis of the d fields, simplifies intothe ordinary Kondo physics away from the half-fillingonce the basis of the w fields is used.Let us now discuss the consequences in more detail.In the basis of the w fields the particle-hole symmetryat ϕ = π manifests itself via the appearance of an ordi-nary Kondo resonance at the Fermi energy in the spec-tral function A WN ( ω ). Additionally, two satellite Hubbardpeaks emerge at approximately ± U/ ϕ = π case, Josephson and other superconducting effectsdo not disrupt Kondo correlations, which is a trivial re-sult already visible from the effective Hamiltonian (35).However, decreasing the phase difference ϕ , while keep-ing parameters U , ∆, Γ N and Γ S constant, drives theeffective system away from its particle-hole symmetricpoint with the asymmetry factor δ of the underlyingasymmetric SIAM increasing due to δ = Γ S cos( ϕ/ ϕ (from right to left), we observethe central peak corresponding to the Kondo resonanceat ϕ = π shifting gradually away form the Fermi en-ergy which is accompanied by its gradual broadening andmild decrease of its height [46]. When the particle-holeasymmetry in the basis of the w fields is relatively small,corresponding to δ = 2 . . . . . . . . . − − − −− − − − − − − − − − π Γ N A W N ( ω ) δ = 5 . δ = 4 . δ = 2 . δ = 0 . π Γ N A D N ( ω ) ϕ = 0 ϕ = π/ ϕ = 2 π/ ϕ = π π Γ N A D A ( ω ) ω/ Γ N ϕ = 0 ω/ Γ N ϕ = π/ ω/ Γ N ϕ = 2 π/ ω/ Γ N ϕ = π FIG. 2. The top row of panels shows the normal spectral functions for the asymmetric SIAM with U/ Γ N = 5 in the directionof decreasing asymmetry parameter δ . According to Eqs. (40) and (41), the top row spectral functions also correspond tothe spectral functions A WN ( ω ) of the ∆ → ∞ model expressed in the basis of the w fields when U and Γ N are the same and δ = Γ s cos( ϕ/ A DN ( ω ) of the ∆ → ∞ model in the basis ofthe d fields obtained by symmetrization (33) of A WN ( ω ) corresponding to the top row. The bottom row shows the anomalousspectral function A DA ( ω ) of the ∆ → ∞ model in the basis of the d fields obtained by antisymmetrization (34) of the top row. overcome its broadening. Consequently, performing sym-metrization operation (33) to obtain A DN ( ω ) makes thecentral peak broader but still singly-peaked at the Fermienergy. However, decreasing the angle ϕ further increasesthe asymmetry parameter of the underlying asymmetricSIAM even more, which results in a very strong decen-tralization of the Kondo-like peak. The broadening isthen not sufficient to overcome the shift and symmetriza-tion (33) leads to a split central peak in the d basis (seethe δ = 4 .
33 and δ = 5 . ϕ ∗ ≈ π/ w fields are suppressed which translates into the d basisas a split central peak in the normal spectral function.Let us also stress out, that the phase-induced splitting ofthe central peak may be completely avoided in parameterregimes with very small hybridization Γ S which is in com-plete analogy to SCIAM. The crucial point is to keep thevalue of Γ S so small that for given interaction strength U the asymmetry parameter δ is sufficiently small to avoidsignificant disruption of the Kondo peak.Let us now follow the number and positions of all off-center peaks visible in the corresponding spectra of the∆ → ∞ model with increasing the angle ϕ as seen in thebasis of the d fields (the middle row of panels in Fig. 2 seen form left to right). At ϕ = 0 (the first column), wesee a pair of off-center peaks at ω ≈ ± N and a secondpair of highly suppressed peaks at ω ≈ ± N . The twoinner peaks have already been connected to the remnantsof the Kondo resonance as seen in the basis of the w fields.The outer pair of peaks relates to the suppressed chargeexcitations of the underlying asymmetric SIAM in thebasis of the w fields. Increasing the angle ϕ causes the δ parameter of the underlying asymmetric SIAM to de-crease. Consequently, the two inner peaks move towardseach other until they finally merge together for ϕ ∗ ≈ π/ δ . There-fore, at around ϕ ∗ , when the Kondo-like central peakforms, two pairs of symmetrically placed off-center peaksbecome recognizable (at ϕ = 2 π/ ± U/ ϕ = π . Here, they corre-spond to the ordinary Hubbard peaks of the symmetricSIAM.For the ∆ → ∞ model, we may thus encountertwo distinctive spectral function shapes which cross oversmoothly into each other. For ϕ < ϕ ∗ we observe one pairof symmetrically positioned peaks that move towards theFermi energy with increasing ϕ → ϕ ∗ while the otherpair of outer peaks is mostly unnoticeable. At ϕ ≈ ϕ ∗ the inner peaks cross each other and continue to movetowards ± U/ ϕ .The pair of the outer peaks becomes more pronouncedat around ϕ ∗ where they also start to move noticeablytowards the Fermi energy. At ϕ = π , they finally mergetogether with the inner pair of peaks at approximately ± U/
2. Taking together, the movement of the variouspeaks is highly reminiscent of the behavior of the ABSstates of the SCIAM. For ϕ < ϕ ∗ the spectral function re-sembles the 0 phase while for ϕ > ϕ ∗ the behavior of thepeaks shows analogy to the π phase of the SCIAM. Nev-ertheless, in both cases the presence of the normal leadcauses the ∆ → ∞ model to acquire a singlet many-bodyground state for all values of ϕ and no quantum phasetransition is present in the three-terminal setup. We willthus refer to the ϕ < ϕ ∗ scenario as 0-like phase whilefor ϕ > ϕ ∗ the term π -like phase will be used. More-over, for 0-like phase we have previously also discussed asecond pair of peaks in the spectral function, which liesfurther the Fermi energy. Although, they can easily beinterpreted within the effective asymmetric SIAM theycannot be explained as broadened analogs of ABS statesof the SCIAM. Crucially, this outer pair of peaks appearsbecause the corresponding many-body states can be ex-cited due to the singlet nature of the ground state whilein the SCIAM corresponding energy eigenstates do alsoexist, but they cannot be excited for ϕ < ϕ ∗ . Additionaldifference appears also in the π -like phase for ϕ > ϕ ∗ ,where the presence of non-zero DOS around the Fermienergy due to the normal electrode always gives rise tothe Kondo screening not present in the SCIAM at all. Inthe π -like phase we thus observe the coexistence of fourbroadened ABS states typical for the SCIAM with theKondo resonance at the Fermi energy. IV. FINITE-GAP MODELA. NRG calculations
As discussed in Sec. II C, the transformation T mapsthe finite-gap three-terminal setup at the half-filling ontoan NRG-tractable one-channel problem. In particular,in the basis of the w fields, the Hamiltonian describes anAnderson impurity coupled to a continuum of bath stateswith modified tunneling DOS which except of ϕ = π isparticle-hole asymmetric. This can be shown by eval-uating the hybridization function Γ W ( ω ) correspondingto Σ W ( ω ) as shown in the Appendix B. In the limit ofinfinitely wide band one obtainsΓ W ( ω ) = Γ N S | ω | Θ( ω − ∆ ) √ ω − ∆ (cid:18) ω cos ϕ (cid:19) , (43) − − − . . . ω/ ∆ ϕ/π FIG. 3. Phase evolution of the tunneling DOS Γ W ( ω ) in the w basis for Γ N / ∆ = 1 and Γ S / ∆ = 2. At ϕ = 0, only theleft BCS singularity does appear and the resulting tunnellingDOS is highly asymmetric. Increasing ϕ diminishes the asym-metry while BCS singularities develop at both gap edges. Thesymmetry is fully restored only at ϕ = π . where Θ is the Heaviside step function. The phase evolu-tion of hybridization function Γ W ( ω ) is shown in Fig. 3for selected parameters. Note that only for ϕ = π thehybridization function is symmetric with asymmetry in-creasing towards ϕ = 0. Consequently, in the w ba-sis, the particle-hole symmetry of the finite-gap three-terminal model is broken for ϕ (cid:54) = π . Moreover, asym-metry decreases with ϕ in analogy to the ∆ → ∞ case.Since Γ W ( ω ) is diagonal, the basis of the w fields allowsus to employ the standard one-channel NRG method ofRef. [38]. To this end, we have utilized NRG Ljubl-jana code [40] in one-channel mode with leads givenby the spin-unpolarized hybridization function Γ W ( ω ).To suppress numerical oscillatory artefacts of NRG, wehave used the intertwined z -discretization according tothe scheme of ˇZitko et al. [45]. Explicitly, z = n/ n ∈ { , . . . } . To achieve greater accuracy in thespectral function, the so-called self-energy trick has beenalso employed. The resulting spectral functions are thensmoother and do also correctly reproduce discontinuitiesat the BCS gap edges which are otherwise smeared out.It is also important to stress that in the main body of thearticle we concentrate on the wide band with the half-bandwidth set to B/ ∆ = 2000. The corrections requiredfor the case of a narrow band are discussed separately inthe Appendix B.The experimentally accessible spectral functions in the d basis have been obtained by means of Eqs. (33) and (34)and the results are discussed in Sec. IV B. On-dot inducedpairing ν = (cid:104) d ↓ d ↑ (cid:105) is trivially connected to the filling n w in the basis of the w fields and can therefore be measureddirectly by the expectation value of the number operatoras shown in Sec. IV C. On the contrary, the operator forthe Josephson current requires both operators of the lo-cal d as well as the lead electrons. Thus, to evaluate theJosephson current we use the integral formula of Ref. [47]which involves the anomalous component of Green func-tion in the d basis. The results on the on-dot induced0pairing and the Josephson current are also discussed inSec. IV C.Let us now briefly discuss the character of many-bodyspectrum for the present three-terminal setup with finitegap. Because of the normal metallic lead, the spectrum isexpected to by continuous. Nevertheless, applying NRGmeans that the model is first discretized and mappedonto a semi-infinite Wilson chain. At each successivelength of Wilson chain a low energy discrete spectrumis obtained while energies are typically rescaled by thelogarithmic discretization factor Λ. In a fixed RG point,one observes a steady separation of energy levels whichthen comprise the energy spectrum of the correspondinguniversal Hamiltonian.To disclose the level structure of the logarithmicallydiscretized finite-gap three-terminal setup we have em-ployed NRG Ljubljana by gradually increasing ϕ at thefixed model parameters while setting the length of thecorresponding Wilson chain to reach the low tempera-ture RG fixed point. The phase evolution of the resultingspectra is presented in Fig. 4. Although we have trans-formed the finite-gap three-terminal setup into the basisof the w fields, the transformation T is unitary and theobtained spectra are thus independent of the use of theunderlying basis of d or w fields. We also stress out thatthe parameter regime in Fig. 4 involves a sign reversal ofthe local pairing at ϕ ∗ pair ≈ . π as well as the Joseph-son current reversal at ϕ ∗ j ≈ . π (see the discussion inSec. IV C).At ϕ = π , the NRG energy spectrum consists of a sin-glet ground state, followed by the quadruplet of higherlying energy levels and an even higher lying sextet of en-ergy levels. Such a discrete spectrum is typical for ordi-nary symmetric SIAM. When ϕ is decreased the quadru-plet splits into two doublets and the sextet splits into aquadruplet with two additional singlets placed symmet-rically around it. Such outcome corresponds, however,to the ordinary asymmetric SIAM. In other words, NRGestablishes that the phase-dependent low energy many-body spectra of the logarithmically discretized finite-gapmodel are essentially indistinguishable from those for theordinary particle-hole asymmetric SIAM at ϕ (cid:54) = π andthe particle-hole symmetric SIAM at ϕ = π . The quali-tative behavior of the spectral functions of the finite-gapthree-terminal model is thus expected to essentially fol-low the results for the ∆ → ∞ case discussed in Sec. III,where the exact mapping onto the particle-hole asym-metric SIAM was shown explicitly. This way, NRG givesalso an important insight into the physical consequencesof the particle-hole asymmetry of the underlying model inthe w basis which is completely given by its hybridizationfunction Γ W ( ω ). In other words, the universal Hamilto-nian describing the low energy properties of the presentmodel corresponds to the ordinary asymmetric SIAMfor ϕ (cid:54) = π and goes smoothly over into the symmetricSIAM for ϕ → π . Therefore, we may expect that thespectral properties of the finite-gap case are only quan-titatively altered compared to the infinite-gap scenario. . . . . . . . GS, Q =0 , M z =1 Q = − , M z =2 Q = − , M z =1 Q =+1 , M z =2 Q =(0 , , M z =(1 , Q =+2 , M z =1 E / ∆ ϕ/π singletdoubletquartet FIG. 4. Low energy eigenvalues of the logarithmically dis-cretized Hamiltonian obtained using NRG Ljubljana for thehalf-bandwidth B/ ∆ = 2000, U/ ∆ = 3, Γ S / ∆ = 1 andΓ N / ∆ = 0 . Q is defined as the total charge of theWilson chain measured with respect to the half-filling and M z ≡ S z with S z being the overall magnetization of theWilson chain. Thus, quantitative changes in quantities like on-dot in-duced pairing and Josephson current are not related tothe change of the ground state or the structure of thehigher lying states, but to the fine details of the cor-responding spectral functions as discussed in Secs. IV Band IV C.
B. Spectral properties and the Kondo scale
The phase evolution of the normal spectral function inthe basis of the d fields is shown in Fig. 5 for two valuesof U/ ∆. Let us first discuss the U/ ∆ = 3 case which in-corporates features of 0-like as well as π -like phase as ob-served in the ∆ → ∞ case (panels ( a ) and ( b ) in Fig. 5).Since qualitatively finite-gap case does not differ muchfrom the ∆ → ∞ case, we will only briefly sum up themost important points. At ϕ = 0, two broadened ABS-like peaks placed symmetrically around the Fermi en-ergy are distinguishable. With increasing ϕ , both peaksmove towards the Fermi energy and merge at a certainvalue ϕ ∗ which depends non-trivially on U and ∆. Sub-sequently, for all ϕ > ϕ ∗ the central Kondo-like peak ispresent. Moreover, four side-peaks (corresponding to thetwo symmetrized Hubbard satellites in the w basis) alsoemerge. Increasing ϕ further shifts the two peaks on eachside of the spectra together, until at ϕ = π they coalesceat approximately ± U/
2. At this point, the tunnelingDOS is symmetric and the ϕ = π spectrum resemblesthe three-peak structure obtained for symmetric SIAMwhere two Hubbard satellites and one central Kondo peakare present. The second case with U/ ∆ = 6 is shown inpanels ( c ) and ( d ) of Fig. 5. Here, the ratio Γ S / ∆ is in-sufficient to generate particle-hole asymmetry leading to1the emergence of the 0-like phase. Such regimes are anal-ogous to the observations made for SCIAM at sufficientlylarge values of U/ ∆.To make the movement of the in-gap peaks explicitlymanifest, we visualized the phase-dependent positionsof their maxima via heatmaps as shown in Fig. 5b,d.We clearly observe in Fig. 5b their crossing at angle ϕ ∗ .When ϕ is further increased, the in-gap peaks move apartagain. However, for all ϕ > ϕ ∗ two additional in-gapstates emerge relatively close to the band edges. Thetwo peaks for ϕ < ϕ ∗ can thus be related to the twoABS states of SCIAM in the 0 phase, while the fouroff-central peaks are in one-to-one correspondence withthe four ABS states observed in the π phase of SCIAM.However, unlike in SCIAM the non-zero DOS around theFermi energy gives rise also to the central Kondo-like res-onance for all ϕ > ϕ ∗ . So in total, there are five in-gappeaks in the analog of π phase for the present model.Taking together, the obtained spectra qualitativelycorrespond to the ∆ → ∞ model and the physical inter-pretation in terms of the particle-hole asymmetry of theunderlying model in the basis of the w fields holds anal-ogously. Nevertheless, there are quantitative differencesto the ∆ → ∞ case which appear once integral quanti-ties such as the filling n w in the basis of the w fields areconsidered. For the ∆ → ∞ model, the filling n w mono-tonically decreases from the n w = 1 value obtained at ϕ = π for all parameter regimes. In the finite-gap three-terminal case, there are parameter regimes where n w firstincreases to values larger than 1 (positive effective chem-ical potential) and then starts to monotonically decreaseto values smaller than 1 (negative effective chemical po-tential). Such integral properties are shown in Sec. IV Cto be crucial for the system to exhibit effects such as pair-ing or Josephson current reversal which are typical of 0- π transition observed in SCIAM. This means that the pre-cise shape of spectral functions plays an important rolewhen analyzing the finite-gap case.Let us now move to the quantitative behavior in the 0-like phase. The pair of strongly pronounced in-gap peaksfor ϕ < ϕ ∗ can easily be interpreted as the broadenedABS states of the purely superconducting case. To thisend, it can be verified that their structure can be fittedwell by Lorentzian line shapes centered around the fre-quency ω ABS of the broadened ABS maxima with Γ N used as the broadening parameter. Consequently, for ϕ < ϕ ∗ the spectra may be not only qualitatively butalso quantitatively interpreted as broadened ABS statesof the pure SCIAM model. However, once ϕ approaches ϕ ∗ the broadening of the in-gap states becomes muchlarger and for ϕ > ϕ ∗ it does not match with Γ N . Onthe other hand, the number of the in-gap peaks is stillfour and they correspond well to the ABS states of the π phase of the SCIAM.On the other hand, in the π -like phase, one observesan unusual co-existence of the Kondo-like peak with fourbroadened ABS states ( ϕ > ϕ ∗ ). Since the prerequisiteof an effective Kondo screening is the non-zero tunnelling DOS around the Fermi energy brought in by the normalelectrode, the central peak forms for arbitrary values of U at ϕ = π . In Sec. IV A, we have already establishedthe correspondence of the low energy many-body NRGspectra to that of the particle-hole asymmetric SIAM for ϕ < π with particle-hole asymmetry monotonically de-creasing towards ϕ = π where it completely vanishes.In the basis of the w fields, one therefore observes thatstarting at the particle-hole symmetric case for ϕ = π and then decreasing ϕ , causes a gradual movement of theoriginal Kondo peak away from the Fermi energy which isinduced by the increasing particle-hole asymmetry. How-ever, as long as ϕ > ϕ ∗ the increasingly large broadeningdoes overcome this shift and symmetrization (33) stillleads to a well defined central peak in the spectral func-tion A WσN ( ω ). Only when ϕ is decreased further, does thebroadening of the central peak stop compensating for therapid movement of the peak, so that the symmetrization(33) results in a doubly peaked spectral function A WσN ( ω )in the basis of the w fields.In the d basis, the lack or presence of a single centralpeak can be thus understood as a sign of the Kondo-likeinteraction-screening efficiency. To capture this quanti-tatively, we have extracted the phase-dependent Kondotemperature T K as the half-width at half maximum(HWHM) value of the zero-energy peak in the spectra ofthe π -like phase. Comparison to T NK , when both super-conducting leads are decoupled, is shown in Fig. 6. Weobserve that, in general, T K is bigger in the hybrid casecompared to the ordinary SIAM. However, this cannotbe naively interpreted as the enhancement of the Kondoscreening in the hybrid reservoirs since once we relate T NK not to just the decoupled BCS electrodes but to the threenormal electrodes obtained by taking ∆ → T K would be smaller. How-ever, approaching ϕ ∗ from the π -like phase clearly en-hances T K as shown in Fig. 6. Moreover, as hypothesizedalready in Ref. [32], T K is proportional to cos ( ϕ/
2) withonly very small deviations appearing when approaching ϕ → ϕ ∗ . In the 0-like phase the central Kondo-like peakvanishes and the definition of T K as HWHM of the cor-responding peak becomes meaningless. C. Pairing and Josephson current
Let us now address also the transport properties in thehybrid three-terminal structure. Since they have been al-ready obtained in Ref. [32] using QMC, we will review theNRG results only briefly and concentrate mostly on howtransport properties may be extracted from the trans-formed basis of the w fields. Although the w basis allowsto obtain a simpler Hamiltonian, it requires redefinitionof the reservoir and thus also the geometry of the setup.The transport properties are however measured in theoriginal d basis. In particular, superconducting effectsare connected to the off-diagonal terms in the Hamil-tonian expressed in the d basis which are then rotated2 − − . . . . . . − − . . . . . . ω/ ∆ ϕ/π (a) ω/ ∆ ϕ/π (c)0 0 . . . . ϕ/π − − . . ω / ∆ (b) 0 0 . . . . ϕ/π − − . . ω / ∆ (d) FIG. 5. ( a ) Phase evolution of the spectral function A DN ( ω ) in the sub-gap region of the finite-gap three-terminal setup for U/ ∆ = 3, Γ N / ∆ = 0 .
1, Γ S / ∆ = 1 and half-bandwidth B/ ∆ = 2000. For ϕ < ϕ ∗ ≈ . π one observes a pair of broadenedABS states while for ϕ > ϕ ∗ an additional pair of broadened ABS states and an additional central Kondo-like peak do appear.In an analogy to the SCIAM, for ϕ < ϕ ∗ the regime is referred to as the 0-like phase while for ϕ > ϕ ∗ as the π -like phase. ( b )Heatmap corresponding to panel ( a ) highlights the phase-dependent position of the maxima of the in-gap peaks (dashed lines).( c ) Phase evolution of the spectral function A DN ( ω ) for the finite-gap three-terminal setup. All parameters apart from U/ ∆ = 6are the same as in panel ( a ). The Kondo correlations dominate the system which remains in the π -like phase for all values of ϕ . For sufficiently large U , such a scenario does occur also in the SCIAM. ( d ) Heatmap corresponding to panel ( c ) shows thephase-dependent position of the maxima of the in-gap peaks (dashed lines). to the diagonal elements of the Green function in the w basis. All effects related to the off-diagonal Green func-tions are therefore to be extracted by suitable backwardstransformations. Let us first discuss the on-dot inducedpairing ν which is indirectly related to the Josephson cur-rent. While the on-dot induced pairing n w in the w basisis by definition zero as G W ( ω ) has no off-diagonal en-tries, it turns out that the filling n w in w basis is directlyrelated to the pairing observed in the original basis asfollows from simple application of the transformation T to the definition of ν . In particular ν ≡ (cid:104) d ↓ d ↑ (cid:105) = 12 − n w . (44)The above equation connects the change of the filling inthe basis of the w fields from n w < n w > ν in the ba-sis of the d fields from positive to negative values. In thecase of ∆ → ∞ model the mapping onto the asymmetricSIAM fixed the possible values of filling to n w < N / Γ S is shown in Fig. 7. Theobserved dependencies are in a close relation to the well-known result obtained for SCIAM. However, there is nojump in the on-dot induced pairing that would indicate atrue phase transition. Instead, the on-dot induced pair-ings for all parameter regimes are continuous functionsand a crossover region of significant drop is only visiblefor small values of Γ N / Γ S . For small enough hybridiza-tion Γ N , see the phase dependencies for Γ N / Γ S = 0 . N / Γ S = 0 . ϕ ∗ pair . Such pairing reversalis however not directly related to the change of spectralproperties at ϕ ∗ and the two values are generally differ-ent, i. e. ϕ ∗ (cid:54) = ϕ ∗ pair . Nevertheless, for Γ N → ϕ ∗ pair approaches ϕ ∗ . The comparison of the presentresults to QMC is shown in the the Appendix C, where3 . . . . l og ( T K / T N K ) cos ( ϕ/ U/ ∆ = 346 FIG. 6. Phase dependence of the Kondo temperature T K obtained by NRG for the half-bandwidth B/ ∆ = 2000, U/ ∆ = 3, Γ N / ∆ = 0 .
1, Γ S / ∆ = 1. The Kondo temperatureis determined as the half-width at half maximum value of thecentral Kondo-like peak observed in the π -like phase. also the effect of the finite bandwidth on the phase de-pendency in the crossover region is discussed in detail. Atthis place, it is sufficient to note that deep in the 0-likeor π -like phase QMC and NRG agree well within theirnumerical accuracy.While the spectra as well as the on-dot induced pair-ing could be directly transformed back into the d basis,the operator for the Josephson current J involves also thelead c electrons and we have to resort to an expression in-volving the anomalous component of the Green function[47] J
2∆ = sin ϕ (cid:90) + ∞−∞ dωπ f ( ω )Im (cid:2) G DA ( ω )Σ DA ( ω ) (cid:3) , (45)where f ( ω ) is the Fermi-Dirac distribution, G DA is theanomalous Green function in the d basis and Σ DA ( ω ) (26)in the limit of the infinite bandwidth readsΣ DA ( ω ) = Γ S √ ∆ − ω Θ(∆ − ω )+ i Γ S sgn( ω ) √ ω − ∆ Θ( ω − ∆ ) . (46)The integral (45) involving anomalous components of theGreen function requires high frequency resolution andreliable broadening procedure when applying NRG.The results of phase-dependent Josephson current areshown in Fig. 7 for U/ ∆ = 3 at different ratios of Γ N / Γ S .The sign reversal of the Josephson current occurs onlyfor small hybridization strengths Γ N which clearly showsthat Kondo correlations are important in the system andmay overcome the superconducting correlations. Onceagain, the ϕ ∗ j at which sign reversal occurs does notmatch ϕ ∗ from crossing of the ABS states nor ϕ ∗ pair , butall three values tend to be the same for Γ N → V. CONCLUSIONS
We have investigated a general finite-gap model of QDwith an arbitrary Coulomb repulsion attached to the hy-brid reservoir composed of one normal lead and two BCSleads with an arbitrary phase difference. To obtain re-liable and method-unbiased results on phase-dependentspectral functions of the problem, the standard NRG wasemployed. However, the full problem with three types ofleads requires within the standard NRG approach to im-plement three-channel calculations which poses severalnon-trivial challenges [41]. To circumvent the numericallimitations we have thus introduced a unitary transfor-mation T of the local dot electrons d , which allows to re-formulate the finite-gap model with hybrid phase-biasedreservoirs at the half-filling, despite the general belief[38], as a one channel problem. Moreover, the approachis by no means limited to the non-zero values of Γ N andholds also in the case of the decoupled normal electrodewhere the three-terminal setup goes over into the ordi-nary SCIAM. Consequently, at the half-filling the trans-formation T transforms SCIAM onto the SIAM with thetunneling DOS that incorporates a hard gap. The vanish-ing tunelling DOS around the Fermi energy consequentlyhampers the formation of the Kondo resonance. How-ever, the quantitative understanding requires to applyNRG techniques with either finite-sized Wilson chains,when the standard logarithmic discretization around theFermi energy is used, or to modify the discretization toavoid the energy range of the hard gap.In this regard, the present three-terminal reservoir hasalways a non-zero tunneling DOS around the Fermi en-ergy and a standard logarithmic discretization in thetransformed basis of the fields w σ with σ ∈ {↑ , ↓} canbe employed. Thus, the open-source NRG Ljubljanacode could be employed unaltered. The obtained phase-dependent spectral functions showed behavior similar tothat observed in the 0 and π phases of the SCIAM.Thus, two regimes, referred here as the 0-like and the π -like phase, have been identified, see Figs. 5 and 7 formore details. However, the two regimes do not consti-tute separate phases, because even though the normalelectrode may be only weakly hybridized with the re-maining system, it always provides a non-zero tunnellingDOS around the Fermi energy which enables the forma-tion of the Kondo many-body ground state of the singletnature. In other words, 0-like and the π -like phases repre-sent separate regimes which are smoothly interconnectedwith each other when the phase bias ϕ is changed con-tinuously. The width of the resulting crossover region,as shown in Figs. 5 and 7, is roughly proportional to thehybridization Γ N .In more detail, when the hybridization Γ N is suffi-ciently small, see Fig. 5, one obtains quite a narrowcrossover region in which the off-center in-gap peaks ofthe corresponding spectral function do cross at ϕ ∗ . More-over, one also observes sign changes in the transportproperties in the crossover region with the on-dot induced4 . . . . . . . . (a) − . . . . . .
25 0 0 . . . . (b) ν ϕ/π Γ N / Γ S = 0 . . . . . J / J ϕ/π Γ N / Γ S = 0 . . . . . FIG. 7. ( a ) Phase evolution of the on-dot induced pairing ν ≡ (cid:104) d ↓ d ↑ (cid:105) of the finite-gap three-terminal setup at T = 0 obtainedusing NRG Ljubljana in the basis of the w fields with subsequent use of Eq. (44) for the half-bandwidth B/ ∆ = 2000 (solidlines) and B/ ∆ = 100 (dashed lines). Parameters of the model are U/ ∆ = 3 and Γ S / ∆ = 1 while Γ N varies. ( b ) Phaseevolution of the Josephson current J obtained using NRG Ljubljana in the basis of the w fields with subsequent use of Eq. (45)for the same parameters as in panel ( a ). pairing changing sign at ϕ ∗ pair and the Josephson currentat ϕ ∗ j , see Fig. 7 for more details. Although generally thevalues of ϕ ∗ , ϕ ∗ pair and ϕ ∗ j do not equal, they do so in thelimit of Γ N →
0, where they reach the SCIAM value ofthe critical phase bias of the 0- π transition. Moreover,at large values of Γ N pairing and Josephson current arealways positive, ϕ ∗ pair and ϕ ∗ j are thus not even definedwhile the distinct spectral properties of the 0-like andthe π -like phase are still present. The identification ofthe 0-like and the π -like phases is therefore always basedonly on their spectral properties. For ϕ > ϕ ∗ the defin-ing spectral property of the π -like phase is the presenceof the two pairs of the in-gap peaks which merge to-gether at ϕ = π and the presence of a central Kondo-likeresonance. The two pairs of the in-gap peaks show ananalogous phase-dependent behavior as the ABS statesof the SCIAM (see the heatmaps in Fig. 5). The twopairs of the in-gap peaks can therefore be understood asbroadened analogs of the ABS states of the SCIAM. How-ever, unlike in the SCIAM, the tunneling DOS due to thepresence of the normal lead is non-zero around the Fermienergy which causes an effective screening of the spin ofthe QD in the π -like phase as manifested by the presencethe Kondo resonance in the spectral function. Thus, un-like in the SCIAM, in the finite-gap three-terminal setupthe broadened ABS states do co-exist with the Kondoresonance in the π -like phase, see panels a ) and c ) inFig. 5 for more details. Nevertheless, at ϕ ∗ the cen-tral peak at the Fermi energy splits and can no longerbe attributed to Kondo-like correlations since the effec-tive underlying model is already far from the half-fillingso that charge excitations dominate. For ϕ < ϕ ∗ thespectral weight at the Fermi energy tends to zero withfurther decreasing the phase difference ϕ , a region corre-sponding to the 0-like phase. Here, the split peak can be interpreted in terms of two broadened ABS states withtheir phase-dependent behavior resembling the SCIAM.However, such a 0-like phase has an additional pair oflow-intensity peaks at higher frequencies, which (unlikein SCIAM) can be excited in the one-particle mannerdue to admixtures of the doublet state induced by thecoupling to the normal lead.Such complex behavior is qualitatively explained viathe transformation T which is thus not merely a technicaltool for the implementation of NRG. The transformation T maps the Hamiltonian of the three-terminal model atthe half-filling onto the SIAM with a structured tunnelingDOS shown in Fig. 3. The tunneling DOS in the basis ofthe w fields is at ϕ = 0 highly particle-hole asymmetric.The asymmetry is then continuously diminished with in-creasing ϕ and vanishes completely at ϕ = π . As shownin Fig. 4 the particle-hole asymmetry of the tunnelingDOS effectively acts as the particle-hole asymmetry inan ordinary SIAM with the concomitant increase of theKondo temperature upon increasing the asymmetry fol-lowed by entering the mixed valence regime and eventu-ally complete destruction of the Kondo resonance. TheKondo temperature T K can be quantitatively assessedvia the phase-dependent HWHM of the central Kondo-like peak. The analysis in Sec. IV C (see also Fig. 6)showed that log T K follows (almost up to the value of ϕ ∗ ) the cos ( ϕ/
2) trend hypothesized already in Ref. [32]based on the infinite-gap limit of the present model andthe second-order perturbation theory.Moreover, using the transformation T we have also ob-tained the phase-dependent on-dot induced pairing andthe phase-dependent Josephson current in various para-metric ranges. Results in the limit of the infinitely wideband are presented in Fig. 7, where also effects of finitewidth of the band are shown (with a detailed deriva-5tion given in the Appendix B). Incorporating these cor-rections allowed also for comparison with another nu-merically exact method, the continuous-time hybridiza-tion expansion (CT-HYB) QMC, which demonstrates anagreement between both methods in the regions outsideof the crossover while the crossover itself shows largetemperature dependence, see Appendix C. As alreadydiscussed, the resulting pairing and supercurrent rever-sals do not occur exactly at ϕ ∗ identified from the spec-tral functions. Moreover, sign reversals appear only atsufficiently low ratios Γ N / Γ S . Once a given thresholdis exceeded the Kondo screening dominates the system,superconducting correlations are essentially suppressedand only modify the phase-dependent transport. Thisbehavior is always enhanced by increasing the interac-tion strength, reducing the gap size, or increasing thehybridization of the normal lead. The observation of the0- π -like crossover is thus possible only in a fairly smallpart of the parameter space corresponding to the weakcoupling of the normal lead to the QD.The mapping T not only significantly reduces the com-plexity of the hybrid normal-superconductor reservoirsduring the NRG implementation, but it also allows forconceptual understanding of the competing Kondo andJosephson effects via the particle-hole asymmetric SIAM.In this regard, it is worth mentioning that the transfor-mation T applies in the same form also to the SCIAM,i.e., an interacting quantum dot coupled to two supercon-ducting leads with hard gap in the spectrum, leading toits mapping onto the problem of normal Anderson impu-rity coupled to an insulator-like electronic reservoir witha hard spectral gap around the Fermi level. The origi-nal 2 × ACKNOWLEDGMENTS
We acknowledge discussions with Rok ˇZitko, MartinˇZonda, and V´aclav Janiˇs. This work was supported byGrant No. 19-13525S of the Czech Science Foundation(PZ, TN), by grant INTER-COST LTC19045 (VP), bythe COST Action NANOCOHYBRI (CA16218) (TN),the National Science Centre (NCN, Poland) via GrantNo. UMO-2017/27/B/ST3/01911 (TN) and by The Min-istry of Education, Youth and Sports from the LargeInfrastructures for Research, Experimental Developmentand Innovations project “IT4Innovations National Super- computing Center LM2015070” (VP). Access to comput-ing and storage facilities owned by parties and projectscontributing to the National Grid Infrastructure Meta-Centrum provided under the programme “Projects ofLarge Research, Development, and Innovations Infras-tructures” (CESNET LM2015042) is also appreciated.
Appendix A: Interaction term in the basis of the wfields
To calculate the interaction part H U in the basis of the w fields, let us first apply the transformation T − to thefollowing quantities d †↑ d ↑ = 1 + w †↑ w ↑ − w †↓ w ↓ − w †↑ w †↓ + w ↓ w ↑ , (A1) d †↓ d ↓ = 1 − w †↑ w ↑ + w †↓ w ↓ − w †↑ w †↓ + w ↓ w ↑ . (A2)We also define µ = w †↑ w ↑ − w †↓ w ↓ and ξ = w †↑ w †↓ + w ↓ w ↑ which satisfy1 − µ = ξ = − w †↑ w ↑ + 2 n w ↑ n w ↓ + w ↓ w †↓ , (A3) µξ = ξµ = 0 , (A4)where n w ↑ = w †↑ w ↑ and n w ↓ = w †↓ w ↓ . Using these rela-tions we obtain d †↑ d ↑ d †↓ d ↓ = 1 − µ − µξ + ξµ − ξ + ξ n w ↑ n w ↓ − W † ( σ x + σ z ) W . (A5)Consequently, applying the T − transformation to the in-teraction term H U gives H U = U d †↑ d ↑ d †↓ d ↓ − U D † ( σ x + σ z ) D = U w †↑ w ↑ w †↓ w ↓ , (A6)which in the basis of the w fields obtains the form ofthe ordinary Hubbard term. Notice that the additionalquadratic term in H U in the basis of the d fields cancelsexactly the quadratic term of Eq. (A5) when correspond-ingly transformed. Appendix B: Corrections due to the finitebandwidth
It is now our aim to give a detailed derivation of theself-energy contribution Σ D for the hybrid three-terminalsetup under the study. The discussion simplifies when Σ D is divided right from the beginning into its contribution Σ DN due to just the normal lead and its superconductingcontribution Σ DS due to the two phase-biased supercon-ducting leads. Explicitly: Σ D ( z ) = Σ DN ( z ) + Σ DS ( z ) (B1)6with Σ DN ( z ) = (cid:88) k V N k ( z + − E N k ) − V N k Σ DS ( z ) = (cid:88) α ∈{ L,R } , k V α k ( z + − E α k ) − V α k (B2)with z being an arbitrary complex number. Later, onlythe functional form only infinitesimally close to the realaxis needs to be resolved. For that we set z = ω + ≡ ω + iη with ω being a real frequency and η being aninfinitesimally small positive number, thus taking the cutslightly above the real axis. The self-energy contributionfrom the normal lead for the constant density of statesin the band ρ is well-known and reads as Σ DN ( ω + ) = i Γ N (B3)with Γ N ≡ πρ | V N | . The superconducting part on theother hand is non-trivial and shall be treated here in moredetail. It is defined as Σ DS ( ω + ) = (cid:88) α ∈{ L,R } k V α k (cid:0) ω + − E α k (cid:1) − V α k (B4)with notation following that of the main article. The in-verse matrix (the resolvent) appearing in Eq. (B4) yields (cid:0) ω + − E α k (cid:1) − = ω + + ∆ α C α σ x − ∆ α S α σ y + ε k σ z ( ω + iη ) − ∆ − ε k . (B5)Since η is infinitesimally small the denominator can besimplified to (cid:0) ω + − E α k (cid:1) − = ω + ∆ α C α σ x − ∆ α S α σ y + ε k σ z ω − ∆ − ε k + iη sgn( ω ) . (B6)Considering now V α k ≡ V S with V S being real,we may use the resolvent expression (B6) and simplifyEq. (B4) to obtain Σ DS ( ω + ) as Σ DS ( ω + ) = 2 | V S | (cid:88) k ω − ∆ cos (cid:0) ϕ (cid:1) σ x + ε k σ z ω − ∆ − ε k + iη sgn( ω ) . (B7)Moreover, assuming a constant density of states ρ inthe band the k summations may be replaced by the en-ergy integration: Σ DS ( ω + ) = 2Γ S (cid:104) ω − ∆ cos (cid:16) ϕ (cid:17) σ x (cid:105) F ( ω + ) , (B8)where Γ S ≡ πρ V S and F ( ω + ) = 1 π (cid:90) B − B dεω − ∆ − ε + iη sgn( ω ) (B9)with B being the half-bandwidth and the term propor-tional to σ z vanished due to the integrand being an odd function of ε . To resolve the integrals we need to distin-guish the in-gap region ∆ − ω > η → limit can be taken directly) with F ( ω + ) = − π √ ∆ − ω arctan (cid:18) B √ ∆ − ω (cid:19) . (B10)The out-of-the gap region ω − ∆ > F ( ω + ) = − i sgn( ω ) √ ω − ∆ + ln (cid:16) B + √ ω − ∆ B −√ ω − ∆ (cid:17) π √ ω − ∆ . (B11)The resulting Σ D ( ω + ) has thus a non-zero imaginarypart only outside of the in-gap region while all effects ofthe finite-sized band are to be found in its real part whichis non-zero in the whole band. However, once the limit B → ∞ is taken the real part out of the in-gap regionvanishes also.Taking together, the self-energy contribution Σ D ( ω + )takes the form of (26), whereΣ DN ( ω ) = i Γ N + 2Γ S F ( ω + ) (B12)and Σ DA ( ω ) = − S ∆ cos (cid:16) ϕ (cid:17) F ( ω + ) . (B13) Appendix C: Comparison of the NRG results withQMC
In order to assess the ability of the presented NRGscheme to provide reliable results on the integral quan-tities like the on-dot induced pairing and the Josephsoncurrent, we compare the results with a numerically ex-act continuous-time, hybridization-expansion (CT-HYB)QMC, as this method was already successfully used tostudy both the two-terminal [48] and three-terminal [32]setups and agrees with standard NRG results well withinthe QMC error bars.The CT-HYB calculation is performed in the origi-nal basis of d fields employing off-diagonal elements ofthe hybridization function using the TRIQS/CTHYBsolver [49]. The total Hamiltonian of the system does notconserve particle number, therefore the superconduct-ing pairing is introduced to the method using a canon-ical particle-hole transformation in the spin-down sec-tor, mapping the system to an impurity Anderson modelwith attractive interaction [50, 51]. As CT-HYB is an in-herently finite-temperature method, all calculations wereperformed at k B T / ∆ = 0 . B = 100∆ and a cutoff in Matsubarafrequencies ω maxn ≈ . . . . . . . . (a) T = 0 00 . . . . . . . . (c) k B T = 0 . − . . . . . . . . (b) T = 0 − . . . . . . . . (d) k B T = 0 . ν ϕ/π NRG, Γ N / Γ S = 0 . . . . . ν ϕ/π CT-HYB, Γ N / Γ S = 0 . . . . . J / J ϕ/π NRG, Γ N / Γ S = 0 . . . . . J / J ϕ/π NRG, T = 0, Γ N = 0CT-HYB, Γ N / Γ S = 0 . . . . . FIG. 8. Left column: Phase-dependent on-dot pairing ν ≡ (cid:104) d ↓ d ↑ (cid:105) (panel a ) and the Josephson current J (panel b) calculatedusing NRG at zero temperature in the basis of the w fields with subsequent use of Eqs. (44) and (45) for the parameters B/ ∆ = 100, U/ ∆ = 3, Γ S / ∆ = 1 and various values of Γ N . Right column: Equivalent results calculated using CT-HYB QMCin the basis of the d fields for small finite temperature k B T / ∆ = 0 .
025 (symbols with error bars): on-dot pairing ν (panel c )and Josephson current J (panel d) with J = 2 e ∆ / (cid:126) . Lines are splines of QMC data and serve only as a guide to the eye. for the induced pairing ν as a function of phase difference ϕ for B/ ∆ = 100, U/ ∆ = 3, Γ S = ∆ and various val-ues of Γ N at zero temperature are plotted in panel a (topleft). The importance of the finite-bandwidth correctionswere already discussed in Fig. 7. The equivalent resultsof CT-HYB for small finite temperature are plotted inpanel c (top right). The curves match within QMC errorbars for small and large values of ϕ . In the crossover re-gion, the results slightly differ as the finite temperatureis a source of additional smearing, having a similar effectas Γ N , as already discussed in Ref. [32]. In panel b (bot- tom left) we plotted the NRG results for the Josephsoncurrent for the same set of parameters as in panel a . Weadded a Γ N = 0 result from Ref. [32] to mark the posi-tion of the QPT in a case of detached normal electrode.The relevant CT-HYB result is again plotted in panel d (bottom right). Comparison again shows good agree-ment up to the finite-temperature effects, proving thatthe presented NRG scheme is able to provide correct re-sults on the Josephson current even though it cannot bemeasured directly and must be reconstructed from thecalculated anomalous spectral function. [1] S. J. Tans, M. H. Devoret, H. Dai, A. Thess, R. E. Smal-ley, L. J. Geerligs, and C. Dekker, Nature , 474(1997).[2] 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).[3] P. Jarillo-Herrero, J. A. van Dam, and L. P. Kouwen-hoven, Nature , 953 (2006).[4] H. I. Jørgensen, K. Grove-Rasmussen, T. Novotn´y,K. Flensberg, and P. E. Lindelof, Phys. Rev. Lett. , , 53(2006).[6] A. Eichler, R. Deblock, M. Weiss, C. Karrasch, V. Meden,C. Sch¨onenberger, and H. Bouchiat, Phys. Rev. B ,161407 (2009).[7] J.-D. Pillet, C. H. L. Quay, P. Morfin, C. Bena, A. L.Yeyati, and P. Joyez, Nat. Phys. , 965 (2010).[8] R. Maurand, T. Meng, E. Bonet, S. Florens, L. Marty,and W. Wernsdorfer, Phys. Rev. X , 011009 (2012).[9] J. D. Pillet, P. Joyez, R. ˇZitko, and M. F. Goffman,Phys. Rev. B , 045101 (2013).[10] J. A. van Dam, Y. V. Nazarov, E. P. A. M. Bakkers,S. De Franceschi, and L. P. Kouwenhoven, Nature ,667 (2006).[11] E. J. H. Lee, X. Jiang, R. Aguado, G. Katsaros, C. M.Lieber, and S. De Franceschi, Phys. Rev. Lett. ,186802 (2012).[12] E. J. H. Lee, X. Jiang, R. Zitko, R. Aguado, C. M. Lieber,and S. D. Franceschi, “Scaling of sub-gap excitations in asuperconductor-semiconductor nanowire quantum dot,”(2016), arXiv:1609.07582.[13] S. Li, N. Kang, P. Caroff, and H. Q. Xu, Phys. Rev. B , 014515 (2017).[14] A. C. Hewson, The Kondo Problem to Heavy Fermions ,Cambridge Studies in Magnetism (Cambridge UniversityPress, 1993).[15] K. G. Wilson, Rev. Mod. Phys. , 773 (1975).[16] P. Kopietz, L. Bartosch, and F. Sch¨utz, Introduction tothe functional renormalization group , Lect. Notes Phys.(Springer-Verlag Berlin Heidelberg, 2010).[17] S. Streib, A. Isidori, and P. Kopietz, Phys. Rev. B ,201107 (2013).[18] K. Edwards, A. C. Hewson, and V. Pandis, Phys. Rev.B , 165128 (2013).[19] V. Janiˇs and P. Augustinsk´y, Phys. Rev. B , 165108(2007).[20] M. R. Buitelaar, T. Nussbaumer, and C. Sch¨onenberger,Phys. Rev. Lett. , 256801 (2002).[21] A. Eichler, M. Weiss, S. Oberholzer, C. Sch¨onenberger,A. Levy Yeyati, J. C. Cuevas, and A. Mart´ın-Rodero,Phys. Rev. Lett. , 126602 (2007).[22] T. Sand-Jespersen, J. Paaske, B. M. Andersen, K. Grove-Rasmussen, H. I. Jørgensen, M. Aagesen, C. B. Sørensen,P. E. Lindelof, K. Flensberg, and J. Nyg˚ard, Phys. Rev.Lett. , 126603 (2007).[23] C. Buizert, A. Oiwa, K. Shibata, K. Hirakawa, andS. Tarucha, Phys. Rev. Lett. , 136806 (2007).[24] K. Grove-Rasmussen, H. I. Jørgensen, and P. E. Lindelof,New J. Phys. , 124 (2007).[25] D. J. Luitz, F. F. Assaad, T. Novotn´y, C. Karrasch, andV. Meden, Phys. Rev. Lett. , 227001 (2012).[26] A. Mart´ın-Rodero and A. Levy Yeyati, Adv. Phys. ,899 (2011).[27] V. Meden, Journal of Physics: Condensed Matter ,163001 (2019).[28] S. De Franceschi, L. Kouwenhoven, C. Sch¨onenberger, and W. Wernsdorfer, Nat. Nanotechnol. , 703 (2010).[29] H. I. Jørgensen, T. Novotn´y, K. Grove-Rasmussen,K. Flensberg, and P. E. Lindelof, Nano Lett. , 2441(2007).[30] R. ˇZitko, J. S. Lim, R. L´opez, and R. Aguado, Phys.Rev. B , 045441 (2015).[31] A. Jellinggaard, K. Grove-Rasmussen, M. H. Madsen,and J. Nyg˚ard, Phys. Rev. B , 064520 (2016).[32] T. Doma´nski, M. ˇZonda, V. Pokorn´y, G. G´orski,V. Janiˇs, and T. Novotn´y, Physical Review B , 045104(2017).[33] K. Satori, H. Shiba, O. Sakai, and Y. Shimizu, J. Phys.Soc. Japan. , 3239 (1992).[34] T. Yoshioka and Y. Ohashi, J. Phys. Soc. Jpn. , 1812(2000).[35] J.-G. Liu, D. Wang, and Q.-H. Wang, Phys. Rev. B ,035102 (2016).[36] Y. Tanaka, A. Oguri, and A. C. Hewson, New J. Phys. , 115 (2007).[37] A. Oguri and Y. Tanaka, J. Phys.: Conf. Ser. , 012146(2012).[38] R. Bulla, J. Keller, and T. Pruschke, Z. Phys. B Con-dens. Matter , 195 (1994).[39] T. Hecht, A. Weichselbaum, J. von Delft, and R. Bulla,J. Phys.: Cond. Mat. , 275213 (2008).[40] R. ˇZitko, “NRG Ljubljana - open source numerical renor-malization group code,” (2014), nrgljubljana.ijs.si.[41] A. K. Mitchell, M. R. Galpin, S. Wilson-Fletcher, D. E.Logan, and R. Bulla, Phys. Rev. B , 121105 (2014).[42] T. Novotn´y, A. Rossini, and K. Flensberg, Phys. Rev.B , 224502 (2005).[43] A. Kadlecov´a, M. ˇZonda, and T. Novotn´y, Phys. Rev. B , 195114 (2017).[44] H. R. Krishna-murthy, J. W. Wilkins, and K. G. Wilson,Phys. Rev. B , 1044 (1980).[45] R. ˇZitko and T. Pruschke, Phys. Rev. B , 085106(2009).[46] The broadening has been also predicted analytically inRefs. [32, 52] via the Schrieffer-Wolff transformation. Us-ing the transformation T , the ∆ → ∞ can be mappedonto the particle-hole asymmetric SIAM, from which theenhancement of the exchange coupling J compared to theparticle-hole symmetric case follows trivially. Thus, as aconsequence of the locally induced SC pairing, the T K isenhanced and the central peak is broader compared tothe single-channel normal SIAM.[47] M. ˇZonda, V. Pokorn´y, V. Janiˇs, and T. Novotn´y, Sci.Rep. , 8821 (2015).[48] A. Kadlecov´a, M. ˇZonda, V. Pokorn´y, and T. Novotn´y,Phys. Rev. Applied , 044094 (2019).[49] P. Seth, I. Krivenko, M. Ferrero, and O. Parcollet, Com-put. Phys. Commun. , 274 (2016).[50] D. J. Luitz and F. F. Assaad, Phys. Rev. B , 024509(2010).[51] V. Pokorn´y and M. ˇZonda, Physica B , 488 (2018).[52] T. Doma´nski, I. Weymann, M. Bara´nska, and G. G´orski,Sci. Rep.6