{'fraction_non_alphanumeric': 0.03803915327073054, 'fraction_numerical': 0.01238261976762693, 'mean_word_length': 4.548569410102437, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 8, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Thermal transport through ac-driven transparent Josephson weak links

We discuss how phase coherence manifests in the heat transport through superconducting single and multichannel Josephson junctions in time dependent situations. We consider the heat current driven through the junction by a temperature difference in dc voltage and ac phase biased situations. At low bias, the electromagnetic driving mainly modifies the form of the coherent resonance that transports a large part of the heat current, rather than simply dissipating energy in the junction. We find a description for the heat current in terms of quasiparticle n-photon absorption and emission rates, and discuss analytical and numerical results concerning them. In addition to the ensemble average heat transport, we describe also its fluctuations.


I. INTRODUCTION
Time-dependent effects in superconducting junctions have been under extensive study since the discovery of the Josephson effect. As far as charge transport is considered, the physics is well understood both theoretically and experimentally and has been studied in detail -both for tunnel junctions 1 and for finite junction transparencies. [2][3][4][5][6][7][8] In contrast, time-dependent heat transport through Josephson junctions has attracted interest only quite recently, 9 in the wake of theoretical [9][10][11][12][13][14][15] and experimental 16,17 advances in understanding of how the heat transport is modulated by phase differences of the superconducting order parameter.
In a stationary situation, electronic heat current between two superconductors connected by a weak link is mainly mediated by the transport of above-gap quasiparticles. The heat transport is determined by their energy dependent transmission rates Y (E), [10][11][12][13][14] which in the simplest picture can be found by solving the Bogoliubovde Gennes scattering problem associated with the weak link. 15,18 Due to the presence of phase-coherent Andreev reflections, the transport depends on the phase difference between the superconducting order parameters of the two sides, and moreover exhibits transmission resonances which are of importance especially in low-transparency junctions. [9][10][11][12][13][14][15] Presence of an electromagnetic field in the weak link enables additional possibilities for quasiparticle transfer, as energy can be absorbed from or emitted to the field [see Fig. 1]. Transport is no longer elastic, and processes such as photon assisted tunneling 19 and multiple Andreev reflections 20,21 appear. The energy absorption physics is also related to the quasiparticle dynamics of the Andreev bound states (ABS) residing in the junction. This physics has been studied by a number of recent experiments using quantum point contacts with a small number of transport channels, [22][23][24] diffusive metal junctions, 25,26 and phase dependent absorption has been FIG. 1: (a) Heat current J flows through a Josephson junction between identical superconductors R and L, driven by a temperature difference TR > TL. The simultaneously applied time-dependent voltage V (t) not only dissipates power in the junction, but also modulates the phase coherent effects related to quasiparticle transport, thereby changing J even when the dissipated power is negligible compared to J . (b) Heat is carried by quasiparticles traversing the junction. They can absorb or emit energy by interacting with the electromagnetic field in the junction. For monochromatic ac bias V (t) ∝ sin(ωt), energy absorbed is an integer k multiple of the frequency. Heat transport is fully described by the total rates Y k (E) of the corresponding processes. ∆L and ∆R denote the temperature-dependent energy gaps of the two superconductors. studied in tunnel junctions 27 .
To describe the effects of time-dependent drive on the heat current carried by quasiparticles, the energy dependent transmission rates Y (E) of the stationary state need to be generalized to depend on both the initial E and final E quasiparticle energies Y (E, E ). If the junction is driven with a periodic signal at frequency ω 0 (or dc voltage eV = ω 0 ), the final energy is E = E − kω 0 , and the results can be described in terms of a set of rates {Y k (E)} related to emission/absorption of k photons to the field, as illustrated in Fig. 1(b).
Heat transport and its fluctuations in the presence of time-dependent driving have been recently discussed in terms of Floquet scattering matrices in Refs. 28,29 for normal-state systems. In this description, the rates {Y k (E)} can be directly related to the Floquet scattering matrices. Although these works consider normal-state systems, the formal description of electronic heat current in time-dependent superconducting junctions turns out to be similar. This should be contrasted to charge current: the superconducting condensate transports charge without requiring quasiparticle transfer, which implies that transport of charge and heat are not as closely connected in superconductors as in normal-state systems.
In this work, we study theoretically the heat transport in superconducting junctions of general transparency driven by classical fields, as shown in Fig. 1. We derive results for the rates {Y k (E)} and for the heat transport and its fluctuations using quasiclassical Green function theory. Previously, a similar approach has been used in Ref. 30 for studying time-dependent heat transport in NIS junctions, and it has been extensively applied for examining the charge transport (see e.g. Refs. 2,6,8,31). For tunnel junctions, the approach here is equivalent to that of Ref. 9, but it can also describe the time-dependent heat transport in junctions with a finite transparency away from the tunneling limit.
We consider situations with dc and ac bias, and discuss how the electromagnetic driving modulates the temperature-driven heat current. In particular, for weak driving, we find the effect is mainly dominated by a modulation and phase-averaging of the coherent resonances appearing in the rates {Y k (E)}. We find analytical descriptions in certain limiting cases, and discuss representative features in terms of numerically exact results.
This manuscript is structured as follows. Section II introduces the theoretical and numerical approaches used. Section III A discusses the effect of a dc voltage bias on the heat transport, and III B concentrates on the ac bias case. Section IV summarizes the main points. We postpone details of the derivation of the analytical results to a number of appendixes, referring to the salient points in the main text.

II. MODEL
We consider the situation depicted in Fig. 1: heat transport between two identical superconducting terminals that are connected by a generic contact described by the transmission eigenvalues {τ j } of the normal-state scattering matrix of the junction. The terminals are kept at different constant temperatures, T L and T R , and a time dependent voltage or phase bias is applied between them. Moreover, they are assumed to have sufficient impurity concentration to lie in the dirty limit. To study the transport in this system, we make use of the Keldysh-Nambu Green function formulation for transport in superconducting structures, [32][33][34] and the quasiclassical boundary condition description of a weak link between bulk superconductors in the diffusive limit. 2,35,36 At equilibrium, electrons inside a superconducting terminal at temperature T with superconducting gap ∆ are described by the quasiclassical equilibrium Green functionǧ 0 (E), Rτ 3 , whereτ 3 is the third spin matrix in the Nambu space. We assume the energy gap ∆(T ) depends on the temperature according to the BCS gap relation. 20 The Green function describing a superconductor at a nonzero electric potential V is found by a gauge transform where . We consider the phase as a classical field, neglecting its quantum fluctuations, which simplifies its Keldysh structure.
For cases where the frequencies in the drive φ and the dc bias voltage are commensurate, with a base harmonic ω 0 , we can use a (Floquet) matrix representation 38 , In this representation, time convolution ∞ −∞ dt 1 A(t, t 1 )B(t 1 , t) is equivalent to a matrix product, k A nk ( 0 )B km ( 0 ). Below, products between two-time quantities with time or energy arguments omitted imply time convolutions.
We also denote φ(t, t ) ≡ δ(t − t )φ(t), so that Eq. (3) is written aš g = e iτ3φǧ 0 e −iτ3φ . The charge and heat currents flowing between two connected superconductors can be found from boundary conditions applicable to quasiclassical Green functions. 2,35 Transport properties of the junction are described by a matrix current 35 whereǧ L andǧ R are the Green functions inside the terminals at the two sides of the junctions. In particular, the charge current, and the heat current entering terminals L and R read (cf. also Ref. 30), whereσ are spin matrices in Keldysh space, and (t, t ) = i∂ t δ(t − t ) and µ(t, t ) = δ(t − t )∂ t φ(t). Moreover, we fix the potential of the right terminal to zero. The heat currents as defined above satisfy conservation of energy at each instant of time, µ(t)I(t) =Q L (t) +Q R (t).
To make a connection with Ref. 9, it is advantageous to rewrite the heat currents in a different form, using e iφ = ( + µ)e iφ , where,ġ R = g R,0 andġ L = e iφτ3 g L,0 e −iφτ3 . The quantity appearing in the second terms is The total time derivative ∂ t Q J does not contribute to time-averaged heat currents. Below, as far as heat currents are concerned, we mostly discuss the time-averaged quantities,Q dc,L/R , as these are experimentally more easily accessible.
In the case of low-transparency tunnel junctions, we haveŽ = 1/R T , where R T is the tunnel resistance. The approach then coincides with that discussed in Ref. 9. Comparison to the tunneling Hamiltonian calculation of Ref. 9 identifies first terms of Eqs. (9), (10) as the heat carried by quasiparticles into the bulk of the terminals. Moreover, the term Q J ∼ H T (t) , the component of the internal energy associated with the part of the Hamiltonian coupling the two terminals together. For slow phase variations,φ → 0, one can verify that the expression is related to the Josephson energy of the junction at equilibrium, ∂ t Q J (t) → I J (t)φ(t) = ∂ t E J (φ(t)).
One can compute the sums over τ n for different transparency distributions. A relevant case is the eigenvalue distribution corresponding to a dirty interface 39,40 which has been experimentally found to agree with hightransparency oxide junctions. 41 From this, The result for diffusive junctions in the short junction limit (dwell time τ D /∆) can be found either by solving the Usadel equation, 34,42,43 or averaging over the corresponding transparency distribution 44 . The former approach readily yields: As discussed in Ref. 15, the resulting heat current in diffusive junctions in the absence of voltage bias turns out to be independent of the phase difference ϕ over the junction.
The above approach can also be used to compute statistics of heat transport. The problem can be formulated as follows: 45,46 the internal energy Q R of one terminal is probed two times, first at t = 0 and a second time at t = t 0 . The cumulant generating function can in the long-time limit (t 0 → ∞) be written as [45][46][47] The resulting statistics in the absence of driving are discussed in Ref. 15, and results including driving are for the tunnel limit discussed in Ref. 9. Below, as far as fluctuations are concerned, we mainly concentrate on the zero frequency heat noise The above formulation is directly accessible to numerical calculations in the matrix representation (4). The numerical results below are found by keeping only a finite number of the harmonics k and truncating the energy range to energies within several |∆|, T from the Fermi level. Matrix inverses, exponentials, square roots, and logarithms can be computed efficiently numerically. 48,49

III. DRIVEN QUASIPARTICLE TRANSPORT
The electronic dc heat current is carried by quasiparticles, which are either transferred from one side to the other, or created or destroyed in pairs by breaking or creating Cooper pairs in the condensate. During their traversal through the junction area, the quasiparticles can either gain or lose energy via interaction with the electromagnetic field.
Within the formalism described here, the dc heat current entering terminal α = L/R can in general be written in terms of the quasiparticle distributions in the two terminals: Expressions for Y αβk (E) can be found with some algebraic manipulations discussed in Appendix B. The functions Y α,β,k (E) can be understood to be proportional to the total rates of processes in which a quasiparticle starting at energy E + kω 0 in terminal β ends up in terminal α at energy E while losing or gaining energy kω 0 by interacting with the electromagnetic field in the junction. One should note that the rates Y do not correspond to elementary energy transfer processes in the counting statistics sense. Rather, they are simply proportional to time-averaged total rates. Analogous expressions as Eq. (17) for the heat current can in normal-state systems be written in terms of the Floquet scattering matrices. 28,29 The total power dissipated in the junction is In general, it is directly connected to the charge current in the junction, by conservation of work done by the bias, P = IV . For the discussion below, we also define the temperature-driven heat current as The result is a well-defined temperature-driven current provided the dissipation is low, P |J |, which requires a large enough temperature difference between the two sides of the junction.

A. Dc voltage bias
In dc voltage biased junctions with a finite transparency, the well-known multiple Andreev reflection processes (MAR) 20,21 play a role in the quasiparticle transport. They lead to diffusion of quasiparticles upward in energy inside the junction, and in general give rise to several side bands Y LR,k =0 , corresponding to quasiparticles Andreev reflected several times inside the junction.
Similarly as for the charge transport 21,31,50 , the calculation of the heat transport and its fluctuations in the dc voltage biased case reduces to solving a relatively simple recurrence equation, describing the propagation of quasiparticles upwards in energy. We defer the derivation of the equations to Appendix C. They can be solved analytically in certain limiting cases, and are in general straightforward to solve numerically.
The tunnel junction limit is one of the analytically tractable cases, and for it one finds 9,11 where N L/R are the densities of states of the superconductors. The result is shown in Fig. 2 for reference. The resulting expression for the electronic heat current is similar to that in normal-state junctions, with the superconductivity only modifying the quasiparticle densities of states of the two terminals.
In junctions with finite transparency, Andreev reflections as expected start to facilitate quasiparticle transport across the energy gap region. The corresponding factors Y k are plotted in Fig. 3 for a fully transparent junction (τ = 1). One can note that in addition to the direct quasiparticle transmissions, one obtains MAR processes for |k| > 2∆/V .
Results for the heat current under the simultaneous application of a temperature difference and bias voltage are shown in Fig. 4, corresponding to a quantum point contact with partial transparency. One can note that features of multiple Andreev reflections do not appear to play a large role in the temperature-driven heat current as long as P J . This can be seen explicitly by considering the contributions in Eq. (17) separately for each sideband k (not shown), which indicates direct quasiparticle transfer (Y k=±1 ) is responsible for the low-bias features. At higher bias when P J , especially in the hightemperature curves for which ∆ R (T R ) < ∆ L (T L ), MAR processes do not contribute in the same way toQ R anḋ Q L . This results to sharp changes in J = (Q R −Q L )/2 at voltages close to thresholds of the different possible MAR processes.
The nonmonotonic behavior at low voltages in arises from modification of the resonance in Y LR,±1 at |E| > ∆. The evolution of the resonance in Y with increasing bias voltage is shown in Fig. 5. Initially, transport at energies immediately above the gap is reduced. A somewhat similar effect also occurs in the tunnel junctions, as described by Eq. (20), where the peak in Y LR,1 near E = ∆ is reduced from N (E) 2 to N (E)N (E + V ). The behavior of partially transparent junctions differs from tunnel junctions in that the amplitude of the resonance increases again at larger voltages, resulting to a non-monotonic behavior in the heat current. The effect is the largest at low temperatures, as visible in Fig. 4, as transport close to gap edges has the most importance in that case.  The above nonmonotonic behavior becomes more prominent as the junction transparency increases and the suppression of heat current is lifted. As expected, MAR processes also become more significant as the transparency increases: the higher the transparency, the lower the bias voltage V at which the dissipated power P overwhelms the heat current J .
How multiple Andreev processes contribute to the heat transport is also visible in the heat transport noise, shown in Fig. 7. At low bias, the voltage dependence of the noise follows that of the heat current itself, a feature similar to what occurs in a stationary situation without bias. 15 At high bias, the noise increases rapidly as the MAR processes activate. The MAR features are best visible in the ratio S E /Q, which has sharp features at the activa- tion thresholds V = 2∆/n. For a Poissonian process of absorbing energy packets E 0 , the ratio S E /Q = E 0 indicates directly the energy transferred by a single elementary process. Although the situation in the present case is more complicated, the fact that order n = ceil(2∆/V ) MAR process is associated with absorption of energy E 0 ∼ nV ∼ ceil(2∆/V )V results to oscillatory behavior similar as seen in Fig. 7 (right panel).
Finally, Fig. 8 shows the bias dependence of the temperature-driven heat current for the dirty interface and diffusive channel distributions. The behavior of the heat current in these cases is dominated by the large number of low-transparency channels, and is mainly similar to low-transparency junctions in that the heat current decreases with increasing bias without significant oscillations. The peaks appearing in the dirty interface case occur at V = ∆ L (T L )−∆ R (T R ), and are due to matching the peaks in the above-gap DOS in the superconductors. A similar feature is present also in tunnel junctions. 9 In contrast, the results for the diffusive distribution contain less low-transparency channels, and the gap-difference features are not clearly resolved.

B. Ac phase bias
Consider now the case of a junction under monochromatic ac phase bias, φ(t) = ϕ 0 /2 + 2s cos(ωt), where s = V ac /(2 ω). Such an excitation in general has two effects -first, the microwave photons ω can be absorbed, which can result to dissipation in the form of photonassisted tunneling and breaking of Cooper pairs. Second, the time dependent phase modulates the coherent transmission resonance responsible for carrying a significant portion of the heat current.
Let us first discuss the dissipation in the junction, for which Cooper-pair breaking transitions from negative to positive energies are important. In particular, the Andreev bound states are crucial for this, as they reduce the frequency threshold for Cooper-pair breaking. 4,5,51 An ABS at energy E A < ∆ enables pair breaking by 1-photon processes to occur via processes involving only the bound states at ω = 2E A and via processes involving also the continuum at ω > ∆ + E A , as also seen e.g. in Refs. 5,6. These frequency thresholds occur prominently in the dissipated power shown in Fig. 9. 58 Analytical results for the dissipated power in the representative case ∆ L = ∆ R = ∆, T L = T R = T in the limit s → 0 can be deduced from results in Ref. 5. Figure 9 also displays the heat transport noise, computed numerically via Eq. (15). At high frequencies, the result converges towards S E → ωQ R , and one can check that also the third cumulant ∂ 3 u S R (u)| u=0 ≈ it 0 ω 2Q R in this limit. That the cumulants coincide suggests that photon absorption follows Poisson-like statistics. The situation at lower frequencies ω < ∆ + E A is somewhat more complicated, as multiple photons are required for each absorption event. For Poissonian multiphoton processes, we expect S E = nωQ R for (E A + ∆)/n < ω < (E A +∆)/(n−1). This is in line with the results in Fig. 9, apart from features that come from the second possible transition process involving only the bound states.
At frequencies below the transition thresholds, dissipation by photon absorption via processes involving the ABS become less important. Rather, the main physics comes from quasiparticles in the continuum and from the modification of their phase dependent transport resonances. Analytical results can be obtained by taking the drive amplitude s as a small parameter. Working to order s 2 and considering only processes in the continuum (E > ∆), we find for a fully transparent QPC (τ = 1), A ω = 2(E + ω) 2 − (1 + cos ϕ 0 )∆ 2 , The divergence in Y RL,0 at |E| = ω + ∆| cos(ϕ 0 /2)| is to be understood in the principal value sense -the above result does not describe the dissipative processes involving the Andreev bound states. Details of the calculation and results for τ = 1 are in Appendix D. In the tunneling limit, τ → 0, on the other hand, The results for τ = 1 are shown in Fig. 10. The energy dependence of the driven response reflects the behavior in the stationary state: there is a resonant contribution at energies immediately above the superconducting gap, which is for transparent junctions reduced when the Andreev bound states at E A = 1 − τ sin 2 (ϕ/2) separate from the gap edges.
The resulting change δG th = G th (s) − G th (s = 0) in the heat conductivity G th = J /(T L − T R ) in the limit s → 0 is shown in Fig. 11 for a transparent junction and in Fig. 12 for a low-transparency junction. For the transparent junction case, the effect of the ac excitation is the largest at phase differences ϕ 0 ≈ 0, as expected on the basis of Fig. 10. Moreover, the excitation tends to suppress the heat flow more than it enhances it. The dependence on ϕ 0 in the tunneling limit is opposite to the transparent limit, reflecting the similarly reversed phase dependence in the stationary case 13 . In the tunneling limit, the resonance has the largest effect at ϕ 0 ≈ π, and the ac drive reduces its contribution to the heat current.
The amount of power dissipated relative to the change in heat current for a 0.1∆ 0 temperature difference is also indicated in Figs. 11-12 with contour lines. In the regime shown, power dissipation is mainly due to processes at E > |∆|, which are well described by Eqs. (22)- (24). The features seen in the left panels mainly come from the fact that the change δG th reverses its sign around ϕ 0 ≈ π/4 both for tunnel and transparent junctionsas seen in Figs temperature and frequency dependent. The main feature of the dissipation in the right panel is that if the phase excitation amplitude is kept fixed, the dissipated power grows with increasing frequency, whereas the effect of the drive on the temperature-driven heat current saturates towards higher frequencies. Moreover, at sufficiently high frequency, the Cooper pair breaking processes activate, and start to dominate the heat flows.
Phase oscillation of heat conductivity as a function of increasing drive amplitude is shown in Fig. 13. Initially, the response is well described by the leading-order theory, but starts to deviate from the s → 0 limit when the drive amplitude exceeds s 0.1. The result moreover stays close to the adiabatic limit G th (ϕ(t)), obtained by averaging the stationary heat-current-phase relation over the drive cycle.
The heat current in diffusive short junctions is independent of the phase difference ϕ 0 . 15 Nevertheless, the response of the heat current to a small ac excitation is phase dependent. This is illustrated in Fig. 14, which shows the temperature-driven heat current J = (Q R −Q L )/2 as a function of ϕ 0 . The inset of Fig. 14 displays the comparison to results from Eqs. (D5)-(D12), which neglect transitions involving the sub-gap Andreev bound states in the junction, but still capture the main features. The changes in J are however not very large, and are of a similar order of magnitude or smaller compared to the power dissipated in the junction.
The power dissipated in a diffusive junction grows rapidly at phase differences away from zero. Such behavior is expected due to (multi)photon absorption processes involving the sub-gap bound states, as these processes are limited by the presence of a phase-dependent energy gap 2E g = 2|∆|| cos(ϕ/2)| in the density of states inside the junction. 5 As ϕ increases, the energy gap decreases, and n-photon processes activate after corresponding thresholds are crossed: The onsets of dissipation in the numerical results in Fig. 14 align with these thresholds, even though the phase oscillation amplitude 4s is not very small compared to π. Note that in contrast to the above-gap heat transport, the sub-gap absorption physics can be sensitive to electron-phonon and other inelastic relaxation processes, 52,53 as excited sub-gap quasiparticles cannot easily escape to the leads. Such mechanisms are neglected in the above results, and the sub-gap populations are determined via the drive-induced multiphoton coupling to the above-gap states.

IV. DISCUSSION
The effect of electromagnetic driving on heat transport in superconducting junctions is twofold: first, quasiparticles can absorb energy from the field and the field can break Cooper pairs, both of which leads to dissipation. Second, a nonzero electromagnetic field inside the junction implies a time dependence in the the phase difference in superconducting order parameters. This results to modulation and averaging of phase coherent transport effects, which does not directly dissipate power in the junction.
In general, we find that at low drive frequencies or voltages, there is a finite parameter range in which the response of the system is non-adiabatic while the power dissipated by the electromagnetic fields is still small compared to the temperature-driven heat currents. That is, time-dependent effects associated with modification of the transport resonances of the continuum quasiparticles can be seen before the dissipation heats up the system significantly.
The steady-state effects on the temperature-driven heat current J can in principle be measured with the help of established experimental thermometry techniques 54 for studying heat in mesoscopic structures, similar to those applicable to a stationary situation. In addition, the energy absorption P is experimentally accessible via certain spectroscopic approaches, used for quantum point contacts in Ref. 24 and it can also be deduced from susceptibility measurements, used in Ref. 25 for diffusive junctions.
In summary, we consider heat current driven by a temperature difference across superconducting junctions of different transparencies. We obtain analytical and numerical results for the transport rates describing transport combined with energy absorption of k photons from the electromagnetic field. The results are applicable to quantum point contacts and to generic quantum coherent multichannel junctions.
From the above discussion, and the symmetries which is a form useful for analytic considerations. In particular, this can be rewritten in the form used in Sec. III, with The results corresponding toQ L,dc are found by exchanging φ ↔ −φ and L ↔ R. A sum rule applies, where the second step follows from ( RĝR Rτ 3 , and the fact thatĝ R is diagonal in the energy space. Finally, let us point out how the heat current in diffusive junctions is phase-independent in the absence of bias. In the absence of spectral broadening, quasiparticle heat current can only originate from region |E| > |∆ L |, |∆ R |. From the form of the equilibrium Green functions (1) it follows that for |E| > |∆ L |, |∆ R |, .

(B10)
Therefore, for this energy range, Integrating over thej distribution yields which indeed has no phase modulation. It is also clear that phase modulation in general prevails if the distributionj(ζ) is nontrivial.

Appendix C: DC voltage bias recurrence
We proceed now to evaluate the cumulant generating function of energy exchange 15,45,47 Tr ln[q n e iφτ3/2ǧ L,0 e −iφτ3/2 (C2) + e −iφτ3/2 e iu σ1/2ǧ R,0 e −iu σ1/2 e iφτ3/2 ] in the dc voltage-biased case. One can simplify the problem by reusing the approach of Ref. 31, where a similar generating function for the charge rather than the heat current is discussed. The key idea is to make use of the fact that the sum of the Green functions is a blocktridiagonal matrix in energy space, so that one can apply a determinant recursion formula In the case here, it is useful to work in a gauge where the superconducting terminals are biased to ±V /2, so that the frequency components read where E n = E 0 + nV , and we work in Nambu⊗Keldysh instead of Keldysh⊗Nambu: . Expressions foȓ f are obtained by replacing cosh → sinh.
In the recursion formula for the generating function, we therefore have Clearly, a j stays diagonal in the recursion. Finding the generating function now reduces to solving a pair of 2 × 2 matrix recurrences v(E + V, u) = q ngL (E + V ) +g R (E, u) (C13a) where a normalization was added to make the result convergent; this does not affect observables. Differentiation yields directly the average heat current, ∂ u S R | u=0 = iQ R,dc .
The recursion can be solved numerically in a straightforward way, asf → 0 for |E| → ∞. Analytical results can be obtained in some particular limiting cases.

Tunnel limit
We can consider the solution of the recursion in the tunnel limit q → ∞. There, to order O(q 0 ), it reads where we used the normalization conditionsg 2 −f 2 =1, gf =fg. From this, This yields the result for heat transport by tunneling quasiparticles, also obtained in Ref. 9, where N = Re cosh θ are the densities of states on both sides, and F L/R Fermi functions.

Zero bias
In the zero-bias limit V → 0, the results coincide with averaging the corresponding d.c. formulas in Refs. 13,15 over the phase difference ϕ = 0 . . . 2π. However, in this case it is also possible to solve the recurrence relation exactly.
We first observe that a factorization applies: and similarly forf . For L → R, theh R matrix obtains the additional u-dependent unitary transformation (C9).
This V → 0 heat transport statistics deviates from the stationary Levitov-Lesovik form 56,57 due to the averaging over the slowly varying phase ϕ(t) = 2V t in the long time limit. Expanding in u, we find so that the heat current iṡ We can verify that this coincides with the average of the d.c. result in Ref. 13: For comparison with the inset of Fig. 5, let us also state the result for τ = 1/2: which coincides with numerics.