Non-equilibrium quasiparticles in superconducting circuits: photons vs. phonons

We study the effect of non-equilibrium quasiparticles on the operation of a superconducting device (a qubit or a resonator), including heating of the quasiparticles by the device operation. Focusing on the competition between heating via low-frequency photon absorption and cooling via photon and phonon emission, we obtain a remarkably simple non-thermal stationary solution of the kinetic equation for the quasiparticle distribution function. We estimate the influence of quasiparticles on relaxation and excitation rates for transmon qubits, and relate our findings to recent experiments.


Introduction
Superconducting devices have been under development for several years for applications in detection [1], nonlinear microwave amplification [2][3][4], quantum information processing [5], and quantum metrology [6,7]. Intrinsic excitations in the superconductor, known as Bogoliubov quasiparticles, can be detrimental to such devices, for example limiting the sensitivity of detectors and causing decoherence in qubits. The devices are routinely operated at temperatures so low that no quasiparticles should be present in thermal equilibrium. However, a significant number of residual non-equilibrium quasiparticles is typically detected, and it is well established that their density (normalized to the Cooper pair density) can be x qp ∼ 10 −5 -10 −8 [8][9][10]. Much less is known about the energy distribution of these quasiparticles.
Many experiments involving residual quasiparticles in qubits [11][12][13] are successfully described by the theory [14,15], based on the assumption of a fixed average quasiparticle distribution which perturbs the device operation; the resulting net effect of the quasiparticles is equivalent to that of a frequency-dependent resistance included in the circuit. The fixed distribution assumption is valid in the weak signal regime, when the back-action of the superconducting condensate excitations (resonator photons or qubit excitations, hereafter all referred to as photons for the sake of brevity) on the quasiparticles can be neglected, or for studying observables which are not sensitive to the quasiparticle energy distribution function, but only to the quasiparticle density. This assumption must be reconsidered when the probing signal is strong enough to modify the quasiparticle distribution and the latter can affect the quantities which are measured.
Some efforts in this direction have been made. In the numerical work [16], the external circuit was treated classically, so it tended to heat the quasiparticles to infinite temperature, and the distribution was stabilized by phonon emission only (see also [17] and references therein for related experiments). In [18], the quasiparticle distribution was assumed to be entirely determined by the external circuit (both heating by the drive and cooling by photon emission were included at the quantum level), while phonon emission was neglected.
Here we study analytically the competition between quasiparticle heating via absorption of low-energy photons (i. e., with energy ω 0 2∆, twice the superconducting gap; we use units with ħ h = 1, k B = 1 throughout the paper) and cooling via both photon emission in the external circuit and phonon emission in the material, as schematically shown in Fig. 1. We obtain a nonthermal stationary solution of the kinetic equation for the quasiparticle distribution function. Using our solution we estimate the influence of non-equilibrium quasiparticles on relaxation and excitation rates for transmon qubits, and relate our findings to recent experiments.
The paper is organized as follows: we first introduce the kinetic equation and discuss the main properties of its solution for different regimes (conditions for applicability of the kinetic equation and detailed derivation of the solution are given in appendices). We then discuss our results in light of recent experiments with superconducting qubits.

Kinetic equation
The quasiparticle distribution function f (ε) satisfies the kinetic equation that we write as where we measure the energy ε from the superconducting gap and assume 0 < ε ∆. The two terms on the right-hand side represent collision integrals due to absorption/emission of  Figure 1: A schematic representation of the model under consideration. The quasiparticle distribution function f (ε) (shaded areas) is affected by photon exchange with e.g. a qubit and by the emission of phonons. As explained in the text, for "cold" quasiparticles the width T * * of the distribution function [Eq. (11)] is small compared to the photon frequency ω 0 , and the distribution function has a second, smaller peak at that frequency (replicas at higher multiples of ω 0 can be neglected in the cold quasiparticle limit).
phoTons and phoNons, respectively. The photon term can be written as [18] Here the square roots approximate the density of states near the gap edge, ε ∆; for ε < ω 0 the first term should be set to zero.n is the average number of excitations in the superconducting subsystem with transition frequency ω 0 ; for a qubit 0 ≤n < 1, while for a harmonic system such as a resonator anyn ≥ 0 is allowed; this is in contrast to the classical approach of Ref. [16], valid only forn 1. The effective temperature T 0 is defined by the relation e ω 0 /T 0 = (1 ∓n)/n, with the upper (lower) sign for a two-level qubit (harmonic) subsystem. Γ 0 characterizes the rate of photon absorption/emission by the quasiparticle. For a weaklyanharmonic, single-junction qubit such as the transmon it is given by with δ being the mean level spacing of the superconducting islands. For a resonator, one should use the mean level spacing in the whole resonator volume, occupied by the quasiparticles, while 2π∆ should be replaced by Qω 0 with Q being the resonator quality factor if the material resistivity were the same as in the normal state. We do not address photon dynamics here, assuming the photon state to be stationary and fixed by the external circuit. In principle, one can write coupled equations for the photon density matrix and f (ε), as in Ref. [18]; their solution in the presence of phonons remains an open question. Also, we neglect multiphoton absorption/emission and assume the photon system to be either strictly two-level or strictly harmonic; then each act of absorbtion/emission involves only energy ω 0 , and the photon state enters only vian. Going beyond these assumptions would allow quasiparticles to absorb/emit multiples of ω 0 or energies slightly different from ω 0 due to weak anharmonicity in the resonator or broadening of the qubit transition (we remind that the ratio between anharmonicity and broadening determines if the photon system should be treated as two-level or harmonic [18]).
We assume the phonons to be kept at constant low temperature T ph ω 0 . Then, the quasiparticle-phonon scattering collision integral can be written as [19][20][21]: where we assume the quasiparticles to be non-degenerate, f (ε) 1, so all terms quadratic in f (ε) (in particular, quasiparticle recombination) are neglected. This approximation is discussed in detail in Appendix A. The function F (ω) is defined as Here ζ(x) is the Riemann zeta function, Ξ is the deformation potential, v F and v s are the Fermi velocity and the speed of sound, ρ 0 is the mass density of the material, ν n is the density of states at the Fermi level for the material in the normal state, taken per unit volume and for both spin projections. These material parameters are conveniently wrapped into the coefficient Σ, which controls energy exchange between electrons and phonons for the material in the normal state: the power per unit volume transferred from electrons to phonons, kept at temperatures T e and T ph , respectively, is given by Σ(T 5 e − T 5 ph ) [22]. This relationship as well as Eq. (5) are appropriate for a clean metal, when the electron elastic mean free path due to static impurities is longer than the mean free path due to the electron-phonon scattering. In the opposite limit, F (ω) is proportional to a different power of ω: F (ω) ∝ ω for impurities with fixed positions [23] and F (ω) ∝ ω 3 for impurities which move together with the phonon lattice deformation [24,25]. Although at low energies, discussed here, the system is expected to be in the diffusive limit, the clean-limit expressions usually agree better with experiments (see Ref. [26] for a review). Thus, here we use the clean-limit formula; the diffusive limit is discussed in Appendix D where we show that the results are qualitatively similar.
We assume the phonon temperature to be very low, much smaller than the typical quasiparticle energy. Then the last term in Eq. (4) determines the quasiparticle relaxation rate by phonon emission 1 : When we neglect the absorption of phonons (i.e., setting T ph = 0), the phonon collision integral simplifies: Some consequences of T ph > 0 are explored in Appendix C. From now on, with kinetic equation we will mean Eq. (1) with the right hand side given by the sum of this simplified expression plus Eq. (2). In the steady state, ∂ f /∂ t = 0, the kinetic equation reduces to an integral equation for f (ε). In the next section we study the solution of this equation and we identify two main regimes (see Fig. 2). 1 To avoid confusion, we note that the symbol τ ph used here corresponds to τ ph (∆) of Ref. [18]. It is related to Figure 2: A schematic representation of the quasiparticle distribution function f (ε) (green shaded areas) in the two regimes, cold and hot (left and right panels), described in Secs. 3.1 and 3.2, respectively. For cold quasiparticles, most population is concentrated at energies ε T * * ω 0 ; above T * * , the distribution function decays as f (ε) ∼ ε −4 ; f (ε) has a second, smaller peak at ε = ω 0 (replicas at higher multiples of ω 0 can be neglected). In the hot quasiparticle regime, f (ε) has a sawtooth dependence with the period ω 0 , on top of a smooth envelope (black dashed curve), characterized by a scale ε ∼ T * ω 0 . (Note the different scales for ω 0 in the two panels.)

Cold quasiparticle regime
The first regime occurs when there are few excitations to absorb and/or the electron-phonon relaxation is sufficiently strong, so the quasiparticles are more likely found at low energies, ε ω 0 . From time to time, a quasiparticle absorbs a quantum ω 0 and is promoted to energies above ω 0 , and then quickly relaxes back to lower energies by emitting either a quantum ω 0 , or a phonon. Thus, we can find the distribution at energies just above ω 0 by perturbation theory from the stationary kinetic equation. At those energies, we can neglect the first term in the right hand side of Eq. (7) because the ε integral is dominated by ε − ε ω 0 , while in the second term ε ≈ ω 0 (the smallness of the first term can be checked self-consistently after the solution is found). Then, neglecting also the occupation at energies larger than 2ω 0 , we obtain: where we introduced a dimensionless parameter The condition e ω 0 /T 0 1 or Λ 1 ensures that the perturbative approach is valid and thus defines the cold quasiparticle regime. 2 2 One may worry that if the condition T 0 ω 0 is not satisfied, f (ε + ω 0 ) ≈ f (ε) for small ε, so the perturbation theory does not work. Note, however, that if we try to find f (ε+2ω 0 ) from the stationary kinetic equation focusing Substituting f (ε + ω 0 ) from Eq. (8) into the kinetic equation at low energies ε ω 0 , we obtain On the right-hand side, the first term in the brackets represents the out-scattering to the higher levels at ε > ω 0 (minus the photon emission contribution), while the second term is due to phonon emission, when the quasiparticle is transferred to lower energies. This second term dominates at energies and the cold quasiparticle condition can be expressed as T * * ω 0 . The last term on the lefthand side of Eq. (10) represents the uniform incoming flux from the energies between ω 0 and 2ω 0 (higher energies are neglected, as discussed above). This term is smaller than the first one (incoming flux from energies below ω 0 ) as long as ε ω 0 (see Appendix B). In both integrals, we can push the upper integration limit to infinity if the integrand decreases quickly enough. Then, we are left with the equation which has a remarkably simple exact solution f (ε) = C/ε 4 with an arbitrary constant C. This power-law form justifies the change in integration limit and is valid for energies T * * ε ω 0 ; it resembles the Kolmogorov spectrum of turbulence [27], and describes a similar physical situation: the flow of probability, which is injected at high energies, ε ω 0 , and flows to lower energies, down to ε ∼ T * * where it encounters an effective sink, represented by the first term on the right-hand side of Eq. (10). In our case, the probability is reinjected back at higher energies by absorbing a quantum ω 0 .
The behavior of the solution at energies ε T * * depends on the relation between the two dimensionless parameters, e ω 0 /T 0 and Λ. While we are unable to find an analytical solution in this region, it is possible to establish the qualitative character of the solution. Let us drop the second term on the left-hand side of Eq. (10) and introduce new variables where r = Λ 6 /e 7ω 0 /T 0 . This gives an equation for F r (x): on energies close to 2ω 0 and neglecting f (ε + 3ω 0 ) as well as the first term in Eq. (7), we obtain so having just one of the conditions T 0 ω 0 or T * ω 0 suffices to neglect f (ε + 2ω 0 ) and higher.
The 1/ε 4 asymptotics of f (ε) translates into F r (x) ∼ 1/ x at x x * * ≈ min{r 2/7 , r 1/3 }, so the integral converges at the upper limit. F r (0) is finite, since at low energies the integral is also well-behaved. The contribution of y x * * is suppressed by the numerator, so F r (0) and F r (x * * ) are of the same order. This qualitative analysis is supported by numerical findings which show that to good accuracy the function F r (x), upon rescaling x by x * * , takes a form independent of r (Appendix B).
Note that because of the square root divergence of f (ε) as ε → 0, the normalization coefficient depends weakly (logarithmically) on the phonon temperature when T ph T * * . On the other hand, Eq. (8) and the power-law decay at large energy remain valid even when T ph > T * * , although the latter starts from an energy large compared to T ph , as discussed in Appendix C.

Hot quasiparticle regime
Ifn (and hence T 0 ) is sufficiently large, the width T * * of the distribution function may exceed ω 0 . When both ω 0 /T 0 , Λ 1, the quasiparticles can absorb many excitations and their typical energy becomes large compared to ω 0 , which defines the hot quasiparticle regime. Still, this does not mean that f (ε) becomes smooth on the scale ε ∼ ω 0 : the singularity in the density of states at ε → 0 is imprinted in f (ε → ω 0 ), f (ε → 2ω 0 ), etc., by photon absorption, giving rise to a series of peaks in f (ε) at integer multiples of ω 0 (cf. Fig. 2), which were observed in the numerical results of Ref. [16]. Below we focus on large-scale features of the distribution function, its fine structure lying beyond the scope of the present paper. Then, if f (ε) is understood as the smooth envelope, we can approximate the photon collision integral in Eq. (2) by a diffusion operator: whose form ensures the correct steady-state solution f ∼ e −ε/T 0 when phonon emission is neglected [18]. In this approximation, it is convenient to introduce the temperature scale For the quasiparticle to be hot, we need T 0 , T * ω 0 . Let us first assume ω 0 T * T 0 (more precisely, we are considering the limits ω 0 → 0, T 0 → ∞, while keeping T * = const). Then we can replace the exponential terms in the curly brackets of Eq. (15) with unity and introduce dimensionless variables x = ε/T * and y = ε /T * . As a result, the steady-state equation acquires a parameter-free form: We can find approximate solutions to Eq. (17) in the low (x 1) and high (x 1) energy limits; the solutions will contain unknown coefficients, which we will fix by solving the equation numerically. In the high-energy regime, we can drop the last term, and with the change of variable x = 4 1/6 z we obtain the Airy equation: Thus, the high-energy part of the solution is expressed in terms of the Airy function: with some coefficient f ∞ . At low energies x → 0, in Eq. (17) we can drop the term proportional to x 7/2 and set the lower integration limit to zero. We can then find a solution in the form of an expansion: with a = 21 32 and some coefficient f 0 , whose relation to f ∞ is in principle determined by matching the solutions at x ∼ 1; here we find the relation accurately by comparing to numerical result. We solve Eq. (17) numerically as follows: we discretize variables x and y, so that the equation becomes a matrix equation and we find the vector with the smallest eigenvalue in absolute value. We repeat this process on finer grids and then extrapolate to zero spacing; in the extrapolation, the smallest eigenvalue tends to zero. By fitting the numerical solution with Eq. (19) we find f ∞ / f 0 3.00, while using the definitions in Eq. (21), numerical integration gives a 0.564 and b 0.119. With these values both the high-and low-energy formulas, Eqs. (19) and (20), fit well the numerical solution, see Fig 3. Finally, to find the proper normalization, we can use the numerical result Let us briefly discuss the opposite limit, T * T 0 ω 0 . In this regime we can use Eq. (15), but not the simpler form that led to the first term in the right hand side of Eq. (17). In this limit, the bulk of the distribution is f (ε) ∝ e −ε/T 0 , and only at very high energies, ε T ∞ ≡ T 3 * /(2T 0 ), the tail is suppressed by the phonon emission: We see that when T 0 T * , the Boltzmann-like exponential decay f (ε) ∼ e −ε/T 0 is valid up to the very high energy T * T * /2T 0 , above which the distribution function decays faster than exponentially. On the other hand, when T 0 T * , the quasiparticle distribution function is governed by T * rather than T 0 : it changes little for ε T * and decays faster than exponentially for ε > T * .

Discussion: relevance to experiments
The above results for the quasiparticle distribution function are valid both for a harmonic and a two-level (qubit) superconducting system, if the appropriate relationship betweenn and T 0 is used. We focus henceforth on the qubit case to investigate to what extent a qubit can heat the quasiparticles. Since T 0 /ω 0 = 1/ ln (1/n − 1), the condition T 0 > ω 0 translates intō n > 1/(1 + e) ≈ 0.27. Noting that the excited state occupationn of an undriven qubit is generally below 0.3 (in most cases being between a few percent and not much above 10%, with some qubits being much colder,n ∼ 0.1% [28]), we conclude that an undriven qubit cannot make the quasiparticles hot. For a strongly driven qubit, we can haven → 1/2 − and hence T 0 → ∞ [18]. In this limit, we should check if T * ω 0 . The largest values of T * /ω 0 are reached for low-frequency qubits made of small islands, since the rate Γ 0 is inversely proportional to the volume. Considering a typical aluminum island of volume ∼ 0.02 µm 3 , we have Γ 0 ∼ 10 5 Hz [18]. For the time τ ph , we use the experimental results for thin films in Refs. [29] and [30] to estimate τ ph ∼ 10 ns. Then for typical qubit frequencies ω 0 = 4 to 8 GHz, we find that T * /ω 0 varies between 1.2 and 0.8. Hence we conclude that even under strong driving a qubit cannot significantly heat the quasiparticles to energies much above its frequency. This conclusion is not too sensitive to the used parameters: for instance, increasing τ ph by two orders of magnitude only doubles the calculated T * because of the very weak power Λ −1/6 in Eq. (16).
We can extend the considerations above to include the possibility that the quasiparticles are not heated directly by the qubit, but indirectly due to the interaction of the qubit with another mode. This mode could be a resonator mode, or a spurious (harmonic) mode of the chip. For instance, take a resonator, coupled to the qubit with coupling strength g; typically, the coupling strength is of order 100 MHz and the resonator-qubit detuning δω is of order 1 GHz. This implies that the quasiparticle rate Γ 0 is suppressed by about two orders of magnitude, Γ 0 → Γ 0 (g/δω) 2 . Therefore, even for a small qubit, we do not expect significant heating unless the mode is populated with hundreds of photons.
Let us now consider an undriven qubit; as discussed above, for the experimentally observed range ofn the quasiparticles are cold. However, for small-island qubits the distribution function width T * * is given by ω 0 e −ω 0 /(3T 0 ) , while for larger qubits (3D transmon [9,11], Xmon [31]) with electrode volume of order 10 3 -10 4 µm 3 , it is given by ω 0 Λ −2/7 , see Eq. (11). The different regimes for T * * lead to different behaviors for the quasiparticle-induced excitation rate. In both cases, since T * * ω 0 the relaxation rate is approximately given by [15] where is the quasiparticle volume density n qp normalized by the Cooper pair density. For the excitation rate, we use Eq. (8) to write This ratio has the detailed balance form; note, however, that it is determined by the qubit occupation rather than the energy scale of the quasiparticles, which are not in thermal equilibrium. For large qubits, e −ω 0 /(3T 0 ) > Λ −2/7 , we cannot neglect the term proportional to Λ in the denominator, and therefore we conclude Γ qp 01 /Γ qp 10 n/(1−n). We should point out that for large qubits we estimate T * * T ph : not surprisingly, in a large device the quasiparticles should be close to thermal equilibrium with the phonon bath. However, the found inequality for the ratio of quasiparticle rates still holds, since it is based on the use of Eq. (8). This inequality, together with Eq. (24), validates the use of the density x qp as the relevant dynamical variable affecting qubit relaxation. The dynamics of x qp in the presence of traps has been studied in recent works [32,33], where possible ways to improve qubits performance are analyzed.
Some recent experiments [12,34] with large qubits report that at low temperatures the main relaxation mechanism is a "residual" (non-quasiparticle) one, Γ r 10 Γ qp 10 . Then, if we assume that the qubit is the main source of quasiparticle non-equilibrium, the above considerations imply Γ Γ r 01 , whereas experiments indicate that the opposite inequality holds [28,34]. This discrepancy could indicate that one should account for a (currently unknown) energy-dependent source of pair breaking photons and/or phonons generating hot quasiparticles; in other words, the source of the quasiparticles should be identified, possibly via experiments with resonators [35]. Alternatively, a more detailed modeling of the systems under investigation might be necessary, for example by including details about the geometry of the superconducting electrodes or properties of the substrate and interfaces that have been neglected so far.

Conclusion
In summary, we have studied the distribution function for quasiparticles interacting with photons of frequency ω 0 and a low-temperature phonon bath. We have identified the dimensionless parameters that control whether the quasiparticles are cold (i.e., mostly occupying states with energy below ω 0 ) or hot, see Eqs. (9) and (16) and the text following them; those parameters are determined by the average number of photon excitationsn, the product Γ 0 τ ph characterizing the relative strengths of quasiparticle coupling to photons and phonons, and the ratio between photon frequency and superconducting gap ∆. Applying our results to transmon qubits, we conclude that the qubit itself cannot overheat the quasiparticles. The extension of our findings to other systems, such as resonators and multi-junction circuits (arrays [36,37], fluxonium qubit [13]), will be presented elsewhere.

A Neglected terms in the kinetic equation
Besides quasiparticle-phonon scattering, phonon absorption can result in production of a pair of quasiparticles, and a pair of quasiparticles can recombine by emitting a phonon. In the absence of extrinsic processes, the balance between these two processes determines the quasiparticle density (that is, the normalization of the distribution function), since the scattering processes considered in the main text redistribute the quasiparticle energy but do not change their number. However, if the generation/recombination rates are small compared to the rates of the scattering processes, the latter determine the shape of the distribution function. The generation rate is proportional to e −2∆/T ph and is therefore negligibly small. To estimate the importance of recombination, we note that it can be included in the kinetic equation (1) by adding a term where the recombination rate can be expressed using Eq. (25): The combination 2n qp /ν n has the meaning of the normal-state level spacing in a volume occupied by one quasiparticle. For quasiparticles excited by photons to energy of order ω 0 , we can consider the condition Γ rec Γ ph (ω 0 ), which translates into For typical values ω 0 5 GHz, ∆ ∼ 50 GHz, and x qp 10 −5 , this condition is satisfied. For quasiparticles near the gap edge, the relevant process is that of photon absorption with ratē nΓ 0 , see the last term in Eq. (2); then the conditionnΓ 0 Γ rec sets a lower bound onn. For qubits with small aluminum islands, assuming at most a few quasiparticle in the islands we find x qp ∼ 10 −5 and hencen > 10 −2 , which is typically satisfied, as discussed in the main text; of course if there is only one quasiparticle in an island, recombination cannot take place, and there is evidence that in systems with several islands the average number of quasiparticles in each island is less than one [10] or it can be suppressed to less than one [38]. For large volume qubits, even assuming x qp ∼ 10 −8 , we would find the requirementn > 1, which cannot be met. On one hand, this means that in large qubits we cannot neglect the effect of generation/recombination in determining the energy dependence of the distribution function; on the other hand, since Eq. (30) is satisfied, Eq. (8) is still valid and our considerations about transition rates are unaffected.

B Solution to Eq. (14)
As discussed in the main text, at x x * * ≈ min{r 2/7 , r 1/3 }, the asymptotic solution to Eq. (14) takes the form F r (x) ∼ 1/ x for any r. That equation can be solved numerically by iteration; that is, starting with a trial function with the correct asymptotic behavior, we calculate the (n + 1)th iteration by inserting the nth iteration into the right hand side of Eq. (14). Numerically, the integral is evaluated by splitting it into a "low" energy region, extending from x to about 10x * * , and a high energy one; the contribution of the low-energy region is then obtained  ) and (32). The dashed line is the analytic prediction given by the next-to-leading asymptotic terms [see text after Eq. (32)]; the agreement between the two curves at large x further validates our numerical procedure. by discretizing the integral on an equally-spaced grid with steps of order 0.02x * * , while the high-energy part is estimated analytically by using the leading asymptotic form of the solution. The calculation is stopped when the desired accuracy is reached; typically, the maximum relative change in the numerical solution becomes less than 10 −6 after a small number ( 8) of iterations. We have checked that the solution thus found is only weakly sensitive (relative deviations at most of order 10 −3 ) to e.g. doubling the value of the splitting point of the integral or the resolution. The exact form of the initial trial function is unimportant: as long as we set it proportional to 1/ x for x > x * * , it can be represented for x < x * * by an array of random numbers between 0 and 1. Interestingly, as we show next, the dependence of F r on r is, to a good degree, accounted for with an appropriate rescaling.
Let us consider the limiting cases r → 0 and r → ∞; by further rescaling variables by r 1/3 and r 2/7 , respectively, we find the equations and Their asymptotic solutions at large x are, up to an overall coefficient, At arbitrary x, we solve these equation numerically, as explained at the beginning of this section, and normalize the results so that F 0 (x)/F ∞ (x) → 1 as x → ∞. Then the ratio between the two functions deviates from unity by less than 1 % at all x, as shown in Fig. 4. This result suggests that it should be possible to express F r for any r in terms of a single function with good accuracy. Indeed, let us define the average functionF (x) = [F 0 (x) + F ∞ (x)]/2, and normalize it so thatF (0) = 1; our numerical result for this function is shown in Fig. 5 and it is accurately fit (within 0.1 %) by the Padé-like expression For arbitrary value of r, we then find (within ∼ 1 %) with α(r) ≈ 1/r 1/3 + 1/r 2/7 − 0.7345/r 1/6+1/7 .
Here the power of the last term in the right hand side is arbitrarily set as the average between the powers of the first two, asymptotic terms, and the numerical factor is obtained by comparison with numerical solutions in the range r from 10 −3 to 10 3 . A more precise definition of the energy scale T * * of Eq. (11) can be given as The above considerations were based on neglecting the second term in Eq. (10) in comparison to the first. To check this assumption, we note that using Eq. (12), at ε T * * using the definitions in Eq. (13) we estimate the first term to be (e ω 0 /T 0 F r (0)/r α(r)Λ 2 ) ω 0 /ε. For the second term, using again those definitions and introducing the change of variable x = t 2 /α(r), we find the approximate upper bound 2(e ω 0 /T 0 F r (0)/r α(r)Λ 2 ). Therefore we can indeed neglect the second term when ε ω 0 .

C Effect of finite phonon temperature
Most of the arguments given in this paper lead to finite results when the phonon temperature T ph → 0, so that only phonon emission is allowed. However, a problem arises when one tries to normalize the distribution function in Eq. (13) according to Eq. (25): the integral diverges logarithmically at low energies. This divergence should be cut off at ε ∼ T ph T * * : where the dimensionless slope β depends in general on the electronic mean free path and ω D is the Debye frequency. Assuming again low phonon temperature, we find for the relaxation rate Γ ph (ε) = 16π 5 where the last expression redefines the characteristic electron-phonon time τ ph in the diffusive regime [cf. Eq (6)]. The kinetic equation (neglecting phonon absorption) has now the form Proceeding as in the main text, we find that for cold quasiparticles Eq. (10) becomes where Λ is now defined as The equation can be simplified for energies [cf. Eq. (11)] in which case it reduces to The solution to this equation is f (ε) = C/ε 3 : although the power of this tail is different from that found for the clean metal case, the qualitative behavior of the distribution function is the same. Next, we show that the qualitative similarities between clean and diffusive cases persist also for hot quasiparticles. For hot quasiparticles we approximate again the photon collision integral via a diffusion operator, Eq. (15), and the kinetic equation becomes The temperature scale T * should be now defined as T * = ω 0 Λ −1/5 , so for T * T 0 , using dimensionless variables, we obtain the steady-state equation as At high energies, x 1 (but still x T 0 /T * ), we can neglect the last term, and with the change of variable x = 4 1/5 z we arrive at the generalized Airy equation whose solution can be written in terms of the modified Bessel function of the second kind to give with asymptotic behavior Note that proper normalization can be found if b is known. As in Sec 3.2, using the numerical solution to Eq. (49) in these definitions we find a ≈ 0.468 and a ≈ 0.121. Fitting the numerical solution we also obtain f ∞ ≈ 0.53 f 0 . The numerical solution and the analytical approximations are plotted in Fig. 6.
In the regime T * T 0 , we find again that the distribution function takes the Boltzmann form up to high energies: where T ∞ = 2 1/3 T * (T * /2T 0 ) 2/3 .
In the hot regime, the generalized Airy equation obtained at high energies after the substitution whose solution is expressed in terms of the modified Bessel function K 2/7 .