# Josephson Current through a Quantum Dot Connected with Superconducting Leads

Authors: Satoshi Kawaguchi
Journal: Advances in Condensed Matter Physics
 (2019)

---

## Abstract

In this study, we consider the Josephson current in a system composed of a superconductor/quantum dot/superconductor junction. In the model, the Coulomb interaction in the quantum dot is taken into consideration, and the Lacroix approximation is applied to study the electron correlation. We derive Green’s function of the quantum dot by applying the Lacroix truncation. Although the Andreev bound state does not occur in our formulations, theπ-junction occurs for a restricted parameter range. On comparing the Kondo temperature with that estimated by another method, it is found that our Lacroix approximation does not capture well the Kondo physics in the superconductor/quantum dot/superconductor junction.


---

## Body

## 1. Introduction

In the last few decades, the Josephson current in a system composed of a superconductor/quantum dot/superconductor (S/QD/S) junction has been extensively studied [1–4]. When identical superconducting leads are separated by a thin layer of insulator, the Josephson current can flow because of the coherent tunneling of Cooper pairs across the insulator in the absence of a potential difference. When the tunneling amplitude across the barrier is small and the spin is conserved, the current depends on the superconductor phase difference θ between the left and right leads. The current is expressed as J=Jcsinθ, where Jc denotes the critical current, which is proportional to the normal conductance through the barrier [5]. When θ = 0, the Josephson current is zero and the junction is in the ground state. When θ = π, the current becomes zero; however, in this case, the junction energy is maximum and it is in an unstable state. Very recently, carbon-nanotube Josephson junction systems have been studied intensively [6–8].When we consider physics at low temperatures under the Coulomb interaction in the QD, the physical behavior of the system depends on the relative magnitude of the Kondo temperatureTK and the BCS gap Δ [9]. When TK≫Δ, the Kondo effect is sufficiently strong to break the Cooper pair at the Fermi level, and the localized spin in the QD is screened; it is expected that a Kondo singlet will be formed. This results in a positive critical current (0-junction). By contrast, when TK≪Δ, the Cooper pair is strongly coupled, and the Kondo screening is essentially negligible. In this case, the Cooper pair is subjected to a localized magnetic moment in the QD. When the Coulomb repulsion in the QD is large, the ground state of the QD is a magnetic doublet, and the electrons in a Cooper pair can tunnel one by one via virtual processes. The spin ordering of the Cooper pair is reversed, resulting in a π-junction [3, 10, 11]. The 0-π transition is expected to occur around TK∼Δ. The tuning of the π-junction in S/QD/S systems has been studied extensively [3, 12–14].The current density in the S/QD/S system is composed of two parts: continuous and discrete spectrums. The former (latter) arises from outside (inside) the BCS gap. The current from the continuous current density is calculated by the usual numerical integral. By contrast, the current from the discrete current density is calculated by applying the complex function theory. The discrete current density arises from the Andreev bound states (ABSs) inside the BCS gap. The ABSs are determined as poles of the QD Green’s function, and there is a pair of ABSsω±0 in the absence of Coulomb interaction. In particular, when the energy level of the QD coincides with the Fermi level ωF=0, ω±0=0. In this situation, the current jumps discontinuously when the BCS phase difference is ±π [15, 16]. Usually, the current from the discrete current density is much larger than that from the continuous current density.The Coulomb interaction in the QD is studied by many methods: the numerical renormalization group (NRG) method [17–19], noncrossing approximation [20], quantum Monte Carlo (QMC) method [21–23], and so on. Although the NRG and QMC methods are very precise methods, they are computationally expensive. The simplest method is the Hartree approximation, which corresponds to the zeroth-order approximation. In this method, a two-particle Green’s function is truncated by decoupling at the mean-field level. Beyond the Hartree approximation, the Hartree–Fock (HF) approximation is proposed, where up to the first order of the tunneling amplitude is considered. This method was applied to the level-crossing quantum phase transition between the BCS-singlet and the magnetic doublet states [24–27]. The above two approximations can be applied to describe single-particle physics. The higher-order Green’s functions are not taken into consideration; therefore, these approximations are not sufficient when we consider Kondo physics [28]. To overcome this defect, the Lacroix approximation has been proposed [29]. In this approximation, a greater higher-order correlation effect is included in the QD Green’s function by truncation in the second order. Although the Lacroix approximation suffers from several defects, the mathematical procedures to derive the QD Green’s functions are a simple application of the equation of motion. Although there have been many studies on QD systems that employ the Lacroix approximation [28, 30–34], only a few studies have been conducted on the current in S/QD/S systems [33, 35].From these standpoints, we examine the current in a system composed of an S/QD/S junction with Coulomb interaction by applying the Lacroix approximation. Under second-order truncation and simplification, Green’s function of the QD is obtained. Using Green’s functions, we calculate the electron occupation number in the QD and the Josephson current. We can observe theπ-junction in a restricted parameter range, but our Lacroix approximation does not capture well the competition between the Kondo effect and superconductivity.


## 2. Model and Formulation

We first introduce the setup of the system and give its Hamiltonian. For the system, Green’s functions of the QD and the Josephson current are derived by employing the equation of motion.


### 2.1. Model and Green’s Functions of the QD

We consider a system composed of an S/QD/S junction, where Coulomb interaction exists in the QD. The geometry of the setup is shown in Figure1. The total Hamiltonian of the system is written as(1)H=HL+HR+HD+HT,where(2)Hα=∑k,σεkcα,k,σ†cα,k,σ+∑kΔαcα,k,↑†cα,−k,↓†+H.c.,HD=∑σεddσ†dσ+Und↑nd↓,HT=t∑k,σcL,k,σ†dσ+cR,k,σ†dσ+H.c.,where α=L,R and σ=↑,↓. Hα represents the superconducting lead α; cα,k,σ denotes the annihilation operator of an electron with energy εk, wave number k, and spin σ in the lead; the order parameter Δα=Δeiθα with the BCS gap Δ and the BCS phase θα; HD represents the QD; dσ denotes the annihilation operator of an electron with spin σ; and εd denotes the QD energy level. The occupation number of an electron in the QD with spin σ is defined by ndσ=dσ†dσ, and U represents the Coulomb repulsion between electrons with up- and downspins. HT represents the electron tunneling between the leads and the QD. The coupling strength between electrons in the QD and the leads is defined by Γ0=πνt2, where the tunneling amplitude t is real and v denotes the normal density of states (DOS) at the Fermi level.Figure 1
Schematic diagram of a superconductor/quantum dot/superconductor (S/QD/S) junction.t denotes the tunneling amplitude between the leads and the QD.To describe the above system, we introduce the 2 × 2 Nambu representation. We set the spinor field operators as(3)Ψαkσ=cα,k,σcα,−k,σ¯†,Ψαkσ†=cα,k,σ†cα,−k,σ¯,Ψdσ=dσdσ¯†,Ψdσ†=dσ†dσ¯,where σ¯=↓,↑. Hereinafter, Green’s function written with a hat symbol represents a 2 × 2 matrix. Using the above operators, the retarded (advanced) 2 × 2 matrix Green’s function is defined as G^αkσ,βk′σ′rat,t′=∓iθ±t∓t′Ψαkσt,Ψβk′σ†t′. The lesser Green’s functions are defined as follows:(4)G^αkσ,α′k′σ′<t,t′=icα′,k′,σ′†t′cα,k,σtcα′,−k′,σ¯′t′cα,k,σtcα′,k′,σ′†t′cα,−k,σ¯†tcα′,−k′,σ¯′t′cα,−k,σ¯†t,G^dσ,αkσ<t,t′=icα,k,σ†t′dσtcα,−k,σ¯t′dσtcα,k,σ†t′dσ¯†tcα,−k,σ¯t′dσ¯†t.Employing the equation of motion, we derive all Green’s functions. We denoteg^dr as the retarded Green’s function of the QD in the absence of coupling with the leads, and g^ασr is denoted as the retarded Green’s function of the lead α with spin σ in the absence of coupling with the QD. In the Nambu representation, these functions are written as(5)g^drω−1=ω−εd+i0+00ω+εd+i0+,g^ασrω=−iπνρω1σΔαωσΔ¯αω1,respectively, where a symbol with an overbar represents the complex conjugate of the symbol, and the factor ρω is defined as(6)ρω=ωω2−Δ2,ω>Δ,ωiΔ2−ω2,ω<Δ.It is noteworthy thatρ denotes the ordinary BCS DOS for ω>Δ; however, it is imaginary for ω<Δ.We define Green’s functions of the QD with spinσ using the Zubarev notation as(7)G^dσrω=dσdσ¯†;dσ†dσ¯=dσ;dσ†dσ;dσ¯dσ¯†;dσ†dσ¯†;dσ¯.By using the Dyson equation, we obtain the retarded Green’s functions of the QD as(8)G^dσrω−1=G^dσ0rω−1−Σ^σUω,where Σ^σUω is the self-energy due to the on-site Coulomb interaction. Following references [19, 36], Σ^σUω=UF^dσωG^dσrω−1, where F^dσω is defined as(9)F^dσω=dσndσ¯;dσ†dσndσ¯;dσ¯−dσ¯†ndσ;dσ†−dσ¯†ndσ;dσ¯.In the absence of Coulomb interaction,G^dσ0rω is given by(10)G^dσ0rω−1=g^drω−1−Σ^0ω,where Σ^0 is the noninteracting self-energy, which is calculated by(11)Σ^0ω=T^g^Lσrω+g^RσrωT^,where T^=tσ^z, and σ^z denotes the third component of the 2 × 2 Pauli matrices. Each element of Σ^0 is given in Appendix.In the presence of Coulomb interaction, it is necessary to calculateΣ^σU self-consistently. In this study, we derive Green’s functions, G^dσr, under the Lacroix approximation. Although the Lacroix truncation was originally proposed for the Anderson model with normal conducting leads [29], we extend the method to superconducting leads. The detailed derivations are given in Appendix. The final expression of the (11)-component of G^dσr is(12)G^dσrω11=1Rσω−UQσω+Tσω/Pσω,where(13)Pσω=ω−εd−U1−ndσ¯+3iΓ0+UA1σω−A2σω,Qσω=ω−εd−Undσ¯A1σω−A2σω+Undσ¯1−ndσ¯−dσ¯dσdσ†dσ¯†−2iΓ0ndσ¯−B1σω+B2σω,Rσω=ω−εd−Σ110ω−Undσ¯−Σ120ω+Udσ¯dσΣ210ω+Udσ†dσ¯†ω+εd−Σ220ω+Undσ,Tσω=2εd+U1+ndσ−ndσ¯×dσ¯dσΣ210ω+Udσ†dσ¯†ω+εd−Σ220ω+Undσ,and the (21)-component of G^dσr is given by(14)G^dσrω21=Σ210ω+Udσ†dσ¯†ω+εd−Σ220ω+UndσG^dσrω11.We can obtainG^dσrω12 and G^dσrω22 in a similar way, and the final expressions are given in Appendix.


### 2.2. Current and Electron Occupation Number in the QD

The current from the leadL to the lead R is given by the time derivative of the electron occupation number in the lead L. Applying the equation of motion [37, 38] and the definition NL=∑k,σcL,k,σ†cL,k,σ, we obtain the current IL as follows [39]:(15)IL=−edNLdt=ieℏNL,H=ieℏt∑k,σcL,k,σ†dσ−dσ†cL,k,σ=eℏt∑k,σGdσ,Lkσ<t,t−GLkσ,dσ<t,t,where the matrix elements of the lesser Green’s functions are denoted as follows: Gdσ,Lkσ<t,t′=icL,k,σ†t′dσt. Similarly, IR, which is the current from the lead R to the lead L, is obtained. Using the 2 × 2 matrix Green’s functions and their Fourier-transformed forms, we obtain the Josephson current I as(16)I=12IL−IR=12et2πℏ∑k,σ∫dωG^dσ,Lkσ<ω−G^Lkσ,dσ<ω−G^dσ,Rkσ<ω+G^Rkσ,dσ<ω11,where G^ij is the i,j-th element of G^. The occupation number of an electron with spin σ in the QD is given by(17)nd↑=12πi∫dωG^d↑<ω11,nd↓=12πi∫dωG^d↓<ω11.The averaged occupation number in the QD isnd=nd↑+nd↓.To calculate the current (16), the retarded and lesser Green’s functions between the QD and the leads are necessary. We can derive these retarded Green’s functions by employing the equation of motion, as follows:(18)G^dαr=G^dσrT^g^α,σr,G^αdr=g^α,σrT^G^dσr.In the above expressions, we used a simple notation for summation over a wave vectork: ∑k,k′G^αk,βk′≡G^αβ. The Langreth theorem states that the lesser Green’s function of a matrix A^B^ is given as A^B^<=A^rB^<+A^<B^a [40]. For equilibrium, by employing the fluctuation dissipation theorem, we obtain G^<ω=fωG^aω−G^rω, where fω denotes the Fermi distribution function. This function is given as fω=1/eωβ0+1 with β0=1/kBT, and the Fermi energy is set as ωF=0. The advanced Green’s function, Gijaω, is calculated by using the retarded Green’s function as G^ijaω=G^jirω†. Thus, we can construct all Green’s functions on the basis of G^dσrω.The Josephson current consists of two components: the current due to the continuous spectrum forω>Δ and that due to the discrete spectrum for ω<Δ:(19)I=Idis+Icon,Idis=e2πℏ∫−ΔΔdωfωjdω,Icon=e2πℏ∫−∞−Δ+∫Δ∞dωfωjcω,where jdcω denotes the discrete (continuous) current density. Under the HF approximation, using the definitions of Σ^0 and Σ^dσHF given in equations (11) and (A.10), respectively, we find that G^dσrω has singular points (poles). We write the determinant of G^dσrω−1 as Dω, which becomes(20)Dω;εd,θL,θR=detG^dσrω−1=ad−bc,where(21)a=ω−εd+2iΓ0ρω−Undσ¯,b=−σiΓ0ΔL+ΔRρ0ω−Udσ¯dσ,c=−σiΓ0Δ¯L+Δ¯Rρ0ω−Udσ†dσ¯†,d=ω+εd+2iΓ0ρω+Undσ,where ρ0ω=ρω/ω. In our previous study [27], we considered the spin-flip effects on the current in an S/QD/S junction with a Josephson junction in the range ω<Δ. For the system, the HF approximation was applied by neglecting the pairing correlation functions dσ†dσ¯† and dσ¯dσ. In the absence of direct tunneling between leads and spin-flip effects in the QD, Dω;εd,θ=detG^ddrω−1 in that paper coincides with equation (20) without pairing correlation functions. Dω given by equation (20) is a quadratic function in terms of ω, and consequently, two poles can occur in G^dσrω. When ρω is given, Dω has a finite imaginary part for ω>Δ. On the contrary, Dω has a real part with an infinitesimal imaginary part for ω<Δ. The ABSs exist inside the gap ω<Δ, and their positions correspond to the poles of Green’s function G^dσrω. The discrete current Idis originates from the ABSs, and to calculate Idis, we employ the Sokhotski–Plemelj formula, limω⟶ω0∫dω/gω=−iπ/g′ω0, where ω0 is a real solution of gω0=0 denoting the position of the ABS. On the contrary, the continuous spectrum outside the gap ω>Δ contributes to the continuous current Icon, which is calculated by means of the usual numerical integral.In this study, we consider the system to be at zero temperature, and all the energy quantities are scaled by the BCS gapΔ. Without loss of generality, we put θR=−θL, and the BCS phase difference between the left and right superconducting leads is set as θ=θL−θR. We set the DOS in the leads in the normal state v=1/2W in the range ω≤W, where we choose the half bandwidth W=20 such that W is much larger than all other energy scales.


## 2.1. Model and Green’s Functions of the QD

We consider a system composed of an S/QD/S junction, where Coulomb interaction exists in the QD. The geometry of the setup is shown in Figure1. The total Hamiltonian of the system is written as(1)H=HL+HR+HD+HT,where(2)Hα=∑k,σεkcα,k,σ†cα,k,σ+∑kΔαcα,k,↑†cα,−k,↓†+H.c.,HD=∑σεddσ†dσ+Und↑nd↓,HT=t∑k,σcL,k,σ†dσ+cR,k,σ†dσ+H.c.,where α=L,R and σ=↑,↓. Hα represents the superconducting lead α; cα,k,σ denotes the annihilation operator of an electron with energy εk, wave number k, and spin σ in the lead; the order parameter Δα=Δeiθα with the BCS gap Δ and the BCS phase θα; HD represents the QD; dσ denotes the annihilation operator of an electron with spin σ; and εd denotes the QD energy level. The occupation number of an electron in the QD with spin σ is defined by ndσ=dσ†dσ, and U represents the Coulomb repulsion between electrons with up- and downspins. HT represents the electron tunneling between the leads and the QD. The coupling strength between electrons in the QD and the leads is defined by Γ0=πνt2, where the tunneling amplitude t is real and v denotes the normal density of states (DOS) at the Fermi level.Figure 1
Schematic diagram of a superconductor/quantum dot/superconductor (S/QD/S) junction.t denotes the tunneling amplitude between the leads and the QD.To describe the above system, we introduce the 2 × 2 Nambu representation. We set the spinor field operators as(3)Ψαkσ=cα,k,σcα,−k,σ¯†,Ψαkσ†=cα,k,σ†cα,−k,σ¯,Ψdσ=dσdσ¯†,Ψdσ†=dσ†dσ¯,where σ¯=↓,↑. Hereinafter, Green’s function written with a hat symbol represents a 2 × 2 matrix. Using the above operators, the retarded (advanced) 2 × 2 matrix Green’s function is defined as G^αkσ,βk′σ′rat,t′=∓iθ±t∓t′Ψαkσt,Ψβk′σ†t′. The lesser Green’s functions are defined as follows:(4)G^αkσ,α′k′σ′<t,t′=icα′,k′,σ′†t′cα,k,σtcα′,−k′,σ¯′t′cα,k,σtcα′,k′,σ′†t′cα,−k,σ¯†tcα′,−k′,σ¯′t′cα,−k,σ¯†t,G^dσ,αkσ<t,t′=icα,k,σ†t′dσtcα,−k,σ¯t′dσtcα,k,σ†t′dσ¯†tcα,−k,σ¯t′dσ¯†t.Employing the equation of motion, we derive all Green’s functions. We denoteg^dr as the retarded Green’s function of the QD in the absence of coupling with the leads, and g^ασr is denoted as the retarded Green’s function of the lead α with spin σ in the absence of coupling with the QD. In the Nambu representation, these functions are written as(5)g^drω−1=ω−εd+i0+00ω+εd+i0+,g^ασrω=−iπνρω1σΔαωσΔ¯αω1,respectively, where a symbol with an overbar represents the complex conjugate of the symbol, and the factor ρω is defined as(6)ρω=ωω2−Δ2,ω>Δ,ωiΔ2−ω2,ω<Δ.It is noteworthy thatρ denotes the ordinary BCS DOS for ω>Δ; however, it is imaginary for ω<Δ.We define Green’s functions of the QD with spinσ using the Zubarev notation as(7)G^dσrω=dσdσ¯†;dσ†dσ¯=dσ;dσ†dσ;dσ¯dσ¯†;dσ†dσ¯†;dσ¯.By using the Dyson equation, we obtain the retarded Green’s functions of the QD as(8)G^dσrω−1=G^dσ0rω−1−Σ^σUω,where Σ^σUω is the self-energy due to the on-site Coulomb interaction. Following references [19, 36], Σ^σUω=UF^dσωG^dσrω−1, where F^dσω is defined as(9)F^dσω=dσndσ¯;dσ†dσndσ¯;dσ¯−dσ¯†ndσ;dσ†−dσ¯†ndσ;dσ¯.In the absence of Coulomb interaction,G^dσ0rω is given by(10)G^dσ0rω−1=g^drω−1−Σ^0ω,where Σ^0 is the noninteracting self-energy, which is calculated by(11)Σ^0ω=T^g^Lσrω+g^RσrωT^,where T^=tσ^z, and σ^z denotes the third component of the 2 × 2 Pauli matrices. Each element of Σ^0 is given in Appendix.In the presence of Coulomb interaction, it is necessary to calculateΣ^σU self-consistently. In this study, we derive Green’s functions, G^dσr, under the Lacroix approximation. Although the Lacroix truncation was originally proposed for the Anderson model with normal conducting leads [29], we extend the method to superconducting leads. The detailed derivations are given in Appendix. The final expression of the (11)-component of G^dσr is(12)G^dσrω11=1Rσω−UQσω+Tσω/Pσω,where(13)Pσω=ω−εd−U1−ndσ¯+3iΓ0+UA1σω−A2σω,Qσω=ω−εd−Undσ¯A1σω−A2σω+Undσ¯1−ndσ¯−dσ¯dσdσ†dσ¯†−2iΓ0ndσ¯−B1σω+B2σω,Rσω=ω−εd−Σ110ω−Undσ¯−Σ120ω+Udσ¯dσΣ210ω+Udσ†dσ¯†ω+εd−Σ220ω+Undσ,Tσω=2εd+U1+ndσ−ndσ¯×dσ¯dσΣ210ω+Udσ†dσ¯†ω+εd−Σ220ω+Undσ,and the (21)-component of G^dσr is given by(14)G^dσrω21=Σ210ω+Udσ†dσ¯†ω+εd−Σ220ω+UndσG^dσrω11.We can obtainG^dσrω12 and G^dσrω22 in a similar way, and the final expressions are given in Appendix.


## 2.2. Current and Electron Occupation Number in the QD

The current from the leadL to the lead R is given by the time derivative of the electron occupation number in the lead L. Applying the equation of motion [37, 38] and the definition NL=∑k,σcL,k,σ†cL,k,σ, we obtain the current IL as follows [39]:(15)IL=−edNLdt=ieℏNL,H=ieℏt∑k,σcL,k,σ†dσ−dσ†cL,k,σ=eℏt∑k,σGdσ,Lkσ<t,t−GLkσ,dσ<t,t,where the matrix elements of the lesser Green’s functions are denoted as follows: Gdσ,Lkσ<t,t′=icL,k,σ†t′dσt. Similarly, IR, which is the current from the lead R to the lead L, is obtained. Using the 2 × 2 matrix Green’s functions and their Fourier-transformed forms, we obtain the Josephson current I as(16)I=12IL−IR=12et2πℏ∑k,σ∫dωG^dσ,Lkσ<ω−G^Lkσ,dσ<ω−G^dσ,Rkσ<ω+G^Rkσ,dσ<ω11,where G^ij is the i,j-th element of G^. The occupation number of an electron with spin σ in the QD is given by(17)nd↑=12πi∫dωG^d↑<ω11,nd↓=12πi∫dωG^d↓<ω11.The averaged occupation number in the QD isnd=nd↑+nd↓.To calculate the current (16), the retarded and lesser Green’s functions between the QD and the leads are necessary. We can derive these retarded Green’s functions by employing the equation of motion, as follows:(18)G^dαr=G^dσrT^g^α,σr,G^αdr=g^α,σrT^G^dσr.In the above expressions, we used a simple notation for summation over a wave vectork: ∑k,k′G^αk,βk′≡G^αβ. The Langreth theorem states that the lesser Green’s function of a matrix A^B^ is given as A^B^<=A^rB^<+A^<B^a [40]. For equilibrium, by employing the fluctuation dissipation theorem, we obtain G^<ω=fωG^aω−G^rω, where fω denotes the Fermi distribution function. This function is given as fω=1/eωβ0+1 with β0=1/kBT, and the Fermi energy is set as ωF=0. The advanced Green’s function, Gijaω, is calculated by using the retarded Green’s function as G^ijaω=G^jirω†. Thus, we can construct all Green’s functions on the basis of G^dσrω.The Josephson current consists of two components: the current due to the continuous spectrum forω>Δ and that due to the discrete spectrum for ω<Δ:(19)I=Idis+Icon,Idis=e2πℏ∫−ΔΔdωfωjdω,Icon=e2πℏ∫−∞−Δ+∫Δ∞dωfωjcω,where jdcω denotes the discrete (continuous) current density. Under the HF approximation, using the definitions of Σ^0 and Σ^dσHF given in equations (11) and (A.10), respectively, we find that G^dσrω has singular points (poles). We write the determinant of G^dσrω−1 as Dω, which becomes(20)Dω;εd,θL,θR=detG^dσrω−1=ad−bc,where(21)a=ω−εd+2iΓ0ρω−Undσ¯,b=−σiΓ0ΔL+ΔRρ0ω−Udσ¯dσ,c=−σiΓ0Δ¯L+Δ¯Rρ0ω−Udσ†dσ¯†,d=ω+εd+2iΓ0ρω+Undσ,where ρ0ω=ρω/ω. In our previous study [27], we considered the spin-flip effects on the current in an S/QD/S junction with a Josephson junction in the range ω<Δ. For the system, the HF approximation was applied by neglecting the pairing correlation functions dσ†dσ¯† and dσ¯dσ. In the absence of direct tunneling between leads and spin-flip effects in the QD, Dω;εd,θ=detG^ddrω−1 in that paper coincides with equation (20) without pairing correlation functions. Dω given by equation (20) is a quadratic function in terms of ω, and consequently, two poles can occur in G^dσrω. When ρω is given, Dω has a finite imaginary part for ω>Δ. On the contrary, Dω has a real part with an infinitesimal imaginary part for ω<Δ. The ABSs exist inside the gap ω<Δ, and their positions correspond to the poles of Green’s function G^dσrω. The discrete current Idis originates from the ABSs, and to calculate Idis, we employ the Sokhotski–Plemelj formula, limω⟶ω0∫dω/gω=−iπ/g′ω0, where ω0 is a real solution of gω0=0 denoting the position of the ABS. On the contrary, the continuous spectrum outside the gap ω>Δ contributes to the continuous current Icon, which is calculated by means of the usual numerical integral.In this study, we consider the system to be at zero temperature, and all the energy quantities are scaled by the BCS gapΔ. Without loss of generality, we put θR=−θL, and the BCS phase difference between the left and right superconducting leads is set as θ=θL−θR. We set the DOS in the leads in the normal state v=1/2W in the range ω≤W, where we choose the half bandwidth W=20 such that W is much larger than all other energy scales.


## 3. Numerical Calculation Results

Here, we show the numerical calculation results obtained by the Lacroix approximation. We define the total density of states in the QD asρdω=−1/π∑σImG^dσrω11, and ρd is shown in Figure 2. In our calculation, the multiple integral over ω is divided into two regions: ω>Δ and ω<Δ. ρω defined by equation (6) diverges at ω=Δ; therefore, a sharp peak and dip are observed around ω=1 in ρd. In the case of the particle-hole symmetric case, εd=−U/2, we observe the Coulomb peaks around ω=±U/2 (Figure 2(a)). However, it is known that the Kondo peak at the Fermi level (ω=0) does not appear in the Lacroix approximation. On the contrary, in the case of the particle-hole asymmetric case, we observe the Coulomb peaks around ω=εd and εd+U, and the antiresonant dip is observed around ω=U+2εd (Figure 2(b)). This unphysical resonance is a specific characteristic of the Lacroix approximation, which disappears in the limit U⟶∞ [29, 41]. In Figure 2(c), the Kondo peak is observed around the Fermi level.Figure 2
Total local density of states in the QD forΓ0=0.5 and U = 12. (a)εd=−6. (b)εd=−5. (c)εd=−4.
(a)
(b)
(c)The dependence of the current on the BCS phase difference is shown in Figure3. In the case of the HF approximation, the current I is expressed in the form I=Icsinθ, where Ic is called the critical current and θ is the BCS phase difference. In the HF approximation, the amplitude Ic depends slightly on θ, and the phase shift occurs when εd is chosen in the range −U<εd<0 (Figure 3(a)). That is, the current-phase relations change from I=Icsinθ to I=−Icsinθ with changing εd, which is called the 0-π transition [42]. A similar property is observed in the Lacroix approximation (Figure 3(b)). However, a difference is that the current amplitude Ic depends strongly on the phase θ. This is caused by different calculation schemes employed by the two approximations; the HF approximation is a type of mean-field theory, taking only the nonconnected part in F^dσ into consideration, and the macroscopic parameters ndσ and dσ†dσ¯† are determined self-consistently. On the contrary, in the Lacroix approximation, by taking the connected term F^dσc into consideration, these macroscopic parameters are self-consistently determined for each θ. As a result, these parameters are dependent on θ. Although the current satisfies I=0 at θ=0,±π, it deviates from the relation I=Icsinθ with a fixed Ic.Figure 3
Dependence of the current on the BCS phase difference forΓ0=0.1 and U = 5. (a) HF approximation: the solid, dashed, and dotted curves correspond to cases εd = 1, −3, and −7, respectively. (b) Lacroix approximation: the solid, dashed, and dotted curves correspond to cases εd = 0, −4, and −7.5, respectively. (c) Lacroix approximation: the solid and dashed curves correspond to cases εd = −5.5 and −6.0, respectively.
(a)
(b)
(c)Let us consider the 0-π transition in terms of the Kondo temperature and the BCS gap under a large Coulomb repulsion. The Kondo temperature TK in the normal conductor connected to a QD is given by several calculation schemes. TK calculated by scaling theory is TK=ΓU/2Γexpπεd1+εd/U/2Γ, where Γ=2πνt2 [43]. For the SUN Anderson impurity model, TK calculated by the renormalization group scaling theory is TK=expπεd/NΓ [44]. However, for the same model, TK calculated by the Lacroix approximation is TK=Wexpπεd/N−1Γ [29, 45]. Although the formulation by the Lacroix truncation scheme is rather crude, it can qualitatively capture the 0-π transition in the S/QD/S system, where the Kondo effect and superconductivity compete against each other [35]. The relative magnitude of Δ and TK determines the characteristics of the dependence of the current on the BCS phase difference. In the weak coupling case, Δ ≫ TK, the π-junction is observed, and the ground state of the QD is doublet (dashed curve in Figure 3(b)). On the contrary, in the strong coupling case, TK≫Δ, the 0-junction is observed, and the ground state of the QD is the Kondo singlet (solid curve in Figure 3(b)). Comparing the results obtained by our Lacroix approximation with those obtained by the HF approximation (Figure 3(a)), the dependence of the current on the BCS phase difference θ is highly nonsinusoidal. For the intermediate coupling case, TK∼Δ, and the dependence of the current breaks into three different regions. In the central region around θ=0, the behavior of the current resembles that in the ballistic short junction: the current changes linearly with θ. In the surrounding two regions, the dependence of the current on the phase is similar to that in a π-junction (Figure 3(c)). At the transition point TK∼Δ, we obtain εd∼N−1ΓlnΔ/W/π by using the Lacroix expression. In our parameters with N = 2, εd∼−0.2. This value does not agree with our numerical calculation result, εd∼−5.5. In spite of this discrepancy, the above-mentioned properties of dependence of current on the BCS phase difference are qualitatively consistent with the results obtained by the NRG method [13]. In our calculation of HF approximation, a complex dependence of the current on the phase, such as that in Figure 3(c), did not occur around the transition point TK∼Δ.The dependence of the averaged occupation number of the electrons in the QD on the QD energy level is shown in Figure4. For a much smaller εd such as εd≪−U, nd∼2. Intermediate values of εd such as −U<εd<0, where nd∼1, define a magnetic region in which a π-junction occurs. For larger values of εd satisfying εd>0, nd∼0. This property is identical for both the HF and the Lacroix approximations.Figure 4
Dependence of the averaged occupation number in the QD on the QD energy level forΓ0 = 0.1 and θ = π/2. The solid, dashed, and dotted curves correspond to cases U = 5, 8, and 10, respectively. (a) HF approximation. (b) Lacroix approximation.
(a)
(b)The current density distribution under the HF approximation is shown in Figure5(a). The total current density is defined as jω=jdω+jcω. When εd is chosen such that εd∼−U, we can find sharp peaks within the BCS gap Δ in the current density distribution j. The positions of these peaks correspond to the zero of D. There is an ABS below the Fermi level, and it carries a negative current. Although jc>0, it is generally small and satisfies jc≪jd so that j < 0 (π-junction) [46]. The dependence of the pole positions of Green’s function G^dσrω11 on εd is shown in Figure 5(b). We defined ω±0 as the real solutions of Dω=0 within the BCS gap, ω±0<Δ. Here, ω±0 corresponds to the position of the ABS, and we find that ω−0≃−ω+0 is satisfied [47]. On the contrary, in the Lacroix approximation, G^dσrω11 is given by equation (12), which is described by the variables Pσ, Qσ, Rσ, and Tσ. If Pσ, Qσ, and Tσ are neglected, G^dσrω11=1/Rσω, which agrees with the result obtained from the HF approximation [27]. When the connected term F^dσc is considered, Pσ, Qσ, and Tσ appear in the expression of G^dσrω11. We see that Pσ and Qσ have finite imaginary parts due to the terms proportional to iΓ0, and Rσ and Tσ have a similar property: for ω>Δ, Rσ and Tσ have finite imaginary parts. On the contrary, for ω<Δ, Rσ and Tσ have real parts with an infinitesimal imaginary part. In our Lacroix approximation, we define Dω as D=PσRσ−UQσ+Tσ. Then, unlike the HF approximation, D has a finite imaginary part for any value of ω, and there is no real solution for D = 0 in the range ω<Δ.Figure 5
Current density and dependence of the pole position on the QD energy level. The data are obtained by the HF approximation forΓ0 = 0.1 and θ = π/2. (a) Current density. The solid and dashed curves correspond to cases U = 3 with εd=−3 and εd=−4, respectively. (b) Dependence of the pole positions of G^dσ0r on the QD energy level. The solid curves with symbols •, ∘, and ∗ correspond to cases U = 0, 3, and 5, respectively.
(a)
(b)The dependence of the critical current on the QD energy level is shown in Figure6. In the case U = 0, Ic is positive, and it diverges at εd=0 (Figure 6(a)). For a finite U, we apply the HF or Lacroix approximations. Generally, for U<Γ0, nd↑ and nd↓ are almost equal, and the current is positive with a maximum at εd=−U/2. By contrast, for U>Γ0, nd↑ and nd↓ are not equal, and therefore, the QD becomes magnetic. In the HF approximation, the π-junction occurs in the range −U<εd<0 (Figure 6(b)). By contrast, in our Lacroix approximation, although the π-junction occurs around εd∼−U or 0, it does not occur in the entire range −U<εd<0 (Figure 6(c)). Current I is composed of Idis and Icon, and the dependence of Idis and Icon on εd is shown in Figure 6(d). Because there is no real solution for D = 0, Idis and Icon were calculated by the usual integral. We see that Idis and Icon have similar properties; Idis and Icon are negative around εd∼−U and Δ, but Idis and Icon are positive with a sharp peak around εd∼−U+Δ and 0. Thus, in our calculations, even if there is no real solution for D = 0 in the range −U<εd<0, a negative current (critical current Ic<0) occurs around εd∼−U and Δ.Figure 6
Dependence of the critical current on the QD energy level forΓ0=0.1 and θ=π/2. (a) U = 0. (b) HF approximation. (c) Lacroix approximation. The solid, dashed, and dotted curves correspond to cases U = 5, 8, and 10, respectively. (d) Component of Ic. U = 8 and the data are obtained by the Lacroix approximation. The solid, dashed, and dotted curves represent Ic, Idis, and Icon, respectively.
(a)
(b)
(c)
(d)


## 4. Conclusions

We have studied the current in an S/QD/S junction system with Coulomb interaction. In the HF approximation, up to the first order of tunneling amplitudet was considered, while the electron correlation was not considered sufficiently; therefore, it was necessary to include a higher-order Green’s function. In this study, the Lacroix approximation was applied, where up to the second order of t was considered. In the HF approximation, the connected term F^dσc was neglected in the cluster expansion, and the macroscopic parameters ndσ, dσ†dσ¯†, and dσ¯dσ were self-consistently determined. It turned out that this simplification corresponded to the neglect of dσndσ¯;dσ†c, cα,k,σndσ¯;dσ†c, dσ¯†cα,k,σ¯dσ;dσ†c, and cα,k,σ¯†dσ¯dσ;dσ†c in the expansion of the higher-order Green’s functions. The lack of these connected terms is a crucial defect when we consider low-temperature physics such as Kondo physics. In our Lacroix approximation, the connected term F^dσc was taken into consideration under several simplifications; the superconducting correlation functions involved in the QD and leads and all the pairing corrections obtained from higher-order Green’s functions were neglected. Under these simplifications, we calculated dσ;dσ† and dσ¯†;dσ¯ self-consistently.In our Lacroix approximation,dσ;dσ† was expressed by the variables Pσ, Qσ, Rσ, and Tσ. When the connected term F^dσc was not considered, Pσ, Qσ, and Tσ vanished, giving dσ;dσ†=1/Rσ, which recovered the expression obtained by the HF approximation [27]. The denominator of dσ;dσ† was defined as Dω. In the HF approximation, when we chose εd∼−U or 0 for a finite U, there were two real solutions for D = 0 in the range ω<Δ. One solution below the Fermi level yielded a negative current, and the π-junction occurred in the range −U<εd<0. On the contrary, in our Lacroix approximation, there was no real solution for D = 0, which results from our calculation scheme based on the Lacroix approximation; the variables Pσ and Qσ were complex for the entire range of ω because of the terms proportional to iΓ0. As a result, although the π-junction did not occur in the entire range −U<εd<0, the negative current occurred only around εd∼−U and Δ.In the HF approximation, the higher-order electron correlation is inherently not taken into consideration, and therefore, a discussion of the characteristics of the current in terms of the ratioΔ/TK is inappropriate. Although higher-order Green’s functions are taken into consideration under simplifications in the Lacroix approximation, it is not a precise calculation compared with the NRG, functional RG, and QMC methods [48]. In fact, our numerical calculation results did not agree well with TK estimated by Lacroix. In spite of this discrepancy, our numerical calculation results captured several aspects of the 0-π transition and the Kondo resonance in restricted parameters. However, our Lacroix approximation could not capture well the competition between Kondo physics and π-junction for all the parameters. The Lacroix approximation employs a truncation at the second order of t and applies cluster expansion for higher-order Green’s functions. Under these simplifications, reliable calculation results are obtained in restricted parameters, and only the qualitative features of Kondo physics in superconductors are obtained.



---

*Source: 5639487-2019-11-18.xml*