Correlation effects in superconducting quantum dot systems

TL;DR

This paper investigates how electron correlations influence a quantum dot coupled to superconductors, using advanced computational methods to compare different theoretical approaches for realistic parameters.

Contribution

It introduces a comprehensive analysis of correlation effects in superconducting quantum dots using quantum Monte Carlo, perturbation theory, and NRG, highlighting their comparative strengths.

Findings

Correlation effects significantly modify the quantum dot's electronic properties.

Quantum Monte Carlo provides accurate results for experimentally relevant parameters.

Comparison reveals strengths and limitations of perturbation theory and NRG.

Abstract

We study the effect of electron correlations on a system consisting of a single-level quantum dot with local Coulomb interaction attached to two superconducting leads. We use the single-impurity Anderson model with BCS superconducting baths to study the interplay between the proximity induced electron pairing and the local Coulomb interaction. We show how to solve the model using the continuous-time hybridization-expansion quantum Monte Carlo method. The results obtained for experimentally relevant parameters are compared with results of self-consistent second order perturbation theory as well as with the numerical renormalization group method.

Equations

$$\mathcal{H}={\mathcal{H}}_{dot}+\sum _{s=L,R}({\mathcal{H}}_{lead}^s+{\mathcal{H}}_c^s)$$

$${\mathcal{H}}_{dot}=\sum _{\sigma }{\epsilon }_{\sigma }^{}{d}_{\sigma }^{†}{d}_{\sigma }^{}+U{d}_{\uparrow }^{†}{d}_{\uparrow }^{}{d}_{\downarrow }^{†}{d}_{\downarrow }^{}.$$

$$\begin{array}{rl}{\mathcal{H}}_{lead}^s& =\sum _{\mathbf{k}\sigma }\epsilon (\mathbf{k})c_{s,\mathbf{k}\sigma }^{†}{c}_{s,\mathbf{k}\sigma }^{}\\ & -\Delta \sum _{\mathbf{k}}(e^{i{\Phi }_s}{c}_{s,\mathbf{k}\uparrow }^{†}{c}_{s,-\mathbf{k}\downarrow }^{†}+\text{H.c.})\end{array}$$

$${\mathcal{H}}_c^s=-\sum _{\mathbf{k}\sigma }{t}_s(c_{s,\mathbf{k}\sigma }^{†}{d}_{\sigma }^{}+\text{H.c.})$$

$$\begin{array}{rl}{d}_{\uparrow }^{†}& \to d_{\uparrow }^{†},d_{\downarrow }^{†}\to d_{\downarrow }^{},d_{\uparrow }^{}\to d_{\uparrow }^{},d_{\downarrow }^{}\to d_{\downarrow }^{†},\\ c_{\mathbf{k}\uparrow }^{†}& \to c_{\mathbf{k}\uparrow }^{†},c_{\mathbf{k}\downarrow }^{†}\to c_{-\mathbf{k}\downarrow }^{},c_{\mathbf{k}\uparrow }^{}\to c_{\mathbf{k}\uparrow }^{},c_{\mathbf{k}\downarrow }^{}\to c_{-\mathbf{k}\downarrow }^{†},\end{array}$$

Full text

††journal: Physica B: Condensed Matter

1 Introduction

Conventional Josephson junctions had become a standard building blocks of various electronics devices including SQUIDs [1], RSFQs [2], and qubits [3] in quantum computing. No wonder that their tunable generalizations, the superconducting quantum dots, gain a lot of attention from both theorist and experimentalist. These hybrids, where a quantum dot is placed between two superconducting leads, promise a great deal of technological advances such as quantum supercurrent transistors [4], monochromatic single-electron sources [5] or single-molecule SQUIDs [6]. Of no less importance is the fact that they are rich and relatively easy to deal with playgrounds for studying various physical phenomena. These include the competition between Kondo effect and superconductivity [7], the Andreev subgap transport [8] as well as quantum phase transitions from impurity spin-singlet to spin-doublet ground state which are in the experiments signaled by the sign reversal of the supercurrent (0-$\pi$ transition) and accompanied by a crossing of the subgap Andreev bound states (ABS) [9, 10, 11].

The superconducting quantum dot can be adequately described by a single impurity Anderson model (SIAM) coupled to BCS leads [12]. Consequently, a lot of different theoretical approaches have been applied to study this system. Many important results have been obtained using various (semi)analytical methods based on different perturbation approaches [8, 13, 14, 15]. Moreover, as was shown in recent studies [10, 16], a surprisingly large portion of the parametric space of the superconducting SIAM can be reliably covered with a properly formulated second-order perturbation theory (2ndPT) in the on-dot electron Coulomb interaction. Unfortunately, this method cannot describe the $\pi$-junction behavior due to the ground-state degeneracy.

None of the mentioned (semi)analytical perturbative methods can cover all experimentally relevant cases. Therefore there is a big demand for “heavy” numerical methods. A very good quantitative agreement with the experiments can be obtained with the numerical renormalization group (NRG) [17, 9] and quantum Monte Carlo (QMC) [18, 12] methods. Although both methods have some disadvantages, including their computational demands, they have, besides parametric universality, one big practical advantage. Namely, the existence of well-tested versatile open-source software packages.

In the present paper we focus on the continuous-time hybridization-expansion quantum Monte Carlo (CT-HYB) [19] implementation for experimentally inspired parameters [9] representing a strong coupling regime, which is beyond the reach for most (semi)analytical techniques.

We show how to include the superconductivity into CT-HYB quantum Monte Carlo solver. Then we study various single-particle quantities as functions of the gate voltage. We discus how they behave near quantum phase transition and show that the CT-HYB can be reliably used to obtain the phase diagram. We also use numerical analytical continuation to obtain the spectral function. We compare all obtained CT-HYB results with either 2ndPT or the NRG.

2 The model Hamiltonian

We describe the system by the single-impurity Anderson model with BCS leads. The Hamiltonian reads

[TABLE]

where $s=L,R$ denotes the left and right leads. The impurity Hamiltonian describes a single-level atom with local Coulomb repulsion $U$ and on-site energy $\varepsilon_{\sigma}=\varepsilon+\sigma B$, where $B$ is the external magnetic field

[TABLE]

The Hamiltonian of the BCS leads reads

[TABLE]

where $\Delta e^{i\Phi_{s}}$ is the complex gap parameter. We assume the same gap size in both leads, $\Delta_{L}=\Delta_{R}=\Delta$, meaning that the leads are made from the same material, as it is usual in the experimental setups. Finally, the coupling part reads

[TABLE]

where $t_{s}$ denotes the tunneling matrix element.

Hamiltonian (1) does not conserve the electron number and therefore cannot be solved directly using standard CT-HYB technique. To circumvent this problem we utilized a canonical particle-hole transformation in the spin-down sector

[TABLE]

previously used by Luitz and Assaad [18] to include superconductivity in the continuous-time interaction-expansion (CT-INT) QMC calculations. The new quasiparticles are identical to electrons in the spin-up sector and to holes in the spin-down sector. This transformation maps our system to SIAM with attractive interaction $-U$ and off-diagonal hybridization of the quantum dot with the leads. The local energy levels transform as $\varepsilon_{\sigma}\rightarrow\sigma\varepsilon_{\sigma}$. Since $\varepsilon_{\sigma}=\varepsilon+\sigma B$ and $\sigma^{2}=1$, this transformation maps the local energy $\varepsilon$ on the magnetic field $B$ and vice versa. The dispersion and tunneling matrix elements transform in the same manner, $\varepsilon(\mathbf{k})\rightarrow\sigma\varepsilon(\mathbf{k})$ and $t_{s}\rightarrow\sigma t_{s}$. The resulting Hamiltonian conserves the total electron number and can be solved using standard CT-HYB implementations.

3 The CT-HYB method

We use the TRIQS/CTHYB Monte Carlo solver [20, 21]. We consider a flat density of states in the leads of finite half-width $D=30\Delta$. The coupling of the quantum dot to the leads is described by tunneling rates $\Gamma_{s}=\pi|t_{s}|^{2}/(2D)$. We denote $\Gamma=\Gamma_{R}+\Gamma_{L}$ and consider only the symmetric coupling $\Gamma_{R}=\Gamma_{L}$. Any asymmetric coupling $\Gamma_{R}\neq\Gamma_{L}$ with the same total $\Gamma$ can be easily gained from the symmetric solution using a simple analytical relation derived in Ref. [22].

Continuous-time quantum Monte Carlo belongs to a family of inherently finite-temperature methods and the calculations are usually restricted to rather high temperatures. However, since the typical energy scale in our setup is the superconducting gap $\Delta\sim 100\mu$eV, it allows us to easily reach experimental range of temperatures $T\sim 10-100$mK.

The biggest disadvantage of CT-HYB in comparison with NRG or 2ndPT is that the calculation is performed on the imaginary-time axis. Obtaining the spectral function from imaginary-time data is a well-known ill-defined problem. Various numerical methods are used to perform the analytic continuation, the maximum entropy method being the most common one [23]. However, this method fails to resolve sharp spectral features like the Andreev bound states. Therefore we use the Mishchenko’s stochastic optimization method (SOM) [24] in its recent implementation [25] which is better suited to our needs.

4 Results

Our calculations are inspired by the experiment of Pillet et. al. [9]. The paper describes the tunneling spectroscopy measurement performed on a carbon nanotube connected to superconducting aluminum leads. Experimental results show the Andreev bound states as functions of gate voltage and are nicely reproduced using the NRG method. The superconducting gap is $\Delta=150\mu$eV, Coulomb interaction $U\approx 2$meV and the phase difference is zero ($\Phi_{L}=\Phi_{R}$). We use these parameters in our calculations and set the magnetic field $B$ to zero. It is worth to note that we did not encounter any fermionic sign problem during the calculation.

In Fig. 1 we plot the diagonal (panel a) and the off-diagonal (panel b) part of the occupation matrix as functions of the shifted local energy level $\varepsilon_{U}=\varepsilon+U/2$ ($\varepsilon_{U}=0$ represents the half-filled dot). From now on we use $\Delta$ as the energy unit. We restrict ourselves to positive values of $\varepsilon_{U}$ as the rest can be determined from symmetry. We chose parameters $U=13.3\Delta$ and $\Gamma_{L}=\Gamma_{R}=0.45\Delta$ which are within the experimental range. The three solid lines are CT-HYB results calculated at inverse temperatures $\beta\Delta=10$ (green), $20$ (blue) and $40$ (red) that correspond to temperatures $T=175$ mK, $77$ mK and $44$ mK, respectively. The diagonal part corresponds to the electron density $n=\langle d^{{\dagger}}d\rangle$. It varies very weakly in $\pi$ phase, then changes abruptly at the phase transition.

The position of the phase transition can be more easily determined from the off-diagonal part, which represents the induced gap $\mu=\langle d^{{\dagger}}d^{{\dagger}}\rangle$. This parameter is negative in the $\pi$-phase and positive in the [math]-phase. We see that the temperature has very little effect on the position of the phase transition which takes place at $\varepsilon_{U}\approx 5.5\Delta$. We included also the results of the 2ndPT method, which is available only in the [math]-phase. It fits very well the CT-HYB results in this phase although it gives the phase transition at $\varepsilon_{U}=6.15\Delta$ (c.a. $12\%$ error). However, this discrepancy is expected as we are investigating a strong coupling regime ($U/\Delta\gg 1$ and $U/\Gamma\gg 1$).

The inset of panel a shows the average perturbation order $\langle k\rangle$ scaled by the inverse temperature $\beta$ for the CT-HYB calculations in the main panels. This quantity is an estimator of the kinetic energy [26]. It scales linearly with $\beta$ and exhibits a maximum just above the phase transition point. Although the zero temperature extrapolation of the position of maxima could be in principle used to estimate the phase transition point, safe determination of its position from $\langle k\rangle$ requires a rather elaborate procedure [27].

Calculation of a spectral function requires much more precise QMC data than the calculation of an expectation value due to the underlying, ill-defined analytic continuation procedure. While the expectation values with reasonable error bars can be obtained within few CPU-hours, calculation of a spectral function, including the SOM procedure, takes usually more than 100, depending strongly on the temperature and the coupling strength $\Gamma$. In Fig. 2 we plot the color map of the spectral function at $\beta\Delta=40$ ($T=44$mK) in the vicinity of the gap region as it can be directly compared to the experimental data. We use the same parameters as in Fig. 1 and compare it with the position of ABS calculated using NRG and 2ndPT methods at $T=0$. NRG results were obtained using NRG Ljubljana code [28]. The ingap maxima of the spectral function in the [math]-phase and around the transition point are in very good agreement with the positions of ABS calculated by NRG. In the $\pi$-phase the position of the maxima tends to shift to higher energies. This is surprising as the electron density and the induced gap values are in good agreement with NRG even in this region. The position of the peaks also does not depend on the temperature and while it does depend on the width of the non-interacting band $D$, this dependence is weak and it effects the position of the maxima equally in [math] and $\pi$ phase, therefore it cannot explain this discrepancy.

In order to get some insight into this problem we plotted in Fig. 3 the spectral functions calculated using NRG and CT-HYB methods. The top panel shows results for $\varepsilon_{U}=7\Delta$ which is in [math]-phase. We see that the arrows that represent ABS from NRG calculation match the maxima of the spectral function obtained using SOM procedure from CT-HYB data. We also see that CT-HYB spectra are missing the structure just above the gap edges at $\pm\Delta$. Bottom panel shows spectral functions for $\varepsilon_{U}=2\Delta$ ($\pi$-phase). The mismatch between the arrows and the maxima is clearly visible and we do not have a satisfactory explanation of this discrepancy.

The local energy level $\varepsilon$ is a parameter that can be easily tuned in the experimental setups by changing the gate voltage. On the other hand, the tunneling rate $\Gamma$ is very hard to measure and it is usually obtained from numerical fits [9]. Studying the relation between these parameters is therefore important for interpretation of the experimental results. In Fig. 4 we plotted the phase diagram in the $\varepsilon_{U}-\Gamma$ plane. Red solid line represents the phase boundary calculated using CT-HYB at inverse temperature $\beta\Delta=20$ that corresponds to $T=77$ mK. This boundary was determined from the positivity of the induced gap $\mu$. This line does not change with further decreasing temperature beyond the resolution of the plot. We also included 2ndPT result for comparison. As this perturbation expansion is performed in $U/\Gamma$ parameter it differs more for small values of $\Gamma$ where it develops a “hump” as already pointed out in Ref. [16]. The two lines then meet for $\Gamma=0$ at the exact result $\varepsilon_{U}=U/2$. The blue arrow marks the cut at $\Gamma=0.9\Delta$ along which the data in Figs. 1 and 2 are plotted.

5 Conclusions

We studied a $0-\pi$ quantum phase transition in a single-level quantum dot connected to two superconducting BCS leads using the continuous-time hybridization-expansion quantum Monte-Carlo method. We used the $U$ and $\Delta$ parameters from experiment described in Ref. [9] in order to stay in a realistic region of the parameter space. Performing an electron-hole transformation in the spin-down sector we mapped the system on a model that can be solved using CT-HYB method as implemented in the TRIQS package. We presented results as functions of the gate voltage $\varepsilon$ as this parameter is easily tunable in the experiment. We showed how the $0-\pi$ quantum phase transition point can be extracted from the behavior of the induced gap and presented the finite-temperature spectral function as well as the phase diagram in the $\varepsilon-\Gamma$ plane that can be used to determine the value of the tunneling rate $\Gamma$.

In summary, we showed that CT-HYB is an effective method for studying superconducting quantum dot systems, where the interaction strength is the dominant energy scale. The present formulation is sign problem free and one can access the low-temperature region using reasonable amount of computational resources. We also showed how the spectral function can be obtained using analytic continuation based on the Mishchenko’s stochastic sampling method in order to study the behavior of the subgap Andreev bound states. Comparing the position of the subgap maxima with ABS frequencies from the NRG calculation shows good agreement in the [math]-phase but a discrepancy in the $\pi$-phase which is of unknown origin. Furthermore, the model can be generalized to include a normal (non superconducting) electrode. As already pointed out in Ref. [7] where such a three-terminal device was studied, NRG and 2ndPT methods fail in this setup except special cases and CT-HYB becomes the method of choice.

Acknowledgments

Research on this problem was supported by Grant No. 15-14259S of the Czech Science Foundation (V.P.) and the Grant No. DEC-2014/13/B/ST3/04451 of the National Science Centre (Poland) (M.Ž.). Access to computing and storage facilities owned by parties and projects contributing to the National Grid Infrastructure MetaCentrum provided under the programme “Projects of Large Research, Development, and Innovations Infrastructures” (CESNET LM2015042), is greatly appreciated. V. P. thanks Roberto Mozara for the help with the stochastic optimization method.