Second Order Perturbation Theory for a Superconducting Double Quantum Dot

We extend our approach based on the second order perturbation theory in the Coulomb interaction recently developed for quantum dots coupled to superconducting leads to the superconducting double quantum dot setups. Using our perturbative method we evaluate several single-particle quantities such as on-dot induced gap and generalized occupations together with the Andreev in-gap spectra and compare them with numerically exact results from the Numerical Renormalization Group and Quantum Monte Carlo finding a very good correspondence for not too strongly correlated regimes. Thus we can offer in a wide parameter range this method as an efficient and reliable alternative to the heavy numerical tools exclusively used so far for the description of such experimentally relevant systems.


Introduction
Large number of experiments involving quantum dots attached to superconductors have been conducted in the past two decades [1]. Functional quantum dots in such setups are realized using a variety of systems (single molecules, carbon nanotubes, semiconducting nanowires etc.) and also various arrangements of several superconducting and/or normal leads exist. Parameters of such systems are often tunable, e.g., single-particle energies by the gate voltage and/or phase-difference in generalized Josephson junctions by the magnetic flux piercing the SQUID loops, which contributes to their versatility. The envisioned applications of such systems range from various sensors and detectors (e.g., single-molecule SQUIDs [2,3]) to building blocks of quantum information technologies [1].
While most of the experimental and theoretical studies have been thus far performed for single quantum dots attached to superconducting leads, in past several years the focus has partly shifted to double quantum dots (DQDs). Analogously to the single quantum dot setups there are many physical realizations of the DQD version using carbon nanotubes [2,4], semiconducting InAs [5][6][7][8] or InSb [9] nanowires, or molecular dimers on surfaces [10][11][12]. Apart from their relevance for fundamental physics, studying DQD systems is also a first step in understanding the behavior of longer nanowires which can host topologically non-trivial (Majorana) end states [9] and help to understand the behavior of the emerging Andreev bands.
Theoretical description of DQD setups uses either idealized limiting cases such as the superconducting atomic limit [13] assuming very large superconducting gap or master equation approach [14] with infinitesimally weak coupling to the leads (both of these assumptions are not realistic as in experiments the smallest energy scale involved is typically the superconducting gap) or heavy numerics such as Numerical Renormalization Group (NRG) for realistic experimental parameters [15]. This is the same situation as in the single dot setups before our discovery that a simple second-order perturbation theory in Coulomb interaction gives excellent results in a vast (and experimentally relevant) part of the parameter space [16,17].
In this work, we extend our second order perturbation theory (2PT) to the DQD setup and critically examine its performance by detailed comparison with two numerically exact techniques, NRG and the hybridization-expansion quantum Monte Carlo (CT-HYB). We find that, analogously to the single dot setups, the 2PT performs very well up to moderately strong interaction and can thus be used as a reliable and highly efficient tool for semi-quantitative description of DQD systems in a large portion of the parameter space. The Hamiltonian for the two-impurity Anderson model with two superconducting leads in the serial configuration sketched in Fig. 1 where the left lead is connected to the left dot and the right lead to the right one reads

The model Hamiltonian
where describes the two QDs and is the Hamiltonian describing a BCS superconducting lead and its coupling to the quantum dot.
Here d † iσ creates an electron on site i = L, R with spin σ and energy ε iσ , t σ is the inter-dot hopping amplitude, U i is the local Coulomb interaction (charging energy) on site i, c † ikσ creates an electron with spin σ and energy ε ikσ in lead i, ∆ i e iΦ i is the complex superconducting order parameter in lead i and V ikσ is the hopping between the lead i and the corresponding quantum dot.
For the sake of simplicity we neglect possible capacitive coupling U n L n R between the dots as well as the Rashba coupling that is present in semiconductor nanowires [13]. From now on we also consider the gap of the same size in both leads, ∆ L = ∆ R ≡ ∆, as the leads in an experiment are usually made from the same material and use ∆ as the energy unit in all results. We also assume spin-independent system with t ↑ = t ↓ ≡ t, V ik↑ = V ik↓ ≡ V ik , and ε ik↑ = ε ik↓ ≡ ε ik . In all our calculations we use constant tunneling densities of states 2 ) within the band of finite half-width D. Unfortunately, we cannot use the symmetry-asymmetry relation from the single quantum dot setup, where both leads are connected to the same dot, to map asymmetric coupling situations Γ L Γ R onto the symmetric one by a suitable change of the phase difference [18]. Yet, as in any Josephson junction also observables in this setup can only depend on the superconducting phase difference Φ = Φ L − Φ R , not on the absolute values of the two phases.
It is often useful to compare the results to the solution in the ∆ → ∞ superconducting atomic limit [19]. The Hamiltonian in this limit reads Its descrete spectrum gives a picture about the structure of the spin multiplets and helps to better understand the results from numerical methods.

Second order perturbation theory
The presented second-order perturbation theory is a direct generalization of the method introduced in Ref. [16,17] for a single superconducting quantum dot. This method can be straightforwardly applied to a non-interacting Green function describing a DQD that is introduced in the Appendix. It is based on a diagrammatic expansion in the powers of the Coulomb interaction U up to the second order. This simple and fast method provides a reliable description of the system in the weak and intermediate interaction regimes provided the ground state is a singlet, as this method fails for degenerate ground states due to the violation of the Gell-Mann-Low theorem. However, the calculations performed in the ∆ → ∞ superconducting atomic limit as well as NRG calculations [15] suggest that the ground state of a double quantum dot system is a singlet in the vast part of the parameter space, making the 2PT method usable in most situations.
The method was implemented in Matsubara (imaginary frequency) formalism using the TRIQS libraries [20]. All calculations were performed at a small finite temperature k B T = 10 −2 ∆ with a cutoff in imaginary frequencies ω max n ≥ 500∆. The spectral functions were obtained by analytic continuation to the real frequency domain, G(iω n ) → G(ω + iη) using the Padé approximation [21]. The continuation was performed from the first 50 Matsubara frequencies as adding more frequencies did not change the result in any way. A small imaginary part η = 10 −3 ∆ was added to the real frequency to guarantee the correct analytic properties of the continued function.

Numerical renormalization group
We compared our results obtained with 2PT against the zero temperature data calculated by the numerical renormalization group [22,23], using the open source package NRG Ljubljana [24,25]. The logarithmic discretization parameter was set to 4, the maximum number of states kept in the truncation was 4000 (times the multiplicity), the cut-off energy was set to 10 in the units of characteristic NRG energy scale and the minimal number of kept states was set to 1000.

Quantum Monte Carlo
We also used the TRIQS/CTHYB continuous-time, hybridization-expansion quantum Monte Carlo solver [26] to cross-check the NRG results and assess the effects of the finite temperature on the system. Hamiltonian (1) does not conserve the electron number and therefore cannot be solved directly using the standard CT-HYB technique. To circumvent this problem we utilized a canonical particle-hole transformation in the spin-down sector as described in Ref. [27]. The continuous-time quantum Monte Carlo is an inherently finite-temperature method and all calculations were performed at k B T = 0.05∆. If we consider a typical gap of an InAs nanowire proximitized to Al, ∆ ≈ 200µeV [5], this temperature corresponds to T ≈ 120mK. Some of the results were recalculated at k B T = 0.025∆, showing very little temperature dependence of the measured quantities.

Results
All presented results were calculated for a symmetric setup U L = U R ≡ U, ε L = ε R ≡ ε and Γ L = Γ R ≡ Γ, for fixed coupling strengths Γ = 2∆, t = ∆ and zero phase difference Φ = 0. The half-bandwidth of the flat tunneling density of states in the leads is set to D = 100∆ in all numerical calculations. The numerical values for λ in almost the whole range are less precise, yet the overall shape of the curve follows closely the exact numerics including the existence and position of its minimum.
In Fig 3a we plotted the in-gap normal spectral function ρ(ω) ≡ − Im G 11 (ω + i0)/π (see Appendix) obtained from the imaginary-frequency 2PT solution by analytic continuation to the real axis using the Padé approximation and compared it to the excitation spectrum calculated by NRG. Only the singlet-doublet transition (red line) is visible in the one-electron spectral function (color map). The singlet-triplet transition (blue dashed) and singlet-other singlet transition (black dashed) violate the ∆S z = 1/2 selection rule and therefore are invisible in the one-electron spectral function. The 2PT correctly predicts the general behavior of the Andreev bound states (e.g. the minimum in their energy as function of U) and their position for small values of U, but naturally fails to describe their exact position in the strong interaction regime.
In order to better understand the behavior of the Andreev bound states, we also calculated the spectrum of the atomic Hamiltonian H ∞ (4) using the atomic solver implemented in TRIQS and plotted it in Fig 3b. We see that the ground state in the wide-band limit is always a singlet and the first excited state for small U is a degenerate pair of doublets. This degeneracy can be lifted by applying gate voltage ε, leading to a splitting of the Andreev bound states (see below). The second pair of doublets lies high in energies and the transition to them falls within the continuous band of the finite-∆ model. Interaction strength U is a material-dependent quantity and cannot be changed easily in the experiment. On the other hand, local energy level ε can be tuned by changing the gate voltage. In Fig. 4 we plotted the behavior of a DQD as a function of the energy level for two values the interaction strength U = 5∆ and U = 10∆ (marked by arrows in Fig 2a). The correlation effects get weaker as we move away from half-filling (ε = −U/2) because the average dot occupation n ≡ d † L↑ d L↑ = d † R↑ d R↑ (panel a) is decreasing. As a result, 2PT data agree with the NRG better as we increase ε.
Panel d shows the in-gap spectral function ρ(ω). The ground state is again a singlet and we see two pairs of Andreev bound states corresponding to the singlet-doublet transitions. This is a result of the splitting of the two degenerate doublets by the applied gate voltage mentioned earlier in the text. As the average dot occupation is decreasing, the Andreev bound states move closer to the gap edges, eventually merging into the continuum, leaving the spectral function free of any in-gap states. The singlet-triplet (blue) and singlet-other singlet (black) transitions are again invisible to the one-particle spectral function as mentioned in the previous text and Fig. 3a.

Conclusions
We presented a fast and simple perturbation method to calculate properties of a DQD system with superconducting leads and benchmarked it against exact numerical techniques as the numerical renormalization group and the quantum Monte Carlo. It is usable in the weak and intermediate correlation regimes, provided that the ground state is a singlet, which is true in most situations. In the strongly correlated regime it can still provide qualitatively correct description of the system. The main limitation of this approach lies in its inability (so far) to describe phases with degenerate ground states (doublet or triplet in DQD) and, consequently, also to provide finite temperature results close to phase transitions where the adjacent phases thermally mix [28]. On the other hand, this method can be straightforwardly generalized to more complicated setups with various geometries, metallic leads and the inter-dot capacitive coupling. where ω n = (2n + 1)πk B T ,Î is a 4 × 4 unit matrix, Hereε describes the local energy levels and hoppings in the isolated DQD andΓ i (iω n ) is the hybridization function describing the coupling between the quantum dot i = L, R and the corresponding superconducting lead with 2/π arctan D/ ∆ 2 + ω 2 n being a correction due to finite bandwidth D.