Nonequilibrium quasiparticle distribution in superconducting resonators: effect of pair-breaking photons

Many superconducting devices rely on the finite gap in the excitation spectrum of a superconductor: thanks to this gap, at temperatures much smaller than the critical one the number of excitations (quasiparticles) that can impact the device's behavior is exponentially small. Nevertheless, experiments at low temperature usually find a finite, non-negligible density of quasiparticles whose origin has been attributed to various non-equilibrium phenomena. Here, we investigate the role of photons with energy exceeding the pair-breaking threshold $2\Delta$ as a possible source for these quasiparticles in superconducting resonators. Modeling the interacting system of quasiparticles, phonons, sub-gap and pair-breaking photons using a kinetic equation approach, we find analytical expressions for the quasiparticles' density and their energy distribution. Applying our theory to measurements of quality factor as function of temperature and for various readout powers, we find they could be explained by assuming a small number of photons above the pair-breaking threshold. We also show that frequency shift data can give evidence of quasiparticle heating.


I. INTRODUCTION
A detector should ideally respond in a well-understood way to the physical quantity under investigation while being unaffected by other processes.However, unwanted effects can contribute to noise that obscures the signal.The noise mechanism can be specific to the type of detector: even within superconducting detectors, different mechanisms are prominent for different designs.For example, in superconducting nanowire single photon detectors [1,2] vortices crossing the nanowire lead to dark counts [3,4], and in kinetic inductance detectors (KIDs) [5] the generation-recombination noise due to the creation and annihilation of quasiparticles increases the noise equivalent power [6][7][8].A KID consists of a resonator whose inductance and hence resonant frequency changes when photons of energy above twice the superconducting gap ∆ break Cooper pairs.The response of the resonator is monitored via a probe tone that maintains a large number n ≫ 1 of photons with frequency ω 0 < 2∆ in the resonator (hereinafter we set ℏ = k B = 1).The device is operated at temperatures much smaller than the critical one, T ≪ T c , where the thermal equilibrium number of quasiparticles is expected to be negligible.However, a variety of experiments with not only resonators [9,10] but also superconducting qubits [11,12] indicates the presence of many quasiparticles whose origin is often unclear.These excess quasiparticles influence basic parameters of the resonator such as the quality factor, which is why the goal of this work is to quantify their effect and propose a possible source.
To model the non-equilibrium state of a superconductor, a system of coupled kinetic equations for the distribution functions of quasiparticles and phonons has been proposed long ago [13]; even in the steady-state and for uniform systems -that is, considering only the energy dependence of the distribution functions, -these equations are usually solved numerically [14,15].Only recently, approximate analytical solutions have been obtained in the parameter regime relevant to KIDs [16,17].The results of these two works can be summarized as follows: a large number of (non-pair-breaking) photons can "heat up" the quasiparticles by pushing them to higher energy as compared to the phonon temperature; these quasiparticles can relax by emitting phonons, so that the latter are also driven out of equilibrium.Quasiparticles pushed to energy above 3∆ can emit a relatively large number(as compared to the thermal equilibrium value of pair-breaking phonons with frequency ω > 2∆; because of these phonons, in the steady state the quasiparticle density is then much larger than the equilibrium value, a situation that is encountered at low temperatures.These findings can quantitatively describe the experimental data for a resonator's internal quality factor Q i of Ref. [18] at intermediate temperatures at which the number of quasiparticles is close to the equilibrium value but their distribution function is not thermal.As temperature decreases, Q i is predicted to saturate, but at values few to several orders of magnitude larger than those measured. The decrease of Q i with readout power observed at low temperature in Ref. [18] is incompatible with the expectation from losses due to two-level systems [19].On the other hand, clear evidence for the presence of pair-breaking photons is provided by measurements of socalled parity switching rates in superconducting qubits: effects initially attributed to "hot" quasiparticles [20] have been explained in terms of tunneling assisted by pair-breaking photons [21], as confirmed by additional experiments [22,23].Motivated in part by these results, here we extend the model [16,17] discussed in the previ-FIG.1.A full model of a thin-film superconducting resonator (dark box) takes into account the resonant mode of frequency ω0, quasiparticles, phonons, and their interactions (with coupling constants c QP P hot and τ0).Additionally, quasiparticles can be generated by pair-breaking photons of frequency ωP B (coupling constant c QP P hot,P B ; this process was not considered in Ref. [17]) and the phonons interact with the substrate, which act as a thermal bath (coupling constant τ l ).
ous paragraph to include the effect of a small number of pair-breaking photons of frequency ω P B > 2∆ and show that they could be responsible for the low-temperature behavior of Q i reported in Ref. [18].
In Sec.II we present the kinetic equation that determines the quasiparticle distribution; it extends the previously used kinetic equations describing the interaction of quasiparticles with photons of energy below the pair breaking threshold 2∆ [17,24,25] by including a contribution from a mode of energy above the threshold.In Sec.III we derive approximate analytical solutions for the case of zero temperature and negligible number of photons below the threshold, and validate the results numerically.The effect of the low-energy photons on the distribution's shape is investigated in Sec.IV.In Sec.V the results of the preceding section are used to calculate quality factor and resonance frequency shift in thin-film resonators, and we analyse experimental data [18] for this quantities.Section VI summarizes our findings.

II. KINETIC EQUATION
The kinetic equation for the quasiparticle distribution function f (E) in a homogeneous superconductor has the form df (E) dt = St P hon {f, n} + St P hot {f, n} + St P hot P B {f, nP B }.
(1) with E the energy measured from the Fermi level.The collision integrals in the right-hand side account for the interaction between quasiparticles and phonons, St P hon {f, n}, non-pair-breaking photons, St P hot {f, n}, and pair-breaking ones, St P hot P B {f, nP B }, respectively.While this structure of the kinetic equation is quite general, we will mostly consider thin-film superconducting resonators, as depicted in Fig. 1.The phonon collision integral St P hon {f, n} and the photon one for non-pair breaking photons St P hot {f, n} can be found for instance in Ref. [17].In this work we extend the model by including photons of energy ω P B > 2∆ via the collision integral: Here, St P hot {f, nP B } is a number-conserving scattering term accounting for the redistribution of quasiparticles in energy due to the absorption of pair-breaking photons.It can be obtained from the photon integral St P hot for nonpair-breaking photons of frequency ω 0 [see diagrams a) and b) in Fig. 2], by replacing photon energy ω 0 , photon number n and coupling constant c QP P hot by the respective values of the pair breaking photons ω P B , nP B and c QP P hot,P B [the appearance of n or nP B in the argument of teh collision integral indicates which of these quantities should be used].We use here the notation of Ref. [17] by defining Ref. [26] for a discussion of coherence factors).This collision integral conserves the number of quasiparticles, ∆ dE ρ(E)St P hot = 0. We note that both for non-and pair-breaking photons we consider a single mode of definite frequency; this is justified for non-pair-breaking photons by the fact that, as mentioned in the Introduction, in applications such as KIDs one high-quality factor mode of the resonator is probed.For pair-breaking photons this is in general a simplification; however, at least in some of the regimes we will consider, the extension to multiple modes is straightforward (see Secs.III and IV).Moreover, we will show that this simplification does not alter qualitatively our interpretation of experimental data for the quality factor in Sec.V.
In addition to the scattering term, there is a term accounting for the generation of new quasiparticles [diagram d) in Fig. 2], and a term describing recombination accompanied by the emission of a photon, In Appendix A we discuss how to estimate the coupling constants c QP P hot and c QP P hot,P B .The phonon collision integral can be obtained from the pair-breaking  3)], c) the emission of a phonon [c.f.first term in the bracket of Eq. ( 6)], d) quasiparticle generation by pair-breaking photons [Eq.( 4)], a process not included in Ref. [17].Diagram e) depicts generation by a pair-breaking phonon [last term of Eq. ( 6)], f) provides the bare recombination rate r0xqp, and g) renormalizes the bare recombination coefficient r0 to r in Eq. (6) [see also Eq. ( 8)].Diagram h) describes an additional pair-breaking mechanism due to phonon emitted by quasiparticles with energy E > 3∆; this process is not included in the low-energy approximation of Eq. ( 6), but if n is sufficiently large (cf.Sec.IV) it can affect the density, see the term proportional to G(T * /∆) in Eq. ( 49) and Ref. [17].In Sec.. III we only consider c), d) and f), in Sec.IV we include diagrams a) and b), and in Secs.IV C and V we consider all diagrams.We do not include the diagram corresponding to photon mediated recombination [similar in structure to f)], as its contribution is generally negligible compared to phonon mediated recombination (see Sec. III A) photon one, Eq (2), by replacing , and integrating over ω > 0; the integration limits are chosen so that the second argument of U ± is larger than ∆.In these replacements, n(ω) is the phonon distribution function and τ 0 is the time scale characterizing the strength of the electronphonon interaction [27].
To complete the description of the system, one should also consider the kinetic equation for the phonon distribution function n(ω); we don't give it here as it will not be needed explicitly, and rather refer the reader to our previous work [17].That equation contains a thermalization term −[n(ω) − n T (ω, T B )]/τ l accounting for the relaxation of the phonons back to equilibrium at temperature T B over the time scale τ l .For the numerical calculations in this work, we either assume fast equilibration to zero temperature, n(ω) = 0, or solve the full coupled system of kinetic equations for quasiparticle and phonons, as will be specified.The numerical solutions are obtained using a straightforward extension of the algorithm described in Ref. [17].
In the next sections, we discuss approximate analytical solutions to Eq. ( 1) in the steady state df /dt = 0.While different approximations will be employed in different regimes, some approximations are common to all cases we will consider: we generally assume small occupation probability at all energies, f ≪ 1, so that we can replace Pauli blocking factors with unity, (1 − f ) → 1.The validity of this assumption, which limits our considerations to temperatures small compared to the critical one and small number of pair-breaking photons, can be checked once the solution to the kinetic equation is found; importantly, number-conserving contributions to the collision integrals [cf.Eq (3)] remain number conserving in this approximation.Also, we focus on quasiparticles of low energy, ∆ < E ≲ 2∆; assuming the phonon bath temperature T B to be sufficiently low (at least compared to the gap, although more stringent conditions will be discussed), phonon absorption can be neglected and the phonon collision integral is approximately given by [cf.
Ref. [16] and diagrams c), e), and f) in Fig. 2] where ϵ = E − ∆ is the energy measured from the gap [with a slight abuse of notation, we also substitute f (E) = f (∆ + ϵ) → f (ϵ)] and St P hon g is obtained from Eq. (4) following the steps given after Eq. ( 5).The last term inside square brackets accounts for quasiparticle recombination accompanied by phonon emission; here we have introduced the quasiparticle density normalized by the Cooper pair density, x qp = N qp /2ρ F ∆, where the quasiparticle density is given by The parameter r = 4(∆/T c ) 3 /τ 0 is the (normalized) recombination coefficient, where depends on the thermalization time τ l of the phonons and τ P B 0 is the lifetime of a phonon at the pair-breaking threshold [the reduction of the recombination coefficient for τ l ̸ = 0 accounts for the processes in diagram g)].This time is proportional to τ 0 times the ratio between ion density over Cooper pair density and times (T c /ω D ) 3 , with ω D the Debye frequency; for aluminum, we use τ 0 /τ P B 0 = 1.7 × 10 3 [17].The two terms inside square bracket in Eq. ( 6) define the energy-dependent spontaneous phonon emission rate and the density-dependent phonon-assisted recombination rate 1 τ qp r = rx qp = RN qp (10) with the recombination coefficient R = 2∆ 2 /ρ F τ0 T 3 c (which, strictly speaking, depends weakly on energy and/or the width of the distribution function above the gap -that is, the "effective" quasiparticle temperature, -see [17,27]; we neglect this factor-of-order-one correction in this work).The fact that the recombination time depends on the density represents the main non-linearity of the kinetic equation; as we will show, however, the density can be calculated without a detailed solution of the kinetic equation in an approach reminiscent of the phenomenological one pioneered by Rothwarf and Taylor [28].In fact, to use that approach one must also take into account the last term in the right-hand side of Eq. ( 6), namely the generation of quasiparticles by phonons; for f ≪ 1, that term depends only on the phonon distribution function n and not on f .In this work we will not need the detailed form of St P hon g (see e.g.Ref. [17]), rather the integral of its product time the quasiparticle density of states; for phonons in thermal equilibrium at temperature T B we have where r 0 is the recombination coefficient r for τ l = 0 and in the last term we introduce the notation for the generation coefficient G T B at the given temperature T B .
Two additional approximations can be introduced for the collision integral with pair-breaking photons, Eq. (2).First, for quasiparticles in the relevant energy range, in St P hot {f, nP B } the terms proportional to U + (E, E − ω P B ) [cf.Eq. (3)] vanish since E − ω P B < ∆; we also assume that f ≪ 1 at all energies and that the distribution function decays quickly with increasing energy, so that f (E + ω P B ) ≪ f (E)n P B /(n P B + 1), and arrive at where in the last equality we have introduced the quasiparticle lifetime against the absorption of pair-breaking photon.Second, we assume that photon-assisted recombination can be neglected in comparison to the phononassisted one, that is [see Eqs. ( 5) and (10)] Clearly, this inequality is violated as E → ω P B −∆ due to the divergence of the superconducting density of states; however, we will show that the violation happens so extremely close to that energy to have no impact on our results.
In closing this introductory section, we note that in all the formulas above the gap value ∆ should be understood as given by the self-consistent equation for the order parameter.In the presence of a small density of quasiparticles, the gap is smaller than its no-quasiparticle value ∆ 0 ; defining δ∆ = ∆ 0 − ∆, at leading order we have δ∆/∆ 0 ≃ x qp [17] and for most purposes the difference between ∆ and ∆ 0 can be ignored.One exception is the evaluation of the resonator's frequency shift, which we discuss in Sec.V A.

III. GENERATION BY PAIR BREAKING PHOTONS
As a first step, in this Section we study the steady-state quasiparticle distribution in the presence of pair-breaking photons only,a situation complementary to that analyzed in Refs.[16,17] in which only non-pair-breaking photons were taken into account.Concretely, we assume that there are no modes below the pair-breaking threshold and we set St P hot {f, n} = 0 (that is, c QP P hot = 0).We also assume that the phonons are at zero temperature, n(ω) = 0, or in other words fast thermalization (τ l → 0) with a T B = 0 bath.Therefore, the kinetic equation reduces to 0 = St P hot {f, nP B } + St P hot P B,g {f, nP B } + St P hon {f, n} with the pair-breaking photon collision integrals of Eqs.(3) [appropriatiely modified for pair-breaking photons] and (4) and the phonon collision integral of Eq. (6) [in which by assumption the last term accounting for generation by phonons vanishes].
Even in this simplified case, we can in principle distinguish two different regimes, depending on which process is dominant near the gap.Since at those energies phononassisted recombination is faster than phonon emission [1/τ qp e,n (ϵ) ≪ 1/τ qp r as ϵ → 0], we should compare the rate for the former process, 1/τ qp r [Eq.(10)], to the rate of absorption of photons, 1/τ qp abs,P B [Eq. (12)].We show below (Sec.III A) that the situation of experimental relevance for Al resonators is In fact, to check when this inequality holds, we need to eliminate the dependence of the right-hand side on the quasiparticle density.To that end, we multiply Eq. ( 14) by ρ(E), integrate over E, and use that St P hot is number-conserving to find, for f (ϵ) ≪ 1, where and the structure factor S − (x) [21,29] has the approximate form S − (x) ≃ π(x − 2)/2 for x − 2 ≪ 2. Using Eq. ( 16), the inequality in Eq. ( 15) can be rewritten as where from now on we assume ξ < ∆ (ω P B < 3∆).When the above condition is true, in most cases quasiparticles created by photon pair breaking relax towards the gap and/or recombine by phonon emission before they can absorb a photon; indeed, the condition can be satisfied if the number of photon is sufficiently small.If a photon is absorbed before a recombination event, this will lead to a second peak in the quasiparticle distribution function at energies above 3∆.At those energies, quasiparticle relaxation by phonon emission is much faster than near the gap, so assuming that to be the dominant process, one can treat the second peak in a perturbative way, similarly to the "cold" regime of Ref. [16].We do not pursue this approach further in this section, since if the condition in Eq. ( 15), or equivalently Eq. ( 18), is satisfied, the effect of the second peak can be neglected and only energies ϵ < ξ are relevant (see, however, Sec.IV B for the expression for the second peak).So long as Eq. ( 15) holds, we can further simplify Eq. ( 14) by ignoring the first term on the right-hand side to get Next, we introduce the dimensionless energy variable and the function where the dimensionless parameter determines whether at the highest energy at which quasiparticles are generated (that is, ϵ = ξ) recombination is faster, γ ′ * > 1, or slower, γ ′ * < 1, than relaxation by phonon emission.With this notation, the steady-state equation (19) takes the form of a Volterra integral equation of the second kind with the kernel In Eq. ( 23) the first term on the right-hand side accounts for quasiparticle out-scattering by spontaneous phonon emission and the second one for the corresponding inscattering process; in that term the upper limit of integration is set to 1 consistently with ignoring absorption of pair-breaking photons and thus occupation at energies ϵ > ξ, as discussed earlier in this section.The term on the left-hand sides originates from the quasiparticle generation by photons; its divergence as γ → 1 originates from the diverging density of states and, as we will see, leads to a divergent distribution function.This unphysical divergence is in fact cut off by the Pauli blocking factors that we have neglected by assuming f ≪ 1; however, the divergence is integrable and to our knowledge it does not lead to unphysical behavior of any observable, so it is not further considered.Depending on the value of γ ′ * we can distinguish two cases.For γ ′ * ≳ 1, the solution to Eq. ( 23) can be given in terms of a Neumann series [30] with over the whole range γ ∈ [0, 1].In fact for γ ′ * ≫ 1 we can ignore (γ ′ ) 7/2 in the denominator of the kernel, Eq. ( 24), and readily see that the terms in the series are suppressed by (γ ′ * ) −7j/2 .The second case, γ ′ * ≪ 1, is relevant to small quasiparticle density.In this case the solution has to be constructed differently.For γ ≫ γ ′ * , the kernel can be simplified to Using this simplified kernel, the solution to Eq. ( 23) can again be written as a Neumann series, Eq. ( 25), but using the approximate kernel Ĩ in Eq. ( 26) instead of I.
The first and second order terms can be calculated analytically and are given in Appendix B. In fact, such a solution is valid under a weaker assumption than γ ≫ γ ′ * : Using the simplified kernel is a good approximation if the dominant contribution to the integral in Eq. ( 23) comes from the interval γ ′ ∈ [γ ′ * , 1].This holds for γ ≫ γ * , with γ * defined by (using the simplified kernel in the left-hand side is a reasonable approximation, while using it in the right-hand side overestimates the value of that integral, since the full kernel is smaller than the simplified one for small γ ′ ; this leads to a conservative estimate for γ * ).While γ ′ * determines whether out-scattering at energy γ is dominated by spontaneous phonon emission (γ > γ ′ * ) or recombination (γ < γ ′ * ), γ * distinguish whether in-scattering mostly originates from states for which spontaneous phonon emission (γ > γ * ) or recombination dominates (γ < γ * ).Equation ( 28) can be solved numerically using a bisection algorithm, the solution is presented in Fig. 3.We see there that γ * is at least an order of magnitude smaller than γ ′ * ; taking the geometrical average between the two quantities, we then estimate that the Neumann series constructed with the simplified kernel is a good approximation for γ ≳ γ c ≡ γ ′ * /3.Except for the zeroth order, the terms in the "simplified" Neumann series diverge as γ → 0; this is a result of the use of the simplified kernel, which has a nonintegrable singularity for γ = 0 as γ ′ → 0, while the exact kernel approaches zero in that limit.Thus, the solution for γ ≲ γ c must be constructed differently.We proceed by considering a series expansion for ϕ up to third order in γ, FIG. 4. Coefficients of the low-energy expansion for ϕ, Eq. ( 29).For better visibility, the products a l (γ ′ * ) l , l = 0 to 3, are plotted.Points are calculated using the numerical result for ϕ in Eq. ( 31), while lines using the approximate Neumann series as described in the text.
Neglecting small factors of order (γ c /γ ′ * ) 7/2 and higher, the lower integration limit in the right-hand side of Eq. ( 23) can be set to zero, and for consistency the lefthand side should be expanded also only up to third order in γ.Then the expansion coefficients for ϕ are given by with In this expression ϕ should be understood as the exact solution, extending over the whole interval [0, 1].The coefficients can be then calculated numerically, see the dots in Fig. 4.However, for an approximate estimate of the coefficients, we can set the lower integration limit in Eq. ( 31) to γ c and use the "simplified" Neumann series solution for ϕ including terms up to j = 2 (see Appendix B).
The approximate coefficients calculated in this way are shown as lines; except for a 3 , we find good agreement between exact and approximate coefficients.For later use, we note that a 0 ≃ 0.073 log 3 (12/γ ′ * ) − 0.35 log 2 (12/γ ′ * ) + 1.23 log(12/γ ′ * ) is a good approximation in the range covered by the dots in Fig. 4 (the dependence on powers of the logarithm originates from the behavior of ϕ (j) , j = 0, 1, 2, at small γ, while the numerical coefficients are obtained by comparison to the numerics.
We exemplify the results obtained in this section in Fig. 5.The black lines show the distribution function as obtained by numerical solution of the kinetic equation, while the colored lines are the analytical approximations.Both the "simplified" Neumann series [cf.Eq. ( 25)] and the low-energy expansion [Eq.(29)] are in excellent agreement with the numerical results in their respective regimes of validity.
The consideration of this Section can be straightforwardly generalized to account for multiple pair-breaking modes: to determine the quasiparticle density x qp , it suffices to replace the right-hand side of Eq. ( 16) with its sum over all modes.Then, Eq. ( 19) [or equivalently Eq. ( 23)] being linear in the distribution function, it can be solved for the contribution f i (ξ i γ i ) of each pairbreaking mode i as for a single mode.The distribution function f is finally given by the sum over all f i .

A. Validity of assumptions
In obtaining the approximate solution for the distribution function, we made a number of assumptions whose validity depend on parameters such as the recombination rate r and the quasiparticle-phonon coupling constant c QP P hot,P B ; for aluminum, as order of magnitude we take r ≃ 10 7 Hz [11], and we estimate in Appendix A that for resonators of this material c QP P hot,P B ≃ 10 Hz.The condition in Eq. ( 18) can be interpreted as a bound on the number of pair-breaking photons, which with these parameters reads nP B ≪ 10 6 ξ/∆.Even for pair-breaking photons near the threshold (ξ → 0), this conditions is in practice always satisfied for nP B ≲ 1: while for subgap photons high quality factors of order 10 6 -10 7 are possible [8,18], the quality factor of above-gap modes is more than two orders of magnitude smaller, see Appendix A; this implies ξ/∆ > 10 −5 , since the frequency of the pair-breaking mode is known with relative precision given by the inverse quality factor.Alternatively, for pair-breaking photons of sufficient energy, ξ/∆ > 0.1, the bound on the photon number becomes nP B < 10 5 .We note that if Eq. ( 18) holds, then we also find from Eq. ( 16) that x qp ≪ 1; this ensures that the suppression of the gap δ∆/∆ ≃ x qp can also be neglected for the purposes of this section.
We now turn to the validity of Eq. ( 13), that is, negligibility of photon-mediated recombination.Setting ε = ω P B − ∆ − E, in the limit ε → 0 (i.e., ϵ → ξ) we can rewrite that inequality in the form where we used Eqs.( 21) and (29).For nP B ≫ 1 and assuming γ ′ * ≳ 1, we have a 0 ≃ 1 and, for the parameter discussed above we have ε/ξ ≫ 10 −14 n2 P B ; even at the upper bound nP B ≃ 10 5 determined in the previous paragraph, the right-had side is small, of order 10 −4 .In the opposite regime nP B , γ ′ * < 1, the first factor on the right-hand side is of order unity (as can be verified by inspection of Fig. 4 or by using the approximate analytical expression for a 0 given previously) and the second factor is ∼ 10 −12 , meaning that the requirement becomes significantly less stringent.Therefore, as mentioned at the end of Sec.II, recombination by photon emission can be ignored.
At the beginning of this section, we assumed that there are no other modes beside that of the pair-breaking pho-tons by setting c QP P hot = 0.If such a mode is present (c QP P hot > 0), even if it is unpopulated (n = 0) it can in principle affect the distribution function, since quasiparticles could relax by emitting a photon of energy ω 0 .However, similarly to the process of photon-mediated recombination considered in the previous paragraph, we now show that we can ignore this relaxation mechanism.Indeed, this mechanism is accounted for by the last line of Eq. ( 3), and the corresponding relaxation rate diverges for ϵ → ω 0 .We compare this rate to that of relaxation by phonon emission, Eq. ( 9) (assuming that 1/τ qp e,n (ω 0 ) > 1/τ qp r , as we expect to be the case at low temperature and hence quasiparticle density); we find the condition where ε = ω 0 + ∆ − E ≪ ω 0 .This inequality is equivalent to ε/∆ ≫ (∆/ω 0 ) 7 (105c QP P hot /16r 0 ) 2 , and using c QP P hot ≃ 1 Hz (see Appendix A), we estimate that the second factor on the right-hand side is small, of order 10 −13 ; while the first factor could in principle be large, since ω 0 ≪ 2∆, in practice the resonators are designed to have modes in the frequency range above at least a couple of GHz, in which case (∆/ω 0 ) 7 ≲ 10 9 .Hence we conclude that the above condition is violated only for energy so close to ∆ + ω 0 as to not affect the quasiparticle distribution.In fact, the numerical solutions displayed in Fig. 5 were obtained with c QP P hot = 1 Hz and ω 0 /2π = 4.84 GHz and show no significant deviation from the analytical approximation derived assuming c QP P hot = 0. Finally, we also assumed that we can ignore phonons by setting T B = 0 and τ l = 0.Even within the assumption of fast thermalization (τ l ≪ τ P B 0 ), thermal phonons with T B > 0 can generate quasiparticles by breaking Cooper pairs; therefore, a necessary condition to ignore them is that this generation mechanism gives a negligible contribution to the quasiparticle density.By comparing the second line in Eq. ( 16) (that is, the quasiparticle density generation rate due to pair-breaking photons) to the thermal phonon generation rate [see text after Eq. ( 10)], we define the crossover (or effective generation) temperature with W with the Lambert (product logarithm) function.For T B > TB thermal phonons are the main source of quasiparticles, while for T B < TB , the dominant generation mechanism is photon pair breaking and TB gives the temperature at which quasiparticles in thermal equilibirum would have the same density as those generated by the photons.For the typical parameters discussed above (r 0 = 10 7 Hz, c QP P hot,P B = 10 Hz), even for photons with frequencies close to the pair-breaking threshold, ξ/∆ = 10 −3 , and a low occupation probability, nP B = 10 −10 , we estimate a relatively high crossover temperature TB > 100 mK.In fact, for the crossover temperature to go below e.g.20 mK the photon occupation probability should be extremely small, nP B < 10 −84 , so ignoring generation by thermal phonons should likely be a good approximation for typical low-temperature experiments; in our examples in the plots in Fig. 5 we have assumed nP B ≥ 10 −7 .
While the condition T B < TB implies that the density is not affected by thermal phonons, it is not sufficient to ensure that they do not alter the shape of the quasiparticle distribution.A more stringent condition is found by requiring that the rate of phonon absorption at the gap is small compared to the recombination rate [as defined in Ref. [17], the former rate is approximately given by 3/τ qp e,n (T B )], or equivalently by requiring that the typical energy gained by absorbing a phonon is small compared to the width of the distribution in the absence of phonons, T B < γ ′ * ξ.These requirements can be expressed as x qp > (T B /∆) 7/2 .With the same parameters used above and assuming ξ > 0.1∆, for T B = 10 mK we estimate, using Eq. ( 16), that the condition is satisfied for nP B > 10 −8 .The bound becomes more stringent as the phonon temperature increases; we do not explore this higher temperature regime further here, as we focus next on a competing mechanism affecting the distribution shape, namely the absorption and emission of non-pairbreaking photons.
The results of this Section show that at low temperature a small number or photons above the pair-breaking threshold can lead to experimentally relevant quasiparticle densities.In fact, using Eqs.( 16) and ( 21) one can show that if both ξ and γ ′ * are not too small compared to unity, then f (0) (cf.Fig 5) is of the same order of magnitude as the normalized density x qp , which in various experiments has been estimated to range between 10 −9 and 10 −5 (see [31] and references therein).Even the width ∼ γ ′ * ξ of the quasiparticle distribution can be determined by the interplay between pair breaking by photons and recombination, rather than by absorption of thermal phonons.This could happen in particular in superconducting qubits, for which pair-breaking photons can contribute directly to qubit transitions as well as to quasiparticle generation [32,33], although the actual shape of the distribution function is likely affected by the fact that typically the films forming the qubit's junction have different gaps.While qubits are by necessity operated in a regime where at most a few non-pair-breaking photons are present, this is not the case for resonators, where the large number of non-pair breaking photons can determined the distribution's shape, as we discuss next.

IV. NON-PAIR-BREAKING PHOTONS
When both modes above and below the pair-breaking threshold are populated, n, nP B > 0, various regimes are possible.For instance, for n sufficiently small compared to nP B we can expect the quasiparticle distribution function to be mainly determined by the pair-breaking photons; then the effect of the non-pair-breaking ones can be treated perturbatively.Still, we can expect the effect to be different depending on the frequency ω 0 < 2∆ satisfying ω 0 < γ c ξ, γ c ξ < ω 0 < ξ, or ω 0 > ξ (we have assumed γ c < 1, otherwise ω 0 should be compared to ξ only).In the opposite regime of sufficiently large n, the distribution function dependence on energy is instead mainly due to the below-threshold photons and the pair-breaking one mostly contribute to the overall quasiparticle density.In what follow, we examine in more detail some (but not all) of these many regimes [34].

A. Perturbative regime
We begin by investigating the effect of a small number of non-pair-breaking photons of frequency ω 0 < 2∆ on the quasiparticle distribution function.We assume that the approximations employed in Sec.III are still applicable, so that the steady-state kinetic equation becomes [cf.Eq. ( 19)] 0 = St P hon {f, n} with St P hot {f, n} given by Eq. ( 3).Simplifying that equation in the regime of low quasiparticle energies (ϵ < ∆) and densities (f ≪ 1) we get where the last line vanishes for ϵ < ω 0 .
In the previous section we have considered already when the effect of this term is negligible in the case n = 0 by focusing on the decay rate associated with the last term on the right-hand side of Eq. ( 36) [a finite n adds the factor (1 + n) to the left-hand side of condition (33)].For finite but small occupation, n ≲ 1, to find a necessary condition enabling us to treat the non-pair-breaking photons as a perturbation we consider the absorption rate arising from the second term in the first square brackets, which is highest at the lowest energies ϵ < ω 0 at which the second line of Eq. ( 36) is zero [at higher energies, the photon absorption accounted for by the first term in the second square brackets should also be considered]; we then compare that rate to the recombination rate, Eq. ( 10), to find the condition n ≪ rx qp where x qp depends on nP B as follows from Eq. ( 16).Assuming the inequality holds, we write the quasiparticle distribution in the form where f 0 is the steady-state distribution in the absence of low-energy photons, n = 0, and the kinetic equation can be rewritten as where in the phonon collision integral we dropped the explicit dependence on the phonon distribution as we assume zero temperature, n = 0 (this implies in particular St P hon g = 0).Although the photon collision integral in Eq. ( 3) accounts only for emission or absorption of a single photon at a time, consecutive absorption processes uninterrupted by phonon emission lead to peaks in the energy distribution function at multiples of the photon energy.Thus, we seek an expression for δf in form of an expansion δf = f 1 + f 2 + ... using an iterative approach to find f m once f m−1 is known: with m = 1, 2, .... To find an explicit expression for f m , m ≥ 1, we ignore the in-scattering part of St P hon {δf } -that is, the first term on the right-hand side of Eq. ( 6); we will comment below on this step.We also assume 1/τ qp e,n (ω 0 ) > 1/τ qp r , so that for ϵ > ω 0 we can ignore the term proportional to x qp in St P hon .Furthermore, we note that in St P hot there is at each iteration a term that diverge as ϵ → mω 0 originating from the last line in Eq. (36).Of the two terms in square brackets there, we expect the first one to be dominant; in other words, we assume With these simplifications, we find peaks for ϵ > mω 0 with approximate shape where [16,17] The value of T * being larger than ω 0 means that the nonpair-breaking photons are effective at heating the quasiparticles [16,17]; similarly, here it can counteract the fast suppression of the amplitude of the peaks caused by the terms in the second line in Eq. ( 42) being approximately 1/m 7/2 [(m − 1)!] 4 for ϵ ≃ mω 0 .
The result for the peaks in Eq. ( 42) is compatible with Eq. ( 41) so long as n > 1/(m + 1) 4 .We can also estimate if a given peak is visible by comparing its amplitude (which we estimate at ϵ = mω 0 + γ ′ * ξ) to the value of f 0 at the peak's position.We thus find the "visibility condition" where we assumed ω 0 > γ ′ * ξ.Using this expression, one can for instance check under which conditions the first peak is visible while at the same time Eq. ( 37) holds.More interestingly, for a given choice of parameters as m increases the left-hand side decreases, so there is a last visible peak which we denote with index m l (we have implicitly assumed that m l < ξ/ω 0 , so that it is consistent to ignore occupation at energies ϵ > ξ as in Sec.III).For m > m l , it is reasonable to set f m = 0; in turn, this makes it possible to understand why Eq. ( 42) gives a good description of the last peak but not necessarily of lower-index peaks.Indeed, in our derivation we have ignored the first term on the right-hand side of Eq. ( 6); considering the contribution of that term to Eq. ( 39) for ϵ > mω 0 , one can check that the integral between ϵ and (m + 1)ω 0 can be ignored compared to the first term in square brackets in Eq. ( 6).At energies just above (m + 1)ω 0 , the integral is dominated by the f m+1 term, which we ignore for m = m l .On the other hand, taking for example m = m l − 1, the m l peak gives a non-negligible contribution to the integral, a contribution that accounts for quasiparticle relaxing to lower energies by emitting phonons; this mechanism then leads to the lower-index peaks to be broader than what Eq. ( 42) predicts.In Fig. 6 we show with points the results of numerical solutions of the kinetic equation simulation for different number of non-pair-breaking photons.Only for the smallest number Eq. ( 37) holds; nonetheless, in all cases the number of visible peaks agrees with the expectation from Eq. ( 44) and the last visible peak is well described by Eq. ( 42), see the solid lines.

B. Heating regime
We next consider the regime in which there are many non-pair-breaking photons, such that they can heat the quasiparticles, T * > ω 0 .Even if this condition is met, the pair-breaking photons could affect the shape of the distribution function; therefore, we additionally require that the generation of quasiparticles can be ignored at lowest order, so that the results of Refs.[16,17] can be used as a starting point.Near energy ξ we can approximate the photon generation collision integral, Eq. ( 4), as which diverges as ϵ → ξ; however, since the width of the peaks in the distribution function is limited by the photon frequency ω 0 , for our estimates we replace ξ − ϵ → ω 0 .If the shape is determined by the balance between phonon emission and the absorption/emission of non-pair-breaking photons, the above contribution from St P hot P B,g should be small in comparison to, for instance, the phonon collision integral, or more precisely the term proportional to (ϵ/∆) 7/2 in Eq. ( 6), where we can use [17] for ϵ = ξ.Assuming that the quasiparticle density is determined by the photon pair-breaking mechanism and hence given by Eq. ( 16), we arrive at the condition For the parameters of Fig. 6 this correspond to T * /ξ ≳ 0.29, while for the highest photon number there we estimate T * /ξ ≃ 0.1.Even doubling the value of T * so that T * /ξ = 0.2 (or c QP P hot n ≃ 8.2 × 10 4 Hz), see Fig. 7, the considerations of the previous subsection are still valid with minor modifications.For example, according to Eq. ( 44) the m = 7 peak should be the last visible one, but since that peak is close in energy to ξ, the right-hand side of that equation should be multiplied by a factor of order ξ/ω 0 ; including this factor, we find that the inequality holds only weakly.Further increasing T * above the threshold value by setting T * /ξ = 0.3 (c QP P hot n ≃ 9.3 × 10 5 Hz), the shape of the distribution function follows the predictions of Ref. [17] up to energies beyond ξ.47).We use the same parameters as in Fig. 6 except for the photon number n, given in the bottom left inset in terms of the ratio between T * of Eq. ( 43) and ξ of Eq. ( 17).The top right inset shows the second peak at energies ϵ > ωP B , with the blue dashed curve corresponding to Eq. ( 48) where f (ϵ) is the numerically calculated distribution function at low energies.
In Fig. 7 it is evident that the distribution function changes qualitatively at a crossover energy ϵ c > ξ, becoming mostly smooth in the range ϵ c ≲ ϵ ≲ ω P B .We ascribe this smooth part to "excess" quasiparticles relaxing by phonon emission (and/or scattering with non-pairbreaking photons), the excess being understood as the second broad peak at energies above ω P B .As mentioned in Sec.III, this second peak is due to the absorption of pair-breaking photons, and can be found by following the approach of Ref. [16] [cf.Eq. ( 8) there]: ; (48) we compare this equation to numerical calculations in the top right inset of Fig. 7.These modifications affect the distribution only at the relatively large energies ϵ > ϵ c and, as explained in Sec .II, do not contribute significantly to quasiparticle creation even when considering a finite phonon thermalization time (see also the next subsection).Furthermore, since usually r 0 ≫ c QP P hot,P B , the second peak is much smaller than the first one and thus does not appreciably affect measurable properties like the quality factor of a resonator.

C. Quasiparticle density and its fluctuations
So far we have assumed that quasiparticle generation by phonons can be ignored.Once we enter the heating regime, however, the distribution function is not significantly affected by the pair-breaking photons, so that the latter can only influence the overall quasiparticle density.Therefore, in this regime we can extend the generalized Rothwarf-Taylor equation determining the quasiparticle density by including contributions from both pairbreaking phonons and photons, with G T B of Eq. ( 11) and G TB , accounting for pairbreaking photons, obtained from G T B by replacing r → r 0 and T B → TB , see Eq. ( 34); this is equivalent to G TB = 2πρ F c QP P hot,P B nP B ξ.The term linear in N qp depends on T * because it describes generation by nonequilibrium phonons emitted by quasiparticles heated by the non-pair-breaking photons to energies above 3∆ [see diagram h) in Fig. 2 and Ref. [17]]: Equation ( 49) can be derived by multiplying Eq. ( 1) by ρ(E) and integrating it over energy; in the collision integrals, the (nonequilibrium) phonon distribution can be expressed in terms of the quasiparticle distribution, while the energy dependence of the latter, but not its normalization, are determined by the balance between phonon scattering processes and non-pair-breaking photon absorption; a detailed derivation is given in Ref. [17].
Clearly, as far as the density is concerned, pairbreaking phonons and photons play the same role; for r ∼ r 0 , as discussed after Eq. ( 34) one of the two generation mechanisms is dominant, depending on which of the two temperatures, T B or TB , is the largest.Then defining T G = max{T B , TB } and T * B = T 3 * /∆ 2 [17], we can distinguish two steady-state regimes, in analogy to those discussed in Ref [17] in the absence of pairbreaking photons.At high heating, T * B > T G , we can neglect the first two terms on the right-hand side and find N qp,0 ≃ G(T * /∆)/R, independent of phonon temperature and number of pair-breaking photons (we use subscript 0 to denote the steady state); for low heating, T * B < T G , we can neglect the term linear in N qp in Eq. ( 49) and the density is approximately N qp,0 ≃ (G T B + G TB )/R, independent of the number of nonpair-breaking photons n.In both regimes, the dependence of the density on bath temperature are similar: at low temperature -T B < TB for low heating and T B < T * B for high heating -the density is approximately constant, while above these crossover temperatures it increases exponentially, being approximately the same as in thermal equilibrium.As an example, in Fig. 8 we plot the density as function of temperature for a few values of TB in the low-heating regime; the behavior resembles that in the high-heating regime, see Ref. [17].
Interestingly, the two regimes display quantitatively different power spectral densities for quasiparticle number fluctuations.In general, as shown in Appendix C the spectral density (per unit volume) can be written in the form where τr = 1/[2RN qp,0 − G(T * /∆)] is the relaxation time of quasiparticle number fluctuations and we assumed fast phonon dynamics, τr ≫ τ P B 0 , τ l .At low heating, we have τr ≃ 1/2RN qp,0 and we recover the result of Ref. [35].At high heating τr ≃ 1/RN qp,0 , which modifies the relation between measurable quantities such as τr and S N (0) and quantities than can be derived from them, such as the steady-state quasiparticle density.
It should be noted that in experiments the fluctuations measured are not directly those of the quasiparticle number, but those of phase and amplitude of the complex transmission.Therefore, comparison to experiments requires knowledge of the responsivities [9,36] where Q is the loaded quality factor, α the kinetic inductance fraction, σ 1(2) the real (imaginary) part of the ac conductivity, Q i,qp the contribution of quasiparticles to the internal quality factor, and δω 0 the frequency shift.For a coupling-limited quality factor, both responsivities are approximately independent of quasiparticle number.However, in the heating regime the responsivities depend on the photon number n via Q i,qp and δω 0 , as investigated in Ref. [17] -see also the next section.
In closing this section, we remark that the considerations presented in this work for the heating regime can be straightforwardly extended to the case of multiple pair-breaking modes (cf.Ref. [37]) by introducing more sources like the term G TB in Eq. ( 49).

V. COMPARISON TO EXPERIMENTS
We now turn our focus to the implications of the results of the previous sections for experiments.We reconsider in particular the measurement of quality factor and resonant frequency as functions of temperature and for various readout powers reported in Ref. [18].We reproduce in Fig. 9 the internal Q-factor data from that reference together with fit lines calculated using the theory of quasiparticle heating of Ref. [17], which does not include the effect of pair-breaking photons, with the addition of an "extrinsic" dissipation mechanism of unknown origin; this second mechanism is needed in order to capture the low temperature saturation of the Q-factor.Concretely, the inverse internal quality factor is assumed to be of the form where, at leading order in T * /∆, we have [17] Q leading to a linear relation between inverse quality factor and quasiparticle density N qp , with a slope that decreases with increasing readout power (the condition T * > ω 0 ensures that the quality factor is not affected by the detailed shape of the peaks on the scale ω 0 , but only on the overall width T 0 of the distribution function).In other words, since at low heating the quasiparticle density is independent of the photon number n, the quality factor is expected to increase with readout power, as n and hence T * increase; indeed, when the coupling quality factor is small, Q c ≪ Q i , we have n = 2Q c P read /ω 2 0 [17].While the increase of Q i with power agrees with observations at sufficiently high temperatures, it is qualitatively inconsistent with measurements at the lowest temperatures (crossover to the highheating regime with decreasing temperature cannot explain the data, see Ref. [17]).We emphasize that Q i,ext is introduced purely in a phenomenological way, and ranges from approximately 2.5 × 10 6 at low power to 0.7 × 10 6 at high power.The decrease at higher power excludes the standard explanation of the low-temperature value of Q i,ext as originating from two-level systems, since the saturation of the latter with increasing power would lead to an increase of the quality factor [19].
To try and explain the low-temperature behavior of the quality factor, we now consider the effect of pair-breaking photons, assuming we remain in the low-heating regime.The quasiparticle density can be strongly affected by a small number of such photons, see the discussion after FIG. 9. Quality factor vs temperature for different readout powers.The experimental points are taken from Ref. [18].The solid lines are based on Eq. ( 54) and are calculated using the following parameters: ∆ = 189 µeV, c QP P hot = 58 mHz, τ0 = 63 ns, ω0/2π = 5.3 GHz, α = 0.13, τ l = 170 ps.More details can be found in Ref. [17].
Eq. (34).In fact, in a series of works [38,39] radiation emitted by a higher temperature (3-4 K) stage of the refrigerator was identified as a source of quasiparticle generation, and at least partially mitigated by using a 'box in a box' design.Furthermore, to reduce the influence of pair-breaking photons propagating through coaxial cables, filters have been used at the cold stage that, extrapolating measurements below ∼ 10 GHz to twice the gap, are expected to suppress the number of pair-breaking photons by almost 4 orders of magnitude.Despite such careful engineering, the "leakage" of pairbreaking photons from higher temperature stages is likely not completely eliminated.However, at low temperatures such leakage would result in a density of quasiparticles independent of n and of T B (up to TB , see Fig. 8), and the quality factor should still increase with n, in contrast with observation.
Since the assumption of a fixed number of pairbreaking photons is inconsistent with experiment, we are lead to amend it by setting where the first term on the right-hand side accounts for the leakage from higher temperature stages discussed above.The second term introduces a dependence of nP B on n and can originate from high-frequency noise produced by the source used to generate the signal probing the resonator.This dependence in turns implies that when n ≳ nP B,0 /µ the quasiparticle density increases with readout power even in the low-heating regime, because the second term on the right-hand side of Eq. ( 49) increases with n.As shown in Fig. 10, with this assumption we can obtain a satisfactory fit to the measurements of Ref. [18] using reasonable parameters: the "leaked" photon number n P B,0 is about one order of magnitude bigger than what one would expect assuming a higherstage temperature of 3.5 K with a 10 −4 attenuation factor for photons of energy ω P B = 2.8∆, but the discrepancy would drop to just a factor of 3 for photons at the pairbreaking threshold [40].As for the value of µ = 10 −9 , it can be explained by a combination of 10 −4 attenuation times a 10 −5 ratio between power of the generated signal and of the high-frequency noise [41].
The assumption in Eq. ( 56) is independent of temperature, so we can employ it also when evaluating the temperature dependence of the quality factor, as shown in Fig. 11.We stress that while in Fig. 9 the value of Q i,ext is chosen independently for each of the six readout powers as to best fit the data, in Fig. 11 just the two parameters appearing in the right-hand side of Eq. ( 56) are used to fit all the data at the same time [42].The agreement between theory and data in Figs. 10 and 11 thus support our proposed explanation for the physical mechanism responsible for the saturation of the quality factor in terms of excess pair-breaking photons.
We also remark that the way pair-breaking photons affect the quality factor is not equivalent to adding an extrinsic mechanism: the pair-breaking photons influence Q i,qp via the quasiparticle density N qp , and in the lowheating regime we can write the latter in the form with N T B = G T B /R and the analogous definition for FIG.11.Quality factor vs temperature for different readout powers.As in Fig. 9, the experimental points are from Ref. [18].Here the solid lines include only the quasiparticle contribution, with no extrinsic loss mechanisms.The lines are calculated using the same parameters as in Fig. 10 but, in contrast to Fig. 9, here we take into account effect of a powerdependent number of pair-breaking photons, see Eq. ( 56) and the discussion that follows it.
N TB .Defining Q i,T B ( TB ) by the replacement N qp → N T B (N TB ) in Eq. ( 55), we find that which should be contrasted with Eq. (54).Finally, we point out that the above analysis of the quality factor is not restricted to a single pair-breaking mode, but is valid in the heating regime if multiple pair-breaking modes (or even a broadband source) are present: the heating regime is by definition that in which the shape but not necessarily the normalization of the distribution function is determined by the interplay between non-pair-breaking photons and phonons, see Sec.IV B; in this regime, the pair-breaking photons influence the quality factor only by contributing to the total quasiparticle density as in Eq. (57).

A. Frequency shift
The approach presented above to calculate the internal quality factor can also be applied to evaluate the shift in resonant frequency as function of the readout power.For thin-film resonators, the relative change in frequency is proportional to the kinetic inductance fraction times the relative change in the imaginary part of the ac conductivity [44].For the latter quantity, we use here the order-of-magnitude estimate discussed in Ref. [17] to find FIG.12. Main panel: Dots: relative shift in resonance frequency from Ref. [18] at low temperature, TB < 150 mK, relative to the lowest readout power, P read = −100 dBm.Solid line: numerical calculation of the shift for TB = 0 and other parameters as in Fig. 10.Dashed line: order-of-magnitude estimate from Eq. (59); both lines calculated as deviations from the value for P read = −100 dBm.Dotted line: prediction from Ref. [43] (corrected for the finite kinetic inductance fraction, see text).Inset: temperature evolution of the shift in linear scale (the colors are the same as in Fig. 11).
for the relative shift The terms in square brackets originate from the suppression of the order parameter by the nonequilibrium quasiparticles, while the last term from the direct effect of the nonequilibrium distribution on the imaginary part of the ac conductivity.The latter term is in fact only a rough estimate, since for its actual calculation knowledge of the shape of the first peak (between energies ∆ and ∆ + ω 0 ) in the distribution function is needed [17].Equation (59) predicts at low temperature a nonmonotonic dependence of the magnitude of the shift on readout power: at low power, x qp is independent of power as the quariparticle generation is dominated by the leakage term n P B,0 in Eq. (56), while the factor in curly brackets decreases slowly with power.As power increases and the last term in Eq. (56) becomes relevant, the magnitude of the shift increases with power as the increase in x qp dominates on the decrease of the curly-bracket factor.In practice, however, one can measure the shift from the lowest readout power, not from zero power, so depending on the parameters the nonmonotonic behavior might not be measurable.In fact, for the parameters used to fit the quality factor data the difference in relative shift between the two lowest powers is about 10 −6 which, given quality factors of order 10 6 , is at the limit of being mea-surable.For higher powers, the magnitude of the shift increases, in qualitative agreement with the experiment, see Fig. 12.In that figure, we plot both the analytical estimate above and the numerically calculated relative shift.Both curves deviate from the experimental data by less than a factor of 3.
We note that in Ref. [43] an alternative explanation for the shift was proposed: the authors argue that, in analogy with the case of a dc supercurrent, the nonpair-breaking photons cause a broadening of the peak in the superconducting density of states.This broadening effect can be characterized in terms of a depairing parameter α (not to be confused with the kinetic inductance fraction α) which is proportional to the readout power and which causes a change in the kinetic inductance: δL k /L k ≃ 4.85 α/∆ [43].The relative shift is in turn proportional to the change in kinetic inductance, δω 0 /ω 0 = −αδL k /2L k ; implicitly assuming α = 1, good agreement was found between the theory of Ref. [43] and the experiment of Ref. [18].However, the kinetic inductance fraction in the experiment can be estimated as α ∼ 0.1 [17], meaning that the theoretical result in fact underestimates the shift by approximately an order of magnitude compared to experiment.Interestingly, adding the contributions to the frequency shift of the two mechanisms improves agreement with the experimental data, but further investigation is needed to understand if in fact the two mechanisms coexist.
The inset in Fig. 12 displays the temperature dependence of the relative shift as extracted from the data of Ref. [18].It is instructive to present the same data as function of normalized density x qp instead, since Eq.(59) predicts a linear relationship with a slope dependent on readout power.For low readout power, the experimental data and the numerical estimates are in reasonable quantitative agreement with each other and with the expectation for thermal equilibrium, see Fig. 13.As power increases, the qualitative trend of decreasing magnitude of the slope in the frequency vs density plot is evident in both data and numerics; although there is quantitatve discrepancy at the higher powers, we can conclude that a slope smaller in magnitude than the thermal equilibrium one is indicative of quasiparticle heating.Recently, such a reduced slope has been observed in a transmon qubit in which excess quasiparticles are generated due to the presence of an IR laser beam focused onto the qubit or the substrate next to it [45].

VI. CONCLUSION AND OUTLOOK
In this work, we extend previous analytical studies of the quasiparticle distribution in superconducting resonators in the presence of (sub-gap) photons and phonons [16,17] to include the effect of a small number of photons of energy above the pair breaking threshold 2∆.We derive expressions for the distribution's shape in the limit of low bath temperatures and low numbers of sub-FIG.13.Dots: Experimental shift of resonance frequency from Ref. [18] (see inset in Fig. 12) vs quasiparticle density; the latter is obtained from the numerical solution to the kinetic equation with the same parameters as in the fits in Figs. 10 and 11.Solid lines: numerical results for the frequency shift calculated using the same parameters as for the density; readout powers increases from bottom to top (colors as in Fig. 11).Yellow dotted line: frequency shift assuming a thermal equilibrium distribution; for reference, xqp = 0.002 corresponds to TB = 350 mK.
gap photons, see Secs.III and IV A, complementing the results obtained before in the opposite limit of large numbers of sub-gap photons (see also Sec.IV B).In the latter regime, the quasiparticle density can be determined using a generalized Rothwarf-Taylor equation, Eq. ( 49); at typical temperatures of the phonon bath in the tens of millikelvin, we find that absorption of pair-breaking photons can be the dominant quasiparticle generation mechanism already at extremely low occupation nP B ≪ 1.As a consequence, the saturation of the quality factor of superconducting resonators at low temperatures observed experimentally in Ref. [18] can be explained by assuming a low number of pair-breaking photons, if this number has two contributions: one originating from higher temperature stages of the cryogenic setup and the second from the microwave signal generator, see Eq. (56) and Figs. 10 and 11.In addition to contributing to losses, quasiparticles also cause a shift in the resonator frequency, analogously to the Kramers-Kronig relation between real and imaginary parts of response functions (cf.Ref. [46] for the case of qubits).Interestingly, plotting the frequency shift as function of the quasiparticle density can provide information on the distribution in energy of the quasiparticles, see Eq. (59) and Fig. 13.
Our results show that resonators can be used to probe the high-frequency (that is, above the gap) electromagnetic environment of superconducting circuits, in this respect complementing the use of qubits to probe lowfrequency modes they are dispersively coupled to [47,48].The ability to investigate this environment can guide further efforts in reducing the detrimental impact of pair-breaking photons on the coherence time of superconducting qubits [21][22][23].The change of resonator response with power and/or temperature is also used to explore non-quasiparticle loss mechanisms such as dielectric and two-level system losses [49][50][51], so our results can contribute to improving the characterization of superconducting materials.
The similar calculation for the second order term gives ϕ (2)  and the dilogarithm defined as Li 2 (x) = x 1 dt ln(t) 1−t [note that sometimes a different convention is used, in which the integral is defined as giving Li 2 (1 − x) instead].

Appendix C: Power spectral density
To calculate the power spectrum of quasiparticle number fluctuations, we follow the approach of Ref. [35]; we summarize here the necessary steps solely to make our work self-contained, with no claim of novelty.We assume short phonon lifetimes compared to the quasiparticle lifetimes to neglect fluctuations in the phonon number; the applicability of this assumption to aluminum resonators is discussed also in Ref. [17].Using that each recombination/generation process involves two quasiparticles, we can identify their respective rates in Eq. ( 49), which we then use to set up a master equation for the probability P (N, t|k, 0) of finding N quasiparticles at time t if there were k quasiparticles at time 0, One can verify that Eq. (C4) agrees with Eq. ( 49) if fluctuations are small compared to the average quasiparticle number.In the steady state, the system approaches a probability distribution for which ⟨g(N )⟩ = ⟨r(N )⟩; for small deviations from the steady-state expectation value N 0 = ⟨N ⟩ we can write N = N 0 + ∆N and expand the generation/recombination rates up to second order in ∆N to find The autocorrelation function for the number can be written in the form If quasiparticle creation is dominated by photons or thermal phonons, one can neglect G, and the variance of the quasiparticle number is equal to its expectation value.
In the regime in which phonon trapping dominates, we instead have G(T * /∆) = RN 0 , and the variance is twice the quasiparticle number.Fourier transforming the autocorrelation function (C8) leads to equation (51) for the number spectral density.

FIG. 2 .
FIG.2.Schematic depiction of the main processes entering the collision integrals in Eq.(1).Straight lines correspond to quasiparticles, wavy lines to phonons, and dashed lines to photons.Diagrams a) and b) represent the absorption and emission of non-pair-breaking photons [cf.Eq. (3)], c) the emission of a phonon [c.f.first term in the bracket of Eq. (6)], d) quasiparticle generation by pair-breaking photons [Eq.(4)], a process not included in Ref.[17].Diagram e) depicts generation by a pair-breaking phonon [last term of Eq. (6)], f) provides the bare recombination rate r0xqp, and g) renormalizes the bare recombination coefficient r0 to r in Eq.(6) [see also Eq. (8)].Diagram h) describes an additional pair-breaking mechanism due to phonon emitted by quasiparticles with energy E > 3∆; this process is not included in the low-energy approximation of Eq. (6), but if n is sufficiently large (cf.Sec.IV) it can affect the density, see the term proportional to G(T * /∆) in Eq. (49) and Ref.[17].In Sec.. III we only consider c), d) and f), in Sec.IV we include diagrams a) and b), and in Secs.IV C and V we consider all diagrams.We do not include the diagram corresponding to photon mediated recombination [similar in structure to f)], as its contribution is generally negligible compared to phonon mediated recombination (see Sec. III A)

FIG. 3 .
FIG. 3. Numerically calculated γ * as function of γ ′ * , determining the range of validity γ > γ * for the validity of the Neumann series constructed using the simplified kernel in Eq. (27).

FIG. 6 .
FIG. 6. Distribution function in the presence of non-pairbreaking photons in the perturbative regime.The dots are obtained by solving the kinetic equation numerically with parameters ωP B = 2.8∆, τ l = 0, τ0 = 63 ns, ∆ = 189 µeV, c QP P hot,P B nP B = 10 −6 Hz, and the values of c QP P hot n given in the inset.The solid lines are calculated using Eq.(42).

FIG. 7 .
FIG.7.Distribution function in the crossover from perturbative to heating regimes as identified by Eq.(47).We use the same parameters as in Fig.6except for the photon number n, given in the bottom left inset in terms of the ratio between T * of Eq. (43) and ξ of Eq. (17).The top right inset shows the second peak at energies ϵ > ωP B , with the blue dashed curve corresponding to Eq. (48) where f (ϵ) is the numerically calculated distribution function at low energies.

FIG. 8 .
FIG.8.Quasiparticle density for fixed numbers of pairbreaking photons (see parameter in the inset) as function of temperature in the low-heating regime in which non-pairbreaking photons have negligible influence on the density.Results calculated from numerical solutions to the kinetic equation are displayed with full lines and the corresponding analytical ones obtained by solving Eq. (49) in the steady state with dashed lines.

FIG. 10 .
FIG.10.Low-temperature Q-factor vs readout power.Experimental data points are obtained by averaging the low temperature (TB < 150 mK) data in Ref.[18] and each error bar shows the corresponding standard deviation; the value of the data points coincide with those of Qi,ext [cf.Eq. (54)] used in Fig.9.The solid line is obtained from numerical solutions to the kinetic equation ignoring thermal phonons (that is, setting TB = 0) and using parameters ωP B = 2.8∆, c QP P hot,P B = 27c QP P hot , τ0 = 20 ns, nP B,0 = 2 × 10 −4 , and µ = 10 −9 [see Eq. (56)]; other parameters as in Fig.9.The dashed line shows the analytical estimate, Eq. (55).