Fluctuation of heat current in Josephson junctions

We discuss the statistics of heat current between two superconductors at different temperatures connected by a generic weak link. As the electronic heat in superconductors is carried by Bogoliubov quasiparticles, the heat transport fluctuations follow the Levitov--Lesovik relation. We identify the energy-dependent quasiparticle transmission probabilities and discuss the resulting probability density and fluctuation relations of the heat current. We consider multichannel junctions, and find that heat transport in diffusive junctions is unique in that its statistics is independent of the phase difference between the superconductors. Curiously, phase dependence reappears if phase coherence is partially broken.


I. INTRODUCTION
Heat transport through junctions between superconductors is significantly affected by superconductivity. 1 In tunnel Josephson junctions, superconducting phase coherence manifests as a component of the thermal conductance that oscillates with the phase difference between the superconductors, 2-5 a prediction which was confirmed by recent experiments. 6,7 In general, both the sign and the magnitude of the oscillations depend on the transparency of the junction in question. 4 Studies of heat transport in Josephson junctions have largely concentrated on the ensemble average value of the heat currents. However, the heat transferred in a given time is in general not constant, but has fluctuations. 8,9 When mesoscopic systems are considered, fluctuations in heat transport can lead to fluctuations in other quantities -such as the energy stored on a small metal island -and eventually, in more accessible observables, such as charge current [10][11][12] or temperature measured by a generic temperature probe. 13 In addition to the theoretical question on how the coherent physics of phase differences in superconducting order parameters manifest in statistical properties of heat transport, questions on fluctuations can also be of interest in systems that utilize mesoscopic superconductors in a nonequilibrium settings for example for radiation detection. 14 In this work, we find the statistics of temperature-driven quasiparticle heat transport in a Josephson junction of arbitrary transparency, as illustrated in Fig. 1. Previously, the energy transport statistics has been discussed in the tunneling limit, 15 and we recover these results as a special case. The energy transport separates to elementary quasiparticle transport events. In agreement with the Levitov-Lesovik formula, 16,17 the elementary quasiparticle transmission probabilities are governed by the total transmission eigenvalues T n± (E) [Eq. (13)] of the Bogoliubov-de Gennes scattering problem of the interface. We compute the heat current noise and discuss its parameter dependence in single-channel junctions. We derive results corresponding to dirty-interface 18 1. (a) Heat current I E flows from superconductor L to another superconductor R coupled to it, driven by a difference in the temperatures T L > T R . The heat current is modulated by the phase difference ϕ between the order parameters ∆ L , ∆ R of the superconductors. The Josephson junction connecting the two is described by the set of (spin-independent) transmission eigenvalues {τ n } of the normal-state scattering matrix of the interface. (b) The excitations that carry heat in superconductors are Bogoliubov quasiparticles, whose density of states N (E) is illustrated. Heat current and its fluctuation statistics is fully determined by their energy-dependent transmission probabilities T n± (E). Each normal-state quantum channel (τ n ) splits into inequivalent particle-hole transmission channels (±) due to superconductivity.

II. MODEL
We consider heat transport between two superconducting terminals that are connected by a generic contact described by the transmission eigenvalues {τ j }. The terminals are assumed to lie at different temperatures. We make use of the Keldysh-Nambu Green function formulation for transport in superconducting structures, [20][21][22] and the quasiclassical boundary condition description of a weak link between bulk superconductors in the diffusive limit. 23,24 As in Refs. 2-5, we assume the superconducting order parameter phase difference ϕ across the junction is kept constant by a suitable biasing circuit, similarly as in the previous experiments. 6,7 At equilibrium, electrons inside a superconducting terminal at temperature T with superconducting gap ∆ are described by the quasiclassical equilibrium Green functionǧ(E), in Keldysh(ˇ)-Nambu(ˆ) space. Here,ĝ a = −τ 3 (ĝ r ) †τ 3 , whereτ 3 is the third spin matrix in the Nambu space. Temperature enters in the equilibrium distribution function h 0 = tanh E 2T . To model inelastic effects in the superconductors, we introduce a phenomenological Dynes parameter by replacing E  → E ± iΓ/2 in g r /a , where Γ is a relaxation rate modeling possible microscopic mechanisms. 25,26 The statistics of heat transferred in time t 0 in a two-point measurement protocol can be described via the generating functions 8,9,27 which correspond to changes in the internal energy of terminals α = L, R. Here, H α are the BCS Hamiltonians of the superconducting terminals, u is a scalar counting field, and U(t) the time evolution operator. How to calculate this generating function using a Keldysh approach is discussed in Refs. [27][28][29][30]. Main features of the approach can be seen by first differentiating Eq. (3) to obtain, where ⟨X⟩ u,α = Tr[XU +,α ρ(0)U † −,α ]/Tr[U +,α ρ(0)U † −,α ] is an expectation value computed with modified time evolution operators U ±,α = e ±iu H α /2 Ue ∓iu H α /2 including the counting field with differing signs on Keldysh branches during time [0,t 0 ]. As noted in Ref. 27, this corresponds to time shifts on observables corresponding to lead α, similarly to charge counting where the counting field also shifts the conjugate variable of the measured observable. In the Fourier transformed Green function representation 25 above, such time shifts on the Keldysh contour are represented by 27 The expectation value in Eq. (4) can be evaluated via the quasiclassical boundary condition approach, 23,24 and integrating in u produces the action of superconducting contacts [27][28][29][30] ln where C is a normalization constant, and Tr includes energy integration in addition to Keldysh-Nambu matrix trace. Note that the counting field enters as time shifts in a reservoir correlation function, 27 similarly as in the master equation formulation of Ref. 9. Here and below, we fix α = R and describe the statistics of energy transferred into the right terminal.

III. GENERATING FUNCTION
An important difference in energy transport compared to charge statistics follows from the fact that an Andreev reflection does not transfer energy. In the present formulation, this is visible in the fact that ( There is no energy transfer at sub-gap energies, where there are no quasiparticles, assuming no broadening in the spectrum of the superconductors. The generating function can be found by direct substitutions into Eq. (6). It is however useful to make use ofǧ 2 L/R = 1 and rewrite where q n = −1 + 2/τ n + 2 √ 1 − τ n /τ n are the eigenvalues of the hermitian square of the corresponding transfer matrix, 28,31 where C ′ is a normalization constant. The product structureǧ =ĝ r ⊗V of Eq. (7) implies that the Keldysh and Nambu components can be diagonalized separately. This yields where t 0 is the measurement time, and {λ, 1/λ} and { µ, 1/µ} are the eigenvalues ofĝ r Lĝ r R anď V LVR (u), respectively.
Solving the eigenvalue problems, Eq. (10) can be rewritten in the form of a characteristic function of a multinomial distribution: with the event probabilities Here, f L/R are the electron Fermi distribution functions in the left and right terminals, and Here, ϕ is the phase difference between the order parameters of the two superconductors. The generating function describes events where a quasiparticle at energy E > ∆ attempts to move from the left to the right (k = +1) or from the right to the left (k = −1). The probability that such a transmission event succeeds depends both on the transparency of the transmission channel (q n , τ n ), and on superconductivity (λ). The generating function is non-periodic in u, and the corresponding probability density, which is obtained as a Fourier transform of W R (u), 8,9 describes a continuous distribution. 27 This is in contrast to charge transport, where the fact that transferred electronic charge is quantized to multiples of e is reflected as a periodicity in the counting field, S( χ) = S( χ + 2π), also in the superconducting case. 29 Note however that the calculation for statistics of charge current between superconductors leads to "event probabilities" not satisfying 0 ≤ p ≤ 1, attributed to a problem in interpreting the supercurrent in terms of classical transport events. 29 This problem does not arise in the results discussed here, as the supercurrent flow does not contribute directly to energy transfer.
The probabilities T n,± are the transmission eigenvalues of a Bogoliubov-de Gennes (BdG) scattering problem. The corresponding scattering matrix (for ∆ L = ∆ R ) can be found in Ref. 32. Direct evaluation of the transmission eigenvalues using the results there gives which coincides with Eq. (13) above. A similar connection can be made to results in N/S junctions, 33 (∆ L = 0, ∆ R = ∆ > 0), after identifying the barrier reflectivity in Ref. 33 as Z = (q − 1)/  4q. That the energy statistics is related to these scattering matrices stems from the fact that the counting statistics of Bogoliubov quasiparticles, which carry all of the electronic energy current, must follow the Levitov-Lesovik result. 16,17 Superconductivity has a significant impact on the transmission eigenvalues. In the normal state, they are simply T n β = 4q n /(1 + q n ) 2 = τ n independent of the particle-hole channel index β = ±1. In the superconducting state, they deviate significantly from this, as illustrated in Fig. 2. In particular, there is a transmission resonance in one of the channels: The above result applies for ∆ L = ∆ R , but the resonance appears also for ∆ L ∆ R . This has a large effect especially for junctions whose normal-state transparency is small, in which a large part of the total heat current is carried by the resonance. 4

A. Including broadening of density of states
We can extend the above results to include broadening of the density of states in the superconductors, by adding a Dynes parameter Γ > 0. Also in this case, the final generating function obtains the form (11). As the factorization (7) does not apply, the expressions for the transmission eigenvalues T n β (E) need to be extracted from the determinant in Eq. (9). Straightforward calculation yields for T n,± (E) the indirect expressions T n,+ T n,− = 16q 2 n cos 2 (Im θ L ) cos 2 (Im θ R ) where θ L/R = arctanh ∆ L/R E+i(Γ/2) . The above formulas reduce to Eq. (13) for Γ → 0 + . For Γ > 0, the results cannot however be written in terms of a single λ as in Eq. (13).
The eigenvalues T n,± (E) are plotted for Γ > 0 in Fig. 3. In general, a finite Γ induces a small number of states inside the gaps of the superconductors, so that the junction is transparent for heat transport also at sub-gap energies. Similarly as for above-gap transport, the transmission through the two channels ± can differ significantly. Moreover, we can observe that sharp sub-gap resonances appear in one of the two channels -these are associated with the Andreev bound states and are most prominent in transparent junctions. Moreover, we note that the phase dependence of sub-gap heat transport in high-transparency channels (left panel) can be substantial, whereas it is not as prominent at low transparencies (right panel).

IV. HEAT FLOW STATISTICS
Let us now compute the first moments of the heat statistics. Direct differentiation of Eq. (11) with respect to u yields the average heat current, where r n = 1 − τ n . This result coincides with that found in Ref. 4. Note, however, that the coefficients D ee , D he defined in Ref. 4 do not coincide with T n± , even though the sums do.
Taking the second derivative, we find the heat current zero-frequency noise In the tunneling limit, this result coincides with that found in Ref. 15. The behavior of the heat current noise away from equilibrium is shown in Fig. 4, for a fixed value of T R . These results take into account the temperature dependences ∆ L (T L ), ∆ R (T R ) of the gap on the two sides of the junctions, as a difference ∆ L ∆ R reduces contributions from the resonant transmission. The noise decreases as the temperature T L decreases, and saturates at T L < T R to a constant value that depends on T R . The result becomes independent of the phase ϕ at T L > T c , when the left terminal switches to the normal state. We emphasize the difference of S E between the opaque (left panel) and transparent junction limit (right panel). In particular, in the former case, the heat current noise is minimized for ϕ = 0 whereas in the latter junction S E is minimized for ϕ = π. This behavior resembles the one of the thermal conductance of a temperature-biased Josephson weak-link, as predicted in Refs. 4 and 5. Moreover, when normalized by conductance, the heat noise is smaller for transparent junctions, due to the last negative term in Eq. (19), which is second order in T n β . At equilibrium (T L = T R = T), the heat current noise obeys the fluctuation relation S E = 2G t h T 2 that connects it to the thermal conductance G t h = dI E /dT R | T L =T R . This relation can be used to define an effective noise temperature T * (T L ,T R , {τ n }, ϕ) such that S E = 2T 2 * G th (T * ), which can be useful for characterizing the behavior of the noise, as much of the variation with phase and transparency is contained in the thermal conductance. Such a temperature plot is shown in Fig. 5. The result is fairly insensitive to the values of τ n and ϕ, and the results fall nearly on the same FIG. 4. Heat current noise S E as a function of temperature T L , for different values of ϕ and τ n (single channel), normalized to its equilibrium normal-state value at T = T c . The temperature T R = T c /2 is kept fixed, and Γ = 0. Here, T c is the BCS critical temperature, and we take the temperature dependence of the energy gaps ∆ L (T L ), ∆ R (T R ) into account via the BCS relation, taking ∆ L (0) = ∆ R (0).
curves for different values of these parameters. At low temperatures, the result converges towards T * = max(T L ,T R ).
More generally, the generating function obeys the fluctuation relation 9 For the probability distribution of transferring energy ε into terminal R in time t 0 this implies This fluctuation relation is independent of the channel transmissions T n β , and therefore does not contain information about superconductivity.

A. Diffusive junctions
It is possible to average the above generating functions over known distributions of transmission eigenvalues {τ n }, to obtain results for certain types of multichannel junctions. Let us in particular consider the universal 31,34 transmission eigenvalue distribution of a short diffusive junction, 19 Changing the integration variable to q for β = + and q −1 for β = − in Eq. (11): (22) cancels the τ(q) dependent prefactors. The integral can then be evaluated: Note that the result does not depend on λ. The heat transport statistics of a diffusive junction is independent of the phase difference between the superconductors. Moreover, the only difference to the normal-state result is the presence of a gap max(|∆ L |, |∆ R |) in the energy integration. The absence of phase oscillations arises as coherent contributions from different channels cancel each other. The sign of the phase oscillation is different for large and small τ n , and diffusive junctions contain the exact balance of low and high-transparency channels necessary for the sum to cancel. The diffusive limit distribution is unique in the sense that log q n are uniformly distributed, 31 which is crucial for the above result.
Moreover, it should be noted that the phase dependence reappears in the diffusive limit if inelastic scattering is present (Γ > 0). The technical reason for this is that it is only in the limit Γ → 0 that superconductivity appears solely as a scale factor λ in the square transfer matrix eigenvalues q. We can illustrate the reappearance of phase oscillations via the normalized value M(E) =  n,± T n,± (E)/[2  n τ n ] that can be understood as a transparency for the heat current (17). The result is shown in Fig. 6. At sub-gap energies, the emergence of the diffusive junction minigap of size E g = |∆∥ cos ϕ 2 | is evident. Moreover, as pointed out above for single transparent channels, in contrast with low-transparency tunnel junctions, the relative sub-gap change in heat conductivity can be several orders of magnitude in transparent junctions (left panel).

B. Dirty interfaces
A second universal, potentially experimentally interesting eigenvalue distribution is the "dirty interface" distribution 18,35  Experimentally, this has been found to match results existing in high-transparency tunnel junctions. 36 The generating function can be evaluated also in this case, starting from Eq. (23): where z(E,u) = arccosh[1 + 2F(E,u)]. The associated heat transparency entering the heat current is illustrated in Fig. 7. We can note that the above-gap transport resembles that in tunnel junctions, and the sub-gap part that of diffusive junctions. However, in contrast to the tunnel junction result, 3,4 the energy integral in Eq. (26) is convergent, and does not require cutoffs in the above-gap resonance.

V. DISCUSSION
Electronic transport of heat in Josephson junctions is facilitated by transfer of quasiparticles from one side to the other. Andreev reflections transfer no heat, and so only quasiparticles contribute. Consequently, the heat transport statistics follows directly from the counting statistics of these excitations. Exactly as in normal-state junctions, 16,17 this is determined by the transmission eigenvalues of an appropriate scattering problem. In the superconducting state, this statistics is described by Bogoliubov-de Gennes transmission eigenvalues. We find their analytical expressions, Eqs. (13) and (16).
We considered both above-gap and sub-gap transport of heat, the latter by using a toy model for the broadening of the density of states in the superconductors. The general picture that emerges is that sub-gap heat transport is significantly more sensitive to the phase difference in transparent junctions than in tunnel junctions. Interestingly, in diffusive junctions it is in fact only the sub-gap transport that has any phase dependence at all.
The fluctuation statistics of heat transport can be measured using similar approaches as previously used for studying the heat currents. 6,7 For example, one can make the heat capacity of one of the terminals small, and then observe the fluctuation of its total energy via temperature measurements using established experimental techniques. 14 Setups utilizing several metal islands 10,13 equipped with fast NIS electron thermometry 37 could also be a viable approach for probing the statistics experimentally.
In summary, we consider heat transport driven by a temperature difference across Josephson junctions of varying transparency. We obtain the generating function of fluctuation statistics in closed form. In addition to describing the fluctuations in heat transport, the results provide a way to understand the average heat current in terms of elementary transmission events. Open questions remain in taking into account the effects of voltage, temperature, and order parameter fluctuations in the superconductors, and modeling inelastic effects starting from microscopic rather than phenomenological models.