Josephson Current through a Quantum Dot Connected with Superconducting Leads

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.


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][2][3][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. e current is expressed as J � J c sin θ, where J c 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][7][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 temperature T K and the BCS gap Δ [9]. When T K ≫ Δ, 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. is results in a positive critical current (0-junction). By contrast, when T K ≪ Δ, 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. e spin ordering of the Cooper pair is reversed, resulting in a π-junction [3,10,11]. e 0-π transition is expected to occur around T K ∼ Δ. e tuning of the π-junction in S/QD/S systems has been studied extensively [3,[12][13][14].
e current density in the S/QD/S system is composed of two parts: continuous and discrete spectrums. e former (latter) arises from outside (inside) the BCS gap. e 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. e discrete current density arises from the Andreev bound states (ABSs) inside the BCS gap. e 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.
e Coulomb interaction in the QD is studied by many methods: the numerical renormalization group (NRG) method [17][18][19], noncrossing approximation [20], quantum Monte Carlo (QMC) method [21][22][23], and so on. Although the NRG and QMC methods are very precise methods, they are computationally expensive. e simplest method is the Hartree approximation, which corresponds to the zerothorder 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. is method was applied to the level-crossing quantum phase transition between the BCS-singlet and the magnetic doublet states [24][25][26][27]. e above two approximations can be applied to describe singleparticle physics. e 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][31][32][33][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.

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.

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. e geometry of the setup is shown in Figure 1. e total Hamiltonian of the system is written as where 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 Δ α � Δe iθ α with the BCS gap Δ and the BCS phase θ α ; H D represents the QD; d σ denotes the annihilation operator of an electron with spin σ; and ε d denotes the QD energy level. e occupation number of an electron in the QD with spin σ is defined by n dσ � d † σ d σ , and U represents the Coulomb repulsion between electrons with up-and downspins. H T represents the electron tunneling between the leads and the QD. e coupling strength between electrons in the QD and the leads is defined by Γ 0 � π]t 2 , where the tunneling amplitude t is real and v denotes the normal density of states (DOS) at the Fermi level.
To describe the above system, we introduce the 2 × 2 Nambu representation. We set the spinor field operators as 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 e lesser Green's functions are defined as follows: 2 Advances in Condensed Matter Physics Employing the equation of motion, we derive all Green's functions. We denote g r d 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 respectively, where a symbol with an overbar represents the complex conjugate of the symbol, and the factor ρ(ω) is defined as 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 By using the Dyson equation, we obtain the retarded Green's functions of the QD as where Σ U σ (ω) is the self-energy due to the on-site Coulomb interaction. Following references [19,36] In the absence of Coulomb interaction, G (0)r dσ (ω) is given by where Σ 0 is the noninteracting self-energy, which is calculated by 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 r dσ , 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. e detailed derivations are given in Appendix. e final expression of the (11) where and the (21)-component of G r dσ is given by We can obtain [G r dσ (ω)] 12 and [G r dσ (ω)] 22 in a similar way, and the final expressions are given in Appendix.

Current and Electron Occupation Number in the QD.
e current from the lead L 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 N L � k,σ c † L,k,σ c L,k,σ , we obtain the current I L as follows [39]: where the matrix elements of the lesser Green's functions are denoted as follows: Similarly, I R , 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 where [G] ij is the i, j-th element of G. e occupation number of an electron with spin σ in the QD is given by e averaged occupation number in the QD is 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: In the above expressions, we used a simple notation for summation over a wave vector k: k,k′ G αk,βk′ ≡ G αβ . e Langreth theorem states that the lesser Green's function of a matrix AB is given as [40]. For equilibrium, by employing the fluctuation dissipation theorem, we obtain G is function is given as f(ω) � 1/(e ωβ 0 + 1) with β 0 � 1/k B T, and the Fermi energy is set as ω F � 0. e advanced Green's function, G a ij (ω), is calculated by using the retarded Green's function as G a ij (ω) � (G r ji (ω)) † . us, we can construct all Green's functions on the basis of G r dσ (ω).
e Josephson current consists of two components: the current due to the continuous spectrum for |ω| > Δ and that due to the discrete spectrum for |ω| < Δ: where j d(c) (ω) denotes the discrete (continuous) current density. Under the HF approximation, using the definitions of Σ 0 and Σ HF dσ given in equations (11) and (A.10), respectively, we find that G r dσ (ω) has singular points (poles). We write the determinant of (G r dσ (ω)) − 1 as D(ω), which becomes 4 Advances in Condensed Matter Physics where 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 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 r dσ (ω). When ρ(ω) is given, D(ω) has a finite imaginary part for |ω| > Δ. On the contrary, D(ω) has a real part with an infinitesimal imaginary part for |ω| < Δ.
e ABSs exist inside the gap |ω| < Δ, and their positions correspond to the poles of Green's function G r dσ (ω). e discrete current I dis originates from the ABSs, and to calculate I dis , 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 I con , 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.

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/π) σ Im[G r dσ (ω)] 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 particlehole 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)). is 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.
e dependence of the current on the BCS phase difference is shown in Figure 3. In the case of the HF approximation, the current I is expressed in the form I � I c sin θ, where I c is called the critical current and θ is the BCS phase difference. In the HF approximation, the amplitude I c depends slightly on θ, and the phase shift occurs when ε d is chosen in the range − U < ε d < 0 (Figure 3(a)). at is, the current-phase relations change from I � |I c |sin θ to I � − |I c |sin θ 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 I c depends strongly on the phase θ.
is 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 〈n dσ 〉 and 〈d † σ d † σ 〉 are determined self-consistently. On the contrary, in the Lacroix approximation, by taking the connected term F c dσ 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 � I c sin θ with a fixed I c .
Let us consider the 0-π transition in terms of the Kondo temperature and the BCS gap under a large Coulomb repulsion. e Kondo temperature T K in the normal conductor connected to a QD is given by several calculation schemes. T K calculated by scaling theory is [43]. For the SU(N) Anderson impurity model, T K calculated by the renormalization group scaling theory is T K � exp[πε d /(NΓ)] [44]. However, for the same model, T K calculated by the Lacroix approximation is T K � W exp[πε 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]. e relative magnitude of Δ and T K determines the characteristics of the dependence of the current on the BCS phase difference. In the weak coupling case, Δ ≫ T K , 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, T K ≫ Δ, 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, T K ∼ Δ, 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 T K ∼ Δ, we obtain ε d ∼ [(N − 1)Γ ln(Δ/W)]/π by using the Lacroix expression. In our parameters with N � 2, ε d ∼ − 0.2.
is value does not agree with our numerical calculation result, ε d ∼ − 5.5. In spite of this discrepancy, the abovementioned 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 T K ∼ Δ. e dependence of the averaged occupation number of the electrons in the QD on the QD energy level is shown in Figure 4. For a much smaller ε d such as Intermediate values of ε d such as − U < ε d < 0, where 〈n d 〉 ∼ 1, define a magnetic region in which a π-junction occurs. For larger values of ε d satisfying ε d > 0, 〈n d 〉 ∼ 0. is property is identical for both the HF and the Lacroix approximations. e current density distribution under the HF approximation is shown in Figure 5(a). e total current density is defined as j(ω) � j d (ω) + j c (ω). When ε d is chosen such that ε d ∼ − U, we can find sharp peaks within the BCS gap Δ in the current density distribution j. e positions of these peaks correspond to the zero of D. ere is an ABS below the Fermi level, and it carries a negative current. Although j c > 0, it is generally small and satisfies j c ≪ |j d | so that j < 0 (π-junction) [46]. e dependence of the pole positions of Green's function [G r dσ (ω)] 11 on ε d is shown in Figure 5(b). We defined ω 0 ± as the real solutions  Advances in Condensed Matter Physics 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 r dσ (ω)] 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 r dσ (ω)] 11 � 1/R σ (ω), which agrees with the result obtained from the HF approximation [27]. When the connected term F c dσ is considered, P σ , Q σ , and T σ appear in the expression of [G r dσ (ω)] 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 σ − U(Q σ + T σ ).
en, 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 |ω| < Δ.
e dependence of the critical current on the QD energy level is shown in Figure 6. In the case U � 0, I c 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 , n d↑ and n d↓ are almost equal, and the current is positive with a maximum at ε d � − U/2. By contrast, for U > Γ 0 , n d↑ and n d↓

Advances in Condensed Matter Physics
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 I dis and I con , and the dependence of I dis and I con on ε d is shown in Figure 6(d). Because there is no real solution for D 0, I dis and I con were calculated by the usual integral. We see that I dis and I con have similar properties; I dis and I con are negative around ε d ∼ − U and Δ, but I dis and I con are positive with a sharp peak around ε d ∼ − (U + Δ) and 0. us, in our calculations, even if there is

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 amplitude t 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 c dσ was neglected in the cluster expansion, and the macroscopic parameters 〈n dσ 〉, 〈d † σ d † σ 〉, and 〈d σ d σ 〉 were self-consistently determined. It turned out that this simplification corresponded to the neglect of 〈〈d σ n dσ ; d † σ 〉〉 c , 〈〈c α,k,σ n dσ ; 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. e 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  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 c dσ was not considered, P σ , Q σ , and T σ vanished, giving 〈〈d σ ; d † σ 〉〉 � 1/R σ , which recovered the expression obtained by the HF approximation [27]. e 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 Δ/T K 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 T K 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. e 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.

A. Derivation of Green's Function of the QD
In this section, we present the derivation of G r dσ (ω) following reference [49], in which the normal conductor/QD/S junction was considered. We use the Zubarev notation for the retarded Green's function, where [A, B] ± � AB ± BA. We first derive the equations of motion of each element of the retarded Green's function G r dσ (ω) given by equation (7): For the symbols σ � ↑, ↓, + and − are, respectively, assigned for the up-and downspins in the mathematical calculations. Replacing the sum over the wavenumber k by an integral, we obtain Here, θ(ω) and sgn(ω) are the Heaviside step function and the sign function, respectively.