Time-dependent Landauer-Büttiker formalism for superconducting junctions at arbitrary temperatures
Riku Tuovinen, Robert van Leeuwen, Enrico Perfetto, Gianluca Stefanucci
TTime-dependent Landauer–Büttiker formalism forsuperconducting junctions at arbitrary temperatures
Riku Tuovinen , Robert van Leeuwen , , Enrico Perfetto andGianluca Stefanucci , , Department of Physics, Nanoscience Center, FIN 40014, University of Jyväskylä, Jyväskylä,Finland Dipartimento di Fisica, Università di Roma Tor Vergata, Via della Ricerca Scientifica 1,I-00133 Rome, Italy Laboratori Nazionali di Frascati, Istituto Nazionale di Fisica Nucleare, Via E. Fermi 40,00044 Frascati, Italy European Theoretical Spectroscopy Facility (ETSF)E-mail: [email protected]
Abstract.
We discuss an extension of our earlier work on the time-dependent Landauer–Büttiker formalism for noninteracting electronic transport. The formalism can withoutcomplication be extended to superconducting central regions since the Green’s functions in theNambu representation satisfy the same equations of motion which, in turn, leads to the sameclosed expression for the equal-time lesser Green’s function, i.e., for the time-dependent reducedone-particle density matrix. We further write the finite-temperature frequency integrals in termsof known special functions thereby considerably speeding up the computation. Simulations insimple normal metal – superconductor – normal metal junctions are also presented.
1. Introduction
The process of Andreev reflection[1] (AR) occurring at the interface between a normal metal (N)and a superconductor (S) is of great importance with applications in spintronics and quantumcomputing. An incoming electron from N to S produces a Cooper pair in S and a reflected hole inN[2–4]. In an NSN junction normal metal electrodes are spatially separated by a superconductingcentral region and an entangled electron–hole pair can be transported. This can be seen whenthe junction separation is of the order of the superconducting coherence length for the studiedmaterial, and when the incident electron energies are less than the superconducting gap for theAR process to occur[5].The quantum transport problems are typically time dependent; there is no guarantee thatthe system would in an instant relax to a steady-state configuration once the junction is“switched on” (as in connecting different devices or driving them out of equilibrium by anexternal perturbation). In contrast, there are transient effects depending on, e.g., the system’sgeometry[6–9], its predisposition to external perturbations[10–13], the physical properties of thetransported quanta and their mutual interactions[14–19]. Even if the transport mechanismswere discussed in an idealized noninteracting setting, it is therefore important to consider a fullytime-dependent description of the studied processes. a r X i v : . [ c ond - m a t . m e s - h a ll ] A p r he Landauer–Büttiker formalism is simple to understand as it relates to an intuitive physicalpicture of charge transport in a multiterminal junction[20, 21]. Including the transient descriptionto the formalism by studying the nonequilibrium Green’s function approach does not complicatethe final result[22–27]; the physical picture is still clear and intuitive as different features of thetransport setup can be directly linked to the time-dependence[28–30]. In this paper, we presentan extension to earlier results for both superconducting junctions and arbitrary temperatures(Sec. 2). We present a formula for the time-dependent one-particle reduced density matrix(TD1RDM) as such since it is a closed, analytic expression which can readily be implemented fornumerical model systems. Further details of the derivation are to be found in another work[31].In Sec. 4 we illustrate the features of the formula by studying transients in simple NSN junctions.
2. Background and Nambu representation
We consider a quantum transport setup similar to one studied in the previous volume of thisconference series[29]. In this setup, a noninteracting central region is connected between metallicleads, and the Hamiltonian takes the form ˆ H = ˆ H leads + ˆ H central + ˆ H coupling = (cid:88) kασ [ (cid:15) kα + θ ( t ) V α ] ˆ c † kασ ˆ c kασ + (cid:88) mnσ T mn ˆ c † mσ ˆ c nσ + (cid:88) mkασ (cid:16) T mkα ˆ c † mσ ˆ c kασ + h.c. (cid:17) . (1)The operators ˆ c ( † ) annihilate (create) electrons from (to) a region specified by the subscriptindices: kα is the k -th basis element of the α -th lead, m, n label the basis elements of thecentral region, and σ ∈ {↑ , ↓} is a spin- index. These operators moreover obey the fermionicanticommutation relations { ˆ c xσ , ˆ c † yσ (cid:48) } = δ xy δ σσ (cid:48) . The Hamiltonian structure is determined bythe single-particle levels in the leads (cid:15) kα and the tunneling matrices T between the states of thecentral region ( T mn ) and between the states of the central region and the leads ( T mkα ). Thesystem is driven out of equilibrium for times t > by a sudden shift of the lead energy levelsby V α . In addition to a sudden bias, it is also possible to include time-dependent bias profileswithout complicating the following derivations[32]. In order to describe a superconducting island,we will now add a pairing field operator ˆ ∆ to the Hamiltonian of the central region by[33, 34] ˆ H central → (cid:88) mnσ T mn ˆ c † mσ ˆ c nσ + (cid:88) m ∆ m ˆ c † m ↑ ˆ c † m ↓ + (cid:88) m ∆ ∗ m ˆ c m ↓ ˆ c m ↑ . (2)The one-electron Green’s function in the above setup can be defined via the Nambu spinor ˆ Φ m = ( ˆ Φ m , ˆ Φ m ) T = (ˆ c m ↑ , ˆ c † m ↓ ) T , which obeys the anticommutation relation in a tensor productsense { ˆ Φ µm , ˆ Φ νn } = δ mn δ µν , as a contour-ordered product[34, 35] G rs ( z, z (cid:48) ) = − i (cid:104)T γ [ ˆ Φ r ( z ) ⊗ ˆ Φ † s ( z (cid:48) )] (cid:105) (3)where the contour-ordering operator T γ is taken for the variables z, z (cid:48) on the Keldysh contour γ [28]. When the product in Eq. (3) is expanded, the elements in the resulting × Nambumatrix are the normal and the anomalous components of the Green’s function[36, 37]. As alreadydenoted above, we put an underline for quantities in the Nambu space. The matrix elements inthe Green’s function label the transport setup in the following block form [notice that the indices r, s in Eq. (3) may belong to any block and that the number of leads is arbitrary] h = h · · · h C h · · · h C ... ... . . . ... h C h C · · · h CC ; G = G G · · · G C G G · · · G C ... ... . . . ... G C G C · · · G CC (4)ith ( h αα (cid:48) ) kk (cid:48) ( t ) = (cid:18) [ (cid:15) kα + θ ( t ) V α ] δ αα (cid:48) δ kk (cid:48) − [ (cid:15) kα + θ ( t ) V α ] δ αα (cid:48) δ kk (cid:48) (cid:19) , (5) ( h CC ) mn = (cid:18) T mn ∆ m δ mn ∆ ∗ m δ mn − T mn (cid:19) , (6) ( h Cα ) mkα = (cid:18) T mkα − T mkα (cid:19) , ( h αC ) kαm = (cid:18) T kαm − T kαm (cid:19) (7)for the leads, central region and couplings, respectively. It is important to notice that eventhough only the central region is superconducting, all the blocks in the Hamiltonian are writtenin the form of Bogoliubov–de Gennes[38, 39]. Including the pairing field in the Hamiltonianof the central region adds no extra complication to the evolution of the Green’s function[31];the only difference, compared to the earlier work in Refs. [28–30], is in the interpretation of thematrices in Nambu space. It is also possible to include non-local pairing field ∆ mn with arbitraryspin-coupling leading to -component Nambu spinors which includes the possibility to study alsoMajorana fermions[31, 40]. The Hamiltonian and the Green’s function are connected via theequation of motion (with the boundary condition that the Green’s function is antiperiodic alongthe contour, i.e., the Kubo–Martin–Schwinger boundary conditions[28, 41, 42]) (cid:20) i dd z − h ( z ) (cid:21) G ( z, z (cid:48) ) = δ ( z, z (cid:48) ) (8)and the corresponding adjoint one. We describe the leads within the wide-band approximation(WBA), where the electronic levels of the central region are in a narrow range compared to thelead bandwidth which gives for the retarded embedding self-energy Σ R α,mn ( ω ) = (cid:88) k ( h Cα ) mkα ω − (cid:15) kα − V α + i η ( h αC ) kαn ≈ − i Γ α,mn / , (9)where the bandwidth matrices satisfy Γ = (cid:80) α Γ α . For the lead Green’s function between thecoupling Hamiltonians in Eq. (9) the structure is similar to that of Eq. (5). Approximating theembedding self-energy this way, as a purely imaginary constant, closes the equation of motion (8),and it can then be solved analytically[28–30]. As the solution we get the lesser Green’s function(in region CC ) in the equal-time limit and the TD1RDM by ( ρ CC ) mn ( t ) = − i( G
3. An introductory example
Let us motivate the discussion for the NSN setup by means of a simple example. Consider asingle dot connected to two leads for which the Hamiltonian can be separated in parts for leads,tunneling and dot, respectively as ˆ H = (cid:88) kασ (cid:15) kα ˆ c † kασ ˆ c kασ + (cid:88) kασ t kα ˆ c † kασ ˆ c σ + (cid:88) kασ t ∗ kα ˆ c † σ ˆ c kασ + (cid:15) (cid:88) σ ˆ c † σ ˆ c σ + ∆ ˆ c † ↑ ˆ c † ↓ + ∆ ∗ ˆ c ↓ ˆ c ↑ (16)with (cid:15) kα giving the level structure of the leads α ∈ { L, R } , t kα corresponding to the tunnelingstrength between the leads and the dot, and (cid:15) , ∆ being the energy and the pairing strength inthe dot, respectively. Let us introduce a new set of operators ˆ˜ c xσ = ˆ c † xσ obeying the fermionicanticommutation relation. The Hamiltonian can now be rewritten in terms of the new and oldoperators as ˆ H = (cid:88) kα (cid:15) kα ˆ c † kα ↑ ˆ c kα ↑ + (cid:88) kα ( − (cid:15) kα )ˆ˜ c † kα ↓ ˆ˜ c kα ↓ + (cid:88) kα (cid:15) kα + (cid:88) kα (cid:16) t kα ˆ c † kα ↑ ˆ c ↑ + t ∗ kα ˆ c † ↑ ˆ c kα ↑ (cid:17) + (cid:88) kα (cid:104) ( − t kα )ˆ˜ c † ↓ ˆ˜ c kα ↓ + ( − t ∗ kα )ˆ˜ c † kα ↓ ˆ˜ c ↓ (cid:105) + (cid:15) ˆ c † ↑ ˆ c ↑ + ( − (cid:15) )ˆ˜ c † ↓ ˆ˜ c ↓ + (cid:15) + ∆ ˆ c † ↑ ˆ˜ c ↓ + ∆ ∗ ˆ˜ c † ↓ ˆ c ↑ (17)where two constant shifts (cid:80) kα (cid:15) kα and (cid:15) occur due to the anticommutation relations. Eachterm in Eq. (17) has a similar structure, ˆ c † ˆ c , and we may model the dot part as in Fig. 1 whereigure 1: The dot viewed as a two-level system for different spins. Figure 2: Couplings between thelead and the dot in a transportsetup.the matrix is of the form of Eq. (6) and the corresponding eigenvalues are (cid:15) ± = ± (cid:112) (cid:15) + | ∆ | .The transport setup corresponding to Eq. (17) can then be viewed through the energy diagramin Fig. 2 where we notice the nature of the constant shifts in Eq. (17); they could be regardedas a chemical potential. The energy levels for the lead sector are raised so that the energylevel continuum for the spin-up particles goes up from (cid:80) kα (cid:15) kα and the energy level continuumfor the spin-down particles goes down from (cid:80) kα (cid:15) kα . Similarly, for the dot sector we have theenergy levels raised by (cid:15) . The coupling terms t kα connect separately the spin-up and spin-downparticles between the leads and the dot, and the pairing strength term ∆ acts as a hopping termflipping the spins within the dot.According to this picture for the NSN setup and the Nambu structure in Eq. (10) we willevaluate the local bond currents for the more general structure in Eq. (6) by J mn ( t ) = − (cid:88) σ (cid:104) T mn ( G
4. TD response of a superconducting junction
As a first study we show a numerical confirmation of the presented formula; in the limit whenthe superconducting gap ∆ and the temperature /β vanish we should recover equal resultswith the formula in Ref. [30]. For the sake of simple interpretation of the transients let us take,as an example, a -site tight-binding dimer coupled to two semi-infinite one-dimensional leads.Let the hopping parameter between the sites be equal everywhere: t α = t αC = t C =: (cid:15) (for α = L, R ) leading to the tunneling rate Γ α = 2 t αC / | t α | = 2 (cid:15) . This parameter is the strengthof the bandwidth matrix elements in the form of Eq. (9). In fact, WBA is not a very goodapproximation in this case, as the resonances are comparatively rather wide, but we are only J ( t ) / Γ (a) NS, ∆ = 0 S, ∆ = (cid:15) / S, ∆ = (cid:15) / tΓ β → ∞ β = 10 /(cid:15) β = 2 /(cid:15) β = 1 /(cid:15) Figure 3: (Color online) Transient current through a -site dimer; comparison between differentformulae and parameters. (a) Normal vs. superconducting central region at zero temperature;(b) normal central region at varying temperatures.comparing different formulae within the WBA, so the results should only be taken as comparative.Let us also bias the leads symmetrically to V L = − V R = (cid:15) with respect to the chemical potential µ = 0 . We calculate, from the TD1RDM, the local bond current between the two sites in thedimer using Eq. (18). (Due to equal hopping parameters through the setup this is equal to thecurrent through the lead interfaces modulo a minor time delay.)In Fig. 3(a) we compare normal central region to a superconducting one by evaluating theTD1RDM from a normal Hamiltonian without the pairing field (non-Nambu), and from a NambuHamiltonian with varying pairing field strength at zero temperature. We see how the N and S( ∆ = 0 ) cases are on top of each other, and increasing the value for ∆ decreases the absolutevalue of the current through the central region as the energy levels of the central region get raisedby (cid:112) ∆ + (cid:15) . In Fig. 3(b) we compare normal central regions at varying temperatures. In thiscomparative benchmark of varying temperature, we do not consider superconducting centralregions as the temperature effects for the pairing field ∆ ( T ) should also be taken into accountaccording to the self-consistent gap equation[4]. (In further simulations, also these parametersare considered in more detail.) The zero-temperature limit, β → ∞ , is evaluated from the resultsin Ref. [30] which roughly agrees with an evaluation with β = 10 /(cid:15) . (Increasing β even morewould naturally bring the curves exactly on top of each other.) Because the level structure ofthe studied system is symmetric around the chemical potential, increasing the temperature /β decreases the current due to broadening of the distribution function close to the Fermi level. Ingeneral, however, there is a possibility of enhancing the current by increasing the temperatureif, for instance, the electronic levels were all above the Fermi level. In that case, the lead-stateswith energy higher than the energy of the levels get occupied, leading to an enhanced current.After the comparisons presented above we can confidently conclude that the formulation ofthe NSN transport setup and the implementation of the formula for the TD1RDM is workingproperly. Now, we turn to a more concrete and physically relevant example, and we analyzethe transient features in more detail. We consider a superconducting island made of a benzene-like molecule belonging to the class of quasi-one-dimensional polyacene chains[45]. Transportin simple island setups has been studied, e.g., in Refs. [46–50], in a single-electron-tunnelinglevel, where Coulomb blockade region is explored, and it is shown how the superconducting gap ∆ strongly and non-trivially affects the tunneling process. In polyacene samples (and in othercarbon based materials, such as graphene) the superconductivity could be induced, e.g., bycharge injection, chemical doping or using the proximity effect leading to critical temperaturesranging from to K[51–54]. Our setup is shown schematically in Fig. 4. We model thebenzene molecule in a single π -orbital tight-binding framework with the hopping parameter t C = − . eV[55], and relate other energies to this scale. We also saturate the molecule’s edges S N Figure 4: (Color online) Transport setup in an NSN junction. Normal metal leads of continuumstates undergo a level shift due to the bias voltage V L/R with respect to the chemical potential µ = E F (Fermi level). The discrete level structure of the central region is determined by thetight-binding and gap ∆ parameters, and it is also broadened due to the coupling ( Γ ). Possibletransition mechanisms are shown as CT, AR and CAR; see text for description.(longitudinally, in the transport direction) by hydrogen with modified tight-binding parametersfor hydrogen on-site energies and hydrogen–carbon hopping[56], respectively, so that there is noband gap in equilibrium. This condition is set because we want to isolate the effects from thesuperconducting gap ∆ without complicating the spectrum with the semiconducting gap. Thecoupling strength between the molecule and the leads and the lead hopping are chosen so thatwe are in weak coupling regime Γ = 0 . eV.Looking at the setup in Fig. 4 more closely suggests different transition mechanisms dependingon the parameters. The simplest case is when the bias voltage V α is larger than thesuperconducting gap ∆ . In this case, all the levels inside the bias window act as transportchannels, and transitions through the superconducting states are disrupted since the energy forthe incoming electrons is high enough to break possible Cooper pairs (CP); this is referred to asnormal tunneling (NT). If the bias voltage is smaller than the superconducting gap, this opens apossibility for the formation of a CP in the central region. In this case, it is possible to observeAndreev reflection (AR) between an electron and a hole in the source (or drain) lead forming theCP in the center, or to observe a crossed Andreev reflection (CAR) where an electron from thesource (drain) lead is coupled to a hole in the drain (source) lead through the CP in the center.Also, direct tunneling of an electron via the CP, referred to as cotunneling (CT), is a possibletransmission channel. Next, we simulate these different processes by a suitable parameter choice.We start with a familiar example by simulating NT: The condition is such that the biaswindow is larger than the gap. We will, in addition, fix the temperature, β = 100 / | t C | , wellbelow the critical temperature so that the gap can be approximated as the (constant) value atzero temperature ∆ ( T = 0) . The sample is a benzene molecule consisting of atomic sites ( carbon and hydrogen), coupled to the leads from four sites overall, see Fig. 4. The bias voltageis symmetrically set to V L = − V R = 3 | t C | / and the gap ∆ is varied but kept smaller than orequal to this value. The transient currents through the sample [calculated by summing the bondcurrents from Eq. (18) transversally in the middle of the molecule], the corresponding Fourier J ( t ) / Γ Fixed voltage V = 3 | t C | / ∆ = 0 ∆ = | t C | / ∆ = 3 | t C | / tΓ | P ( t ) | ω/ | tC ||F{ J ( t ) } ( ω ) | Figure 5: (Color online) Transient currents(top panel) and pair densities (bottompanel) in the molecule when varying ∆ .The inset shows the absolute value of theFourier-transformed current. -0.10-0.050.000.050.100.150.20 J ( t ) / Γ Fixed gap ∆ = | t C | / V = | t C | / V = | t C | / V = | t C | tΓ | P ( t ) | ω/ | tC ||F{ J ( t ) } ( ω ) | Figure 6: (Color online) Transient currents(top panel) and pair densities (bottompanel) in the molecule when varying V .The inset shows the absolute value of theFourier-transformed current.transforms and the pair densities [calculated by summing the pair densities within the molecule: P ( t ) = (cid:80) m P m ( t ) from Eq. (19)] can be seen in Fig. 5. Increasing the gap ∆ decreases theoverall current as the conducting states are being pushed away from the bias window. This alsoleads to shifts in the transient frequencies seen in the Fourier spectrum. By looking also at thespectral function A ( ω ) = − π Im Tr[ G R ( ω )] , where G R is the normal Nambu component of theretarded Green’s function, plotted in Fig. 8 we can further identify the transitions.In Fig. 5, when ∆ = 0 , we see two intramolecular transitions at frequencies ω = | t C | and ω = 2 | t C | which move a little when the gap is increased to ∆ = | t C | / corresponding to theshifted energy levels in the spectral function. We also observe two lead–molecule transitions at ω = | t C | / and ω = 5 | t C | / when ∆ = 0 ; these frequencies shift with the peaks in the spectralfunction corresponding to the fixed bias window at V = 3 | t C | / . The reason why we do not seea lead–molecule transition at ω = 3 | t C | / when ∆ = 0 , even though there is a zero-energy statein the molecule, is due to the fact that this state corresponds to the wavefunction’s nodal planesbeing located exactly at the lead interface, and therefore it is an inert state not taking part to thetransient dynamics[29, 30]. Also, as the conditions are for NT, we observe the pair density withinthe molecule going to zero from its equilibrium value when ∆ < V ; this means that there are noout-of-equilibrium CPs forming in the central region, and we see no AR or CAR processes. Whenwe set the gap equal to the bias window, we notice, first of all, that the steady-state current goesto zero since there are no transport channels within the bias window. Some transient oscillationsare still present due to the states in the vicinity of the resonant window, which is seen as anintramolecular transition at ω ∼ | t C | / .Next, we will adjust the transport conditions to visualize the AR and CAR processes. Wehave the same benzene molecule as the central region at the same temperature, β = 100 / | t C | . J ( t ) / Γ Fixed gap ∆ = 3 | t C | / V = | t C | / V = 3 | t C | / V = 5 | t C | / tΓ | P ( t ) | ω/ | tC ||F{ J ( t ) } ( ω ) | Figure 7: (Color online) Transient currents(top panel) and pair densities (bottompanel) in the molecule when varying V .The inset shows the absolute value of theFourier-transformed current. -3 -2 -1 0 1 2 3 ω/ | t C | − − A ( ω ) ∆ = 0 ∆ = | t C | / ∆ = 3 | t C | / Figure 8: (Color online) Spectral func-tions of the coupled benzene moleculewhen varying ∆ .In Fig. 6 we choose the gap as ∆ = | t C | / , and in Fig. 7 we, on the other hand, choose thegap as ∆ = 3 | t C | / . When V ≤ ∆ we observe, in both cases, the transient current oscillatingtowards a zero steady-state current. The oscillation frequencies seen in the Fourier spectrum canbe interpreted from the spectral function in Fig. 8 mainly as intramolecular transitions (around ω = | t C | and ω = 2 | t C | ). Interestingly, in Fig. 6 also, when we increase the bias voltage above thesuperconducting gap V = | t C | there still are no other states within the bias window except theinert state split by ∆ . This resonant window does not add anything to the transient dynamics(due to the inert state), but the static part of the density matrix is modified leading to nonzerosteady-state current. This is interpreted as CT where, in addition to AR and CAR, also anelectron is transferred from one lead to another via the CP. This is also confirmed by lookingat the pair density as it remains nonzero also for V = | t C | . With the larger gap in Fig. 7 andfor V ≤ ∆ we observe CP formation within the molecule but for V > ∆ the pair density goesto zero. For the smaller voltages we mainly find the first intramolecular transition at around ω = 7 | t C | / . For larger voltages we also see the lead–molecule transitions at lower frequencies,and we recover again the NT regime as the bias voltage is high enough for breaking the CPswithin the molecule.In this transport setup, it is not, in general, easy to distinguish between AR and CAR processesas they involve multiple steps. One possible case is when an electron hops from the lead tothe central region (molecule–lead transition), then “transfers” from the spin-up sector in thecentral region to the spin-down sector for which the probability is given by the pairing strength(“intramolecular” transition between the two branches of states split by ∆ ), and then finally hopsback to the lead as a hole (molecule–lead transition). As all these transitions are visible in thetransient oscillations, we may only conclude whether AR and CAR processes are present or not. . Conclusions and outlook We presented an extension to the time-dependent Landauer–Büttiker formalism, discussedin Refs. [28–30], to include superconducting central region in the transport setup, and toevaluate the TD1RDM at arbitrary temperatures. The derived formulae are analytic and closedexpressions involving known special functions, and they can readily be implemented to studyvarious quantum transport problems very efficiently and also at large temporal and spatial scales.As an application of the presented formalism we simulated transport in a superconductingbenzene-like molecule attached to two-dimensional normal metal leads. Assigning a properparameter set for the transport window and the superconducting gap, we observed formation ofCooper pairs within the central molecule leading to Andreev reflection processes.
Acknowledgments
R.T. thanks the Väisälä Foundation of The Finnish Academy of Science and Letters for financialsupport. R.v.L. thanks the Academy of Finland for support. E.P. and G.S. acknowledge fundingby MIUR FIRB Grant No. RBFR12SW0J.
References [1] Andreev A 1964
Sov. Phys. JETP Phys. Rev.
Phys. Rev.
Phys. Rev.
Phys. Rev. Lett. Phys. Chem. Chem. Phys. Phys. Rev. B Phys. Chem. Chem. Phys. Nanoscale Phys. Rev. Lett.
Phys. Rev. B Phys. Rev. B Phys. Rev. B Phys. Rev. Lett. Phys. Rev. Lett. EPL Phys. Rev. B Phys. Rev.B Phys. Rev. B IBM J. Res. Dev. Phys. Rev. Lett. J. Phys. C J. Phys. C Phys. Rev. B Phys. Rev. Lett. Phys. Rev. B. Phys. Rev. B Nonequilibrium Many-Body Theory of Quantum systems:A Modern Introduction (Cambridge Univerisity Press, 2013).[29] Tuovinen R, Perfetto E, van Leeuwen R and Stefanucci G 2013
J. Phys.: Conf. Ser.
Phys. Rev. B Phys. Rev. B Progress inNonequilibrium Green’s Functions ed Bonitz M, Nareyka R and Semkat D (World Scientific)[34] Stefanucci G, Perfetto E and Cini M 2010
Phys. Rev. B. Phys. Rev. B. Phys. Rev.
J. Phys.: Conf. Ser.
Sov. Phys. JETP Superconductivity of Metals and Alloys (New York: Benjamin)[40] Leijnse M and Flensberg K 2012
Semicond. Sci. Technol. J. Phys. Soc. Jpn. Phys. Rev. http://mathworld.wolfram.com/DigammaFunction.html [44] Weisstein E W Hypergeometric function. From MathWorld—A Wolfram Web Resource URL http://mathworld.wolfram.com/HypergeometricFunction.html [45] Kivelson S and Chapman O L 1983
Phys. Rev. B Phys. Rev. Lett. Phys.Rev. Lett.
Phys. Rev.B Phys. Rev. Lett. Phys. Rev.B Phys. Rev. B Phys. Rev. B Phys. Rev. B Nat Commun Phys. Rev. B Phys. Rev. Lett.101