Irreducible Green’s functions method for a quantum dot coupled to metallic and superconducting leads

Physica E 89 (2017) 21–28

Contents lists available at ScienceDirect

Physica E journal homepage: www.elsevier.com/locate/physe

Irreducible Green’s functions method for a quantum dot coupled to metallic and superconducting leads

MARK

⁎

Grzegorz Górski , Krzysztof Kucab Faculty of Mathematics and Natural Sciences, University of Rzeszów, ul. Pigonia 1, 35-310 Rzeszów, Poland

A R T I C L E I N F O

A BS T RAC T

Keywords: Quantum dot Equation of motion Green function Andreev transmittance

Using irreducible Green’s functions (IGF) method we analyse the Coulomb interaction dependence of the spectral functions and the transport properties of a quantum dot coupled to isotropic superconductor and metallic leads (SC-QD-N). The irreducible Green’s functions method is the modification of classical equation of motion technique. The IGF scheme is based on differentiation of double-time Green’s functions, both over the primary and secondary times. The IGF method allows to obtain the spectral functions for equilibrium and nonequilibrium impurity Anderson model used for SC-QD-N system. By the numerical computations, we show the change of spectral and the anomalous densities under the influence of the Coulomb interactions. The observed sign change of the anomalous spectral density can be used as the criterion of the SC singlet-Kondo singlet transition.

1. Introduction

state is superconducting-like. For a weak coupling ΓS and strong values of U , the dominant state is the Kondo singlet state [9]. The transition between these two states can be determined experimentally by means of transport spectroscopy. The measurements of the subgap differential conductance give the evidence for the Andreev bound states. Zero bias anomalies mark the quantum phase transition between the doublet and singlet state. The properties of SC-QD-N systems have been studied using various theoretical approaches: Hartree-Fock theory (HF) [3,17], numerical renormalization group (NRG) [11,13,18], iterative perturbation technique (IPT) [8,9], equation of motion (EOM) [19,20] and noncrossing approximation (NCA) [21]. All of these methods have some restrictions. To describe the systems where Coulomb correlations are not very important one can use the Hartree-Fock theory; e.g. Lee and co-workers [3] used this approximation to describe the influence of external magnetic field on the value of the differential conductance dI / dV . The Coulomb correlations can be described by the use of the NRG, IPT and EOM approaches. The NRG method can be used only for equilibrium state, but IPT and EOM approaches can be used to study the local Coulomb interaction effect on the non-equilibrium transport properties. In this work we use the irreducible Green’s functions (IGF) technique [22] to describe the properties of SC-QD-N system. This technique involves the use of modified equation of motion scheme for the double-time temperature Green’s functions. This scheme permits us to construct the relevant dynamic solutions in a self-consistent way,

The Kondo effect in nanoscopic structures (quantum dots based on semiconductors, carbon nanotubes, nanowires) is actually widely studied both theoretically and experimentally. One of the interesting systems is a quantum dot coupled to one superconducting electrode and one metallic electrode (SC-QD-N). This system is intensively studied experimentally [1–5], and theoretically [6–15]. The connection of the superconducting electrode to a quantum dot through proximity effect induces electron pairing in quantum dot and it causes that the ground state of a quantum dot is the superconducting singlet state. The characteristic feature of the QD-SC connection is the appearance of subgap excitations in the quantum dot spectrum, so-called Andreev or Yu-Shiba-Rusinov bound states. These states play a crucial role in the transport properties of mesoscopic superconducting devices, especially in the subgap regime ( ε < Δ). The strong Coulomb repulsion, between opposite spin electrons, opposes double occupancy of the quantum dot and prefers the doublet ground state. In SC-QD-N heterojunction the Kondo-type correlations have to compete with the proximity induced electron pairing. As a result of this competition, the system exhibits a quantum phase transition between the doublet and the singlet states [8,9,11,16]. The transition between the doublet and the singlet states is determined by different energy scales: the Coulomb interaction (U ), the superconducting gap ( Δ) and the coupling to the superconducting lead (ΓS ). In the large superconducting gap limit and strong coupling ΓS , the singlet ⁎

Corresponding author. E-mail address: [email protected] (G. Górski).

http://dx.doi.org/10.1016/j.physe.2017.01.026 Received 27 September 2016; Received in revised form 12 January 2017; Accepted 27 January 2017 Available online 28 January 2017 1386-9477/ © 2017 Elsevier B.V. All rights reserved.

Physica E 89 (2017) 21–28

G. Górski, K. Kucab

⎡ 〈〈nˆ d ; d +〉〉 〈〈nˆ d ; d 〉〉 ⎤ (1) ↑ ε ↓ ε d↓ ↑ d↓ ↑ ⎥. Γˆ d = ⎢ ⎢⎣〈〈nˆd ↑d↓+; d↑+〉〉ε 〈〈nˆd ↑d↓+; d↓〉〉ε ⎥⎦

without decoupling the chain of equations of motion. This technique was used in many issues of solid-state physics, e.g. itinerant antiferromagnetism [23], itinerant [24] and Heisenberg [25] ferromagnetism, superconductivity in disordered transition metal alloys [26]. We also used this technique to describe the quantum dot coupled to two metallic leads (N-QD-N) [27]. The obtained results show that IGFEOM method gives the results for density of states and the differential conductance comparable to the numerical methods results (NRG) both for the particle-hole symmetric case and also for the asymmetric cases. In this work, to describe a quantum dot coupled to one superconducting electrode and one metallic electrode (SC-QD-N) we will use the IGF-EOM method. In Section 2 we analyse the single impurity Anderson model using the IGF-EOM approach. In this approach, we can describe the SC-QD-N system. We also compute the expressions for the self-energy and the Green’s function in the presence of the Coulomb repulsion. In Section 3 we present numerical results for the spectral density and the Andreev transmittance. We analyse the influence of the Coulomb interaction U on the spectral density and the Andreev transmittance. Additionally we study the influence of parameter ΓN , which characterizes the coupling between the quantum dot and the normal lead, on the spectral density. The results of non-equilibrium transport are also shown in Section 3. Final conclusions are given in Section 4.

(6)

(1) Γˆ d

we will apply the technique of irreducible Green’s For the functions [22,26]. These functions are defined as ir

〈〈[A, H ]− ; B〉〉ε = 〈〈[A, H ]− − zA; B〉〉ε ,

(7)

where the z constant is given by z = 〈[[A, H ]− ; B]± 〉/〈[A; B]± 〉 and it represents the self-energy in the Hartree-Fock-Bogoliubov approximation. Using the irreducible Green’s function technique we can express Eq. (6) as the sum of irreducible function and the mean-field solution. As a result the last part of Eq. (4) can be written as (1) (1) ˆ , ˆ HF G U τˆ 3Γˆ d = U τˆ 3ir Γˆ d + Σ U d

(8)

where the interaction part of self-energy,

ˆ HF , Σ U

is given by

⎡ U 〈nd ↓〉 U 〈d↓d↑〉 ⎤ ˆ HF = ⎢ ⎥. Σ U ⎢⎣U 〈d↑+d↓+〉 −U 〈nd ↑〉⎥⎦

(9)

Inserting Eqs. (8), (A.5) and (A.8) into Eq. (4) we obtain the following relations

ˆ = Iˆ + U τˆ ir Γˆ (1), gˆ dHF −1G d 3 d gˆ dHF −1

(10)

ˆ HF . Σ U

gˆ −1 d

ir

(1) Γˆ d

= − where Neglecting the irreducible function in Eq. (10) we obtain the well-known Hartree-Fock-Bogoliubov approximation which is widely applied to the systems with quantum dot, for which the Coulomb correlations are less important (e.g. for metallic dots, see [15,28]). Because we are interested in the influence of the Coulomb correlations on the transport properties, we have to calculate ir (1) the irreducible function Γˆ .

2. The model Using the Anderson-type Hamiltonian we analyse the system which is built out of the quantum dot connected to one metallic lead and one superconducting lead. The Hamiltonian of this model has the following form

d

H=

In the previous papers where the classic EOM approach in the SIAM model was used (e.g. [29–31]), the higher order Green’s (1) functions Γˆ was calculated by reusing of Eq. (3). Such approach

∑ εd ndσ + Und ↑nd ↓ + ∑ (εkα − μα )nkασ + ∑ (Vkαdσ+ckασ + h. c.) σ

kσα

− Δ∑

+ + (ckS ↑c−kS ↓

kσα

+ c−kS ↓ckS ↑) ,

d

(1)

where dσ+(dσ ) are creation (annihilation) operators for the dot electron with spin σ , ck+ασ (ckασ ), α = N , S are creation (annihilation) operators for the electron in the normal (N) and superconducting (S) lead, εkα is the energy dispersion of α lead, μα is the chemical potential of α lead, εd is the dot energy, U is the on-site Coulomb interaction between electrons on the dot, and Vkα is the coupling between the α lead and the dot. We are looking for expressions for the matrix Green’s function in the Nambu space defined as

allows to obtain Abrikosov-Suhl resonance outside the particle-hole symmetric system. For nd = 1 the Abrikosov-Suhl peak disappears. The classic EOM approach does not fulfil the unitary limit for conductance in the particle–hole symmetry case. Out of the particle–hole symmetric case one obtain a narrow Abrikosov-Suhl resonance peak, whose height and width are small, resulting in an underestimation of Kondo temperature. In order to correct this defect of the classic EOM approach, by ir (1) calculating the Γˆ function we will use the extended equation of

⎡ 〈〈d (t ); d +(t′)〉〉 〈〈d (t ); d (t′)〉〉 ⎤ ↑ ↑ ↓ ˆ (t , t′) = ⎢ ↑ ⎥. G d ⎢⎣〈〈d↓+(t ); d↑+(t′)〉〉 〈〈d↓+(t ); d↓(t′)〉〉⎥⎦

motion approach based on differentiating Green’s function over the second time variable (t′). In the energy representation this leads to the equation [22]

k

d

(2)

−ε⟨⟨A; B⟩⟩ε = − [A, B]+ +

In our analysis, we will use the Green’s function method and the equation of motion (EOM) technique. In general, the EOM for Green’s functions is obtained by differentiation with respect to primary time (t ). After taking the Fourier transform, we obtain the following equation

ε⟨⟨A; B⟩⟩ε = [A, B]+ + ⟨⟨[A, H ]− ; B⟩⟩ε .

A; [B, H ]−

ε.

(11)

Using this equation we obtain the following relation for the irreducible Green’s function

−εir ⟨⟨A; B⟩⟩ε =ir ⟨⟨A; [B, H ]− ⟩⟩ε ,

(3)

(12) ir

(1) which gives the following expression for the Γˆ d function

Applying Eq. (3) to the Hamiltonian (1) we obtain the relation (see Appendix A)

ir

(1) Γˆ d gˆ 0−1 = d

∑

ir

(1) ˆ * τˆ + U ir Γˆ (2)τˆ , Γˆ d kαV d 3 kα 3

kσα

ˆ = Iˆ + U τˆ Γˆ (1), gˆ d −1G d 3 d

(4)

where

ir

⎡ε − ε − Σ 0 ⎤ −Σd012 d d11 ⎥ gˆ d −1 = ⎢ 0 0 ⎢⎣ −Σd 21 ε + εd − Σd 22 ⎥⎦

(13)

where

(5)

⎡ ir〈〈nˆ d ; c + 〉〉 (1) d ↓ ↑ kα ↑ ε Γˆ d kα = ⎢ ir ⎢⎣ 〈〈nˆd ↑d↓+; ck+α ↑〉〉ε

〈〈nˆd ↓d↑; ckα ↓〉〉ε ⎤ ⎥ 〈〈nˆd ↑d↓+; ckα ↓〉〉ε ⎥⎦

(14)

〈〈nˆd ↓d↑; nˆd ↑d↓〉〉ε ⎤ ⎥. ir 〈〈nˆd ↑d↓+; nˆd ↑d↓〉〉ε ⎥⎦

(15)

ir ir

and

ˆ0 Σ d

given by Eq. (A.8). These self-energies come from with self-energies the coupling between dot and the normal or superconducting lead. The (1) two-particle, higher order Green’s functions Γˆ are given by

ir

d

22

⎡ ir〈〈nˆ d ; n d +〉〉 (2) d↓ ↑ d↓ ↑ ε Γˆ d = ⎢ ir ⎢⎣ 〈〈nˆd ↑d↓+; nd ↓d↑+〉〉ε

ir

Physica E 89 (2017) 21–28

G. Górski, K. Kucab

Applying our extended version of EOM (Eq. (11)) to the function we obtain ir

ir

(1) Γˆ dkα

1

0.75

1

(1) ir (1) ˆ gˆ . Γˆ d kα = Γˆ d τˆ 3V kα kα

(16)

0.75

(1) ir (1) ˆ gˆ V ˆ * ˆ = U ir Γˆ (2)τˆ . Γˆ d gˆ 0−1 − Γˆ d ∑ τˆ 3V kα kα kατ 3 d 3 d

πΓNρd/2

Inserting Eq. (16) into Eq. (13) we have ir

1

Δ=5

(17)

kσα

0.5

0.5

0

0.25 0

0.5

-0.2 -0.1 0 0.1 0.2

Using the irreducible Green’s functions method again we obtain ir

(1) Γˆ d

ir (2)ir = U Γˆ d τˆ 3gˆ dHF

0.25

(18)

and

ˆ = Iˆ + U 2 τˆ ir Γˆ (2)ir τˆ gˆ HF , gˆ dHF −1G d 3 d 3 d

0

(19)

-6

-4

-2

0 ε

where ir

⎡ ir〈〈nˆ d ; n d +〉〉ir (2)ir d↓ ↑ d↓ ↑ ε Γˆ d = ⎢ ir ⎢⎣ 〈〈nˆd ↑d↓+; nd ↓d↑+〉〉irε

2

4

6

Fig. 1. The normalized spectral density ρd (ε ) of correlated quantum dot as a function of

〈〈nˆd ↓d↑; nˆd ↑d↓〉〉irε ⎤ ⎥. ir 〈〈nˆd ↑d↓+; nˆd ↑d↓〉〉irε ⎥⎦ ir

energy ε for different values of superconducting energy gap Δ obtained for ΓS = ΓN = 1 and U = 4 .

(20)

(2)ir

The second order self-energy, U 2 τˆ ir3 Γˆ d τˆ 3, appearing in Eq. (19), ˆ (2). will be replaced further by the symbol Σ

(see [9]) and we repeat steps 3-6 until the convergence of the selfenergy and the expectation values is achieved.

U

Next, we will use the Dyson relations

ˆ = gˆ 0 + gˆ 0 Σ ˆ . ˆ G G d d d U d

3. Numerical results

(21)

ˆ =Σ ˆ HF + Σ ˆ ′ , into the first order part We split the self-energy Σ U U U ′ ˆ ). Next, we use Eqs. (19) and (21) and the higher order part (Σ U obtaining

In our numerical analysis, we will concentrate on the equilibrium system where μS = μN = 0 . The energy unit will be the coupling parameter between the dot and the superconducting lead, ΓS = 1. We start with the symmetric case, where the energy of quantum dot level is equal to εd = − U /2 . In Fig. 1 we present the energy dependence of the spectral density of correlated quantum dot, ρd (ε ) = − ImG11(ε )/ π , for different values of superconducting energy gap Δ in the symmetric coupling case, ΓS = ΓN = 1. For Δ = 0 the superconducting lead becomes the normal lead and we have the system N−QD−N . In this case, at the energy ε = 0 , there should be the Abrikosov-Suhl peak for which ρd (0) is equal to the corresponding DOS of the uncorrelated quantum dot. This is the Kondo singlet state. At Δ ≠ 0 there is an increase of the DOS at ε = 0 and the Abrikosov-Suhl peak is narrowing. The narrowing of the Abrikosov-Suhl resonance gives rise to a lowering of the Kondo temperature with respect to the normal case [7]. At ε = ± Δ there are formed the gap edge singularities. Because ΓS = ΓN = 1 and U > ΓS the Kondo singlet state still dominates (see [9]). In Fig. 2 we present the spectral density ρd (ε ) = − ImG11(ε )/ π for different values of the Coulomb interaction U . To avoid the edge effects we analyse the subgap energy regime ( Δ > > 0 ). At small Coulomb interactions (U < ΓS ) we obtain two Andreev quasiparticle peaks (see dashed line in Fig. 2). The increase of the Coulomb interaction causes the Andreev resonances get closer to the Fermi level and eventually

ˆ HF ) (Σ U

−1

ˆ ′ = [Iˆ + Σ ˆ (2)gˆ HF ] Σ ˆ (2). Σ U U d U

(22)

ˆ ′ , we have to know the To calculate the higher order self-energy, Σ U ir (2)ir functions (see Appendix B). After equating the form of the Γˆ d −1

ˆ (2)gˆ HF ] [Iˆ + Σ U d

term to unit matrix we obtain the second order (2)

ˆ′ ≈ Σ ˆ , which does not describe the perturbation theory result Σ U U strong Coulomb correlation and concentrations away from the halffilling. It should be noted here that using this term will allow to obtain a correct behaviour of the self-energy both in the weak and strong coupling limits. The denominator of a similar type (as in Eq. (22)) was used in the MPT for the quantum dot in the SC-QD-N system by Martin-Rodero et al. [6–8], and Yamada et al. [9,10]. In MPT the denominator appearing in self-energy expression is introduced intuitively, to fit the MPT results to the correct results of the weak and strong coupling limit. In our approach, it comes out directly from the EOM method improved by us. ˆ ′ we use the following algorithm: To calculate the self-energy Σ U 1. Using the mean-field approximation we compute the mean-field selfˆ HF . We calculate the expected values of the parameters 〈nd↓〉 energy Σ U and 〈d↓d↑〉. ˆ eff = Σ ˆ HF . 2. We introduce the effective self-energy Σ U

πΓNρd/2 1

U

0.8

(2)ir

3. To calculate the τˆ ir3 Γˆ d τˆ 3 function given by Eq. (B.17) we use Green’s ˆ eff . functions g<,>(ε ) with effective self-energy Σ

0.6

U

4. Using the results obtained in the step No. 3 we calculate from Eq. ˆ ′ , where the denominator (22) the higher order self-energy Σ U −1

0.4 0.2

−1

ˆ (2)B ˆ ˆ (2)gˆ HF ] [Iˆ + Σ is replaced by [Iˆ + Σ with parameter U ] U d HF ˆ = Re gˆ . B d 5. The higher order self-energies are then used to calculate selfconsistently the first-order contributions to the self-energies (9). ˆ =Σ ˆ HF + Σ ˆ ′ we calculate the 6. For the obtained self-energy Σ U U U expectation values of parameters 〈nd↓〉 and 〈d↓d↑〉. ˆ eff = Re[Σ ˆ (μ)] 7. We introduce the new value of effective self-energy Σ U

0

1

2 U

3

4

5

-4

-2

0

2

4

ε

Fig. 2. The normalized spectral density ρd (ε ) of correlated quantum dot as a function of energy ε for different values of the Coulomb interaction U , obtained for ΓS = 1, ΓN = 1 and

U

Δ = 15.

23

Physica E 89 (2017) 21–28

G. Górski, K. Kucab

πΓNρoff /2 0.4 0.2 0 -0.2 -0.4 0

1

2 U

3

4

5

-4

0

-2

2 ε

4

Fig. 3. The anomalous spectral density ρoff (ε ) of correlated quantum dot as a function of energy ε for different values of the Coulomb interaction U , obtained for ΓS = 1, ΓN = 1 and

Δ = 15. Fig. 4. The values of renormalized parameters as a function of the Coulomb interaction U , obtained for ΓS = 1, ΓN = 1 and Δ = 15.

they become single Kondo resonance state. At large values of U the narrow Kondo state localized at ε = 0 is dominant. There are also two broad states localized at ±U /2 [32]. In Fig. 3 we show the anomalous spectral density ρoff (ε ) = − ImG21(ε )/ π for different values of U . With increasing the Coulomb repulsion the value of ρoff (ε ) is decreasing. Further increase of U leads to a sign change of ρoff (ε ) in the energy gap region. This sign change can be used as the singlet-doublet transition criterion [11]. As Yamada et al. showed [9], for the SC-QD-N system with ΓN ≠ 0 and finite SC gap, the sign change of an anomalous spectral density can be interpreted as the criterion of the SC singlet-Kondo singlet transition. For the parameters used in Fig. 3 the SC singlet-Kondo singlet transition occurs for U ≈ 4 . For U < 4 and Δ > ε > 0 the anomalous spectral density ρoff (ε ) is negative and the system is in the SC singlet state. For U > 4 and Δ > ε > 0 the anomalous spectral density ρoff (ε ) changes its sign to positive what is characteristic for a Kondo singlet state. For large values of ΓN we observe the broad crossover region of SC singlet-Kondo singlet transition. In this region the anomalous spectral density disappears but the spectral density ρd (ε ) rapidly grows. The similar behaviour of ρd (ε ) and ρoff (ε ) in the transition region was

coupling parameter between the dot and the superconducting lead). For large U we obtain ΓSR → 0 . Decreasing to zero the value of ΓSR causes that the spectral density for the Fermi level tends to the constant value ρd (μN ) → 2/ πΓN . In experimental studies of the subgap regime, the differential conductance of tunnelling current is measured. It was shown by Krawiec and Wysokiński [12] that in this regime the current due to the Andreev tunnelling is important. In this case the transmittance is given by

TA(ε ) = ΓN2 Gd ,21(ε ) 2 .

In Fig. 5 we present the dependence of Andreev transmittance TA(ε ) for different values of the Coulomb interaction. At small values of U < 4 , which corresponds to the superconducting singlet state, the broad transmittance peak with the maximum value close to 1 is observed. With increasing of the Coulomb repulsion the Andreev transmittance is decreasing. At large values of U one can observe the weak peaks localized at ε ≈ ± U /2 . This behaviour can be compared with decreasing of anomalous spectral density ρoff (ε ) (see Fig. 3) for large U . In the real quantum dot systems we have U / Δ > > 1 and also the large value of ΓS / ΓN ratio [1,2]. In Fig. 6 we present the normalized spectral density, ρd (ε ), and the normalized anomalous spectral density, ρoff (ε ), for different values of ΓN and for U / Δ = 2.5. The results show that for ΓN = ΓS the coupling with the lead in the normal state is so strong that the Andreev resonances disappear and there is one broad resonance peak with the maximum at ε = 0 . When ΓN < < ΓS there are two Andreev resonances localized near ±Δ and there is also the strong central resonance peak. In Fig. 6(b) we can see the sign change

observed by Žitko et al. [11]. The increase of ρd (ε ) in a crossover region can be visualized by introducing the renormalized coupling parameters [13]

ΓNR = zΓN

(23)

and

ΓSR = z(ΓS − 2Σ12(0)),

(24)

where

⎛ Γ dΣ11(ε ) z = ⎜⎜1 + S − 2Δ dε ⎝

⎞−1 ⎟⎟ . ε →0 ⎠

TA 1

(25)

Using these parameters, the DOS near Fermi level can be written as

ρd (ε ) ≈

⎡ ⎤ ΓNR /2 ΓNR /2 z ⎢ ⎥ + 2 2 2π ⎢⎣ (ε + Γ R /2)2 + (Γ R /2)2 (ε − ΓSR /2) + (ΓNR /2) ⎥⎦ S N

0.75 (26)

0.5

and the anomalous spectral density ρoff (ε ) of correlated quantum dot is now given by

0.25

ρoff (ε ) ≈

⎡ ⎤ ΓNR /2 ΓNR /2 z ⎢ ⎥. − 2 2 2 2 2π ⎢⎣ (ε + Γ R /2) + (Γ R /2) (ε − ΓSR /2) + (ΓNR /2) ⎥⎦ S N ΓNR

(28)

0 (27)

1

2 U

ΓSR

In Fig. 4 we show the dependence of and on the Coulomb repulsion. With increasing the Coulomb repulsion we observe the reduction of renormalized coupling parameters. This reduction is especially strong for the ΓSR parameter (which is the renormalized

3

4

5

-4

-2

0

2

4

ε

Fig. 5. The Andreev transmittance obtained at T = 0 as a function of energy ε for different values of the Coulomb interaction U . The other parameters are the same as in Fig. 2.

24

Physica E 89 (2017) 21–28

G. Górski, K. Kucab

1 0.8

πρd/2

small values of U (U < ΓS ) we obtain two broad peaks of dI / dV . Increasing the value of U we observe that these two peaks evolve into one, strong central peak. For U ≈ ΓS we obtain the crossover region, where the differential conductance reaches its maximum value equal to 4e2 / h at zero bias voltage. Further growth of the Coulomb interaction (U > ΓS ) causes the decrease of the magnitude of the zero-bias conductance peak [33]. This behaviour is related to the interplay between the Andreev reflection and the Kondo effect. For the values of U even greater, besides the central peak of dI / dV , one can observe two additional peaks [32]. The maximum value of zero-bias conductance peak can be used to determine the conditions for superconducting singlet-Kondo singlet transition. Comparing the results for differential conductance dI / dV (eV ) (Fig. 7) and for the normalized spectral density ρd (ε ) in the equilibrium state (Fig. 2) one can see some differences between obtained profiles of these dependencies. For U < ΓS we observe similar profiles for dI / dV (eV ) and ρd (ε ) with two maxima, which are “moving” toward central point with increasing the value of U . For U > ΓS the magnitude of spectral density of zero-energy Abrikosov-Suhl peak, ρd (ε = 0), grows with increasing value of U , while for dI / dV we observe the magnitude decrease of the zero-bias conductance peak.

ΓN=2 1 0.5 0.2

(a)

0.6 0.4 0.2 0 (b)

0.2 πρoff /2

0.1 0 -0.1 -0.2 -2

-1.5

-1

-0.5

0 ε

0.5

1

1.5

2

4. Conclusions In this paper, we analysed the Coulomb interaction dependence of the spectral density and the anomalous spectral density of quantum dot for SC-QD-N system. In this analysis, we used the irreducible Green’s functions technique. This technique is based on the equation of motion method but in contrast to the classical EOM (see [29–31]) the doubletime Green’s functions are differentiated with respect to both time variables. Using the IGF technique one can obtain the higher order part −1 (2) ˆ] Σ ˆ (2)A ˆ . As the results of modified of self-energy in the form [Iˆ + Σ

Fig. 6. The normalized spectral density ρd (ε ) (a), and the normalized anomalous spectral density ρoff (ε ) (b) for different values of ΓN and the parameters U = 2.5, Δ = 1 and ΓS = 1.

of the anomalous spectral density during superconducting singletKondo singlet transition. The similar dependence of anomalous spectral density on the coupling parameter between the dot and the metallic lead was reported by Žitko et al. [11] where the NRG method was used. The EOM technique used here allow us to calculate the current flowing through the SC-QD-N system. In the subgap regime eV ≤ Δ the charge transport occurs only across the Andreev scattering. The current flowing through the normal lead can be calculated as

I=−

2eΓN Im h

U

U

occurs in the self-energy Eq. (22). It reproduces the correct results in the atomic limit. In MPT the denominator appearing in self-energy expression is introduced intuitively, to fit the MPT results to the correct results of atomic limit. In the IGF technique this term arises from equations of motion. The relations for Green’s functions and the self-energy obtained by us were used to find the Coulomb interaction dependence of the spectral density and the anomalous spectral density. We showed that increasing of the Coulomb interaction causes the Andreev resonances tends to the single Kondo resonance state. For the anomalous spectral density the increase of the Coulomb interaction causes its sign change. The effect of sign change of the anomalous spectral density can be used as the SC singlet-Kondo singlet transition criterion. We also showed the influence of the Coulomb interaction on the non-equilibrium transport properties. The increase of U leads to the enhancement of the Andreev transport. Such behaviour gives the strong zero-bias peak in the subgap differential conductance. Additionally, for large values of U there are additional peaks in the dI / dV structure, localized at the Andreev bound states [32].

∞

∫−∞ [2f (ε − eV )G11(ε) + G11<(ε)]dε,

(29)

< (ε ) is the lesser Green’s function [9,33]. where G11 Using Eq. (29) we compute the differential conductance dI / dV . In the computations, we use the infinity gap limit [33] and asymmetric couplings for tunnelling, ΓS > ΓN . In Fig. 7 we show the differential conductance for several values of the Coulomb interactions U . For

U=0.5 1 1.5 2 2.5

dI/dV [2e2/h]

1 0.8

U

perturbation theory (MPT) show [8,9,27,34–36], this form of selfenergy improves the results of the second order perturbation theory −1 ˆ ] expression which ˆ (2)A (SOPT). Of particular importance is the [Iˆ + Σ

0.6 0.4

Acknowledgments

0.2 0 -2.5

-2

-1.5

-1

-0.5

0

0.5

1

1.5

2

This work was done due to support from Faculty of Mathematics and Natural Sciences University of Rzeszów within the project no. WMP/GD-12/2016 and due to partial support from Centre for Innovation and Transfer of Natural Sciences and Engineering Knowledge of Rzeszów University.

2.5

eV Fig. 7. Differential conductance as a function of bias voltage for several values of the Coulomb interactions U and for the parameters Δ → ∞, ΓS = 2 and ΓN = 0.5.

25

Physica E 89 (2017) 21–28

G. Górski, K. Kucab

Appendix A Applying Eq. (3) to the Hamiltonian (1) we obtain

ˆ − gˆ d 0−1G d

(1)

∑ τˆ 3VˆkαGˆ kdα = Iˆ + U τˆ 3Γˆ d

,

kσα

(A.1)

⎡ ε − εd 0 ⎤ gˆ d 0−1 = ⎢ ⎥, ε + εd ⎦ ⎣ 0

(A.2)

⎡V 0 ⎤ ˆ = ⎢ kα ⎥ V kα † , ⎢⎣ 0 Vkα ⎥⎦

(A.3)

⎡〈〈c ; d +〉〉 〈〈c ; d 〉〉 ⎤ kα ↑ ↓ ε ˆ (ε ) = ⎢ kα ↑ ↑ ε ⎥, G kdα + + + ⎢⎣〈〈ckα ⎦ ↓; d↑ 〉〉ε 〈〈ckα ↓; d↓〉〉ε ⎥

(A.4)

where

Iˆ is the unit matrix and τˆ i (i = 1, 2, 3) is the Pauli matrix. ˆ (ε ) Using EOM we obtain the equation for the Green’s function G kdα

∑ τˆ 3VˆkαGˆ kdα = ∑ τˆ 3Vˆkαgˆ kαVˆkατˆ 3Gˆ d , kσα

(A.5)

kσα

where gˆ kN is the Green’s function of the normal lead 1 ⎡ ⎢ ε − ε kN + μN =⎢ 0 ⎢⎣

gˆ kN

⎤ ⎥ ⎥, 1 ε + ε kN − μN ⎥ ⎦ 0

(A.6)

and gˆ kS denotes the Green’s function of the superconducting lead

⎤ −Δ 1 ⎡ ε + ε−kS − μS ⎥, ⎢ ε − εkS + μS ⎦ −Δ detA ⎣

gˆ kS =

(A.7) 2

where det A = (ε + εkS − μS )(ε − ε−kS + μS ) − Δ . ˆ gˆ V ˆ ˆ can be written as the non-interacting self-energy Using the approximation for the wide band limit, the term ∑kσα τˆ 3V kα kα kατ 3

ˆ0 = Σ d

∑ τˆ 3Vˆkαgˆ kαVˆkατˆ 3 = − kσα

where

γ (ε ) =

ΓN (S ) ε

is

the

⎡1 Δ⎤ i ˆ i ΓN I − ΓSγ (ε )⎢⎢ Δ ε ⎥⎥ , 2 2 ⎣ ε 1⎦

parameter

Θ( ε − Δ) + 2

ε2 − Δ

ε i Δ2 − ε

which

characterizes

(A.8) the

coupling

between

the

dot

and

the

normal

(superconducting)

lead;

Θ (Δ − ε ) . 2

Appendix B The function ir〈〈nˆd ↓d↑; nˆd ↓d↑+〉〉irε can be written as ir

〈〈nˆd ↓d↑; nˆd ↓d↑+〉〉irε =

∞

∫−∞ d (t − t′)[−Θ(t − t′)]

⎡ ε ⎤ × 〈[ird↓+(t )d↓(t )d↑(t ), d↓+(t′)d↓(t′)d↑+(t′)ir ] 〉exp⎢i (t − t′)⎥ , + ⎣ ℏ ⎦

(B.1)

where Θ(t − t′) is the step function

Θ(t − t′) =

i 2π

∞

−ix (t − t )

′

∫−∞ dx ex + i0+ .

(B.2)

The mean value of the first term of anticommutator in Eq. (B.1) is approximated by ir

〈d↓+(t )d↓(t )d↑(t )d↓+(t′)d↓(t′)d↑+(t′)〉ir ≈〈d↓+(t )d↓(t′)〉〈d↓(t )d↓+(t′)〉〈d↑(t )d↑+(t′)〉 − 〈d↓+(t )d↑+(t′)〉〈d↓(t )d↓(t′)〉〈d↑(t )d↓+(t′)〉.

(B.3)

The irreducible Green’s function cannot be reduced to the low-order Green’s function, therefore the averages of operators with the same time variables ( A(t )B(t ) ) in Eq. (B.3) were neglected. Introducing the definitions

⎡ 〈d +(t′)d (t )〉 〈d (t′)d (t )〉 ⎤ ↑ ↑ ↓ ↑ ⎥ −i gˆ <(t , t′) = ⎢ + ⎢⎣〈d↑ (t′)d↓+(t )〉 〈d↓(t′)d↓+(t )〉⎥⎦

(B.4)

and

26

Physica E 89 (2017) 21–28

G. Górski, K. Kucab

⎡ 〈d (t )d +(t′)〉 〈d (t )d (t′)〉 ⎤ ↑ ↑ ↑ ↓ ⎥ i gˆ >(t , t′) = ⎢ + ⎢⎣〈d↓ (t )d↑+(t′)〉 〈d↓+(t )d↓(t′)〉⎥⎦

(B.5)

we can write ir

〈d↓+(t )d↓(t )d↑(t )d↓+(t′)d↓(t′)d↑+(t′)〉ir ≈ < > > ig22 (t′, t )[g22 (t , t′)g11>(t , t′) − g21 (t , t′)g12>(t , t′)].

(B.6)

Calculating the second average in Eq. (B.1) in a similar way we have ir

〈d↓+(t′)d↓(t′)d↑+(t′)d↓+(t )d↓(t )d↑(t )〉ir ≈ > < < −ig22 (t′, t )[g22 (t , t′)g11<(t , t′) − g21 (t , t′)g12<(t , t′)].

(B.7)

Using the Fourier transform

gij<>(t , t′) =

1 2π

∞

∫−∞ gij<>(x )e−ix(t−t′)dx

(B.8)

and inserting Eqs. (B.6) and (B.7) into Eq. (B.1) we obtain ir

〈〈nˆd ↓d↑; nˆd ↓d↑+〉〉irε =

i (2π )3

∭ ε + y −dxdydz x − z + i 0+

< > > > < < ×{g22 (y )[g22 (x )g11>(z ) − g21 (x )g12>(z )] − g22 (y )[g22 (x )g11<(z ) − g21 (x )g12<(z )]}.

(B.9)

Eq. (B.9) can be written in the following form ir

⎧

∞

i 2π

〈〈nˆd ↓d↑; nˆd ↓d↑+〉〉irε =

⎫

∞

∫−∞ ε − xdx+ i0+ ⎨⎩ 21π ∫−∞ [P>(x + y)g22<(y) − P<(x + y)g22>(y)]dy⎬⎭,

(B.10)

where

P >(x ) =

∞

1 2π

∫−∞ [g22>(y)g11>(x − y) − g21>(y)g12>(x − y)]dy,

1 2π

∫−∞ [g22<(y)g11<(x − y) − g21<(y)g12<(x − y)]dy.

(B.11)

and

P <(x ) =

∞

For the off-diagonal term ir

〈〈nˆd ↑d↓+; nˆd ↓d↑+〉〉irε =

ir

〈〈nˆd ↑d↓+;

nˆd ↓d↑+〉〉irε

(B.12)

we can write

∞

∫−∞ d (t − t′)[−Θ(t − t′)]

⎡ ε ⎤ × 〈[ird↑+(t )d↑(t )d↓+(t ), d↓+(t′)d↓(t′)d↑+(t′)ir ] 〉exp⎢i (t − t′)⎥ . + ⎣ ℏ ⎦

(B.13)

Proceeding in the similar way as before we obtain ir

〈d↑+(t )d↑(t )d↓+(t )d↓+(t′)d↓(t′)d↑+(t′)〉ir < ig21 (t′,

> t )[g22 (t ,

t′)g11>(t ,

≈

> t′) − g21 (t , t′)g12>(t , t′)]

(B.14)

and ir

〈d↓+(t′)d↓(t′)d↑+(t′)d↑+(t )d↑(t )d↓+(t )〉ir ≈ > < < −ig21 (t′, t )[g22 (t , t′)g11<(t , t′) − g21 (t , t′)g12<(t , t′)]

(B.15)

hence ir

〈〈nˆd ↑d↓+; nˆd ↓d↑+〉〉irε =

i 2π

∞

∞

∫−∞ ε − xdx+ i0+ ∫−∞ [P>(x + y)g21<(y) − P<(x + y)g21>(y)]dy.

(B.16)

Using the results obtained above one can write the following relation

ˆ (2) = U 2 i Σ U 2π

∞

∞

∫−∞ ε − xdx+ i0+ ∫−∞ [P>(x + y)ˆτ2[g<(y)]T τˆ 2 − P<(x + y)ˆτ2[g>(y)]T τˆ 2]dy.

(B.17)

[3] E.J.H. Lee, X. Jiang, M. Houzet, R. Aguado, C.M. Lieber, S. De Franceschi, Nat. Nanotechnol. 9 (2014) 79. http://dx.doi.org/10.1038/nnano.2013.267. [4] T. Dirks, T.L. Hughes, S. Lal, B. Uchoa, Y.-F. Chen, C. Chialvo, P.M. Goldbart, N. Mason, Nat. Phys. 7 (2011) 386. http://dx.doi.org/10.1038/nphys1911. [5] V. Mourik, K. Zuo, S.M. Frolov, S.R. Plissard, E.P.A.M. Bakkers, L.P. Kouwenhoven, Science 336 (2012) 1003. http://dx.doi.org/10.1126/ science.1222360. [6] A. Martín-Rodero, A. Levy Yeyati, Adv. Phys. 60 (2011) 899. http://dx.doi.org/ 10.1080/00018732.2011.624266.

References [1] R.S. Deacon, Y. Tanaka, A. Oiwa, R. Sakano, K. Yoshida, K. Shibata, K. Hirakawa, S. Tarucha, Phys. Rev. B 81 (2010) 121308(R). http://dx.doi.org/10.1103/ PhysRevB.81.121308. [2] R.S. Deacon, Y. Tanaka, A. Oiwa, R. Sakano, K. Yoshida, K. Shibata, K. Hirakawa, S. Tarucha, Phys. Rev. Lett. 104 (2010) 076805. http://dx.doi.org/10.1103/ PhysRevLett.104.076805.

27

Physica E 89 (2017) 21–28

G. Górski, K. Kucab

[22] A.L. Kuzemsky, Riv. Nuovo Cim. 25 (2002) 1. [23] A.L. Kuzemsky, Physica A 267 (1999) 131–152. http://dx.doi.org/10.1016/S03784371(98)00665-7. [24] G. Górski, J. Mizia, K. Kucab, Phys. Status Solidi B 253 (6) (2016) 1202–1209. http://dx.doi.org/10.1002/pssb.201552722. [25] S.N. Mitra, K.G. Chakrabort, J. Phys.: Condens. Matter 7 (1995) 379–388. http:// dx.doi.org/10.1088/0953-8984/7/2/016. [26] K.I. Wysokiński, A.L. Kuzemsky, J. Low Temp. Phys. 52 (1983) 81. http:// dx.doi.org/10.1007/BF00681267. [27] G. Górski, J. Mizia, K. Kucab, Physica E 73 (2015) 76–82. http://dx.doi.org/ 10.1016/j.physe.2015.05.021. [28] P.W. Brouwer, A. Lamacraft, K. Flensberg, Phys. Rev. B 72 (2005) 075316. http:// dx.doi.org/10.1103/PhysRevB.72.075316. [29] C. Lacroix, J. Phys. F: Met. Phys. 11 (1981) 2389. http://dx.doi.org/10.1088/03054608/11/11/020. [30] P. Trocha, J. Barnaś, Phys. Rev. B 76 (2007) 165432. http://dx.doi.org/10.1103/ PhysRevB.76.165432. [31] J. Barański, T. Domański, Phys. Rev. B 84 (2011) 195424. http://dx.doi.org/ 10.1103/PhysRevB.84.195424. [32] J. Barański, T. Domański, J. Phys.: Condens. Matter 25 (2013) 435305. http:// dx.doi.org/10.1088/0953-8984/25/43/435305. [33] T. Domański, I. Weymann, M. Barańska, G. Górski, Sci. Rep. 6 (2016) 23336. http://dx.doi.org/10.1038/srep23336. [34] A.A. Aligia, Phys. Rev. B 74 (2006) 155125. http://dx.doi.org/10.1103/ PhysRevB.74.155125. [35] H. Kajueter, G. Kotliar, Phys. Rev. Lett. 77 (1996) 131. http://dx.doi.org/10.1103/ PhysRevLett.77.131. [36] M. Potthoff, T. Wegner, W. Nolting, Phys. Rev. B 55 (1997) 16132. http:// dx.doi.org/10.1103/PhysRevB.55.16132.

[7] J.C. Cuevas, A. Levy Yeyati, A. Martín-Rodero, Phys. Rev. B 63 (2001) 094515. http://dx.doi.org/10.1103/PhysRevB.63.094515. [8] E. Vecino, A. Martín-Rodero, A. Levy Yeyati, Phys. Rev. B 68 (2003) 035105. http://dx.doi.org/10.1103/PhysRevB.68.035105. [9] Y. Yamada, Y. Tanaka, N. Kawakami, Phys. Rev. B 84 (2011) 075484. http:// dx.doi.org/10.1103/PhysRevB.84.075484. [10] Y. Yamada, Y. Tanaka, N. Kawakami, Physica C 470 (2010) S875. http:// dx.doi.org/10.1016/j.physc.2009.10.044. [11] R. Žitko, J.S. Lim, R. López, R. Aguado, Phys. Rev. B 91 (2015) 045441. http:// dx.doi.org/10.1103/PhysRevB.91.045441. [12] M. Krawiec, K.I. Wysokiński, Supercond. Sci. Technol. 17 (2004) 103. http:// dx.doi.org/10.1088/0953-2048/17/1/018. [13] Y. Tanaka, A. Kawakami, A. Oguri, J. Phys. Soc. Jpn. 76 (2007) 074701. http:// dx.doi.org/10.1143/JPSJ.76.074701. [14] A. Oguri, Y. Tanaka, J. Bauer, Phys. Rev. B 87 (2013) 075432. http://dx.doi.org/ 10.1103/PhysRevB.87.075432. [15] S.-Y. Hwang, R. López, D. Sánchez, Phys. Rev. B 91 (2015) 104518. http:// dx.doi.org/10.1103/PhysRevB.91.104518. [16] Y. Tanaka, A. Oguri, A.C. Hewson, New J. Phys. 9 (2007) 115. http://dx.doi.org/ 10.1088/1367-2630/9/5/115. [17] A. Martín-Rodero, A. Levy Yeyati, J. Phys.: Condens. Matter 24 (2012) 385303. http://dx.doi.org/10.1088/0953-8984/24/38/385303. [18] K.P. Wójcik, I. Weymann, Phys. Rev. B 89 (2014) 165303. http://dx.doi.org/ 10.1103/PhysRevB.89.165303. [19] T. Domański, A. Donabidowicz, Phys. Rev. B 78 (2008) 073105. http://dx.doi.org/ 10.1103/PhysRevB.78.073105. [20] T. Domański, A. Donabidowicz, K.I. Wysokiński, Phys. Rev. B 78 (2008) 144515. http://dx.doi.org/10.1103/PhysRevB.78.144515. [21] A.A. Clerk, V. Ambegaokar, S. Hershfield, Phys. Rev. B 61 (2000) 3555. http:// dx.doi.org/10.1103/PhysRevB.61.3555.

28