Dephasing due to quasiparticle tunneling in fluxonium qubits: a phenomenological approach

Among the various types of superconducting qubits [1–3], the recently developed fluxonium qubit [4, 5] has the unique advantage of being protected against both charge and flux noise. This is important since both effects in general limit the performance of the qubits by introducing relaxation and dephasing processes. Indeed, over the last few years considerable effort has been made in order to understand, and consequently to reduce, the causes of relaxation and decoherence in different types of superconducting circuits [1, 6–9].

Initially, the fluxonium qubit has been designed in order to reduce the sensitivity of the Cooper pair box to charge noise [4, 5]. Subsequently, it has been argued that it is also insensitive to flux noise [9]. In order to reach a regime in which these relaxation and decoherence processes are exponentially suppressed, the charging energy ${E}_{\rm{C}}={e}^{2}/2C$, with the capacitance C and the elementary charge e, has to be much larger than the inductive scale ${E}_{\rm{L}}={e}^{2}/4{\alpha }^{2}L$, with the inductance L and the fine-structure constant $\alpha ={e}^{2}/{\hslash }c$. While it is impossible to realize such large 'superinductances', $L\gg {\alpha }^{-2}C\simeq {10}^{4}C$, with conventional media, in the fluxonium qubit it has been realized by an array of (large) Josephson junctions [5]. A downside of this approach is that the device is then potentially plagued by spurious phase slips through the array; however, these have been successfully eliminated [10, 11]. Recently, different detrimental effects in the fluxonium qubit due to nonequilibrium quasiparticles have been addressed theoretically [12–16] and experimentally [17].

Of all the processes limiting the performance of the qubits discussed above quasiparticle tunneling is particularly important as it is intrinsic to the superconducting tunneling junction and as such forms an absolute limit. While relaxation and dephasing mechanisms due to quasiparticle tunneling have been extensively studied in charge qubits [18–23], for the flux qubits most studies so far have concentrated on the relaxation due to quasiparticle processes [12, 13, 17]. On the other hand, first investigations of the dephasing due to quasiparticle tunneling on a perturbative level in the junction transmission have been put forward in [14–16]. It turns out that treating the problem in perturbation theory leads to a diverging result due to the sharp peak of the quasiparticle density of states at the gap. This problem has been addressed by introducing as a cutoff a relaxation rate [15] or a dephasing rate [14, 16] that broaden the density of states, where the latter has been determined self-consistently. It has been discussed that in principle the divergence is lifted by treating the tunneling nonperturbatively as the divergence simply signals the presence of a weakly bound Andreev state close to the gap [14].

Recently, it has been recognised that the presence of different nonequilibrium quasiparticle distributions on the different superconducting islands of the qubit (possibly accidentally arising during operation of the qubit [22] and resulting in an effective temperature gradient of stationary nonequilibrium quasiparticle distributions) can lead to an additional decoherence mechanism for the flux qubits, see [24], where in particular the flux qubit in the Delft-qubit design has been addressed. It has been shown that in this case a limitation of the dephasing time arises caused by heat currents carried by quasiparticles which flow through the Josephson junctions of a superconducting qubit as a response to the (effective) temperature gradient.

The microscopic origin of this dephasing mechanism is the fact that a heat current flowing through a Josephson junction depends on the phase difference of the two superconductors separated by the junction. This fundamental effect has been predicted over 50 years ago, see [25], and has later been studied in more detail in [26–29]. Only very recently, the phase dependence of heat currents through Josephson junctions has been measured experimentally [30–32]. Intriguingly, due to the general phase-dependence of the heat current, heat currents flowing across junctions of superconducting flux qubits can depend on the qubit state. This results in the dephasing of the qubits, since heat currents are dissipative [24].

In this paper, we investigate pure dephasing of the fluxonium qubit due to a nonequilibrium quasiparticle distribution. We in particular analyze the relevance of this intrinsic dephasing mechanism (without population relaxation) slightly away from the sweet spot, where it would be suppressed, see [13–15]. We model the nonequilibrium distribution of quasiparticles by an effective temperature which is different on different superconducting segments of the qubit. We study in details the effect of the heat current flowing both through the Josephson junction of the so-called black-sheep junction with Josephson energy E_J, which is the main source of nonlinearity, as well as through the array of larger junctions constituting the superinductance, which is shunting the black sheep. In particular, we are interested in the linear response regime where the difference in temperatures is small. We show that the sensitivity of the fluxonium qubit to heat transport and the resulting pure dephasing is suppressed by a factor ${E}_{\rm{L}}^{2}{/{E}_{\rm{J}}}^{2}$ (up to logarithmic corrections) rendering the fluxonium qubit rather insensitive to this dephasing source. We point out that even if the effective temperature gradient is vanishingly small, the study of the heat conductance (rather than of the heat currents) is relevant: intriguingly, it can be used as a phenomenological approach for the investigation of dephasing due to nonequilibrium quasiparticles⁵ . Different from our phenomenological approach, prior work on dephasing in flux-based qubits [13–16] was based on a microscopic model. Our approach has the advantage that we are able to consider very small temperature gradients—corresponding to superconducting qubit segments with identical superconducting gaps. In this regime, the dominating effect of the heat-current sensitivity to the qubit state stems from the phase-dependence of a weak bound state originating from Andreev reflection. Due to the fact that our approach is nonperturbative in the tunnelling coupling, we are not plagued by divergencies as the weak Andreev bound state acts as a natural phase-dependent cut-off. Additionally, our approach provides a nice link of this intrinsic dephasing mechanism to thermal transport quantities and their ability to 'measure' a qubit state.

The paper is organized as follows: in section 2, we introduce the generic Hamiltonian of the fluxonium qubit. We recall the results for heat currents flowing through a single Josephson junctions for very small temperature gradients in section 3. These findings are then used in section 4 to investigate on the sensitivity of the heat current on the fluxonium states. In section 5, we show in how far the resulting dephasing can limit the operation of the fluxonium qubit.

The fluxonium is a superconducting qubit which can be thought of as a Cooper-pair box inductively shunted with a superinductance [4]. Its electrical circuit is shown in figure 1(a): it consists of a Josephson junction, which we refer to as the black sheep, with Josephson energy E_J and charging energy E_C. This Josephson junction is shunted by an array of M larger Josephson junctions (with Josephson coupling energy ${E}_{\rm{J}}/\beta$, $\beta \lt 1$). When operated at microwave frequencies well below its self-resonant frequency $\sqrt{{E}_{\rm{J}}{E}_{\rm{C}}/\beta }/{\hslash }$, this array emulates a 'superinductance' for sufficiently large M. It is the presence of this superinductance which renders the fluxonium insensitive to charge noise while yielding a highly anharmonic spectrum, both important points for the realization of a qubit.

Zoom In Zoom Out Reset image size

2.1. Hamiltonian

The Hamiltonian of the fluxonium consists of different parts, stemming from the charging and the Josephson energy of the different junctions [34, 35]. Charging effects on the black sheep junction lead to the charging energy $T=4{E}_{\rm{C}}{\hat{n}}_{\rm{b}}^{2}$ with ${\hat{n}}_{\rm{b}}$ the number of Cooper pairs on the capacitance plate of the black sheep. In the Hamiltonian, the charging energy of the black sheep serves as the kinetic energy, while the sum of the different Josephson energies takes the role of the potential energy. In particular, the potential energy of the black sheep is given by ${U}_{\rm{b}}={E}_{\rm{J}}(1-\mathrm{cos}{\varphi }_{\rm{b}})$ with ${\varphi }_{\rm{b}}$ the superconducting phase difference across the black-sheep junction which is the conjugate variable to ${\hat{n}}_{\rm{b}}=-\rm{i}\partial /\partial {\varphi }_{\rm{b}}$. For simplicity, we assume the junctions in the array constituting the superinductance to be equal with a somewhat larger Josephson energy ${E}_{\rm{J}}/\beta$, $\beta \lt 1$. The total potential energy is thus given by $U={U}_{\rm{b}}+M({E}_{\rm{J}}/\beta )(1-\mathrm{cos}{\varphi }_{M})$ with ${\varphi }_{M}$ the phase difference across each of the junctions constituting the superinductance. We neglect the capacitive energies of the M array junctions since their larger area translates into smaller charging effects. The radius r of the superconducting loop is supposed to be so small that the inductive energy ${\Phi }_{0}^{2}/r$ per flux quantum ${\Phi }_{0}={hc}/2e=\pi e/\alpha$ dominates the Josephson energies E_J. In this regime, the magnetic flux in the superconducting loop is quantized which leads to the condition ${\varphi }_{\rm{b}}+M{\varphi }_{M}+2\pi f=0\ (\mathrm{mod}\;2\pi )$ with the dimensionless parameter $f=\Phi {/\Phi }_{0}$ due to the magnetic flux Φ accounting for the external magnetic field penetrating the loop. Using this relation, we can express the Hamiltonian in terms of one dynamical variable, only, which we choose to be the total phase difference across the superinductance, $\varphi =M{\varphi }_{M}=-{\varphi }_{\rm{b}}-2\pi f$. Note that at fixed φ, the phase difference ${\varphi }_{M}$ across each junction in the array constituting the superinductance becomes small for large M. As a result of the small phase difference across each of the large junctions⁶ , we are allowed to expand the cosine in the potential energy to second order in ${\varphi }_{M}$ with the result

$$\begin{eqnarray}&&U=-{E}_{\rm{J}}\mathrm{cos}{\varphi }_{\rm{b}}+{\displaystyle \frac{1}{2}}{E}_{\rm{L}}{\varphi }^{2};\end{eqnarray} \tag{ 1 }$$

(up to an irrelevant constant shift). Here, we have introduced the inductive energy of the superinductance

$$\begin{eqnarray}&&{E}_{\rm{L}}={\displaystyle \frac{{E}_{\rm{J}}}{\beta M}}={\displaystyle \frac{{({\Phi }_{0}/2\pi )}^{2}}{{L}_{\mathrm{eff}}}},\end{eqnarray} \tag{ 2 }$$

with the effective inductance L_eff, which increases linearly with M. It is the inductive term, which we obtained from expanding the cosine potential of the array to second order in ${\varphi }_{M}$, which breaks the $2\pi$ periodicity of the total potential U as a function of φ, thus rendering the device charge insensitive [4, 35].

In conclusion, the total Hamiltonian $H=T+U$ of the fluxonium is given by

$$\begin{eqnarray}&&H=4{E}_{\rm{C}}{\hat{n}}^{2}+{\displaystyle \frac{1}{2}}{E}_{\rm{L}}{\varphi }^{2}-{E}_{\rm{J}}\mathrm{cos}(\varphi +2\pi f)\end{eqnarray} \tag{ 3 }$$

with $\hat{n}=-\rm{i}\partial /\partial \varphi =-{\hat{n}}_{\rm{b}}$ the operator conjugate to φ. The Hamiltonian (3) for $f\approx 1/2$ may serve as a model-Hamiltonian for different qubit types [34]. In the present manuscript, we mainly discuss the fluxonium qubit which is realized in the regime ${E}_{\rm{J}}\gg {E}_{\rm{C}}\gg {E}_{\rm{L}}$, where the latter inequality is made possible by the presence of the superinductance. We compare the result to the case of the Delft qubit realized in the regime ${E}_{\rm{L}}\simeq {E}_{\rm{J}}\gg {E}_{\rm{C}}$; see table 1 for a discussion of the different energy scales.

Table 1.  Characteristic energy ratios for the fluxonium [5] and Delft qubit [36] when described by the Hamiltonian (3) for the dynamic degree of freedom. The ratio ${E}_{\rm{J}}{/E}_{\rm{C}}$ approximately corresponds to the number of levels in a potential minimum, whereas ${E}_{\rm{J}}{/E}_{\rm{L}}$ is a measure of the number of mimima of the potential which typically contribute to the qubit states [34].

	${E}_{\rm{L}}{/E}_{\rm{J}}$	${E}_{\rm{C}}{/E}_{\rm{J}}$
Fluxonium	$6\times {10}^{-2}$	$3\times {10}^{-1}$
Delft qubit	$5\times {10}^{-1}$	$2\times {10}^{-2}$

2.2. Qubit states

Close to $f=1/2$, the potential landscape is given by a double-well potential, which is an ideal starting ground to encode a qubit, see figure 1(b). For the fluxonium qubit, we find numerically that the potential landscape features two well-localized minima for $| \delta f| \lesssim 0.3$ with $\delta f=f-1/2$. The two minima are situated at ${\varphi }_{\rm{L}/\rm{R}}$. The two lowest-lying states in these minima, representing the qubit states, correspond to a current flowing clockwise/counter-clockwise in the device. However, due to the large inductance the corresponding currents are rather small, even though the two states differ by a large magnetic flux.

The dephasing mechanism which we investigate in the following arises due to the fact that the different semiclassical states correspond to different superconducting phases on the islands which in turn makes the phase-dependent heat currents through the junctions dependent on the qubit state. The resulting decoherence thus projects the qubit on the semiclassical states ${\varphi }_{\rm{L}/\rm{R}}$.

We now exactly look at these semiclassical states. The position of the minima ${\varphi }_{\rm{L}/\rm{R}}$ is given by the solution $\bar{\varphi }$ to the equation

$$\begin{eqnarray}&&0={\displaystyle \frac{\partial U}{\partial \varphi }}={E}_{\rm{L}}\bar{\varphi }+{E}_{\rm{J}}\mathrm{sin}(\bar{\varphi }+2\pi f).\end{eqnarray} \tag{ 4 }$$

At $\delta f=0$, the solutions are situated symmetrically around $\varphi =0$ with ${\varphi }_{\rm{R}}=-{\varphi }_{\rm{L}}={\varphi }^{*}$. In general, this transcendental equation can only be solved numerically. Therefore, in the following, we present all analytical results in the regime ${E}_{\rm{L}}\ll {E}_{\rm{J}}$, whereas the numerical results presented in the figures are obtained taking into account the exact solution to equation (4). In the fluxonium limit ${E}_{\rm{L}}\ll {E}_{\rm{J}}$, we have ${\varphi }^{*}\approx \pi (1-{E}_{\rm{L}}/{E}_{\rm{J}})$. For small $\delta f\ne 0$ the minima shift slightly to the right and are given by ${\varphi }_{\rm{R}}\approx {\varphi }^{*}-2\pi (1-\pi {E}_{\rm{L}}/{E}_{\rm{J}})\delta f$ and ${\varphi }_{\rm{L}}\approx -{\varphi }^{*}-2\pi (1-\pi {E}_{\rm{L}}/{E}_{\rm{J}})\delta f$.

In the vicinity of a local minimum $\bar{\varphi }={\varphi }_{\rm{L}/\rm{R}}$ of the potential shown in figure 1(b), the potential energy can be approximated by a harmonic oscillator potential, $U(\varphi )\approx U(\bar{\varphi })+{\displaystyle \frac{1}{2}}{E}_{\rm{J}}{(\varphi -\bar{\varphi })}^{2}$ with the lowest-lying states corresponding to the ground state wave-functions in each of the potential minima. The spread $\delta \varphi$ of the wavefunction is given by $\delta \varphi \simeq \sqrt{{E}_{\rm{C}}{/E}_{\rm{J}}}$. As both the fluxonium qubit and the Delft qubit are in the semiclassical limit, ${E}_{\rm{C}}\ll {E}_{\rm{J}}$, we neglect the finite extent of the wave functions in the following and assume that they are well localized at the single value ${\varphi }_{\rm{L}/\rm{R}}$ of the phase variable. With that our results become independent of E_C.

In the following, we want to study the effect of a small (effective) temperature gradient on the quantum information encoded in the fluxonium. To this end, we need to calculate the resulting heat current flowing through the junctions interrupting the superconducting loop of the fluxonium qubit. We are interested in the effect of accidental (effective) temperature differences; this means that we need to calculate heat currents in response to small temperature differences $\delta T$, that is ${k}_{\rm{B}}\delta T/\Delta \sim {10}^{-2}$. We hence evaluate the heat current in the linear response regime. In this section, we briefly recall the results of [28, 29] for heat currents in a single Josephson junction in the linear-response regime for small temperature gradients before we continue by investigating the effects of the heat current on the fluxonium qubit in the next sections.

We consider two superconducting reservoirs with gaps ${\Delta }_{1}=\Delta ({T}_{1})$ and ${\Delta }_{2}=\Delta ({T}_{2})$ and phase difference ϕ, interrupted by a junction with transmission probability D. Taking the number of channels contributing to transport to be N, the normal-state resistance of the junction is given by $R=h/(2{e}^{2}{ND})$. We now assume that the two reservoirs are kept at temperatures ${T}_{1}=T$ and ${T}_{2}=T+\delta T$. Being interested in small temperature gradients, $\delta T/T\ll 1$, only, we have⁷ ${\Delta }_{1}={\Delta }_{2}=\Delta$. In this regime, the heat current between the two superconducting electrodes is given by the linear response result [28, 29]

$$\begin{eqnarray}&&\dot{Q}(\phi ,T,\delta T)=-\kappa (\phi ,T)\delta T\end{eqnarray} \tag{ 5 }$$

with the thermal conductance

$$\begin{eqnarray}&&\kappa ={\displaystyle \frac{1}{2{e}^{2}{{Rk}}_{\rm{B}}{T}^{2}}}{{\displaystyle \int }}_{\Delta }^{\infty }\rm{d}\omega \;{\displaystyle \frac{{\omega }^{2}}{{\mathrm{cosh}}^{2}(\omega {/2k}_{\rm{B}}T)}}{\displaystyle \frac{{\omega }^{2}-{\Delta }^{2}}{{({\omega }^{2}-{\omega }_{\rm{c}}^{2})}^{2}}}\left[({\omega }^{2}-{\Delta }^{2}\mathrm{cos}\phi )-D{\Delta }^{2}{\mathrm{sin}}^{2}{\displaystyle \frac{\phi }{2}}\right].\end{eqnarray} \tag{ 6 }$$

Here, we have introduced the energy ${\omega }_{\rm{c}}=\Delta {[1-D{\mathrm{sin}}^{2}(\phi /2)]}^{1/2}$ of the weakly bound Andreev state emerging in the junction [37]. Most importantly, the heat current depends on the phase difference ϕ across the junction. We are interested in the tunnelling regime with $D\ll 1$, such that the last term in (6) can be neglected with the result

$$\begin{eqnarray}&&\kappa ={\displaystyle \frac{1}{2{e}^{2}{{Rk}}_{\rm{B}}{T}^{2}}}{{\displaystyle \int }}_{\Delta }^{\infty }\rm{d}\omega \;{\displaystyle \frac{{\omega }^{2}}{{\mathrm{cosh}}^{2}(\omega {/2k}_{\rm{B}}T)}}{\displaystyle \frac{{\omega }^{2}-{\Delta }^{2}}{{({\omega }^{2}-{\omega }_{\rm{c}}^{2})}^{2}}}({\omega }^{2}-{\Delta }^{2}\mathrm{cos}\phi ).\end{eqnarray} \tag{ 7 }$$

In the following, it is important to make the phase dependence of equation (7) explicit. As outlined in appendix A, we can bring the expression in the form

$$\begin{eqnarray}&&\kappa ={\kappa }_{0}-{\kappa }_{1}{\mathrm{sin}}^{2}{\displaystyle \frac{\phi }{2}}\mathrm{ln}\left({\mathrm{sin}}^{2}{\displaystyle \frac{\phi }{2}}\right)+{\kappa }_{2}{\mathrm{sin}}^{2}{\displaystyle \frac{\phi }{2}}\end{eqnarray} \tag{ 8 }$$

with the coefficients

$$\begin{eqnarray}\begin{array}{ccc}{\kappa }_{0} & = & {\displaystyle \frac{1}{2{e}^{2}{{Rk}}_{\rm{B}}{T}^{2}}}{{\displaystyle \int }}_{\Delta }^{\infty }\rm{d}\omega \;{\displaystyle \frac{{\omega }^{2}}{{\mathrm{cosh}}^{2}(\omega /2{k}_{\rm{B}}T)}},\qquad {\kappa }_{1}={\displaystyle \frac{{\Delta }^{3}}{2{e}^{2}{{Rk}}_{\rm{B}}{T}^{2}{\mathrm{cosh}}^{2}(\Delta /2{k}_{\rm{B}}T)}},\\ {\kappa }_{2} & = & {\kappa }_{1}| \mathrm{ln}D| +{\displaystyle \frac{{\Delta }^{2}}{{e}^{2}{{Rk}}_{\rm{B}}{T}^{2}}}{{\displaystyle \int }}_{\sqrt{3}\Delta }^{\infty }{\displaystyle \frac{{\omega }^{2}}{{\mathrm{cosh}}^{2}(\omega /2{k}_{\rm{B}}T)({\omega }^{2}-{\Delta }^{2})}}\\ & & +{\displaystyle \frac{{\Delta }^{2}}{{e}^{2}{{Rk}}_{\rm{B}}{T}^{2}}}{{\displaystyle \int }}_{\Delta }^{\sqrt{3}\Delta }{\displaystyle \frac{\omega }{{\omega }^{2}-{\Delta }^{2}}}\left({\displaystyle \frac{\omega }{{\mathrm{cosh}}^{2}(\omega /2{k}_{\rm{B}}T)}}-{\displaystyle \frac{\Delta }{{\mathrm{cosh}}^{2}(\Delta /2{k}_{\rm{B}}T)}}\right),\end{array}\end{eqnarray} \tag{ 9 }$$

which are independent of the phase difference. Note that, due to the specific derivation that we chose for these parameters (see appendix A), differences in the functional form of these coefficients occur with respect to the ones given in [28, 29]. For small temperatures compared to the critical one $T\ll {T}_{\rm{c}}$, we observe that ${\kappa }_{\mathrm{0,1,2}}\propto {\rm{e}}^{-\Delta {/k}_{\rm{B}}T}$, i.e., they are exponentially suppressed at low temperatures.

We now proceed to investigate how the phase sensitivity of the heat current through a Josephson junction manifests itself in the fluxonium. The situation we have in mind is some stationary nonequilibrium quasiparticle distribution on the fluxonium qubit. As mentioned above, we model this distribution by an effective temperature, in order to keep the discussion simple⁸ . In particular, we treat the situation where one side of the black-sheep Josephson junction is at an elevated temperature ${T}_{2}=T+\delta T$ with respect to the other which is at temperature ${T}_{1}=T$, see figure 1(a). As a result quasiparticles tunnel from the 'hot' to the cold reservoir. As seen in the last section, the resulting heat current depends on the superconducting phase differences and thus on the state of the qubit. As we have assumed the M Josephson junctions of the array emulating the superinductance to be equivalent, the temperature gradient is distributed among the M elements with the temperature difference $\delta {T}_{M}=\delta T/M$ on a single junction. The heat current $\dot{Q}$ flowing into the cold reservoir, which is held at temperature ${T}_{1}=T$ (or equivalently in the heat current flowing out of the hot reservoir, kept at temperature⁹ ${T}_{2}=T+\delta T$) is given by the sum of two terms

$$\begin{eqnarray}&&\dot{Q}=\dot{Q}(-{\varphi }_{\rm{b}},T,\delta T)+\dot{Q}({\varphi }_{M},T,\delta T/M),\end{eqnarray} \tag{ 10 }$$

where the first term is the contribution of the black sheep and the second term is due to the large Josephson junction connecting the array of the superinductance to the cold reservoir.

The decoherence of the fluxonium qubit is triggered by the difference $\delta \dot{Q}=| {\dot{Q}}_{\rm{R}}-{\dot{Q}}_{\rm{L}}|$ of the heat currents flowing in the device when the qubit is in the state R/L with $\varphi ={\varphi }_{\rm{R}/\rm{L}}$. Following our previous work, [24], we introduce the sensitivity

$$\begin{eqnarray}&&s={\displaystyle \frac{\delta \dot{Q}}{| {\dot{Q}}_{\rm{R}}+{\dot{Q}}_{\rm{L}}| }}\end{eqnarray} \tag{ 11 }$$

as a measure of correlation between the heat current and the qubit state. As a large sensitivity corresponds to a large difference of heat currents for different qubit states, we expect that this in turn leads to fast dephasing of the qubit; an expectation which we confirm below. However, in a first step, we want to calculate the sensitivity for the fluxonium.

4.1. Effect of the number of junctions implementing the superinductance

4.2. Sensitivity of the heat current to the state of the fluxonium qubit

Having found that the sensitivity is independent of the specific realization of the superinductance, we here present an analytical expression for the sensitivity s in the regime $\delta f\ll 1$ and ${E}_{\rm{L}}\ll {E}_{\rm{J}}$, which is relevant for the fluxonium qubit. We have seen in section 2.2 that ${\varphi }_{\rm{R}}={\varphi }^{*}-2\pi (1-\pi {E}_{\rm{L}}{/E}_{\rm{J}})\delta f$ and ${\varphi }_{\rm{L}}=-{\varphi }^{*}-2\pi (1-\pi {E}_{\rm{L}}{/E}_{\rm{J}})\delta f$ with ${\varphi }^{*}=\pi (1-{E}_{\rm{L}}{/E}_{\rm{J}})$. Evaluating s of equation (11) to first order in $\delta f$ and leading order in ${E}_{\rm{L}}{/E}_{\rm{J}}$ yields the final result

$$\begin{eqnarray}&&s\approx {\displaystyle \frac{2{\pi }^{2}{\kappa }_{1}}{{\kappa }_{0}}}\left|\mathrm{ln}({E}_{\rm{L}}{/E}_{\rm{J}})\right|{\displaystyle \frac{{E}_{\rm{L}}^{2}}{{E}_{\rm{J}}^{2}}}\delta f.\end{eqnarray} \tag{ 12 }$$

This expression shows that the sensitivity depends quadratically on the ratio ${E}_{\rm{L}}{/E}_{\rm{J}}$ with logarithmic corrections. Thus, the heat current in the fluxonium qubit is less sensitive to the qubit state than it is the case for the Delft qubit, due to the smaller ratio ${E}_{\rm{L}}{/E}_{\rm{J}}$. This comparison can be easily performed, since a description of both qubits is possible following equation (3), with the effective qubit parameters given in table 1. In what follows, we show that the fluxonium qubit therefore enjoys an increased protection with respect to quasiparticle processes as compared to the Delft qubit. In figure 2, we have plotted the approximate result of the sensitivity of equation (12) exhibiting a quadratic dependence for small ${E}_{\rm{L}}{/E}_{\rm{J}}$; A comparison to the exact expression (11) indicates that the approximation is valid for ${E}_{\rm{L}}{/E}_{\rm{J}}\lesssim 0.05$.

Zoom In Zoom Out Reset image size

In a recent work [24], we have demonstrated that the sensitivity of the heat current in a flux qubit (in the Delft qubit design) leads to a dephasing of the qubit. The reason for this dephasing is the fact that for non-zero sensitivities the tunnelling probabilities of quasiparticles depends on the state of the qubit; hence quasiparticles which tunnel through the Josephson junction, can dephase the qubit. Importantly, the heat current incorporates both the phase-dependent quasiparticle transmission probabilities through the junction as well as the quasiparticle distribution functions. This forms the basis for our argument that the heat current captures the relevant properties leading to qubit dephasing due to quasiparticle tunneling. We therefore propose to investigate this transport property, in order to access quasiparticle dephasing on a phenomenological level. We gain additional confidence in our results for the qubit dephasing obtained from this phenomenological approach by noticing that in the regime of large temperatures (where the cut-off due to the Andreev bound state becomes unimportant) our approach reproduces the perturbative results from the microscopic model [14].

Following [24], see also appendix B, we can derive the expression

$$\begin{eqnarray}&&{\tau }_{\phi }^{-1}={\displaystyle \frac{2{\Delta }^{4}{\pi }^{4}}{{e}^{2}R}}{\left({\displaystyle \frac{{E}_{\rm{L}}\delta f}{{E}_{\rm{J}}}}\right)}^{2}{{\displaystyle \int }}_{\Delta }^{\infty }\rm{d}\omega \;{\displaystyle \frac{{\omega }^{2}-{\Delta }^{2}}{{({\omega }^{2}-{\omega }_{\rm{c}}^{*}{}^{2})}^{2}({\omega }^{2}+{\Delta }^{2}){\mathrm{cosh}}^{2}(\omega {/2k}_{\rm{B}}T)}}\end{eqnarray} \tag{ 13 }$$

for the inverse dephasing time, valid to lowest order ${E}_{\rm{L}}{/E}_{\rm{J}}$ and $\delta f$; here, ${\omega }_{\rm{c}}^{*}=\Delta {(1-D{\mathrm{sin}}^{2}{\displaystyle \frac{{\varphi }^{*}}{2}})}^{1/2}$ is the bound state energy at the phase difference ${\varphi }^{*}$ at $\delta f=0$. Note that at the sweet spot the heat currents are equivalent for the two qubit states which results in the vanishing of the dephasing rate, up to exponentially small corrections in ${E}_{\rm{J}}{/E}_{\rm{C}}$, which are not considered in this paper.

Owing to our phenomenological approach, the dephasing time can be directly brought into contact with the sensitivity of the heat currents flowing in the device to the qubit state. In order to estimate this link, we consider the low temperature regime, $T,\delta T\ll \Delta {/k}_{\rm{B}}$, and write down the product between dephasing time and the difference in heat currents in the two qubit states, $\delta \dot{Q}$. This function gives us an idea about the energy which is transferred by the difference of heat currents in the two qubit states in the time, which the qubit needs to dephase. It turns out that to lowest order in ${E}_{\rm{L}}{/E}_{\rm{J}}$ and $\delta f$, we obtain the simple result

$$\begin{eqnarray}&&{\tau }_{\phi }\delta \dot{Q}\approx {\displaystyle \frac{{\Delta }^{2}\delta T}{4{k}_{\rm{B}}{T}^{2}\delta f}}.\end{eqnarray} \tag{ 14 }$$

The full expressions are presented in appendix C. Equation (3) shows that the value of ${\tau }_{\phi }\delta \dot{Q}$ only depends on the detuning $\delta f$ from the sweet spot and is otherwise independent of any of the qubit parameters. The additional parameters $\Delta ,T,\delta T$ entering the expression describe the heat current due to the flow of the quasiparticles.

Starting from equation (14), we only need to find $\delta \dot{Q}$ in order to obtain the dephasing time ${\tau }_{\phi }$ that limits the performance of the superconducting qubit. In the limit $\delta f$ small, the two states $| {\psi }_{\rm{L}}\rangle$ and $| {\psi }_{\rm{R}}\rangle$ have a similar heat current with ${\dot{Q}}_{\rm{L}}+{\dot{Q}}_{\rm{R}}=2\dot{Q}(\varphi ={\varphi }^{*})$, such that we obtain the expression $\delta \dot{Q}=2s| \dot{Q}(\varphi ={\varphi }^{*})|$. Additionally, in the limit ${E}_{\rm{L}}\ll {E}_{\rm{J}}$, ${\varphi }^{*}$ is close to π and thus $\dot{Q}(\varphi ={\varphi }^{*})=-{\kappa }_{0}\delta T$ which leads to the final result

$$\begin{eqnarray}&&{\tau }_{\phi }\approx {\displaystyle \frac{{\Delta }^{2}}{8{\kappa }_{0}{k}_{\rm{B}}{T}^{2}\delta f}}{\displaystyle \frac{1}{s}}.\end{eqnarray} \tag{ 15 }$$

We see that as expected the dephasing time is inversely proportional to the sensitivity of the heat current that is flowing through the structure to the state of the qubit. As we have seen before, s is proportional to ${({E}_{\rm{L}}/{E}_{\rm{J}})}^{2}$ (up to logarithmic corrections); consequently the fluxonium qubit has the benefit of a larger dephasing time than the Delft qubit. Different from the dephasing caused by charge and flux noise, the inverse dephasing time due to the transport of quasiparticles is not exponentially suppressed in the fluxonium regime ${E}_{\rm{L}}\ll {E}_{\rm{J}}$. However, due to the fact that ${\kappa }_{0}\propto \mathrm{exp}(-\Delta {/k}_{\rm{B}}T)$, an exponential improvement of the dephasing time limited by quasiparticle tunneling can be reached by lowering the temperature. This is expected since it corresponds to decreasing the quasiparticle occupation.

In order to make a quantitive statement relevant for applications, it is useful to introduce the dimensionless number $\varepsilon {\tau }_{\phi }/{\hslash }$ with the level splitting of the qubit, rather than looking at the dephasing time only. In particular, since the inverse level splitting, ${\varepsilon }^{-1}$, yields a measure of the time of a qubit operation, the parameter $\varepsilon {\tau }_{\phi }/{\hslash }$ indicates the typical number of single qubit gates which can be performed before coherence is lost. We numerically determine the level splitting from the difference of the energies of the ground and excited state from the Hamiltonian given in equation (3). In figure 3, we show a plot of the function $\varepsilon {\tau }_{\phi }\;\delta f/{\hslash }$. Since with small changes in f, the level splitting increases linearly with $\delta f$, while ${\tau }_{\phi }$ is proportional to $\delta {f}^{-2}$, we choose to multiply the parameter of interest with $\delta f$ in order to get a function which is essentially independent of $\delta f$. Figure 3 shows a logarithmic plot of $\varepsilon {\tau }_{\phi }\;\delta f/{\hslash }$ for small values of the temperature, $T\ll \Delta {/k}_{\rm{B}}$. The ratio between the two time-scales, ${\tau }_{\phi }$ and ${\hslash }/\varepsilon$, occurs to be approximately exponentially suppressed with increasing temperature. This can be seen by using equations (5) and (12), where we see that ${\tau }_{\phi }\propto {T}^{-2}{\kappa }_{1}{(T)}^{-1}$. Moreover, from the definitions of ${\kappa }_{1}(T)$ given in equation (9), for small temperatures we have ${T}^{2}{\kappa }_{1}(T)\propto \mathrm{exp}(-\Delta {/k}_{\rm{B}}T)$.

Zoom In Zoom Out Reset image size

Finally, figure 3 shows that for small values of the temperature it is possible to perform a great number of operations on the qubit state before they become unreliable due to dephasing of the two-level system. For example for $\delta f={10}^{-2}$ and ${k}_{\rm{B}}T\lt 0.15\Delta$, we have that $\varepsilon {\tau }_{\phi }/{\hslash }\gt {10}^{4}$.

In this paper, we have studied the impact of the phase sensitivity of the quasiparticle transport on the coherence properties of the fluxonium qubit. Using a phenomenological approach, based on the study of heat currents carried by quasiparticles and the associated heat conductance, we have shown that the dephasing time is inversely proportional to the sensitivity, a quantity describing to which extent possible heat currents flowing in the device depend on the state of the qubit. We have shown that the sensitivity of the heat current to the qubit state depends quadratically on the ratio ${E}_{\rm{L}}{/E}_{\rm{J}}$ of the characteristic energies of the fluxonium qubit but not on the number of junctions with which the superinductance, a relevant ingredient of the fluxonium qubit, is realized. The independence of the number of array-junctions can qualitatively be traced back to the fact that, at small temperature gradients, the heat current in the arms of the loop constituting the fluxonium qubit is mainly given by the heat current in the so-called black sheep junction, which does not depend on M. The fact that the sensitivity can be reduced by lowering the ratio ${E}_{\rm{L}}{/E}_{\rm{J}}$ has the important result that the fluxonium qubit is less affected by dephasing due to quasiparticle tunnelling through the Josephson junction when compared to the Delft qubit design. We furthermore find that the dephasing mechanism is exponentially suppressed with temperature due to its origin in quasiparticle tunneling. However, we have shown that at moderately low temperatures the resulting dephasing time is demonstrated to be large enough to easily allow an excess of 10⁴ operations before the qubit dephases.

We thank G Catelani and R Fazio for fruitful discussions. FH acknowledges financial support from the Alexander von Humboldt foundation. JS acknowledges financial support from the Knut and Alice Wallenberg foundation through the Wallenberg Academy Fellows program and from the Swedish VR.

We aim to get an understanding of the phase dependence of the linear response coefficient κ of the heat current in the weak tunnel coupling regime, $D\ll 1$ as given in equation (7)

$$\begin{eqnarray}&&\kappa ={\displaystyle \frac{1}{2{e}^{2}{{Rk}}_{\rm{B}}{T}^{2}}}{{\displaystyle \int }}_{\Delta }^{\infty }\rm{d}\omega \;{\displaystyle \frac{{\omega }^{2}}{{\mathrm{cosh}}^{2}(\omega {/2k}_{\rm{B}}T)}}{\displaystyle \frac{{\omega }^{2}-{\Delta }^{2}}{{({\omega }^{2}-{\omega }_{\rm{c}}^{2})}^{2}}}({\omega }^{2}-{\Delta }^{2}\mathrm{cos}\phi ),\end{eqnarray} \tag{ A.1 }$$

by rewriting it in terms of the expression given in equation (9). In this appendix, we outline the procedure, which we apply to obtain the coefficients given in equation (9).

As a starting point, it is useful to evaluate the logarithmic divergence occuring in the linear response coefficient shown in equation (A.1). Since the divergence stems from values of ω in the vicinity of the superconducting gap Δ, this logarithmic divergence can conveniently be extracted by setting $\omega \approx \Delta$ in all contributions of equation (A.1) with a smooth dependence on ω in the vicinity of Δ. This consideration leads us to determine the integral

$$\begin{eqnarray}\begin{array}{ccc}{\kappa }_{\mathrm{div}} & \equiv & {\displaystyle \frac{{\Delta }^{3}{\mathrm{sin}}^{2}{\displaystyle \frac{\phi }{2}}}{{e}^{2}{{Rk}}_{\rm{B}}{T}^{2}{\mathrm{cosh}}^{2}\left(\Delta {/2k}_{\rm{B}}T\right)}}{{\displaystyle \int }}_{\Delta }^{\sqrt{3}\Delta }\rm{d}\omega \;{\displaystyle \frac{\omega }{{\omega }^{2}-{\Delta }^{2}\left(1-D{\mathrm{sin}}^{2}{\displaystyle \frac{\phi }{2}}\right)}}\\ & = & -{\displaystyle \frac{{\Delta }^{3}{\mathrm{sin}}^{2}{\displaystyle \frac{\phi }{2}}}{2{e}^{2}{{Rk}}_{\rm{B}}{T}^{2}{\mathrm{cosh}}^{2}\left(\Delta {/2k}_{\rm{B}}T\right)}}\mathrm{ln}\left(D{\mathrm{sin}}^{2}{\displaystyle \frac{\phi }{2}}\right).\end{array}\end{eqnarray} \tag{ A.2 }$$

The latter is the only divergent contribution in equation (A.1). With ${R}^{-1}\propto D$, we find that the leading contribution to this term is of the order $D\mathrm{ln}D$. The remaining part of the integral can safely be expanded for small $D{\mathrm{sin}}^{2}{\displaystyle \frac{\phi }{2}}$. We find the remaining part by simply substracting the logarithmic divergence from the full coefficient ${\kappa }_{\mathrm{rem}}\equiv \kappa -{\kappa }_{\mathrm{div}}$, leading to

$$\begin{eqnarray}\begin{array}{ccc}{\kappa }_{\mathrm{rem}} & = & {\displaystyle \frac{1}{2{e}^{2}{{Rk}}_{\rm{B}}{T}^{2}}}{{\displaystyle \int }}_{\Delta }^{\infty }\rm{d}\omega \;{\displaystyle \frac{{\omega }^{2}}{{\mathrm{cosh}}^{2}(\omega {/2k}_{\rm{B}}T)}}{\displaystyle \frac{{\omega }^{2}-{\Delta }^{2}}{{({\omega }^{2}-{\omega }_{\rm{c}}^{2})}^{2}}}({\omega }^{2}-{\Delta }^{2}\mathrm{cos}\phi )\\ & & -{\displaystyle \frac{{\Delta }^{3}{\mathrm{sin}}^{2}{\displaystyle \frac{\phi }{2}}}{{e}^{2}{{Rk}}_{\rm{B}}{T}^{2}{\mathrm{cosh}}^{2}\left(\Delta {/2k}_{\rm{B}}T\right)}}{{\displaystyle \int }}_{\Delta }^{\sqrt{3}\Delta }\rm{d}\omega \;{\displaystyle \frac{\omega }{{\omega }^{2}-{\Delta }^{2}\left(1-D{\mathrm{sin}}^{2}{\displaystyle \frac{\phi }{2}}\right)}}.\end{array}\end{eqnarray} \tag{ A.3 }$$

Indeed, the expansion of equation (A.3) to linear order in D shows no divergent behavior anymore. We find

$$\begin{eqnarray*}\begin{array}{ccc}{\kappa }_{\mathrm{rem}} & = & {\displaystyle \frac{1}{2{e}^{2}{{Rk}}_{\rm{B}}T}}{{\displaystyle \int }}_{\Delta }^{\infty }\rm{d}\omega {\displaystyle \frac{{\omega }^{2}}{{\mathrm{cosh}}^{2}\left(\omega {/2k}_{\rm{B}}T\right)}}\\ & & +{\displaystyle \frac{{\Delta }^{2}}{2{e}^{2}{{Rk}}_{\rm{B}}T}}{{\displaystyle \int }}_{\sqrt{3}\Delta }^{\infty }\rm{d}\omega {\displaystyle \frac{{\omega }^{2}}{{\mathrm{cosh}}^{2}\left(\omega {/2k}_{\rm{B}}T\right)}}{\displaystyle \frac{{\mathrm{sin}}^{2}{\displaystyle \frac{\phi }{2}}}{{\omega }^{2}-{\Delta }^{2}}}\\ & & +{\displaystyle \frac{{\Delta }^{2}}{2{e}^{2}{{Rk}}_{\rm{B}}T}}{{\displaystyle \int }}_{\Delta }^{\sqrt{3}\Delta }\rm{d}\omega {\displaystyle \frac{\omega }{{\omega }^{2}-{\Delta }^{2}}}\left[{\displaystyle \frac{\omega }{{\mathrm{cosh}}^{2}\left(\omega {/2k}_{\rm{B}}T\right)}}-{\displaystyle \frac{\Delta }{{\mathrm{cosh}}^{2}\left(\Delta {/2k}_{\rm{B}}T\right)}}\right]{\mathrm{sin}}^{2}{\displaystyle \frac{\phi }{2}}.\end{array}\end{eqnarray*}$$

This expression, together with the contribution from the logarithmic divergence ${\kappa }_{\mathrm{div}}$, yields the results presented in equations (8) and (9) in the main text.

In order to investigate the impact of possible temperature gradients on the dephasing of the fluxonium qubit, we follow the lines of [24]. We therefore use a model Hamiltonian, ${H}_{\mathrm{mod}}={H}_{0}+{H}_{\rm{I}}$, where H₀ describes the qubit as a two-level system with states $| {\psi }_{\rm{L}}\rangle$ and $| {\psi }_{\rm{R}}\rangle$ and normal quasi-particle reservoirs at different temperatures. The coupling between them is given by H_I. More specifically, we have

$$\begin{eqnarray}&&{H}_{\mathrm{mod}}=-{\displaystyle \frac{\varepsilon }{2}}{\tau }^{3}+{\displaystyle \sum }_{l=1,2}{\displaystyle \sum }_{k,\sigma }({\varepsilon }_{l,k}-{\mu }_{l}){c}_{l,k\sigma }^{\dagger }{c}_{l,k\sigma }^{\phantom{\dagger }}+{\displaystyle \sum }_{k,q,\sigma }\left[({V}_{0}{\tau }^{0}+{V}_{3}{\tau }^{3}){c}_{1,k\sigma }^{\dagger }{c}_{2,q\sigma }^{\phantom{\dagger }}+\rm{h}.\rm{c}.\right]\end{eqnarray} \tag{ B.1 }$$

The matrices ${\tau }^{j},j=0,3$ are Pauli matrices in the qubit space. The level splitting between the qubit states is given by ; coupling between them is supposed to be weak and is neglected here. In the reservoirs, l $(l=1,2)$, the creation (annihilation) operators of particles with momentum k and spin σ are given by ${c}_{l,k\sigma }^{\dagger }({c}_{l,k\sigma })$. It is the state-dependent coupling between qubit and reservoirs occurring in the interaction part of the Hamiltonian, ${V}_{0/3}=({V}_{\rm{R}}\pm {V}_{\rm{L}})/2$, together with the density of states of the reservoirs, which takes account for the phase dependence (and hence for the dependence on the fluxonium states) of the heat current due to the superconductors.

We recover the state-dependent qubit-reservoir coupling and the density of states by comparing the heat current obtained from the model Hamiltonian in the linear response regime

$$\begin{eqnarray}&&{\dot{Q}}_{\rm{R}/\rm{L}}^{\mathrm{mod}}=-{\displaystyle \frac{\pi \delta T}{4{\hslash }{k}_{\rm{B}}{T}^{2}}}{{\displaystyle \int }}_{\Delta }^{\infty }\rm{d}\omega \;{\displaystyle \frac{{\omega }^{2}}{{\mathrm{cosh}}^{2}(\omega /2{k}_{\rm{B}}T)}}{V}_{\rm{R}/\rm{L}}^{2}{n}_{1}^{\rm{R}/\rm{L}}{n}_{2}^{\rm{R}/\rm{L}},\end{eqnarray} \tag{ B.2 }$$

to the one calculated starting from equations (5) and (10) for the black sheep. Note that we here neglect the effect of the heat current through the array junctions, since no relevant dependence on the qubit state occurs there. In order to extract the local density of states of the continuum states above the gap, we use the relation [29]

$$\begin{eqnarray*}\begin{array}{ccc}{n}_{l}^{\rm{R}/\rm{L}}(\omega ) & = & {n}_{l}^{0}\langle {\varphi }_{\rm{R}/\rm{L}}| {\displaystyle \frac{| \omega | {({\omega }^{2}-{\Delta }^{2})}^{1/2}}{{\omega }^{2}-{\omega }_{\rm{c}}^{2}}}| {\varphi }_{\rm{R}/\rm{L}}\rangle \\ & = & {n}_{l}^{0}{\displaystyle \frac{| \omega | {({\omega }^{2}-{\Delta }^{2})}^{1/2}}{{\omega }^{2}-{\Delta }^{2}(1-D{\mathrm{sin}}^{2}{\displaystyle \frac{{\varphi }_{\rm{R}/\rm{L}}}{2}})}}.\end{array}\end{eqnarray*}$$

Here, n⁰_l is the normal conducting density of states including spin. Using this form for the density of states we have the following expression for the tunnelling matrix elements in the tunnelling regime, $D\ll 1$,

$$\begin{eqnarray}&&{V}_{\rm{R}/\rm{L}}^{2}(\omega )={V}_{12}^{2}\left(1-{\displaystyle \frac{{\Delta }^{2}}{{\omega }^{2}}}\mathrm{cos}{\varphi }_{\rm{R}/\rm{L}}\right).\end{eqnarray} \tag{ B.3 }$$

Here, V₁₂ is the tunnelling amplitude of the junction figuring as the black sheep in the fluxonium qubit. It is linked to the normal state resistance by $R={\hslash }/(\pi {e}^{2}{n}_{1}^{0}{n}_{2}^{0}{V}_{12}^{2})$.

With the help of this simplified model, we proceed to study the dynamics of the qubit state. Starting from the density matrix of the full system consisting of the qubit coupled to reservoirs, we trace out the reservoir degrees of freedom and write down a master equation for the reduced density matrix of the qubit, $\rho (t)$. It is helpful to rewrite the interaction part of the Hamiltonian as ${H}_{\rm{I}}={P}_{\rm{R}}{B}_{\rm{R}}+{P}_{\rm{L}}{B}_{\rm{L}}$ with the projectors on the qubit states, $\alpha =\rm{R},\rm{L}$,

$$\begin{eqnarray}&&{P}_{\alpha }=| {\psi }_{\alpha }\rangle \langle {\psi }_{\alpha }| ,\quad \quad {B}_{\alpha }={V}_{\alpha }{\displaystyle \sum }_{k,q,\sigma }{c}_{1,k\sigma }^{\dagger }{c}_{2,q\sigma }+\rm{h}.\rm{c}.\end{eqnarray} \tag{ B.4 }$$

Following a standard procedure, see for example [39], the master equation for the density matrix of the qubit then takes the form

$$\begin{eqnarray*}&&\dot{\rho }(t)=-{\displaystyle \frac{\rm{i}}{{\hslash }}}[{H}_{\rm{S}},\rho (t)]+{\displaystyle \sum }_{\alpha ,\beta =\rm{R},\rm{L}}{\gamma }_{\alpha \beta }\left({P}_{\beta }\rho (t){P}_{\alpha }-{\displaystyle \frac{1}{2}}\{{P}_{\alpha }{P}_{\beta },\rho (t)\}\right)\end{eqnarray*}$$

with the transition rates between qubit states $\alpha ,\beta =\rm{R},\rm{L}$,

$$\begin{eqnarray}&&{\gamma }_{\alpha \beta }={\displaystyle \frac{1}{2}}\int \langle \{{B}_{\alpha }(t),{B}_{\beta }(0)\}\rangle \rm{d}t.\end{eqnarray} \tag{ B.5 }$$

The relaxation behavior of the qubit becomes particularly clear when rewriting the master equation in terms of a Pauli rate equation for the pseudo-spin states of the qubit, $\boldsymbol{S}(t)=\mathrm{Tr}$ $[\rho (t)\boldsymbol{\tau }]$ $\;=\;\left[{\rho }_{\mathrm{LR}}(t)\right.$ $+{\rho }_{\mathrm{RL}}(t),$ $\rm{i}{\rho }_{\mathrm{LR}}$ $(t)-\rm{i}$ ${\rho }_{\mathrm{RL}}(t),$ ${\rho }_{\mathrm{LL}}(t)$ $-{\left.{\rho }_{\mathrm{RR}}(t)\right]}^{{\text{}}T}$. We obtain from equation (B.5),

$$\begin{eqnarray}&&\dot{\boldsymbol{S}}(t)=\boldsymbol{S}(t)\times \boldsymbol{h}-{\tau }_{\phi }^{-1}{({S}_{1}(t),{S}_{2}(t),0)}^{{\text{}}T}\end{eqnarray} \tag{ B.6 }$$

The Pauli rate equation (B.6) contains a precession of the pseudospin around a pseudo-magnetic field, $\boldsymbol{h}={(\mathrm{0,0},\varepsilon /{\hslash })}^{{\text{}}T}$, determined by the level splitting between qubit states. Most importantly, there is also a relaxation of the coherences of the reduced density matrix with the dephasing rate ${\tau }_{\phi }^{-1}$, given by

$$\begin{eqnarray}&&{\tau }_{\phi }^{-1}={\displaystyle \frac{1}{2{\hslash }}}({\gamma }_{\mathrm{RR}}-2{\gamma }_{\mathrm{RL}}+{\gamma }_{\mathrm{LL}}).\end{eqnarray} \tag{ B.7 }$$

It is found to have the explicit form

$$\begin{eqnarray}\begin{array}{ccc}{\tau }_{\phi }^{-1} & = & {\displaystyle \frac{\pi }{2{\hslash }}}{{\displaystyle \int }}_{\Delta }^{\infty }\rm{d}\omega {\left[{V}_{\rm{R}}{({n}_{1}^{\rm{R}}{n}_{2}^{\rm{R}})}^{1/2}-{V}_{\rm{L}}{({n}_{1}^{\rm{L}}{n}_{2}^{\rm{L}})}^{1/2}\right]}^{2}{\mathrm{cosh}}^{-2}(\omega {/2k}_{\rm{B}}T)\\ & = & {\displaystyle \frac{1}{2{e}^{2}R}}{{\displaystyle \int }}_{\Delta }^{\infty }\rm{d}\omega \;{\displaystyle \frac{{\omega }^{2}-{\Delta }^{2}}{{({\omega }^{2}-{\omega }_{\rm{c}}^{*}{}^{2})}^{2}{\mathrm{cosh}}^{2}(\omega {/2k}_{\rm{B}}T)}}{\left(\sqrt{{\omega }^{2}-{\Delta }^{2}\mathrm{cos}{\varphi }_{\rm{L}}}-\sqrt{{\omega }^{2}-{\Delta }^{2}\mathrm{cos}{\varphi }_{\rm{R}}}\right)}^{2},\end{array}\end{eqnarray} \tag{ B.8 }$$

where we introduced ${\omega }_{\rm{c}}^{*}=\Delta {(1-D{\mathrm{sin}}^{2}{\displaystyle \frac{{\varphi }^{*}}{2}})}^{1/2}$.

As discussed in the main text, the dephasing time can be directly brought into connection with the heat current flowing through the qubit due to a finite temperature gradient and its sensitivity to the qubit state. In order to estimate this link, we consider the linear response regime $\delta T\ll T$ and write down the product between dephasing time and the difference in heat currents in the two qubit states, $\delta \dot{Q}$,

$$\begin{eqnarray}&&{\tau }_{\phi }\delta \dot{Q}={\displaystyle \frac{{\Delta }^{2}\delta T{{\displaystyle \int }}_{\Delta }^{\infty }\rm{d}\omega \;{\omega }^{2}{\displaystyle \frac{{\omega }^{2}-{\Delta }^{2}}{{({\omega }^{2}-{\omega }_{\rm{c}}^{*}{}^{2})}^{2}}}\left|\mathrm{cos}{\varphi }_{\rm{L}}-\mathrm{cos}{\varphi }_{\rm{R}}\right|{\mathrm{cosh}}^{-2}{\displaystyle \frac{\omega }{2{k}_{\rm{B}}T}}}{2{k}_{\rm{B}}{T}^{2}{{\displaystyle \int }}_{\Delta }^{\infty }\rm{d}\omega {\displaystyle \frac{{\omega }^{2}-{\Delta }^{2}}{{({\omega }^{2}-{\omega }_{\rm{c}}^{*}{}^{2})}^{2}}}{\left(\sqrt{{\omega }^{2}-{\Delta }^{2}\mathrm{cos}{\varphi }_{\rm{L}}}-\sqrt{{\omega }^{2}-{\Delta }^{2}\mathrm{cos}{\varphi }_{\rm{R}}}\right)}^{2}{\mathrm{cosh}}^{-2}{\displaystyle \frac{\omega }{2{k}_{\rm{B}}T}}}}.\end{eqnarray} \tag{ C.1 }$$

As we noticed before, the main contribution to these integrals stems from contributions of ω close to the superconducting gap Δ. We hence introduce $\omega \approx \Delta$ in those factors which are smooth functions of ω in the vicinity of Δ

$$\begin{eqnarray}&&{\tau }_{\phi }\delta \dot{Q}={\displaystyle \frac{{\Delta }^{2}\delta T{{\displaystyle \int }}_{\Delta }^{\infty }\rm{d}\omega \;{\Delta }^{2}{\displaystyle \frac{{\omega }^{2}-{\Delta }^{2}}{{({\omega }^{2}-{\omega }_{\rm{c}}^{*}{}^{2})}^{2}}}\left|\mathrm{cos}{\varphi }_{\rm{L}}-\mathrm{cos}{\varphi }_{\rm{R}}\right|{\mathrm{cosh}}^{-2}{\displaystyle \frac{\omega }{2{k}_{\rm{B}}T}}}{2{k}_{\rm{B}}{T}^{2}{{\displaystyle \int }}_{\Delta }^{\infty }\rm{d}\omega {\displaystyle \frac{{\omega }^{2}-{\Delta }^{2}}{{({\omega }^{2}-{\omega }_{\rm{c}}^{*}{}^{2})}^{2}}}{\left(\sqrt{{\Delta }^{2}-{\Delta }^{2}\mathrm{cos}{\varphi }_{\rm{L}}}-\sqrt{{\Delta }^{2}-{\Delta }^{2}\mathrm{cos}{\varphi }_{\rm{R}}}\right)}^{2}{\mathrm{cosh}}^{-2}{\displaystyle \frac{\omega }{2{k}_{\rm{B}}T}}}}\end{eqnarray} \tag{ C.2 }$$

and observe that this leads to a cancellation of the integral terms in numerator and denominator of equation (C.2). This underlines the close connection between the heat current sensitivity to the qubit state and the occurring dephasing mechanism. Further simplifying we find

$$\begin{eqnarray}&&{\tau }_{\phi }\delta \dot{Q}={\displaystyle \frac{{\Delta }^{2}\delta T}{4{k}_{\rm{B}}{T}^{2}}}{\displaystyle \frac{\left|\mathrm{cos}{\varphi }_{\rm{L}}-\mathrm{cos}{\varphi }_{\rm{R}}\right|}{{(| \mathrm{sin}{\displaystyle \frac{{\varphi }_{\rm{L}}}{2}}| -| \mathrm{sin}{\displaystyle \frac{{\varphi }_{\rm{R}}}{2}}| )}^{2}}}.\end{eqnarray} \tag{ C.3 }$$

In lowest order in ${E}_{\rm{L}}{/E}_{\rm{J}}$ and $\delta f$, this leads to the result presented in equation (3) in the main text.