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.

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 [14]. 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 , where 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 [68].

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 and the BCS gap [9]. When , 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 , 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 . The tuning of the π-junction in S/QD/S systems has been studied extensively [3, 1214].

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 in the absence of Coulomb interaction. In particular, when the energy level of the QD coincides with the Fermi level , . 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 [1719], noncrossing approximation [20], quantum Monte Carlo (QMC) method [2123], 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 [2427]. 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, 3034], 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 Figure 1. The total Hamiltonian of the system is written aswherewhere and . represents the superconducting lead α; denotes the annihilation operator of an electron with energy , wave number k, and spin σ in the lead; the order parameter with the BCS gap and the BCS phase ; represents the QD; denotes the annihilation operator of an electron with spin σ; and denotes the QD energy level. The occupation number of an electron in the QD with spin σ is defined by , and U represents the Coulomb repulsion between electrons with up- and downspins. represents the electron tunneling between the leads and the QD. The coupling strength between electrons in the QD and the leads is defined by , where the tunneling amplitude t is real and 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 aswhere . 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 . The lesser Green’s functions are defined as follows:

Employing the equation of motion, we derive all Green’s functions. We denote as the retarded Green’s function of the QD in the absence of coupling with the leads, and 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 asrespectively, 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 aswhere is the self-energy due to the on-site Coulomb interaction. Following references [19, 36], , where is defined as

In the absence of Coulomb interaction, is given bywhere is the noninteracting self-energy, which is calculated bywhere , and denotes the third component of the 2 × 2 Pauli matrices. Each element of is given in Appendix.

In the presence of Coulomb interaction, it is necessary to calculate self-consistently. In this study, we derive Green’s functions, , 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 iswhereand the (21)-component of is given by

We can obtain and 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 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 , we obtain the current as follows [39]:where the matrix elements of the lesser Green’s functions are denoted as follows: . Similarly, , 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 aswhere is the -th element of . The occupation number of an electron with spin σ in the QD is given by

The 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: . The Langreth theorem states that the lesser Green’s function of a matrix is given as [40]. For equilibrium, by employing the fluctuation dissipation theorem, we obtain , where denotes the Fermi distribution function. This function is given as with , and the Fermi energy is set as . The advanced Green’s function, , is calculated by using the retarded Green’s function as . Thus, we can construct all Green’s functions on the basis of .

The Josephson current consists of two components: the current due to the continuous spectrum for and that due to the discrete spectrum for :where denotes the discrete (continuous) current density. Under the HF approximation, using the definitions of and given in equations (11) and (A.10), respectively, we find that has singular points (poles). We write the determinant of as , which becomeswherewhere . 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 and . In the absence of direct tunneling between leads and spin-flip effects in the QD, in that paper coincides with equation (20) without pairing correlation functions. given by equation (20) is a quadratic function in terms of ω, and consequently, two poles can occur in . When is given, has a finite imaginary part for . On the contrary, 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 . The discrete current originates from the ABSs, and to calculate , we employ the Sokhotski–Plemelj formula, , where is a real solution of denoting the position of the ABS. On the contrary, the continuous spectrum outside the gap contributes to the continuous current , 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 , and the BCS phase difference between the left and right superconducting leads is set as . We set the DOS in the leads in the normal state in the range , where we choose the half bandwidth such that 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 , and 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 in . In the case of the particle-hole symmetric case, , we observe the Coulomb peaks around (Figure 2(a)). However, it is known that the Kondo peak at the Fermi level () 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 and , and the antiresonant dip is observed around (Figure 2(b)). This unphysical resonance is a specific characteristic of the Lacroix approximation, which disappears in the limit [29, 41]. In Figure 2(c), the Kondo peak is observed around the Fermi level.

The 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 , where is called the critical current and θ is the BCS phase difference. In the HF approximation, the amplitude depends slightly on θ, and the phase shift occurs when is chosen in the range (Figure 3(a)). That is, the current-phase relations change from to with changing , 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 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 into consideration, and the macroscopic parameters and are determined self-consistently. On the contrary, in the Lacroix approximation, by taking the connected term into consideration, these macroscopic parameters are self-consistently determined for each θ. As a result, these parameters are dependent on θ. Although the current satisfies at , it deviates from the relation with a fixed .

Let us consider the 0-π transition in terms of the Kondo temperature and the BCS gap under a large Coulomb repulsion. The Kondo temperature in the normal conductor connected to a QD is given by several calculation schemes. calculated by scaling theory is , where [43]. For the Anderson impurity model, calculated by the renormalization group scaling theory is [44]. However, for the same model, calculated by the Lacroix approximation is [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 determines the characteristics of the dependence of the current on the BCS phase difference. In the weak coupling case, , 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, , 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, , and the dependence of the current breaks into three different regions. In the central region around , 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 , we obtain by using the Lacroix expression. In our parameters with N = 2, . This value does not agree with our numerical calculation result, . 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 .

The 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 such as , . Intermediate values of such as , where , define a magnetic region in which a π-junction occurs. For larger values of satisfying , . This property is identical for both the HF and the Lacroix approximations.

The current density distribution under the HF approximation is shown in Figure 5(a). The total current density is defined as . When is chosen such that , 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 , it is generally small and satisfies so that j < 0 (π-junction) [46]. The dependence of the pole positions of Green’s function on is shown in Figure 5(b). We defined as the real solutions of within the BCS gap, . Here, corresponds to the position of the ABS, and we find that is satisfied [47]. On the contrary, in the Lacroix approximation, is given by equation (12), which is described by the variables , , , and . If , , and are neglected, , which agrees with the result obtained from the HF approximation [27]. When the connected term is considered, , , and appear in the expression of . We see that and have finite imaginary parts due to the terms proportional to , and and have a similar property: for , and have finite imaginary parts. On the contrary, for , and have real parts with an infinitesimal imaginary part. In our Lacroix approximation, we define as . 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 .

The dependence of the critical current on the QD energy level is shown in Figure 6. In the case U = 0, is positive, and it diverges at (Figure 6(a)). For a finite U, we apply the HF or Lacroix approximations. Generally, for , and are almost equal, and the current is positive with a maximum at . By contrast, for , and are not equal, and therefore, the QD becomes magnetic. In the HF approximation, the π-junction occurs in the range (Figure 6(b)). By contrast, in our Lacroix approximation, although the π-junction occurs around or 0, it does not occur in the entire range (Figure 6(c)). Current I is composed of and , and the dependence of and on is shown in Figure 6(d). Because there is no real solution for D = 0, and were calculated by the usual integral. We see that and have similar properties; and are negative around and , but and are positive with a sharp peak around and 0. Thus, in our calculations, even if there is no real solution for D = 0 in the range , a negative current (critical current ) occurs around and .

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 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 was neglected in the cluster expansion, and the macroscopic parameters , , and were self-consistently determined. It turned out that this simplification corresponded to the neglect of , , , and 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 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 and self-consistently.

In our Lacroix approximation, was expressed by the variables , , , and . When the connected term was not considered, , , and vanished, giving , which recovered the expression obtained by the HF approximation [27]. The denominator of was defined as . In the HF approximation, when we chose 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 . 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 and were complex for the entire range of ω because of the terms proportional to . As a result, although the π-junction did not occur in the entire range , the negative current occurred only around 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 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 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.

Appendix

A. Derivation of Green’s Function of the QD

In this section, we present the derivation of following reference [49], in which the normal conductor/QD/S junction was considered. We use the Zubarev notation for the retarded Green’s function, , which is composed of the operators A and B. The equation of motion of is given bywhere . We first derive the equations of motion of each element of the retarded Green’s function given by equation (7):where we put , , , and with . 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 , , and , where , , and . Here, and are the Heaviside step function and the sign function, respectively.

To simplify the description of equations (A.2)–(A.5), we introduce a matrix as

Using these expressions, equations (A.2)–(A.5) are written in the matrix form aswhere is an identity matrix andwhere is a noninteracting self-energy with components , , and , where we used the property . It is known that equation (A.6) is not closed because includes the higher-order terms of . To overcome this difficulty, we employ the truncation approximation. We first expand aswhere is the nonconnected part, which is given asand the connected part is

When the connected part is neglected, satisfies the relation

The solution of (A.12) is obtained self-consistently. The occupation number in the QD and the pairing correlation function are given by the relationswhere is the Fermi distribution function with the Fermi level  0.

Beyond the HF approximation, one needs to take into consideration. For the study of Kondo physics, we need to consider only the contributions from the diagonal components in equation (A.11) [50, 51]. Here, we apply the Lacroix truncation, where the second order of connection is considered. In this approximation, although the diagonal components of are considered, anomalous higher-order Green’s functions such as and are neglected. This yields the relationand equation (A.2) becomes

All we have to do is to calculate self-consistently. The procedure is as follows. The equation of motion of is given by

The higher-order Green’s functions, , , and , are included in the right-hand side of equation (A.16). The equations of motion of these higher-order Green’s functions arewhere and . To obtain , we apply cluster expansion under spin conservation to Green’s functions in equations (A.17)–(A.19). That is,

On the contrary, we neglect all the superconducting correlation functions and pairing correlations obtained from higher-order Green’s functions, such as , , and . Applying these cluster expansions and approximations to (A.16), the final equation composed of the connected part of Green’s functions is

To describe , , and in terms of , we apply the expansionswhere

Substituting the approximations of equations (A.25)–(A.27) into equation (A.24), we obtain . Substituting this result into equation (A.15) with the help of equation (A.14), we finally obtain aswhere

In the above, we defined the variables aswhere we put and used the spectral theorem and . The expression of is given bywhere , in which and denote the retarded and advanced Green’s functions, respectively. Using these expressions, the summation over k in equation (A.31) is given bywhere ε and . Similarly, we defined as

The expression of is approximately given aswith

In the numerical calculation of the multiple integral in and , we use the following formulas:

In this Lacroix approximation, the higher-order Green’s functions in equation (A.16) are expressed in terms of and . Eliminating , we calculate and self-consistently. In the formulation, and are taken into consideration through and , respectively. These terms diverge logarithmically at the Fermi level at zero temperature and play an important role in the Kondo effect [29].

Following similar procedures, we can calculate and , and the final expressions arewherewhere we defined , , and aswhere , , , and . In the numerical calculation of the multiple integral in and , we use the following formulas:

Data Availability

The data used to support the findings of this study are included within the article.

Conflicts of Interest

The author declares that there are no conflicts of interest.