Measurement and dephasing of a flux qubit due to heat currents

The charge current through a superconducting weak link is sensitive to the phase difference of the superconductor order parameters on either side of the link. In the absence of a bias voltage, a dissipationless Josephson current flows through the link, which is proportional to the sine of the phase difference. The origin of this supercurrent can be traced back to the Andreev reflection of incoming electrons, and as such it is an interplay of quasiparticles at the interface and the superconducting condensate on both sides.

In 1965 Maki and Griffin [1] theoretically predicted that the heat current flowing through a temperature-biased Josephson tunnel junction is also a periodic function of the phase difference between the electrodes. Due to the invariance of the heat current under time reversal, it has an even parity with respect to the phase difference. The phase dependence of the heat current—carried by quasiparticles residing at energies outside of the energy gap of the superconductor—comes again from an interplay between these quasiparticles and the superconducting condensate.

This effect has recently been demonstrated experimentally [2–6]. A superconducting ring, namely a dc-SQUID (superconducting quantum interference device) with two Josephson junctions was exposed to a temperature gradient. The measurement of the resulting heat current as a function of the magnetic flux penetrating the SQUID demonstrated the sensitivity of the heat current to the phase differences across the junctions. In this way, the SQUID is operated as a heat modulator.

Heat transport through weak links in superconductors has been theoretically studied in great detail [7–9]; see also [10] for a review on interference in heat transport and thermoelectric effects in superconducting weak links. It has been found that the heat current can be modulated by the applied phase gradient [11], and recent experiments have shown that weak links in superconductors can be used to refrigerate small islands [12] and trap hot quasiparticles [13].

An altogether different application of the phase sensitivity of the supercurrent in superconducting rings is the realization of a persistent current flux qubit where the phase sensitivity of the device is used to implement qubit operations. In particular, the Delft design of the flux qubit consists of a superconducting loop interrupted by three Josephson junctions. It is further characterized by the fact that the Josephson coupling of one of the junctions is smaller by a factor $\alpha \simeq 0.75$ [14], which in actual implementations is made tuneable by replacing this third junction with a split Josephson junction [15]. Another important tuning parameter is the external flux $\Phi$ threading the loop. If the flux is close to half a superconducting flux quantum, $\Phi =h/4e$, the superconducting system emulates a particle in a (shallow) double-well potential, where the state in either well corresponds to a circulating persistent current, either flowing clockwise or counterclockwise around the loop. These two states represent the qubit states of the device.

In what follows, we will combine these two intriguing studies on the phase-sensitivity in superconducting rings. We are in particular interested in the dependence of the heat current on the state of the persistent current qubit. We therefore investigate a superconducting ring with three Josephson junctions subject to a temperature gradient. We use a microscopic description of the Josephson junctions in order to investigate the phase-dependent heat current through them. We will show that the heat current in a temperature-biased Delft qubit is indeed sensitive to the qubit state, with typical sensitivities of 4%.

Beyond this, the state-sensitive heat current has an impact on the qubit state. With the help of a master equation approach, we investigate how the temperature gradient influences the dynamics of the qubit system. We determine the rate of coherence suppression, which is shown to be given by the rate with which the difference in heat currents at the two qubit states accumulates an energy difference of approximately the gap energy. The difference in heat currents due to a thermal gradient depending on the qubit state is hence demonstrated to be a qubit-state measurement. Depending on the temperature gradient, the associated typical dephasing times range from nano- to microseconds, with the proviso that the qubit is detuned from the 'sweet spot' of half a flux quantum, $\Phi ={{\Phi }_{0}}/2$, threading the superconducting loop, to a typical 'operation point' of $\Phi =0.495\;{{\Phi }_{0}}$. This adds an additional contribution to the dephasing, which in general is attributed to non-equilibrium quasiparticles [16–18]. In the appendix, we give a self-contained and detailed derivation of the results of Maki and Griffin for the phase-sensitive heat current in a superconducting weak link.

2.1. Persistent current qubit (Delft qubit)

We investigate a persistent current qubit as it is sketched in figure 1. Such a qubit consists of a superconducting loop with three Josephson junctions, which encloses a flux $\Phi$ supplied by an external magnetic field. The three junctions, $i=\text{a},\text{b},\text{c}$, are in general characterized by different Josephson energies $E_{\text{J}}^{i}$, with $E_{\text{J}}^{i}=I_{\text{crit}}^{i}{{\Phi }_{0}}/2\pi$, where $I_{\text{crit}}^{i}$ is the critical current of the junction and ${{\Phi }_{0}}=h/2e$ the superconducting flux quantum. Following the Delft qubit design, we choose two of the junctions to be equivalent, i.e. having the same Josephson coupling energy $E_{\text{J}}^{\text{a}}=E_{\text{J}}^{\text{b}}\equiv {{E}_{\text{J}}}$, and the third junction with a smaller Josephson energy $E_{\text{J}}^{\text{c}}=\alpha {{E}_{\text{J}}}$, where we have introduced the asymmetry parameter $\alpha \leqslant 1$. The contrasting phase differences, ${{\varphi }_{i}}$, across the junctions (the arrows in figure 1 define the direction for a positive phase difference ${{\varphi }_{i}}$) are related to each other due to the fluxoid quantization around the superconducting loop containing the junctions ⁵ ,

$$\begin{eqnarray}&&{{\varphi }_{a}}-{{\varphi }_{b}}+{{\varphi }_{c}}=-2\pi f,\end{eqnarray} \tag{ 1 }$$

where we have defined $f=\Phi /{{\Phi }_{0}}$. The total Josephson energy of the ring is given by the phase-dependent expression $U={{\mathop{\sum }^{}}_{i}}E_{\text{J}}^{i}\left( 1-\cos {{\varphi }_{i}} \right)$. Combining this relation with the flux quantization condition in (1), the Josephson energy can be written as

$$\begin{eqnarray}&&U={{E}_{\text{J}}}\left[ 2+\alpha -\cos {{\varphi }_{\text{a}}}-\cos {{\varphi }_{\text{b}}}-\alpha \cos \left( 2\pi f+{{\varphi }_{\text{a}}}-{{\varphi }_{\text{b}}} \right) \right]\ .\end{eqnarray} \tag{ 2 }$$

The potential U is plotted in figure 1(b) for $\alpha =0.75$ and f = 0.495, a typical operation point of the Delft qubit. The plot shows a periodic structure of two nearby minima. These two minima, indicated by L and R, fulfill the condition ${{\varphi }_{\text{a}}}=-{{\varphi }_{\text{b}}}\equiv \varphi$ and correspond to situations in which the Josephson current in the loop has opposite signs. Due to the periodicity of the potential, all other minima are equivalent to L and R. If the magnetic flux is tuned to $f=\frac{1}{2}$, the flux point usually called the 'sweet spot', then the two minima are equal, ${{U}_{\min }}=2{{E}_{\text{J}}}\left( 1-\frac{1}{\alpha } \right)$, and they are situated at ${{\varphi }_{\text{L}/\text{R}}}=\mp \arccos \left( 1/2\alpha \right)$. Small deviations $\delta f=f-\frac{1}{2}$ from this point yield a shift of the minima by $\delta \varphi =-2\pi \;\delta f\left( 2{{\alpha }^{2}}-1 \right)\text{/}\left( 4{{\alpha }^{2}}-1 \right)$, such that ${{\varphi }_{\text{L}/\text{R}}}=\mp \arccos \left( 1/2\alpha \right)+\delta \varphi$. Consequently, the potential becomes asymmetric as indicated in figure 1(c). For values $\alpha \leqslant 1/2$, the two minima merge into a single minimum; in the following we will hence always assume $\alpha >1/2$.

Zoom In Zoom Out Reset image size

We are interested in the quantum properties of this system and therefore use a Hamiltonian description taking φ as a general coordinate. The dynamics of the system is provided by the fact that each of the junctions adds a small electrical capacitance C. In fact the conjugate momentum to φ is given by the number of Cooper pairs $N=-i\hbar \partial /\partial \varphi$, which charge the capacitances. We arrive at the Hamiltonian

$$\begin{eqnarray}&&{{H}_{\text{qubit}}}=-4{{E}_{\text{C}}}\frac{{{\partial }^{2}}}{\partial {{\varphi }^{2}}}+{{E}_{\text{J}}}\left[ 2+\alpha -2\cos \left( \varphi \right)-\alpha \cos \left( 2\pi f+2\varphi \right) \right]\end{eqnarray} \tag{ 3 }$$

where the first term takes account of the charging energy ${{E}_{\text{C}}}={{e}^{2}}/2{{C}_{\Sigma }}$, where ${{C}_{\Sigma }}$ combines the capacitive effects of the three junctions, and the second term is the potential energy U as given in (2) with $\varphi ={{\varphi }_{\text{a}}}=-{{\varphi }_{\text{b}}}$. The low-energy physics of this system can be described by the two metastable states $\left| {{\psi }_{\text{L}}} \right\rangle$ and $\left| {{\psi }_{\text{R}}} \right\rangle$, corresponding to the ground states of the local minima of the potential as shown in figure 1(c). They will serve as the two qubit states in the following. In the vicinity of the local minima, the Hamiltonian can be approximated using $U\left( \varphi \right)\approx U\left( {{\varphi }_{\text{L}/\text{R}}} \right)+{{E}_{\text{J}}}\left[ \cos \left( {{\varphi }_{\text{L}/\text{R}}} \right)+2\alpha \cos \left( 2\pi f+2{{\varphi }_{\text{L}/\text{R}}} \right) \right]{{\left( \varphi -{{\varphi }_{\text{L}/\text{R}}} \right)}^{2}}$ for $\varphi \approx {{\varphi }_{\text{L}/\text{R}}}$. The qubit states are then given by the oscillator ground states

$$\begin{eqnarray}&&\left\langle \varphi \left| {} \right.{{\psi }_{\text{L}/\text{R}}}\right\rangle ={{\left( \frac{{{\lambda }_{\text{L}/\text{R}}}}{\pi } \right)}^{1/4}}\exp \left\{ -\frac{{{\lambda }_{\text{L}/\text{R}}}{{\left( \varphi -{{\varphi }_{\text{L}/\text{R}}} \right)}^{2}}}{2} \right\}\end{eqnarray} \tag{ 4 }$$

with the inverse of the variance

$$\begin{eqnarray}\begin{matrix} {{\lambda }_{\text{L}/\text{R}}} & = & \frac{{{E}_{\text{J}}}}{2{{E}_{\text{C}}}}\left[ \cos \left( {{\varphi }_{\text{L}/\text{R}}} \right)+2\alpha \cos \left( 2\pi f+2{{\varphi }_{\text{L}/\text{R}}} \right) \right] \\ {} & {} & \approx \frac{{{E}_{\text{J}}}}{2{{E}_{\text{C}}}}\left[ \frac{4{{\alpha }^{2}}-1}{2\alpha }\mp \frac{\pi \left( 1+2{{\alpha }^{2}} \right)\Delta f}{\alpha \sqrt{4{{\alpha }^{2}}-1}} \right]. \\ &&\end{matrix}\end{eqnarray} \tag{ 5 }$$

These states are shown in figure 2 together with the qubit potential. They are coupled through possible quantum tunnelling through the potential barrier between the two minima, of which the height that depends on the values of the asymmetry parameter, α, is tuneable via the Josephson energy of junction c ⁶ . However, as soon as the flux deviates from the value $\Phi ={{\Phi }_{0}}/2$, then the qubit eigenstates occur, well-localized in the potential wells, the coupling between the two states is negligibly small, and they are hence approximately given by $\left| {{\psi }_{\text{L}}} \right\rangle$ and $\left| {{\psi }_{\text{R}}} \right\rangle$.

Zoom In Zoom Out Reset image size

2.2. Microscopic model of the Josephson junctions

In order to calculate the heat currents flowing through the SQUID, it is important to consider the microscopic model of the three junctions. The microscopic Hamiltonian for two superconducting arms, which we here choose to be l = 1, 2, connected by a tunnel contact is given by

$$\begin{eqnarray}\begin{matrix} {{H}_{\text{junction}}} & = & \underset{l=1,2}{\mathop{\mathop{\sum }^{}}}\,\underset{k,\sigma }{\mathop{\mathop{\sum }^{}}}\,{{\xi }_{l,k\sigma }}c_{l,k\sigma }^{\dagger }c_{l,k\sigma }^{{}}-\underset{l=1,2}{\mathop{\mathop{\sum }^{}}}\,\underset{k}{\mathop{\mathop{\sum }^{}}}\,\left( {{\Delta }_{l}}c_{l,k\uparrow }^{\dagger }c_{l,k\downarrow }^{\dagger }+\text{H}\text{.c}\text{.} \right) \\ {} & {} & +\underset{k,{{k}^{\prime }},\sigma }{\mathop{\mathop{\sum }^{}}}\,\left( V_{k{{k}^{\prime }}}^{12}c_{1,k\sigma }^{\dagger }c_{2,{{k}^{\prime }}\sigma }^{{}}+\text{H}\text{.c}\text{.} \right) \\ &&\end{matrix}\end{eqnarray} \tag{ 6 }$$

where ${{\xi }_{l,k\sigma }}={{\varepsilon }_{l,k\sigma }}-{{\mu }_{0}}$ is the electron energy relative to the chemical potential, which we take to be equal for all electrodes, ${{\mu }_{l}}={{\mu }_{0}},\forall l$. The Hamiltonian for the other two junctions of the SQUID is found equivalently. The creation (annihilation) operators for electrons in reservoir l, with momentum k and spin $\sigma =\uparrow ,\downarrow$ are given by $c_{l,k\sigma }^{\dagger }\left( c_{l,k\sigma }^{{}} \right)$. The two superconductors are kept at different temperatures ${{T}_{l}}$ and have a superconducting gap ${{\Delta }_{l}}$, which we here assume to be independent of k. The gap is characterized by its absolute value $|{{\Delta }_{l}}|$ and the phase ${{\phi }_{l}}$. This phase enters the heat current across the junction only in phase differences ${{\varphi }_{12}}$ across a junction, with ${{\varphi }_{12}}={{\phi }_{1}}-{{\phi }_{2}}$. For the SQUID model considered in this manuscript, we have ${{\varphi }_{12}}={{\varphi }_{\text{a}}}$, etc. The temperature dependence of the magnitude of the superconducting gap is approximately given by $\left| {{\Delta }_{l}}\left( {{T}_{l}} \right) \right|\approx {{\Delta }_{0}}\sqrt{1-{{T}_{l}}/{{T}_{\text{crit}}}}$ with ${{\Delta }_{0}}\simeq {{k}_{\text{B}}}{{T}_{\text{crit}}}$ the gap of the superconductor at zero temperature and ${{T}_{\text{crit}}}$ the critical temperature. Here and in the following, we assume that all the superconductors are built from the same material with equivalent geometries, such that they share ${{T}_{\text{crit}}}$ and ${{\Delta }_{0}}$. Tunnelling between the two superconductors 1 and 2 occurs with the tunnelling amplitude $V_{k{{k}^{\prime }}}^{12}$. The resistance of the junction connecting reservoirs 1 and 2 is related to the normal conducting density of states of the reservoirs at the Fermi level (including spin), $N_{l}^{0}$, and the tunneling amplitude; the inverse resistance is given by $R_{12}^{-1}=\pi {{e}^{2}}N_{1}^{0}N_{2}^{0}{{\left| {{V}^{12}} \right|}^{2}}/\hbar$.

2.3. Heat currents in superconducting links

We are in the following interested in the heat currents flowing through the junctions in the SQUID, when the arms between the junctions are kept at different temperatures, as sketched in figure 1(a). The ring is supposed to be large enough such that the arms between the junctions are larger than the quasiparticle coherence length, and we can therefore model the arms as quasiparticle reservoirs and treat the heat current through the junctions separately. Note that the phase differences across the junctions are related through the superconducting fluxoid quantization given in (1). The heat current in electrode l is defined as the flow of energy with respect to the electrochemical potential of electrode $l=1,2,3$,

$$\begin{eqnarray}&&{{\dot{Q}}^{l}}=\frac{\text{d}}{\text{dt}}\left\langle {{H}_{l}}\right\rangle =-\frac{i}{\hbar }\left\langle \left[ H,{{H}_{l}} \right]\right\rangle ,\end{eqnarray} \tag{ 7 }$$

where ${{H}_{l}}$ is given by the first line of (6). We are subsequently interested in the weak tunnel coupling regime [1, 7]; see the appendix for a detailed derivation of the heat current. The heat current through a junction connecting reservoir l and m due to a difference in temperature, ${{T}_{l}}\ne {{T}_{m}}$, with $l,m=1,2,3$, can then be divided into a pure quasiparticle contribution to the heat current, $\dot{Q}_{\text{qp}}^{l}\left( {{T}_{l}},{{T}_{m}} \right)$, and an interference contribution due to an interplay between quasiparticles and the Cooper pair condensate, $\dot{Q}_{\operatorname{int}}^{l}\left( {{T}_{l}},{{T}_{m}} \right)$, namely

$$\begin{eqnarray}&&{{\dot{Q}}^{l}}\left( {{T}_{l}},{{T}_{m}} \right)=\dot{Q}_{\text{qp}}^{l}\left( {{T}_{l}},{{T}_{m}} \right)-\dot{Q}_{\operatorname{int}}^{l}\left( {{T}_{l}},{{T}_{m}} \right)\cos {{\varphi }_{\operatorname{lm}}}.\end{eqnarray} \tag{ 8 }$$

We find the pure quasiparticle contribution to the heat current to be

$$\begin{eqnarray}&&\dot{Q}_{\text{qp}}^{l}\left( {{T}_{l}},{{T}_{m}} \right)=\frac{2}{{{e}^{2}}{{R}_{\operatorname{lm}}}}\mathop{\int }^{}_{|{{\Delta }_{\max }}|}^{\infty }\text{d}\omega \;{{\omega }^{3}}\frac{{{f}_{l}}\left( \omega \right)-{{f}_{m}}\left( \omega \right)}{\sqrt{{{\omega }^{2}}-{{\left| {{\Delta }_{l}} \right|}^{2}}}\sqrt{{{\omega }^{2}}-{{\left| {{\Delta }_{m}} \right|}^{2}}}}\ ,\end{eqnarray} \tag{ 9 }$$

where ${{f}_{l}}\left( \omega \right)={{\left[ 1+\exp \left( \omega /{{k}_{\text{B}}}{{T}_{l}} \right) \right]}^{-1}}$ is the Fermi function of electrode l and $\left| {{\Delta }_{\max }} \right|=\max \left\{ \left| {{\Delta }_{l}} \right|,\left| {{\Delta }_{m}} \right| \right\}$. The interference contribution to the heat current due to the interplay between quasiparticles and the Cooper pair condensate depends on the phase difference ${{\varphi }_{\operatorname{lm}}}$ of the superconducting condensates and yields

$$\begin{eqnarray}&&\dot{Q}_{\operatorname{int}}^{l}\left( {{T}_{l}},{{T}_{m}} \right)=\frac{2}{{{e}^{2}}{{R}_{\operatorname{lm}}}}\mathop{\int }^{}_{\left| {{\Delta }_{\max }} \right|}^{\infty }\text{d}\omega \;\omega \left| {{\Delta }_{l}}{{\Delta }_{m}} \right|\frac{{{f}_{l}}\left( \omega \right)-{{f}_{m}}\left( \omega \right)}{\sqrt{{{\omega }^{2}}-{{\left| {{\Delta }_{l}} \right|}^{2}}}\sqrt{{{\omega }^{2}}-{{\left| {{\Delta }_{m}} \right|}^{2}}}}.\end{eqnarray} \tag{ 10 }$$

We have ${{T}_{l}},{{T}_{m}}\lesssim \left| {{\Delta }_{0}} \right|/{{k}_{\text{B}}}$, and the square root terms are changing faster than the other factors in the integrals of (9) and (10). The magnitude of the heat currents can then be estimated as

$$\begin{eqnarray}\begin{matrix} \dot{Q}_{qp}^{l}\left( {{T}_{l}},{{T}_{m}} \right) & {} & \simeq \dot{Q}_{\operatorname{int}}^{l}\left( {{T}_{l}},{{T}_{m}} \right)\simeq {{{\dot{Q}}}_{\text{typ}}} \\ {} & {} & =\frac{{{\left| {{\Delta }_{\max }} \right|}^{2}}}{{{e}^{2}}{{R}_{\operatorname{lm}}}}K\left( \left| {{\Delta }_{\min }} \right|/\left| {{\Delta }_{\max }} \right| \right)\left[ {{e}^{-\left| {{\Delta }_{\max }} \right|/{{k}_{\text{B}}}{{T}_{l}}}}-{{e}^{-\left| {{\Delta }_{\max }} \right|/{{k}_{\text{B}}}{{T}_{m}}}} \right]\ , \\ &&\end{matrix}\end{eqnarray} \tag{ 11 }$$

with $K\left( k \right)=\mathop{\int }^{}_{0}^{\pi /2}{{\left( 1-{{k}^{2}}{{\sin }^{2}}\phi \right)}^{-1/2}}\text{d}\phi$ the complete elliptic integral of the first kind, and $\left| {{\Delta }_{\min }} \right|$ the superconducting gap at the larger temperature. The elliptic integral is a monotonously increasing function, which starts at $\pi /2$ for small arguments and has a logarithmic divergence with $K\left( 1-{{\varepsilon }^{2}} \right)\sim \ln 1/\varepsilon$ when k approaches 1. Since the contribution of the integrands of (9) and (10) have a maximum for ω being in the vicinity of the superconducting gap $\left| {{\Delta }_{\max }} \right|$, the quasiparticle and the interference contributions to the heat current are generally of the same order of magnitude.

Before discussing the sensitivities of the heat currents to the qubit state, we here want to briefly give an estimate of the order of magnitude of the heat currents for the limits of small and large temperature differences. In the case of a small temperature difference, $\delta T\equiv \left| {{T}_{l}}-{{T}_{m}} \right|\ll {{T}_{l}},{{T}_{m}}$, we obtain from (11) the typical value ${{\dot{Q}}_{\text{typ}}}\simeq {{\left| {{\Delta }_{\max }} \right|}^{3}}K\left[ 1-\delta T/2\left( {{T}_{\text{crit}}}-T \right) \right]{{e}^{-\left| {{\Delta }_{\max }} \right|/{{k}_{\text{B}}}T}}\delta T/\left( {{e}^{2}}{{R}_{\operatorname{lm}}}{{k}_{\text{B}}}{{T}^{2}} \right)$ of the heat current. Assuming furthermore that $\delta T\ll {{T}_{\text{crit}}}-T$, we obtain the estimate ${{\dot{Q}}_{\text{typ}}}\ \simeq {{\left| {{\Delta }_{\max }} \right|}^{3}}\ln \left[ \left( {{T}_{\text{crit}}}-T \right)/{{T}_{\text{cut}}} \right]{{e}^{-\left| {{\Delta }_{\max }} \right|/{{k}_{\text{B}}}T}}\delta T/\left( {{e}^{2}}{{R}_{\operatorname{lm}}}{{k}_{\text{B}}}{{T}^{2}} \right)$. The tunneling approximation gives a cutoff temperature ${{T}_{\text{cut}}}=\delta T$, which leads to a logarithmic divergence of the heat current for small temperature gradients, as already pointed out in [1]. However, as shown in [7] this is an artifact of the tunneling approximation, which fails to properly take into account a resonance in the density of states due to a weakly bound Andreev state. The resonance introduces a new cutoff at the scale ${{T}_{\text{cut}}}=D{{\Delta }_{0}}{{\sin }^{2}}\left( {{\varphi }_{\operatorname{lm}}}\text{/}2 \right)/{{k}_{B}}$, with D the transparency of the tunneling barrier.

In contrast, in the case of a large temperature difference, we have ${{T}_{\min }}\ll {{T}_{\max }}$. Since in this case ${{T}_{\min }}$ is hence also always much smaller than ${{T}_{\text{crit}}}$, the heat current only depends on ${{T}_{\max }}$ and we obtain the estimate ${{\dot{Q}}_{\text{typ}}}\simeq \Delta _{0}^{2}{{e}^{-{{\Delta }_{0}}/{{k}_{\text{B}}}{{T}_{\max }}}}K\left( \sqrt{1-{{T}_{\max }}/{{T}_{\text{crit}}}} \right)/{{e}^{2}}{{R}_{\operatorname{lm}}}$. If we additionally have ${{T}_{\max }}\lesssim {{T}_{\text{crit}}}$, then the elliptic integral is of order one and the estimate simply reads ${{\dot{Q}}_{\text{typ}}}\simeq \Delta _{0}^{2}{{e}^{-{{\Delta }_{0}}/{{k}_{\text{B}}}{{T}_{\max }}}}/{{e}^{2}}{{R}_{\operatorname{lm}}}$.

Thermal currents in a system similar to the one we study here were measured in the experiment by Giazotto and Martinez-Perez reported in [3]. If we use these same experimental values for an estimate, we have ${{\Delta }_{0}}\simeq 200\;\mu \text{eV}$ and $R\simeq 1\;\text{k}$ $\Omega$. For $T=0.1{{T}_{\text{crit}}}$ and large temperature gradient we obtain ${{\dot{Q}}_{\text{typ}}}\simeq {{10}^{-11}}\;\text{W}$, while for a small temperature gradient we obtain the estimate ${{\dot{Q}}_{\text{typ}}}\simeq \left( \Delta T/T \right){{10}^{-14}}\;\text{W}$ (assuming that the logarithm is of order one).

In the following, we want to investigate the sensitivity of the heat current to the state of the persistent current qubit realized by the three-junction SQUID introduced before. We therefore propose to study the difference between the heat currents compared to the sum of the two currents for the qubit being in the state $\left| {{\psi }_{\text{L}}} \right\rangle$ or $\left| {{\psi }_{\text{R}}} \right\rangle$, characterized by the sensitivity,

$$\begin{eqnarray}&&{{s}_{l}}=\frac{\dot{Q}_{\text{R}}^{l}-\dot{Q}_{\text{L}}^{l}}{\dot{Q}_{\text{R}}^{l}+\dot{Q}_{\text{L}}^{l}}\ \ \ \ \ \text{with}\ \ \ \ \ \dot{Q}_{\text{L}/\text{R}}^{l}=\left\langle {{\psi }_{\text{L}/\text{R}}}\left| {{{\dot{Q}}}^{l}} \right|{{\psi }_{\text{L}/\text{R}}}\right\rangle .\end{eqnarray} \tag{ 12 }$$

The expectation values are obtained from the usual integral over φ of the product of the heat currents given in (8) with the wave functions of (4). We evaluate the heat currents in each electrode due to a temperature gradient induced by ${{T}_{1}}={{T}_{2}}<{{T}_{3}}$. This yields heat currents in electrodes 1 and 2 given by the heat flow through the junction with electrode 3 only, while the heat current in electrode 3 has two contributions. To simplify the notation, we now take as a reference the heat current into electrode 1, with ${{\dot{Q}}_{\operatorname{int}}}\equiv \dot{Q}_{\operatorname{int}}^{1}\left( {{T}_{1}},{{T}_{3}} \right)$ and ${{\dot{Q}}_{\text{qp}}}\equiv\dot{Q}_{\text{qp}}^{1}\left( {{T}_{1}},{{T}_{3}} \right)$. The sensitivities then take the simple form

$$\begin{eqnarray}\begin{matrix} {{s}_{1}} & = & \frac{{{{\dot{Q}}}_{\operatorname{int}}}\left( {{C}_{\text{L}}}-{{C}_{\text{R}}} \right)}{2{{{\dot{Q}}}_{\text{qp}}}-{{{\dot{Q}}}_{\operatorname{int}}}\left( {{C}_{\text{L}}}+{{C}_{\text{R}}} \right)} \\ {{s}_{2}} & = & \frac{{{{\dot{Q}}}_{\operatorname{int}}}\left( {{D}_{\text{L}}}-{{D}_{\text{R}}} \right)}{2{{{\dot{Q}}}_{\text{qp}}}-{{{\dot{Q}}}_{\operatorname{int}}}\left( {{D}_{\text{L}}}+{{D}_{\text{R}}} \right)} \\ {{s}_{3}} & = & \frac{{{{\dot{Q}}}_{\operatorname{int}}}\left( \alpha \left( {{D}_{\text{L}}}-{{D}_{\text{R}}} \right)+\left( {{C}_{\text{L}}}-{{C}_{\text{R}}} \right) \right)}{2{{{\dot{Q}}}_{\text{qp}}}\left( 1+\alpha \right)-{{{\dot{Q}}}_{\operatorname{int}}}\left( \alpha \left( {{D}_{\text{L}}}+{{D}_{\text{R}}} \right)+\left( {{C}_{\text{L}}}+{{C}_{\text{R}}} \right) \right)}\ . \\ &&\end{matrix}\end{eqnarray} \tag{ 13 }$$

For a short notation and assuming the two qubit states to be well localized, we here define the phase-dependent factors

$$\begin{eqnarray}\begin{matrix} {{C}_{\text{L}/\text{R}}} & = & \left\langle {{\psi }_{\text{L}/\text{R}}}\left| \cos \left( {{\varphi }_{\text{a}}} \right) \right|{{\psi }_{\text{L}/\text{R}}}\right\rangle =\cos \left( {{\varphi }_{\text{L}/\text{R}}} \right){{e}^{-1/\left( 4{{\lambda }_{\text{L}/\text{R}}} \right)}}, \\ {{D}_{\text{L}/\text{R}}} & = & \left\langle {{\psi }_{\text{L}/\text{R}}}\left| \cos \left( {{\varphi }_{\text{c}}} \right) \right|{{\psi }_{\text{L}/\text{R}}}\right\rangle =\cos \left( 2{{\varphi }_{\text{L}/\text{R}}}+2\pi f \right){{e}^{-1/{{\lambda }_{\text{L}/\text{R}}}}}. \\ &&\end{matrix}\end{eqnarray} \tag{ 14 }$$

We also used the generalized Ambegaokar–Baratoff relations [19, 20] in order to relate the heat currents ${{\dot{Q}}^{\left( i \right)}}$ through the junctions $i=\text{b},\text{c}$ to each other, when ${{T}_{1}}={{T}_{2}}$. The heat currents through the different junctions are furthermore related to the heat currents ${{\dot{Q}}^{l}}$ into the different reservoirs, $l=1,2,3$, by ${{\dot{Q}}^{i=\text{b}}}\equiv {{\dot{Q}}^{1}}$, ${{\dot{Q}}^{i=\text{c}}}\equiv {{\dot{Q}}^{2}}$ and hence ${{\dot{Q}}^{3}}=-{{\dot{Q}}^{i=\text{b}}}-{{\dot{Q}}^{i=\text{c}}}$. By comparing the separate quasiparticle and interference components of these heat currents, see (9) and (10), we then find

$$\begin{eqnarray}&&\frac{\dot{Q}_{\operatorname{int}}^{\text{c}}}{\dot{Q}_{\operatorname{int}}^{\text{b}}}=\frac{\dot{Q}_{\text{qp}}^{\text{c}}}{\dot{Q}_{\text{qp}}^{\text{b}}}=\frac{I_{\text{crit}}^{\text{b}}}{I_{\text{crit}}^{\text{c}}}=\frac{{{R}_{\text{c}}}}{{{R}_{\text{b}}}}=\frac{{{R}_{23}}}{{{R}_{13}}}=\alpha .\end{eqnarray} \tag{ 15 }$$

This finally leads to the compact expressions in (13). The results for these three sensitivities as a function of the flux, in the vicinity of the sweet spot and the operation point of the Delft qubit, are shown in figure 3. The sensitivity of the heat currents to the qubit state hence yields a possible measure of the latter. The heat currents in electrodes 2 and 3 are most sensitive to the qubit state with a sensitivity of about 2% at the 'operation point', f = 0.495 [14]. The plot in figure 3 shows a dependence of the sensitivities as a function of the magnetic flux penetrating the SQUID, which is very close to a linear function. The slopes of the latter depend on the specific realization of the qubit, namely on the ratio α, on the electrode temperatures ${{T}_{l}}$ and the applied thermal gradient, as well as on the ratio of the Josephson and charging energy. This is shown in the approximate result for the heat currents for small deviations $\delta f$ from the 'sweet spot' $f=\frac{1}{2}$, ${{s}_{l}}\approx {{m}_{l}}\left( \alpha \right)\delta f$ with the respective slope ${{m}_{l}}$. The rather complex explicit analytic form of the slopes of the three sensitivities are given in appendix B, and they are shown in figure 4 as a function of the ratio α, which is tuneable in the experiment. We find that the slopes of the sensitivities do not in general need to have the same sign. The slope with the largest absolute value is the one obtained from the heat currents into reservoir 2. This is related to the fact that the heat currents into this reservoir flow uniquely through the junction with the weakest Josephson coupling, namely junction c, which consequently has the largest phase difference and is most sensitive to the qubit state. While for the working point of the Delft qubit, that is at $\alpha \approx 0.75$, the slope of ${{s}_{2}}$ already has a rather large value, and this value can be improved by lowering α. Note however that with α approaching 0.5, the two valleys of the potential get closer and the qubit states are not well defined any more. Equivalently for $\alpha >1$, the SQUID can not be used as a qubit any longer.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

After having demonstrated the sensitivity of the heat currents to the state of the qubit, the aim of this section is to study the impact of a temperature gradient—and the resulting heat current—on the coherence properties of the qubit. Our interest in this point is twofold: on the one hand we want to find out the behavior of the qubit state under measurement, and on the other we are interested in the impact of accidental temperature gradients on the dephasing of the qubit. We therefore consider the two-level system, defined by the states $\left| {{\psi }_{\text{L}}} \right\rangle$ and $\left| {{\psi }_{\text{R}}} \right\rangle$, namely the qubit states obtained from the low-energy physics of the SQUID, in contact with two heat baths, resulting in the model Hamiltonian

$$\begin{eqnarray}\begin{matrix} {{H}_{\text{toy}}} & = & -\frac{\varepsilon }{2}{{\tau }^{3}}-\frac{w}{2}{{\tau }^{1}}+\underset{l=1,3}{\mathop{\mathop{\sum }^{}}}\,\underset{k,\sigma }{\mathop{\mathop{\sum }^{}}}\,\left( {{\varepsilon }_{l,k}}-{{\mu }_{l}} \right)a_{l,k\sigma }^{\dagger }a_{l,k\sigma }^{{}} \\ {} & {} & +\underset{k,q,\sigma }{\mathop{\mathop{\sum }^{}}}\,\left[ a_{1,k\sigma }^{\dagger }a_{3,q\sigma }^{{}}\left( {{V}_{0}}{{\tau }^{0}}+{{V}_{3}}{{\tau }^{3}} \right)+\text{H}\text{.c}\text{.} \right]. \\ &&\end{matrix}\end{eqnarray} \tag{ 16 }$$

The model is depicted in figure 5. Here, the matrices ${{\tau }^{j}},j=0,1,3$ are Pauli matrices in the qubit space. The level splitting between the qubit states is given by ε, and the weak coupling between them is denoted by w. The creation (annihilation) operators of particles with momentum k and spin σ in lead l are given by $a_{l,k\sigma }^{\dagger }\left( a_{l,k\sigma }^{{}} \right)$.

Zoom In Zoom Out Reset image size

In the simplified model (16), we do not explicitly take into account the three superconducting leads with the heat currents, which depend on all three phase differences, but rather discuss a simplified microscopic model, which involves only two leads. The idea is to set the density of states and tunnelling matrix elements such as to reproduce the correct macroscopic thermal current between the reservoirs at temperature ${{T}_{1}}={{T}_{2}}$ and ${{T}_{3}}$ in the three-lead setup. We expect that such a procedure, while being inaccurate for certain microscopic details, will correctly incorporate the effects of the phase-dependent thermal currents on the qubit. In linear response, the Hamiltonian (16) leads to a heat current,

$$\begin{eqnarray}&&\dot{Q}_{\text{L}/\text{R}}^{\text{toy}}=\frac{\pi }{\hbar }\mathop{\int }^{}_{-\infty }^{\infty }\text{d}\omega \;\omega \left( {{\left| {{V}_{0}} \right|}^{2}}\pm {{\left| {{V}_{3}} \right|}^{2}} \right){{N}_{1}}\left( \omega \right){{N}_{3}}\left( \omega \right)\left[ {{f}_{1}}\left( \omega \right)-{{f}_{3}}\left( \omega \right) \right]\end{eqnarray} \tag{ 17 }$$

with ${{N}_{l}}\left( \omega \right)$ the density of states of the electrons in lead l (including spin). If we set the parameters

$$\begin{eqnarray}\begin{matrix} {{\left| {{V}_{0}} \right|}^{2}} & = & \frac{1}{2}{{\left| {{V}^{13}} \right|}^{2}}\left[ 2\alpha +2-{{C}_{\text{R}}}-{{C}_{\text{L}}}-\alpha \left( {{D}_{\text{R}}}+{{D}_{\text{L}}} \right) \right] \\ {{\left| {{V}_{3}} \right|}^{2}} & = & \frac{1}{2}{{\left| {{V}^{13}} \right|}^{2}}\left[ {{C}_{\text{R}}}-{{C}_{\text{L}}}+\alpha \left( {{D}_{\text{R}}}-{{D}_{\text{L}}} \right) \right] \\ {{N}_{l}}\left( \omega \right) & = & N_{l}^{0}\frac{\left| {{\Delta }_{l}} \right|}{\sqrt{{{\omega }^{2}}-{{\left| {{\Delta }_{l}} \right|}^{2}}}}\theta \left( {{\omega }^{2}}-{{\left| {{\Delta }_{l}} \right|}^{2}} \right), \\ &&\end{matrix}\end{eqnarray} \tag{ 18 }$$

with $\theta \left( x \right)$ the unit-step function, we achieve the goal of reproducing the correct qubit-state dependent heat current with ${{\dot{Q}}^{\text{toy}}}={{\dot{Q}}^{3}}=-{{\dot{Q}}^{1}}-{{\dot{Q}}^{2}}$; here and below, we assume the magnitude of the quasiparticle and interference parts of the heat current to be equal.

We are now interested in the dynamics of the qubit state depending on the qubit-state sensitive heat current induced by the temperature gradient. Starting from the full systemʼs density matrix, we therefore trace out the lead degrees of freedom and write down a master equation for the reduced density matrix of the qubit, $\rho \left( t \right)$. If we write the density matrix of the qubit as $\rho \left( t \right)=\frac{1}{2}\left[ +\tau \cdot \boldsymbol{S}\left( t \right) \right]$ with

$$\begin{eqnarray*}&&\boldsymbol{S}\left( t \right)=\text{tr}\left[ \rho \left( t \right)\tau \right]={{\left( {{\rho }_{\text{LR}}}\left( t \right)+{{\rho }_{\text{RL}}}\left( t \right),i\left( {{\rho }_{\text{LR}}}\left( t \right)-{{\rho }_{\text{RL}}}\left( t \right) \right),{{\rho }_{\text{LL}}}\left( t \right)-{{\rho }_{\text{RR}}}\left( t \right) \right)}^{\text{T}}},\end{eqnarray*}$$

we obtain the Pauli rate equation

$$\begin{eqnarray}&&\dot{\boldsymbol{S}}\left( t \right)=\boldsymbol{S}\left( t \right)\times \boldsymbol{h}-\gamma {{\left( {{S}_{1}}\left( t \right),{{S}_{2}}\left( t \right),0 \right)}^{\text{T}}}.\end{eqnarray} \tag{ 19 }$$

This equation contains a precession around a pseudo-magnetic field, $\boldsymbol{h}={{\left( w,0,\varepsilon \right)}^{\text{T}}}$, determined by the qubit properties, and a relaxation of the coherences of the reduced density matrix with the rate γ, while the diagonal elements, namely the occupations of the qubit states, do not decay. This is also appreciable from the solution of the master equation, which for large detuning $\varepsilon \gg w$ with respect to the weak tunneling between the qubit states, is given by

$$\begin{eqnarray}\begin{matrix} {{\rho }_{\text{LL}}}\left( t \right)\approx {{\rho }_{\text{LL}}}\left( 0 \right), & \qquad & {{\rho }_{\text{RR}}}\left( t \right)\approx {{\rho }_{\text{RR}}}\left( 0 \right), \\ {{\rho }_{\text{LR}}}\left( t \right)\approx {{\rho }_{\text{LR}}}\left( 0 \right){{e}^{-\left( \gamma +i\varepsilon \right)t}}, & \qquad & {{\rho }_{\text{RL}}}\left( t \right)\approx {{\rho }_{\text{RL}}}\left( 0 \right){{e}^{-\left( \gamma -i\varepsilon \right)t}}. \\ &&\end{matrix}\end{eqnarray} \tag{ 20 }$$

The value of the dephasing rate γ reads,

$$\begin{eqnarray}&&\gamma =\frac{4\pi {{\left| {{V}_{3}} \right|}^{2}}N_{1}^{0}N_{3}^{0}}{\hbar }\mathop{\int }^{}_{\left| {{\Delta }_{\max }} \right|}^{\infty }\text{d}\omega \ {{\omega }^{2}}\frac{\left[ 1-{{f}_{1}}\left( \omega \right) \right]{{f}_{3}}\left( \omega \right)+\left[ 1-{{f}_{3}}\left( \omega \right) \right]{{f}_{1}}\left( \omega \right)}{\sqrt{{{\omega }^{2}}-{{\left| {{\Delta }_{1}} \right|}^{2}}}\sqrt{{{\omega }^{2}}-{{\left| {{\Delta }_{3}} \right|}^{2}}}}.\end{eqnarray} \tag{ 21 }$$

Importantly, this rate equals zero if the sensitivity of the heat current to the qubit state vanishes and hence ${{\left| {{V}_{3}} \right|}^{2}}\propto {{C}_{\text{R}}}-{{C}_{\text{L}}}+\alpha \left( {{D}_{\text{R}}}-{{D}_{\text{L}}} \right)=0$. Note that this means that the temperature gradient leads to dephasing only when the qubit is tuned away from the sweet spot. Indeed, it is possible to conclude that the qubit-state sensitivity of the heat current represents a measurement process which is reflected in the time-dependent solution of the master equation given in (20).

The dephasing rate is connected to fluctuations in the electronic subsystem which drive the qubit. In equilibrium, we would expect a fluctuation–dissipation relation to hold, which relates the fluctuations to the response coefficient of the system. Naturally, this is not true in the non-equilibrium situation studied here. It is however interesting to compare the response of the system to the temperature gradient, namely the heat current depending on the qubit states, to the related dephasing rate. We therefore introduce the dimensionless ratio $r=\left| {{\Delta }_{\max }} \right|\gamma \text{/}\left| \dot{Q}_{\text{L}}^{3}-\dot{Q}_{\text{R}}^{3} \right|$. As above, we specialize to the case when ${{T}_{1}},{{T}_{3}}\lesssim \left| {{\Delta }_{\max }} \right|\text{/}{{k}_{\text{B}}}$. With a similar calculation as the one following (10), we obtain the estimate

$$\begin{eqnarray}&&r\simeq \coth \left( \frac{\left| {{\Delta }_{\max }} \right|\left| {{T}_{1}}-{{T}_{3}} \right|}{2{{k}_{\text{B}}}{{T}_{1}}{{T}_{3}}} \right)\ .\end{eqnarray} \tag{ 22 }$$

This means that r is universal with respect to microscopic details like the normal-state resistance ${{R}_{13}}$ or the phase difference ${{\varphi }_{j}}$ across the junctions, and only depends on thermodynamical quantities like the temperatures ${{T}_{1}},{{T}_{3}}$ and the gap ${{\Delta }_{0}}$. We see that for small temperature differences, $\delta T={{T}_{1}}-{{T}_{3}}\ll {{T}_{1}},{{T}_{3}}$, this ratio becomes $r\simeq {{k}_{\text{B}}}{{T}^{2}}/\left( \left| {{\Delta }_{\max }} \right|\delta T \right)$.

The dephasing time ${{\tau }_{\Phi }}={{\gamma }^{-1}}$ is in this case given by

$$\begin{eqnarray}&&{{\tau }_{\phi }}\simeq \frac{\Delta _{0}^{2}\left( 1-T/{{T}_{\text{crit}}} \right)\delta T}{{{k}_{\text{B}}}{{T}^{2}}\left| \dot{Q}_{\text{L}}^{3}-\dot{Q}_{\text{R}}^{3} \right|}\end{eqnarray} \tag{ 23 }$$

$$\begin{eqnarray}&&\simeq \frac{{{e}^{2}}{{R}_{13}}{{e}^{{{\Delta }_{\max }}/{{k}_{\text{B}}}T}}}{{{\Delta }_{\max }}\ln \left( {{T}_{\text{crit}}}/{{T}_{\text{cut}}} \right)\left[ {{C}_{\text{L}}}-{{C}_{\text{R}}}+\alpha \left( {{D}_{\text{L}}}-{{D}_{\text{R}}} \right) \right]}.\end{eqnarray} \tag{ 24 }$$

The dephasing time in units of ${{e}^{2}}{{R}_{13}}{{e}^{{{\Delta }_{\max }}/{{k}_{\text{B}}}{{T}_{1}}}}/{{\Delta }_{0}}$ is shown in figure 6 as a function of the temperature difference $\Delta T$ for different values of the minimum temperature ${{T}_{1}}$. The inset (a) of figure 6 compares the full result (green line) with the approximation of equation (23).

Zoom In Zoom Out Reset image size

In the opposite regime of large temperature bias, we have $r\simeq 1$, and thus the dephasing rate is approximately equal to $\gamma \simeq \left| \dot{Q}_{\text{L}}^{3}-\dot{Q}_{\text{R}}^{3} \right|/\left| {{\Delta }_{\max }} \right|$. In particular, we find in this regime that the dephasing time ${{\tau }_{\phi }}$ of the qubit is given by

$$\begin{eqnarray}&&{{\tau }_{\phi }}\simeq \frac{{{\Delta }_{0}}}{\left| \dot{Q}_{\text{L}}^{3}-\dot{Q}_{\text{R}}^{3} \right|}\simeq \frac{{{e}^{2}}{{R}_{13}}{{e}^{{{\Delta }_{0}}/{{k}_{\text{B}}}{{T}_{\max }}}}}{{{\Delta }_{0}}[{{C}_{\text{L}}}-{{C}_{\text{R}}}+\alpha ({{D}_{\text{L}}}-{{D}_{\text{R}}})]},\end{eqnarray} \tag{ 25 }$$

i.e. the time after which the difference in the energy transported by the heat currents in the two qubit states equals the gap ${{\Delta }_{0}}$ of the superconductor. This is confirmed well by the inset (b) of figure 6, which shows the dramatic decrease in the dephasing time with increasing temperature gradient, which in this regime is approximately given by ${{T}_{\max }}$.

Using the values we applied to estimate the heat currents in this system, we can also estimate the dephasing time. Taking ${{\Delta }_{0}}/{{e}^{2}}{{R}_{13}}\simeq 1$ THz as in [3], we have ${{\tau }_{\phi }}\approx 1$ ns for large temperature gradients, $\delta T/{{T}_{\min }}\gg 1$. For small temperature gradients, it can be shown that a temperature ${{T}_{\min }}$ of less than $0.1{{T}_{\text{crit}}}$ has to be reached in order to avoid a strong limitation of the dephasing due to the thermal current. Indeed, for $\delta T/{{T}_{\min }}\ll 1$, ${{T}_{\min }}\simeq {{T}_{\max }}\lesssim 0.1{{T}_{\text{crit}}}$, and taking the logarithm to be of order one, we have ${{\tau }_{\Phi }}\approx 1\;\mu \text{s}$. The actual dephasing times of the Delft qubit range from a few tens of nanoseconds [21] up to a microsecond [22] and thus are of the same order of magnitude. As the nominal temperatures reached for todayʼs superconducting persistent current qubits are usually smaller than $0.1{{T}_{\text{crit}}}$ [23, 24], it is unlikely that the thermal currents do constitute the dominant source of dephasing for those qubits. However, it is well known that quasiparticles in small superconducting structures badly thermalize, leading to problems in reaching the base temperature in the dilution refrigerator [13, 25], and thus the effects of the phase dependence of the thermal current on the coherence properties of the Delft qubit cannot be excluded.

We have shown that due to the phase sensitivity of the heat current that flows in the weak links of a superconducting loop, the heat current due to a temperature gradient applied to a flux qubit depends on the state of the qubit which is formed when the loop is threaded with a magnetic flux that is close to half a superconducting flux quantum. We found that the sensitivity of the heat current to the qubit state can be up to 4%, when the qubit is tuned away from the 'sweet spot' of exactly half a flux quantum threading the loop. This should allow us to identify the state of the flux qubit in experiments of the type performed in [2, 3].

Moreover, we have found that due to this difference in heat currents at different qubit states, a thermal gradient leads to a dephasing of the qubit. In particular, we found that the ratio of the dephasing rate to the difference in the heat currents is universal with respect to microscopic details, and only depends on the temperature of the reservoirs measured in units of the superconducting gap at zero temperature. For example, in the case of large temperature gradients, the dephasing time of the qubit corresponds to the time when the energy difference transported by the difference in heat currents is of the order of the superconducting gap. We have shown that the dephasing time of the flux qubit in the Delft design due to the phase-sensitive heat current can range from microseconds for small temperature differences to nanoseconds for large temperature differences, thus constituting a potential source of dephasing given the fact that the qubits are driven by microwave pulses, which may lead to an imbalance of heating between the different sections of the superconducting loop.

We thank Gianluigi Catelani, Simone Gasparinetti, Francesco Giazotto, Dmitriy Golubev and Anna Napoli for useful comments on our work. We would like to acknowledge financial support from the Excellence Initiative of the German Federal and State Governments, and from the Ministry of Innovation (SS and JS), as well as from the Alexander von Humboldt foundation (FH).

In this appendix, we derive in detail the analytic formulae for the heat current, which we use in (9) and (10) and which were previously found in [1]. The aim of this section is to describe the heat current in a superconductor–insulator–superconductor Josephson junction, biased with a temperature gradient $\delta T$ across it, while no additional voltage is applied so that we have ${{\mu }_{1}}={{\mu }_{2}}={{\mu }_{0}}$. We assume that each superconductor (the two electrodes are denoted by l = 1, 2) is a particle reservoir in equilibrium at temperature ${{T}_{l}}$ and that it is characterized by the mean-field BCS Hamiltonian

$$\begin{eqnarray}&&{{H}_{l}}=\underset{k,\sigma }{\mathop{\mathop{\sum }^{}}}\,\left( {{\varepsilon }_{l,k}}-{{\mu }_{0}} \right)c_{l,k\sigma }^{\dagger }{{c}_{l,k\sigma }}-\underset{k}{\mathop{\mathop{\sum }^{}}}\,\left( {{\Delta }_{l,k}}c_{l,k\uparrow }^{\dagger }c_{l,-k\downarrow }^{\dagger }+\Delta _{l,k}^{*}{{c}_{l,-k\downarrow }}{{c}_{l,k\uparrow }} \right)\ .\end{eqnarray} \tag{ A.1 }$$

In (A.1), $c_{l,k\sigma }^{\dagger }$ and ${{c}_{l,k\sigma }}$ are single-electron creation and annihilation operators in the momentum k and spin σ representation, and ${{\Delta }_{l,k}}$ is the superconducting energy gap of the lth electrode. Tunnelling between reservoirs is described by the tunnelling Hamiltonian

$$\begin{eqnarray}&&{{H}_{\text{T}}}={{\mathop{\sum }^{}}_{k,q,\sigma }}(V_{kq}^{12}c_{1,k\sigma }^{\dagger }{{c}_{2,q\sigma }}+V_{kq}^{12*}c_{2,q\sigma }^{\dagger }{{c}_{1,k\sigma }}),\end{eqnarray} \tag{ A.2 }$$

where k and q are the momentum quantum numbers and the tunnelling matrix element is denoted by $V_{kq}^{12}$. The total Hamiltonian is then written as ${{H}_{\text{tot}}}={{H}_{1}}+{{H}_{2}}+{{H}_{\text{T}}}$.

Assuming that the system is sufficiently isolated and that in particular phonons are frozen out at very low temperatures, the heat current in electrode 1 is carried by electrons entering or leaving it, accompanied by a change in the overall energy ${{H}_{1}}$, with respect to the electrochemical potential. According to the quantum-mechanical equation of motion, the heat current into the first electrode is

$$\begin{eqnarray}&&\frac{\text{d}{{Q}^{\left( 1 \right)}}}{\text{d}t}=\left\langle \frac{\text{d}}{\text{dt}}{{H}_{1}}\right\rangle =\frac{i}{\hbar }\left\langle \left[ {{H}_{\text{tot}}},{{H}_{1}} \right]\right\rangle .\end{eqnarray} \tag{ A.3 }$$

Of the full commutator $\left[ {{H}_{\text{tot}}},{{H}_{1}} \right]$, only the contribution $\left[ {{H}_{\text{tot}}},{{H}_{1}} \right]=\left[ {{H}_{T}},{{H}_{1}} \right]$ is non-zero. Since we are dealing with fermionic annihilation and creation operators, they must obey the anticommutation rule $\left\{ c_{l,k\sigma }^{\dagger },{{c}_{m{{k}^{\prime }}{{\sigma }^{\prime }}}} \right\}={{\delta }_{\operatorname{lm}}}{{\delta }_{k{{k}^{\prime }}}}{{\delta }_{\sigma {{\sigma }^{\prime }}}}$ from which we can easily derive

$$\begin{eqnarray}&&\left[ c_{l,k\sigma }^{\dagger }{{c}_{l,k\sigma }},c_{l,{{k}^{\prime }}{{\sigma }^{\prime }}}^{\dagger } \right]={{\delta }_{k{{k}^{\prime }}}}{{\delta }_{\sigma {{\sigma }^{\prime }}}}c_{l,k\sigma }^{\dagger },\quad \left[ c_{l,k\sigma }^{\dagger }{{c}_{l,k\sigma }},{{c}_{l,{{k}^{\prime }}{{\sigma }^{\prime }}}} \right]=-{{\delta }_{k{{k}^{\prime }}}}{{\delta }_{\sigma {{\sigma }^{\prime }}}}{{c}_{l,k\sigma }}.\end{eqnarray} \tag{ A.4 }$$

Using these commutation relations, the evaluation of $\left[ {{H}_{T}},{{H}_{1}} \right]$ yields

$$\begin{eqnarray}\begin{matrix} \left[ {{H}_{T}},{{H}_{1}} \right] & = & \underset{k,q,\sigma }{\mathop{\mathop{\sum }^{}}}\,\underset{{{k}^{\prime }},{{\sigma }^{\prime }}}{\mathop{\mathop{\sum }^{}}}\,\left\{ V_{kq}^{12}{{\xi }_{1,{{k}^{\prime }}}}\left[ c_{1,k\sigma }^{\dagger },c_{1,{{k}^{\prime }}{{\sigma }^{\prime }}}^{\dagger }{{c}_{1,{{k}^{\prime }}{{\sigma }^{\prime }}}} \right]{{c}_{2,q\sigma }}+V_{kq}^{12*}{{\xi }_{1,{{k}^{\prime }}}}c_{2,q\sigma }^{\dagger }\left[ {{c}_{1,k\sigma }},c_{1,{{k}^{\prime }}{{\sigma }^{\prime }}}^{\dagger }{{c}_{1,{{k}^{\prime }}{{\sigma }^{\prime }}}} \right] \right. \\ {} & {} & -\left. V_{kq}^{12}\Delta _{1,{{k}^{\prime }}}^{*}\left[ c_{1,k\sigma }^{\dagger },{{c}_{1,-{{k}^{\prime }}\downarrow }}{{c}_{1,{{k}^{\prime }}\uparrow }} \right]{{c}_{2,q\sigma }}-V_{kq}^{12*}{{\Delta }_{1,{{k}^{\prime }}}}c_{2,q\sigma }^{\dagger }\left[ {{c}_{1,k\sigma }},c_{1,{{k}^{\prime }}\uparrow }^{\dagger }c_{1,-{{k}^{\prime }}\downarrow }^{\dagger } \right] \right\} \\ {} & = & -2i\text{Im}\left\{ \underset{k,q,\sigma }{\mathop{\mathop{\sum }^{}}}\,V_{kq}^{12}{{\xi }_{1,k}}c_{1,k\sigma }^{\dagger }{{c}_{2,q\sigma }}+V_{kq}^{12}\left( \Delta _{1,-k}^{*}{{c}_{1,-k\uparrow }}{{c}_{2,q\downarrow }}-\Delta _{1,k}^{*}{{c}_{1,-k\downarrow }}{{c}_{2,q\uparrow }} \right) \right\}. \\ &&\end{matrix}\end{eqnarray} \tag{ A.5 }$$

Substituting this expression in (A.3), the heat current is found to be

$$\begin{eqnarray}&&\frac{\text{d}{{Q}^{(1)}}}{\text{d}t}=\frac{2}{\hbar }\text{Im}\left\{ \underset{k,q,\sigma }{\mathop{\overset{{}}{\mathop{\sum }}\,}}\,\left\langle V_{kq}^{12}{{\xi }_{1,k}}c_{1,k\sigma }^{\dagger }{{c}_{2,q\sigma }}+V_{kq}^{12}\left( \Delta _{1,-k}^{*}{{c}_{1,-k\uparrow }}{{c}_{2,q\downarrow }}-\Delta _{1,k}^{*}{{c}_{1,-k\downarrow }}{{c}_{2,q\uparrow }} \right) \right\rangle \right\}.\end{eqnarray} \tag{ A.6 }$$

The next step is to calculate the expectation values in the general expression for the heat current (A.6), yielding a Kubo formula, when expanding in the small tunnelling matrix elements. In general, to first order in perturbation theory, the expectation value of an operator O(t) is

$$\begin{eqnarray}&&\left\langle O\left( t \right)\right\rangle =-i\mathop{\int }^{}_{-\infty }^{t}\text{d}{{t}^{\prime }}{{\left\langle \left[ O\left( t \right),{{H}_{T}}\left( {{t}^{\prime }} \right) \right]\right\rangle }_{0}}{{e}^{\eta \left( {{t}^{\prime }}-t \right)}}\end{eqnarray} \tag{ A.7 }$$

where the brackets ${\left\langle {\cdot }\right\rangle _{0}}$ denote the equilibrium average with respect to the Hamiltonian ${{H}_{0}}={{H}_{1}}+{{H}_{2}}$ without the perturbation ${{H}_{\text{T}}}$, and η is a small parameter which is eventually taken to zero. Using (A.7), the heat current can be written as

$$\begin{eqnarray}\begin{matrix} \frac{\text{d}{{Q}^{\left( 1 \right)}}}{\text{d}t} & = & -\frac{2}{\hbar }\;\operatorname{Re}\left\{ \mathop{\int }^{}_{-\infty }^{t}d{{t}^{\prime }}{{e}^{\eta \left( {{t}^{\prime }}-t \right)}}\underset{k,q,\sigma }{\mathop{\mathop{\sum }^{}}}\,\left\langle \right. \left[ \left( V_{kq}^{12}{{\xi }_{1,k}}c_{1,k\sigma }^{\dagger }\left( t \right){{c}_{2,q\sigma }}\left( t \right)+ \right. \right. \right. \\ \quad & {} & +\left. {{\left. \left. V_{kq}^{12}\left( \Delta _{1,-k}^{*}{{c}_{1,-k\uparrow }}\left( t \right){{c}_{2,q\downarrow }}\left( t \right)-\Delta _{1,k}^{*}{{c}_{1,-k\downarrow }}\left( t \right){{c}_{2,q\uparrow }}\left( t \right) \right) \right],{{H}_{T}}\left( {{t}^{\prime }} \left. \right\rangle \right) \right)}_{0}} \right\}. \\ &&\end{matrix}\end{eqnarray} \tag{ A.8 }$$

As a first step, we need to again evaluate the commutator expression in the integrand, which assumes the form

$$\begin{eqnarray}&&\begin{matrix} \underset{k,q,\sigma }{\mathop{\mathop{\sum }^{}}}\,\underset{{{k}^{\prime }},{{q}^{\prime }},{{\sigma }^{\prime }}}{\mathop{\mathop{\sum }^{}}}\,\left\{ V_{kq}^{12}V_{{{k}^{\prime }}{{q}^{\prime }}}^{12}\left[ {{\xi }_{1,k}}\left( c_{1,k\sigma }^{\dagger }\left( t \right){{c}_{2,q\sigma }}\left( t \right)c_{1,{{k}^{\prime }}{{\sigma }^{\prime }}}^{\dagger }\left( {{t}^{\prime }} \right){{c}_{2,{{q}^{\prime }}{{\sigma }^{\prime }}}}\left( {{t}^{\prime }} \right)-c_{1,{{k}^{\prime }}{{\sigma }^{\prime }}}^{\dagger }\left( {{t}^{\prime }} \right){{c}_{2,{{q}^{\prime }}{{\sigma }^{\prime }}}}\left( {{t}^{\prime }} \right)c_{1,k\sigma }^{\dagger }\left( t \right){{c}_{2,q\sigma }}\left( t \right) \right) \right. \right. \\ &&\quad -\Delta _{1,k}^{*}\left( {{c}_{1,-k\downarrow }}\left( t \right){{c}_{2,q\uparrow }}\left( t \right)c_{1,{{k}^{\prime }}{{\sigma }^{\prime }}}^{\dagger }\left( {{t}^{\prime }} \right){{c}_{2,{{q}^{\prime }}{{\sigma }^{\prime }}}}\left( {{t}^{\prime }} \right)-c_{1,{{k}^{\prime }}{{\sigma }^{\prime }}}^{\dagger }\left( {{t}^{\prime }} \right){{c}_{2,{{q}^{\prime }}{{\sigma }^{\prime }}}}\left( {{t}^{\prime }} \right){{c}_{1,-k\downarrow }}\left( t \right){{c}_{2,q\uparrow }}\left( t \right) \right) \\ &&\quad +\left. \Delta _{1,-k}^{*}\left( {{c}_{1,-k\uparrow }}\left( t \right){{c}_{2,q\downarrow }}\left( t \right)c_{1,{{k}^{\prime }}{{\sigma }^{\prime }}}^{\dagger }\left( {{t}^{\prime }} \right){{c}_{2,{{q}^{\prime }}{{\sigma }^{\prime }}}}\left( {{t}^{\prime }} \right)-c_{1,{{k}^{\prime }}{{\sigma }^{\prime }}}^{\dagger }\left( {{t}^{\prime }} \right){{c}_{2,{{q}^{\prime }}{{\sigma }^{\prime }}}}\left( {{t}^{\prime }} \right){{c}_{1,-k\uparrow }}\left( t \right){{c}_{2,q\downarrow }}\left( t \right) \right) \right] \\ &&\quad +V_{kq}^{12}V_{{{k}^{\prime }}{{q}^{\prime }}}^{12*}\left[ {{\xi }_{1,k}}\left( c_{1,k\sigma }^{\dagger }\left( t \right){{c}_{2,q\sigma }}\left( t \right)c_{2,{{q}^{\prime }}{{\sigma }^{\prime }}}^{\dagger }\left( {{t}^{\prime }} \right){{c}_{1,{{k}^{\prime }}{{\sigma }^{\prime }}}}\left( {{t}^{\prime }} \right)-c_{2,{{q}^{\prime }}{{\sigma }^{\prime }}}^{\dagger }\left( {{t}^{\prime }} \right){{c}_{1,{{k}^{\prime }}{{\sigma }^{\prime }}}}\left( {{t}^{\prime }} \right)c_{1,k\sigma }^{\dagger }\left( t \right){{c}_{2,q\sigma }}\left( t \right) \right) \right. \\ &&\quad -\Delta _{1,k}^{*}\left( {{c}_{1,-k\downarrow }}\left( t \right){{c}_{2,q\uparrow }}\left( t \right)c_{2,{{q}^{\prime }}{{\sigma }^{\prime }}}^{\dagger }\left( {{t}^{\prime }} \right){{c}_{1,{{k}^{\prime }}{{\sigma }^{\prime }}}}\left( {{t}^{\prime }} \right)-c_{2,{{q}^{\prime }}{{\sigma }^{\prime }}}^{\dagger }\left( {{t}^{\prime }} \right){{c}_{1,{{k}^{\prime }}{{\sigma }^{\prime }}}}\left( {{t}^{\prime }} \right){{c}_{1,-k\downarrow }}\left( t \right){{c}_{2,q\uparrow }}\left( t \right) \right) \\ &&\quad +\left. \left. \Delta _{1,-k}^{*}\left( {{c}_{1,-k\uparrow }}\left( t \right){{c}_{2,q\downarrow }}\left( t \right)c_{2,{{q}^{\prime }}{{\sigma }^{\prime }}}^{\dagger }\left( {{t}^{\prime }} \right){{c}_{1,{{k}^{\prime }}{{\sigma }^{\prime }}}}\left( {{t}^{\prime }} \right)-c_{2,{{q}^{\prime }}{{\sigma }^{\prime }}}^{\dagger }\left( {{t}^{\prime }} \right){{c}_{1,{{k}^{\prime }}{{\sigma }^{\prime }}}}\left( {{t}^{\prime }} \right){{c}_{1,-k\uparrow }}\left( t \right){{c}_{2,q\downarrow }}\left( t \right) \right) \right] \right\}. \\ &&\end{matrix}\end{eqnarray} \tag{ A.9 }$$

In order to take the equilibrium expectation value of this expression, it is useful to employ the Greenʼs functions defined in the following way

$$\begin{eqnarray}\begin{matrix} G_{l,k}^{>}\left( t,{{t}^{\prime }} \right) & = & -i \left\langle {{c}_{l,k}}\left( t \right)c_{l,k}^{\dagger }{{\left( {{t}^{\prime }} \right)}}\right\rangle _{0},\quad G_{l,k}^{<}\left( t,{{t}^{\prime }} \right)=i\left\langle c_{l,k}^{\dagger }\left( {{t}^{\prime }} \right){{c}_{l,k}}{{\left( t \right)}}\right\rangle _{0}, \\ F_{l,k}^{>}\left( t,{{t}^{\prime }} \right) & = & -i\left\langle {{c}_{l,k\uparrow }}\left( t \right){{c}_{l,-k\downarrow }}{{\left( {{t}^{\prime }} \right)}}\right\rangle _{0},\quad F_{l,k}^{<}\left( t,{{t}^{\prime }} \right)=i\left\langle {{c}_{l,-k\downarrow }}\left( {{t}^{\prime }} \right){{c}_{l,k\uparrow }}{{\left( t \right)}}\right\rangle _{0}, \\ F_{l,k}^{\dagger >}\left( t,{{t}^{\prime }} \right) & = & i\left\langle c_{l,-k\downarrow }^{\dagger }\left( {{t}^{\prime }} \right)c_{l,k\uparrow }^{\dagger }{{\left( t \right)}}\right\rangle _{0},\quad F_{l,k}^{\dagger <}\left( t,{{t}^{\prime }} \right)=-i\left\langle c_{l,k\uparrow }^{\dagger }\left( t \right)c_{l,-k\downarrow }^{\dagger }{{\left( {{t}^{\prime }} \right)}}\right\rangle _{0}. \\ &&\end{matrix}\end{eqnarray}$$

The linear-response formula for the heat current, equation (A.8), is then given by

$$\begin{eqnarray}\begin{matrix} \frac{\text{d}{{Q}^{(1)}}}{\text{d}t} & = & \frac{2}{\hbar }\;\operatorname{Re}\underset{k,q}{\mathop{\mathop{\sum }^{}}}\,\mathop{\int }^{}_{-\infty }^{t}\text{d}{{t}^{\prime }}{{e}^{\eta \left( {{t}^{\prime }}-t \right)}}\left\{ V_{kq}^{12}V_{-k-q}^{12}\left[ {{\xi }_{1,k}}\left( F_{1,k}^{\dagger >}(t,t\prime )F_{2,q}^{<}\left( t,t\prime \right)+F_{1,-k}^{\dagger <}\left( {{t}^{\prime }},t \right)F_{2,-q}^{>}\left( t\prime ,t \right) \right. \right. \right. \\ \quad & {} & -\left. F_{1,k}^{\dagger <}\left( t,{{t}^{\prime }} \right)F_{2,q}^{>}\left( t,{{t}^{\prime }} \right)-F_{1,-k}^{\dagger >}\left( {{t}^{\prime }},t \right)F_{2,-q}^{<}\left( {{t}^{\prime }},t \right) \right)+\Delta _{1,k}^{*}\left( G_{1,-k}^{>}\left( t,{{t}^{\prime }} \right)F_{2,q}^{>}\left( t,{{t}^{\prime }} \right) \right. \\ \quad & {} & -\left. \left. G_{1,-k}^{<}\left( t,{{t}^{\prime }} \right)F_{2,q}^{<}\left( t,{{t}^{\prime }} \right) \right)+\Delta _{1,-k}^{*}\left( G_{1,-k}^{>}\left( t,{{t}^{\prime }} \right)F_{2,-q}^{<}\left( {{t}^{\prime }},t \right)-G_{1,-k}^{<}\left( t,{{t}^{\prime }} \right)F_{2,-q}^{>}\left( {{t}^{\prime }},t \right) \right) \right] \\ \quad & {} & -{{\left| V_{kq}^{12} \right|}^{2}}\left[ 2{{\xi }_{1,k}}\left( G_{1,k}^{<}\left( {{t}^{\prime }},t \right)G_{2,q}^{>}\left( t,{{t}^{\prime }} \right)-G_{1,k}^{>}\left( {{t}^{\prime }},t \right)G_{2,q}^{<}\left( t,{{t}^{\prime }} \right) \right) \right. \\ \quad & {} & +\Delta _{1,k}^{*}\left( F_{1,k}^{<}\left( {{t}^{\prime }},t \right)G_{2,q}^{>}\left( t,{{t}^{\prime }} \right)-F_{1,k}^{>}\left( {{t}^{\prime }},t \right)G_{2,q}^{<}\left( t,{{t}^{\prime }} \right) \right) \\ \quad & {} & \left. \left. +\Delta _{1,-k}^{*}\left( F_{1,-k}^{>}\left( t,{{t}^{\prime }} \right)G_{2,q}^{>}\left( t,{{t}^{\prime }} \right)-F_{1,-k}^{<}\left( t,{{t}^{\prime }} \right)G_{2,q}^{<}\left( t,{{t}^{\prime }} \right) \right) \right] \right\}. \\ &&\end{matrix}\end{eqnarray} \tag{ A.10 }$$

The next step is to express the Greenʼs functions by their spectral densities

$$\begin{eqnarray}\begin{matrix} G_{l,k}^{>}\left( t,{{t}^{\prime }} \right) & = & -i\mathop{\int }^{}_{-\infty }^{\infty }\frac{\text{d}\omega }{2\pi }{{e}^{-i\omega \left( t-{{t}^{\prime }} \right)}}\left( 1-{{f}_{l}}\left( \omega \right) \right){{A}_{l,k}}\left( \omega \right), \\ G_{l,k}^{<}\left( t,{{t}^{\prime }} \right) & = & i\mathop{\int }^{}_{-\infty }^{\infty }\frac{\text{d}\omega }{2\pi }{{e}^{-i\omega \left( t-{{t}^{\prime }} \right)}}{{f}_{l}}\left( \omega \right){{A}_{l,k}}\left( \omega \right), \\ F_{l,k}^{>}\left( t,{{t}^{\prime }} \right) & = & i\mathop{\int }^{}_{-\infty }^{\infty }\frac{\text{d}\omega }{2\pi }{{e}^{-i\omega \left( t-{{t}^{\prime }} \right)}}\left( 1-{{f}_{l}}\left( \omega \right) \right){{B}_{l,k}}\left( \omega \right), \\ F_{l,k}^{<}\left( t,{{t}^{\prime }} \right) & = & -i\mathop{\int }^{}_{-\infty }^{\infty }\frac{\text{d}\omega }{2\pi }{{e}^{-i\omega \left( t-{{t}^{\prime }} \right)}}{{f}_{l}}\left( \omega \right){{B}_{l,k}}\left( \omega \right). \\ &&\end{matrix}\end{eqnarray} \tag{ A.11 }$$

Where $V_{kq}^{12}V_{-k-q}^{12}={{\left| V_{kq}^{12} \right|}^{2}}$ is the Fermi function of the lth electrode. Substituting these expressions (A.11) into the equation for the heat current (A.10), the latter simplifies significantly

$$\begin{eqnarray}\begin{matrix} \frac{\text{d}{{Q}^{\left( 1 \right)}}}{\text{d}t} & = & \frac{2}{\hbar }\;\text{Im}\underset{k,q}{\mathop{\mathop{\sum }^{}}}\,\mathop{\int }^{}_{-\infty }^{\infty }\frac{d\omega }{2\pi }\mathop{\int }^{}_{-\infty }^{\infty }\frac{d{{\omega }^{\prime }}}{2\pi }{{\left| V_{kq}^{12} \right|}^{2}}\left( {{f}_{1}}\left( \omega \right)-{{f}_{2}}\left( {{\omega }^{\prime }} \right) \right) \\ {} & {} & \times \left( -\Delta _{1,-k}^{*}\frac{{{A}_{1,-k}}\left( \omega \right){{B}_{2,-q}}\left( {{\omega }^{\prime }} \right)}{{{\omega }^{\prime }}-\omega -i\eta }+\Delta _{1,k}^{*}\frac{{{A}_{1,k}}\left( \omega \right){{B}_{2,q}}\left( -{{\omega }^{\prime }} \right)}{{{\omega }^{\prime }}-\omega +i\eta }+2\Delta _{1,k}^{*}\frac{{{B}_{1,k}}\left( \omega \right){{A}_{2,q}}\left( {{\omega }^{\prime }} \right)}{{{\omega }^{\prime }}-\omega +i\eta } \right. \\ {} & {} & +\left. {{\xi }_{1,k}}\frac{B_{1,k}^{*}\left( \omega \right){{B}_{2,q}}\left( {{\omega }^{\prime }} \right)-{{B}_{1,k}}\left( \omega \right)B_{2,q}^{*}\left( {{\omega }^{\prime }} \right)}{{{\omega }^{\prime }}-\omega -i\eta }+2{{\xi }_{1,k}}\frac{{{A}_{1,k}}\left( \omega \right){{A}_{2,q}}\left( {{\omega }^{\prime }} \right)}{{{\omega }^{\prime }}-\omega -i\eta } \right)\ . \\ &&\end{matrix}\end{eqnarray} \tag{ A.12 }$$

Here, since the tunneling matrix element is invariant under time reversal, we used the relation $V_{kq}^{12}V_{ - k - q}^{12} = {\left| {V_{kq}^{12}} \right|^2}$. According to microscopic BCS theory, the spectral densities are ${{A}_{l,k}}\left( \omega \right)=2\pi \left[ {{\left| {{u}_{l,k}} \right|}^{2}}\delta \left( \omega -{{E}_{l,k}} \right)+{{\left| {{v}_{l,k}} \right|}^{2}}\delta \left( \omega +{{E}_{l,k}} \right) \right]$ and ${{B}_{l,k}}\left( \omega \right)=2\pi {{u}_{l,k}}{{v}_{l,k}}\left[ \delta \left( \omega -{{E}_{l,k}} \right)-\delta \left( \omega +{{E}_{l,k}} \right) \right],$ with ${{\left| {{u}_{l,k}} \right|}^{2}}=1/2\left( 1+{{\xi }_{l,k}}\text{/}{{E}_{l,k}} \right)$, ${{\left| {{v}_{l,k}} \right|}^{2}}=1/2\left( 1-{{\xi }_{l,k}}\text{/}{{E}_{l,k}} \right)$, and the quasiparticle energy-momentum relation ${{E}_{l,k}}=\sqrt{\xi _{l,k}^{2}+{{\left| {{\Delta }_{l,k}} \right|}^{2}}}$. To continue the calculation, it is important to notice that the parameters ${{u}_{l,k}}$, ${{v}_{l,k}}$ and ${{\Delta }_{l,k}}$ are not independent, but that their phases are related by $\Delta _{l,k}^{*}{{v}_{l,k}}/{{u}_{l,k}}={{E}_{l,k}}-{{\xi }_{l,k}}$, such that $\Delta _{l,k}^{*}{{v}_{l,k}}/{{u}_{l,k}}$ must be a real number. That is, the phase of ${{v}_{l,k}}$ relative to ${{u}_{l,k}}$ must be equal to the phase of ${{\Delta }_{l,k}}$. Without loss of generality we can choose ${{u}_{l,k}}$ to be real and positive, so that ${{v}_{l,k}}$ and ${{\Delta }_{l,k}}$ must have the same phase [26]. Finally, we introduce the phase difference φ between the electrodes with the relation $\Delta _{1,k}^{*}{{v}_{2,q}}=\left| \Delta _{1,k}^{*}{{v}_{2,q}} \right|\exp \left( i\varphi \right)$.

The next stage of the calculation is to substitute the spectral densities, ${{A}_{l,k}}\left( \omega \right)$ and ${{B}_{l,k}}\left( \omega \right)$, into the heat current expression (A.12) and to perform the sum over the momenta k and q. To do that, the sum over the momenta is transformed into an integral over the electronic energies ${{\xi }_{l,k}}$ with l = 1, 2, such that [27],

$$\begin{eqnarray}\begin{matrix} \underset{k}{\mathop{\mathop{\sum }^{}}}\,{{A}_{l,k}} & = & \frac{\pi N_{l}^{0}}{2}\mathop{\int }^{}_{-\infty }^{\infty }\text{d}{{\xi }_{l,k}}\left( \delta \left( \omega -{{E}_{l,k}} \right)+\delta \left( \omega +{{E}_{l,k}} \right) \right) \\ {} & = & \pi N_{l}^{0}\mathop{\int }^{}_{\left| {{\delta }_{l,k}} \right|}^{\infty }\text{d}{{E}_{l,k}}\frac{{{E}_{l,k}}}{\sqrt{E_{l,k}^{2}-{{\left| {{\delta }_{l,k}} \right|}^{2}}}}\left( \delta \left( \omega -{{E}_{l,k}} \right)+\delta \left( \omega +{{E}_{l,k}} \right) \right) \\ {} & = & \frac{\pi N_{l}^{0}\left| \omega \right|}{\sqrt{{{\omega }^{2}}-{{\left| {{\Delta }_{l}} \right|}^{2}}}}\theta \left( {{\omega }^{2}}-{{\left| {{\Delta }_{l}} \right|}^{2}} \right). \\ &&\end{matrix}\end{eqnarray} \tag{ A.13 }$$

We denote the normal-state density of states (including spin) of the lth electrode by $N_{l}^{0}$ and finally assumed an isotropic superconductor with an energy-independent gap. Similarly, we find for the other terms

$$\begin{eqnarray}\begin{matrix} \underset{k}{\mathop{\mathop{\sum }^{}}}\,{{B}_{l,k}} & = & {sgn} \left( \omega \right)\frac{\pi N_{l}^{0}\left| {{\Delta }_{l}} \right|}{\sqrt{{{\omega }^{2}}-{{\left| {{\Delta }_{l}} \right|}^{2}}}}\theta \left( {{\omega }^{2}}-{{\left| {{\Delta }_{l}} \right|}^{2}} \right), \\ \underset{k}{\mathop{\mathop{\sum }^{}}}\,{{\xi }_{l,k}}{{A}_{l,k}} & = & {sgn} \left( \omega \right)\frac{\pi N_{l}^{0}\left( {{\omega }^{2}}-{{\left| {{\Delta }_{l}} \right|}^{2}} \right)}{\sqrt{{{\omega }^{2}}-{{\left| {{\Delta }_{l}} \right|}^{2}}}}\theta \left( {{\omega }^{2}}-{{\left| {{\Delta }_{l}} \right|}^{2}} \right), \\ \underset{k}{\mathop{\mathop{\sum }^{}}}\,{{\xi }_{l,k}}{{B}_{l,k}} & = & 0. \\ &&\end{matrix}\end{eqnarray} \tag{ A.14 }$$

In (A.12) the term depending only on ${{B}_{1,k}}{{B}_{2,q}}$, which is related to the sole Cooper pair contribution, vanishes, as it is easy to notice using the last expression in (A.14). If we also use the relation ${{\lim }_{\epsilon \to 0}}\text{Im}\left\{ 1/\left( x-i\varepsilon \right) \right\}=\pi \delta \left( x \right)$ and substitute the integrals over the spectral functions (A.13) and (A.14) into the heat current (A.12), we obtain

$$\begin{eqnarray}\begin{matrix} \frac{\text{d}{{Q}^{\left( 1 \right)}}}{\text{d}t} & = & \frac{\pi }{\hbar }{{\left| V_{kq}^{12} \right|}^{2}}N_{1}^{0}N_{2}^{0}\mathop{\int }^{}_{-\infty }^{\infty }\text{d}\omega \;\omega \frac{\left( {{f}_{1}}\left( \omega \right)-{{f}_{2}}\left( \omega \right) \right)\theta \left( {{\omega }^{2}}-{{\left| {{\Delta }_{1}} \right|}^{2}} \right)\theta \left( {{\omega }^{2}}-{{\left| {{\Delta }_{2}} \right|}^{2}} \right)}{\sqrt{{{\omega }^{2}}-{{\left| {{\Delta }_{1}} \right|}^{2}}}\sqrt{{{\omega }^{2}}-{{\left| {{\Delta }_{2}} \right|}^{2}}}} \\ {} & {} & \qquad \times \left[ {{\omega }^{2}}-\left| {{\Delta }_{1}} \right|\left| {{\Delta }_{2}} \right|\cos \varphi \right]. \\ &&\end{matrix}\end{eqnarray} \tag{ A.15 }$$

In the above derivation we assumed that the normal densities of states and the tunnelling matrix elements are energy-independent. Evaluating the theta-functions we finally obtain

$$\begin{eqnarray}\begin{matrix} \frac{\text{d}{{Q}^{\left( 1 \right)}}\left( {{T}_{1}},{{T}_{2}} \right)}{\text{d}t} & = & \frac{\text{d}Q_{qp}^{\left( 1 \right)}\left( {{T}_{1}},{{T}_{2}} \right)}{\text{d}t}-\frac{\text{d}Q_{\operatorname{int}}^{\left( 1 \right)}\left( {{T}_{1}},{{T}_{2}} \right)}{\text{d}t}\cos \varphi \\ {} & = & \frac{2}{{{e}^{2}}{{R}_{12}}}\mathop{\int }^{}_{\left| {{\Delta }_{\max }} \right|}^{\infty }\text{d}\omega \omega \frac{\left( {{f}_{1}}\left( \omega \right)-{{f}_{2}}\left( \omega \right) \right)}{\sqrt{{{\omega }^{2}}-{{\left| {{\Delta }_{1}} \right|}^{2}}}\sqrt{{{\omega }^{2}}-{{\left| {{\Delta }_{2}} \right|}^{2}}}}\left[ {{\omega }^{2}}-\left| {{\Delta }_{1}} \right|\left| {{\Delta }_{2}} \right|\cos \varphi \right] \\ &&\end{matrix}\end{eqnarray} \tag{ A.16 }$$

where we introduced $\left| {{\Delta }_{\max }} \right|=\max \left\{ \left| {{\Delta }_{1}}\left( {{T}_{1}} \right) \right|,\left| {{\Delta }_{2}}\left( {{T}_{2}} \right) \right| \right\}$ and the normal-state conductance of the Josephson junction, defined via the inverse of the normal-state resistance, $R_{12}^{-1}=\pi {{e}^{2}}N_{1}^{0}N_{2}^{0}{{\left| V_{kq}^{12} \right|}^{2}}/\hbar$.

In (A.16), the total heat current through the junction carried by quasiparticles is $\text{d}Q_{\text{qp}}^{\left( 1 \right)}\left( {{T}_{1}},{{T}_{2}} \right)/\text{d}t$, while $\text{d}Q_{\operatorname{int}}^{\left( 1 \right)}\left( {{T}_{1}},{{T}_{2}} \right)/\text{d}t$ is the interference contribution to the heat current due to an interplay between quasiparticles and the Cooper pair condensate. It is easy to see that $\text{d}Q_{\operatorname{int}}^{\left( 1 \right)}\left( {{T}_{1}},{{T}_{2}} \right)/\text{d}t$, which originates from the Josephson effect and is characteristic of weakly coupled superconductors, vanishes when at least one of the superconductors is in the normal state $\left( \left| {{\Delta }_{l}}\left( {{T}_{l}} \right) \right|=0 \right)$.

In this section of the appendix we provide the analytic formulas for the slopes of the sensitivities ${{s}_{l}}$, for small $\delta f\ll 1$. We find

$$\begin{eqnarray}\begin{matrix} {{m_1}}& = &{\frac{{2\pi \left[ {{{(4{\alpha ^2} - 1)}^2}(2{\alpha ^2} - 1) + \alpha \frac{{{E_{\rm{C}}}}}{{{E_{\rm{J}}}}}(1 + 2{\alpha ^2})} \right]}}{{{{\left( {4{\alpha ^2} - 1} \right)}^{5/2}}\left[ {1 - 2\alpha \frac{{{{\dot Q}_{qp}}}}{{{{\dot Q}_{{\mathop{\rm int}} }}}}{e^{\frac{{{E_{\rm{C}}}}}{{{E_{\rm{J}}}}}\frac{\alpha }{{(4{\alpha ^2} - 1)}}}}} \right]}}} \\ {{m_2}}& = &{\frac{{2\pi \alpha \left[ {(4{\alpha ^2} - 1) + 4\alpha \frac{{{E_{\rm{C}}}}}{{{E_{\rm{J}}}}}} \right]}}{{{{\left( {4{\alpha ^2} - 1} \right)}^{3/2}}\left[ {(2{\alpha ^2} - 1) - 2{\alpha ^2}\frac{{{{\dot Q}_{qp}}}}{{{{\dot Q}_{{\mathop{\rm int}} }}}}{e^{\frac{{{E_{\rm{C}}}}}{{{E_{\rm{J}}}}}\frac{{4\alpha }}{{(4{\alpha ^2} - 1)}}}}} \right]}}} \\ {{m_3}}& = &{\frac{{2\pi \alpha \left[ {{{\left( {4{\alpha ^2} - 1} \right)}^2}\left( {2{\alpha ^2} - 1 + {e^{\frac{{{E_{\rm{C}}}}}{{{E_{\rm{J}}}}}\frac{{3\alpha }}{{(4{\alpha ^2} - 1)}}}}} \right) + \alpha \frac{{{E_{\rm{C}}}}}{{{E_{\rm{J}}}}}\left( {2{\alpha ^2} + 1 + {e^{\frac{{{E_{\rm{C}}}}}{{{E_{\rm{J}}}}}\frac{{4\alpha }}{{(4{\alpha ^2} - 1)}}}}(4{\alpha ^4} - 1)} \right)} \right]}}{{{{\left( {4{\alpha ^2} - 1} \right)}^{5/2}}\left[ {\left( {{e^{\frac{{{E_{\rm{C}}}}}{{{E_{\rm{J}}}}}\frac{{3\alpha }}{{(4{\alpha ^2} - 1)}}}}(2{\alpha ^2} - 1) + 1} \right) - 2\alpha (\alpha + 1)\frac{{{{\dot Q}_{qp}}}}{{{{\dot Q}_{{\mathop{\rm int}} }}}}{e^{ - \frac{{{E_{\rm{C}}}}}{{{E_{\rm{J}}}}}\frac{\alpha }{{(4{\alpha ^2} - 1)}}}}} \right]}}.} \\ &&\end{matrix}\end{eqnarray}$$

The results are plotted in figure 4.