Bogoliubov Quasiparticles in Superconducting Qubits

Extending the qubit coherence times is a crucial task in building quantum information processing devices. In the three-dimensional cavity implementations of circuit QED, the coherence of superconducting qubits was improved dramatically due to cutting the losses associated with the photon emission. Next frontier in improving the coherence includes the mitigation of the adverse effects of superconducting quasiparticles. In these lectures, we review the basics of the quasiparticles dynamics, their interaction with the qubit degree of freedom, their contribution to the qubit relaxation rates, and approaches to control their effect.

1 Superconductivity in an isolated metallic island 1

.1 Electron pairing and condensate
Exposition of the theory of superconductivity conventionally starts with considering electrons freely propagating as plane waves in an ideal, translationally-invariant medium [1]. The electron energy spectrum is then continuous. The number of electron states per unit volume per unit energy, usually referred to as the density of states, is a function of energy, with some finite value ν 0 at the Fermi level. To give a bit different perspective on the subject, let us consider, instead, a medium confined to some large (in units of Fermi wavelength) box and containing finite density of impurities which scatter electrons. Confinement to the box renders electron spectrum discrete, while scattering off randomly-positioned impurities would remove any accidental degeneracy of the levels. Under these conditions, the average density of energy levels ξ α of one-electron states α in the vicinity of the Fermi energy is Vν 0 , with V being the volume of the box. The typical spacing between the adjacent energy levels then is δ = 1/Vν 0 . Taking for a crude estimate ν 0 = 1 (eV·Å 3 ) −1 , we find for an island of volume V = 10 −2 µm 3 the average density of levels 10 10 eV −1 , yielding a tiny level spacing of δ = 10 −10 eV ≈ 1 µK. Hereinafter the term "average" means average over an energy interval which includes many levels, but still very small compared to the Fermi energy E F (typically a few eV in a conventional metal). The Kramers theorem indicates that in the absence of magnetization each discrete electron level in a normal-metal island is doubly-degenerate, forming a Kramers pair. This statement is unaffected by the spin-orbit coupling, as it does not break the time-reversal symmetry. For simplicity, however, we will dispense with the spin-orbit coupling and associate the pair of states with the spin-up and spin-down electrons having the same orbital part of the wave function ψ n (r) (this is an excellent approximation for light elements, such as Al). In terms of these states, the second-quantized form of the Hamiltonian is (for brevity, we do not include the spin-triplet channel for the interaction which does not change the conclusions) Here operators c † nσ and c nσ create and annihilate electrons with energies ξ n (measured from the Fermi level), and H klmn = dr 1 dr 2 V (r 1 − r 2 )ψ k (r 1 )ψ l (r 2 )ψ m (r 2 )ψ n (r 1 ) (2) are the matrix elements of interaction, written in terms of the single-particle eigenfunctions. These are strongly oscillating in space, and there is little correlation between the oscillations of the wavefunctions of different states. As a result, there is a strong hierarchy in the matrix elements H klmn : those with pairwise-equal indices are by far the largest ones. We will illustrate it using an example of a contact interaction, V (r) = (λ/ν 0 )δ(r), characterised by dimensionless interaction constant λ. In this example, the double-integral in the right-hand side of Eq. (2) is reduced to an integral over a single variable r with the integrand (λ/ν 0 )ψ k (r)ψ l (r)ψ m (r)ψ n (r). For generic k, l, m, n the product of wave functions is rapidly oscillating as a function or r thus suppressing the value of the integral (the characteristic length scale for the oscillations is set by the Fermi wavelength) and making it zero on average. Having k = n, l = m or k = m, l = n reduces the product of wavefunctions to |ψ k (r)| 2 |ψ l (r)| 2 which is non-negative, no matter if the one-particle wave functions real-or complex-valued. A non-negative integrand leads to matrix elements H kllk and H klkl which only weakly depend on k and l, having non-zero average ∼ λδ . In the presence of a magnetic flux piercing the island wavefunctions the timereversal symmetry is broken, and ψ k (r) are complex-valued. To the contrary, time-reversal symmetry allows one to choose real-valued eigenfunctions ψ n (r) = ψ n (r). That brings yet another paring, k = l, m = n, yielding a non-negative product [ψ k (r)] 2 [ψ m (r)] 2 . The said three types of parings correspond, respectively, to the Hartree term, Fock term, and the Bardeen-Cooper-Schrieffer (BCS) term. These three interaction types are the leading ones, regardless the details of V (r), including its range and sign. Accounting for the Coulomb long-range component of V (r) generates the charging energy out of the Hartree term, while the Fock term induces exchange interaction (which is safe to ignore in the case of a nonmagnetic material); these two interactions are insensitive to breaking the time-reversal symmetry. The BCS term is responsible for the formation of a superconducting state, once V (r) contains a short-range attraction component. Therefore, neglecting the exchange interaction and levelto-level fluctuations, the interaction term in the island Hamiltonian Eq. (1) takes a universal form, independent on the details of the electron wavefunctions in the island. Herê is the operator of the number of electrons, and accounting for allows one to consider superconductivity in case of the attractive interaction between electrons with energies within some range |ξ n | < ω D (for the phonon mechanism of superconductivity, ω D is of the order of phonon Debye frequency). The superconducting phase transition is associated with the appearance of a macroscopically-large value of Ô †Ô defeating the smallness of the factor λδ in Eq. (3).
The electron number N e is conserved in an isolated island, so the included in Eq. (3) polarization charge N g for now reflects only the level from which all energies are measured. Now we consider fixed even N e and therefore fixed charging energy represented by the first term in Eq. (3), and concentrate on the ground state of an isolated superconducting island described by Hamiltonian The termÔ †Ô in Eqs. (3), (6) is the counterpart of the BCS interaction term conventionally written [1] in the basis of plane waves, Either side of Eq. (7) preserves the total electron number and describes coupling involving large number of singlet pairs: in the case of an island, a pair on a level n is coupled with ∼ ω D /δ other pair states labelled by m. Such type of coupling provides a motivation for applying a mean-field treatment for determining the ground-state energy and thermodynamics of the system. In the mean-field approximation, one introduces the average, to replace the quartic term in Eq. (6) by a bilinear one. After the simplified Hamiltonian, is diagonalized, one evaluates the average . . . in the right-hand side of Eq. (8) in terms of ∆, forming this way a self-consistency equation for ∆. This routine for an island is essentially identical to the one for a bulk superconductor. The bilinear mean-field Hamiltonian is diagonalized by the Bogoliubov transformation, Here γ † nσ , γ nσ are creation and annihilation operators for quasiparticle excitations with spin σ =↑, ↓. The Bogoliubov amplitudes are complex numbers; for convenience, we may define gauge by taking u n = u * n , v n = |v n |e iϕ with ϕ being the phase of the order parameter, ∆ = |∆|e iϕ . To preserve canonical commutation relations, their magnitudes satisfy the constraint where n = ξ 2 n + |∆| 2 is the quasiparticle excitation energy. The ground state |GS of the mean-field Hamiltonian is defined by the condition γ nσ |GS = 0. The order parameter ∆ is found self-consistently as Please note that the right-hand side here remains finite in the macroscopic limit δ → 0. Finally, the Hamiltonian in Eq. (9) is transformed into the Hamiltonian for quasiparticle excitations: Thermal averages present in Eq. (13) are evaluated over the Gibbs ensemble with the Hamiltonian Eq. (14). The self-consistency equation defines the absolute value of the order parameter |∆(T )|, leaving its phase ϕ arbitrary; the ground-state energy, excitation spectrum, and thermodynamic potential of an island are independent of ϕ. The zero-temperature solution of Eq. (13) yields |∆| ≈ 2 ω D exp(−1/λ). It remains finite in the limit δ → 0; for islands of a typical size, δ |∆| ω D . The ground-state wave function of the mean-field Hamiltonian (9) with a given ∆ is where |0 is the vacuum for electronic excitations, c nσ |0 = 0. The subscript ϕ indicates that the phase of ∆ enters this definition via the Bogoliubov amplitudes. One can verify that this expression satisfies the condition γ nσ |ψ ϕ = 0 defining the ground state for quasiparticle excitations (cf. Sec. 2.3). Clearly, the defined by Eq. (15) functions are 2π-periodic: |ψ ϕ = |ψ ϕ+2π . While being an eigenstate of the BCS Hamiltonian, |ψ ϕ is not an eigenfunction of the electron number. It is rather a coherent superposition of states with different numbers of electron pairs, so the number of pairs is not defined. The ground-state energy and the excitations spectrum are independent of ϕ, which provides a relief: out of |ψ ϕ functions, we may form a linear combination corresponding to a definite number of electron pairs N P . The relation (16) is gauge-invariant (i.e., invariant with respect to an arbitrary phase shift, ϕ → ϕ + ϕ 0 ). The wave function (16) is an excellent approximation to the ground state of the Hamiltonian (6), which conserves the electron number. The associated with the finite level spacing δ corrections to ∆ of Eq. (13) and to the corresponding ground-state energy of Hamiltonian (6) scale to zero proportionally to δ with the increase of island volume (see [2] for details and further references). The condensate wavefunctions |ψ ϕ and |ψ N P form two bases in the Hilbert space of many-body paired electron states. We may view the N P and ϕ representations as dual ones, similar tox andp representations in the single-particle quantum mechanics. (The most important difference is that ϕ varies between 0 and 2π, making it a compact variable.) In the N P representation, the operator of number of electron pairsN (measured from some large integer corresponding to the filled Fermi sea in the island) acts as a multiplication operator, N = N ·. Now we establish its form in the ϕ representation: That is,N = −id/dϕ. It is important to remember that the functions |ψ ϕ are 2π-periodic, so the spectrum of −id/dϕ is the set of integers (N = 0 means no extra electron pairs on the island). Conversely, the operatorT increasing the number of pairs by 1 is a multiplication operator in the ϕ-representation: so thatT = e −iϕ . Therefore in the space of states we considered here, variableN is a conjugate to the compact variableφ. The two satisfy the appropriate canonical commutation relation invariant with respect to the basis.

Thermodynamics of a superconducting island
The electron condensate in an isolated island accommodates an even number of particles. If the number of electrons on the island is even, they all reside in the condensate in a T = 0 equilibrium state. Under the same conditions, an odd electron in the island does not have a pair and occupies the lowest-energy quasiparticle state, thus raising the energy of the island by |∆|. At higher temperatures, ionization of the Cooper pairs results in a higher number of equilibrium quasiparticles, diminishing this even-odd effect. To see this we evaluate and compare the partition functions Z 0 and Z 1 for the even and odd numbers of electrons, respectively [3]. In the "even" case the states of the island are parametrized by the number 0, 2, 4, . . . of quasiparticles and their quantum numbers, so we write: (hereinafter we disregard the negligible probability of double-occupancy of any state). The series here is easy to sum up: Similarly in the "odd" case the island contains 1, 3, 5, . . . quasiparticles, and the partition function equals One can easily recognize n qp = 2z(T, δ , ∆)/V as the quasiparticle density in the bulk at T ∆ (the factor of 2 accounting for spin). It is convenient to normalize n qp by the "density of Cooper pairs" n CP , n qp = n CP x qp , n CP = 2ν 0 ∆.
In equilibrium, x qp = 2πT /|∆| exp(−|∆|/T ). The difference between the thermodynamic potentials of the even and odd states, namely T ln(Z 0 /Z 1 ), becomes substantial (order-of-∆) once on average there is less than one thermally excited quasiparticle on an island, i.e., z 1. This happens at T below the scale set by the V-dependent Cooper pair ionization temperature The logarithm in the denominator here is pretty big, it is about 14 for an Al island of a typical volume V = 10 −2 µm 3 . As a result, one expects -on the grounds of thermodynamics -no broken Cooper pairs in such island at T T = 0.14K. In our example, we find x qp ≈ 1.9·10 −7 at T = T and may expect a minuscule x qp ≈ 2.1 · 10 −23 at a typical for qubit experiments temperature T = 40mK. However, numerous measurements find x qp = 10 −7 − 10 −5 at these temperatures. The origin of the excess quasiparticles is not known and remains under scrutiny. Meanwhile, it is worth assessing how harmful they are for the qubits operation and look for ways to mitigate their unwanted effects.
This Hamiltonian captures the coherent, non-dissipative tunneling of pairs of electrons; each term of the sum corresponds to transfer of n Cooper pairs in a single tunneling event. For a conventional tunnel junction, electron transmission coefficient is small, so one may safely keep only the lowest-order term (n = 1) in the sum. Furthermore, time-reversal symmetry for tunneling through a non-magnetic insulator dictates C * 1 = C 1 . Using the phase representation (ϕ L , ϕ R ) for the operators T L and T R and omitting the phase-independent const term, we obtain The connected islands at ϕ = 0 constrain the motion of a Cooper pair less than each island separately, so the ground-state energy of the entire system is reduced by the link; it means that in Eq. (27) the only remaining coefficient C 1 < 0, and therefore E J > 0 in Eq. (28). We note in passing that the Josephson energy E J and the normal-state conductance of the junction G do not have to be small compared, respectively, to ∆ and e 2 /h, as the smallness of transmission coefficient may be compensated by a large area of the junction. Later on, we will evaluate E J microscopically and relate it to G. One more remark is due here: we tacitly assumed that the ground state of the system is non-degenerate. Tunneling via a quantum dot carrying an uncompensated electron spin provides a counter-example, as the Kramers degeneracy is preserved at sufficiently weak tunneling [4]. The presence of the localized spin results in E J < 0, so that the lowest energy of the junction is reached [4] at ϕ = π. Formation of a π-junction in tunneling through a quantum dot was demonstrated, e.g., in Ref. [5].
Transfer of N Cooper pairs across the junction creates a charge dipole between the islands. The corresponding electrostatic energy [cf. Eq. (3)] in terms of the operatorN = (1/i)d/dϕ, reads Here charging energy E C = e 2 /2C takes into account the junction capacitance as well as any capacitance shunting the junction; n g is the static charge (in units of 2e) induced by a biasing gate, background charges, and unpaired electrons. Out of the three contributions only the first one is controllable; the two others fluctuate on some large time scale. The contribution stemming from unpaired electrons is discrete, and changes by ±1/2 upon a quasiparticle tunneling across the junction; the background charge may vary continuously. A single-junction qubit is described by the Hamiltonian H J + H C acting in the space of periodic functions, ψ(ϕ) = ψ(ϕ + 2π). Clearly, its spectrum is discrete, non-equidistant, and depends on n g periodically with period 1. The n g -dependence is detrimental for the qubit coherence, due to the uncontrolled variations of n g . In transmons [6], the unwanted sensitivity to n g is countered by increasing the ratio E J /E C . At E J /E C 1, one may separate the quantum dynamics of the phase difference ϕ into small fluctuations around the minima (ϕ = 2πn with integer n) and discrete phase slips by ±2π. The former correspond to the dynamics of an anharmonic oscillator having non-equidistant levels needed for a qubit operation. The latter brings the unwanted sensitivity of the levels (∝ δε n cos 2πn g ) to the uncontrolled variations of n g . The probability amplitude of a phase slip is exponentially small at E J /E C 1, δε n ∝ exp(− 8E J /E C ). This allows one to effectively suppress the influence of n g without affecting the qubit energy levels (the relative anharmonicity α r scales as a power law of E C /E J , α r E C /8E J , see Ref. [6]). A wide variety of experiments demonstrated the prominence of discrete ±1/2 jumps (commonly referred to as e-jumps) in the spectrum of fluctuations of n g . Its average value, n g , can be controlled by a gate electrode in a properly-designed transmon. There are two special values of the gate voltage for which n g = 1/4 or 3/4, making the transmon energy levels insensitive to the ±1/2 jumps. To see this, let us consider, e.g., n g = 1/4. An e-jump changes this initial value of n g to n g = 3/4. Due to the periodicity of the spectrum with n g , the energy levels of the qubit Hamiltonians with n g = 3/4 and n g = −1/4 in Eq. (29) are identical to each other. Lastly, we may change ϕ → −ϕ, as the Josephson energy Eq. (28) is even in ϕ; this change would return the charging energy Hamiltonian after an e-jump to its initial form prior to the jump (n g = 1/4). The data for the qubit transition frequency, accumulated over a large series of sequential measurements, clearly shows the reduced sensitivity to the e-jumps at the said special values of n g , see Fig. 1(a) and Fig. 1(b).
In a fluxonium qubit [7], the protection is achieved by shunting the junction with a highinductance loop. The loop -superinductor -is actually a chain of n s ∼ 100 Josephson junctions with sufficiently large E s J /E s C so that the phase slips probability amplitude in the superinductor is negligible. This allows one to dispense with the periodicity of ψ(ϕ) function and approximate the Hamiltonian of the superinductor by Figure 1: Panel (a) copied from Ref. [8]. Spectroscopy of a qubit as a function of gate-induced charge n g . For each pixel, a Gaussian pulse (σ = 20 ns, corresponding to a pulse on resonance) is applied at the indicated frequency and the qubit is immediately measured. Each pixel is an average of 5000 repetitions (50 ms). Darker pixels correspond to higher homodyne readout voltages that are proportional to the probability of the qubit in the excited state. An "eye"-shaped pattern indicates charge-e jumps associated with the tunneling of nonequilibrium quasiparticles. Panel (b) copied from Ref. [9]. Normalized two-tone spectroscopy measurements of the 0 → 1 transition versus the offset charge.
The inductive energy E L = E s J /n s = (Φ 0 /2π) 2 /L, with Φ 0 = π c/e being the quantum of flux, accounts for the loop inductance L; we also allowed for an external flux Φ e threading the loop formed by the superinductor and the qubit junction. The form of the Hamiltonian (30) dictates the boundary conditions for the wave function ψ(ϕ → ±∞) = 0. This, in turn, allows one to eliminate n g from the Hamiltonian H J + H I + H C , by a simple gauge transformation. The inclusion of the shunt makes the energy levels of this Hamiltonian independent of n g , while allowing for their control by Φ e .
To summarize, Hamiltonian describes a wide variety of superconducting qubits. In the absence of the superinductor (E L = 0) Hamiltonian (31) acts in the space of periodic functions; if E L = 0 then the proper boudary condtion is ψ(±∞) = 0, which allows one to gauge out the n g -dependence.
Hamiltonian (31) acts in the low-energy subspace, meaning that it is good for describing energy levels well below the quasiparticle continuum; that, in turn, sets the requirement E J , E C , E L ∆. In metallic islands, screening length is very short (it is about the interatomic distance) and the energy of plasmon is extremely high (typically of the order of Fermi energy). As a result, quantum fluctuations of chargeN governed by Eq. (31) lead merely to the fluctuations of the potential of an entire island; the associated with the fluctuations electric fields do not penetrate the bulk of an island. Fluctuations of the potential are benign for the gaugeinvariant Hamiltonian (6) and do not affect its excited states and their occupations. Therefore, as long as quasiparticles do not tunnel and thus are not exposed to the potential difference between the islands, their presence is inconsequential to the dynamics of the qubit degree of

Tunneling Hamiltonian and the normal-state conductance of a junction
Consider two normal-state leads separated by an insulating barrier. Electrons can tunnel through the barrier, and we model this system using the Hamiltonian: where H α , α = L, R are the Hamiltonians of the left/right lead, and is the tunneling Hamiltonian, describing transfer of an electron from a state n R in the right lead to a state n L in the left one, and a transfer in the opposite direction; t n L n R is the corresponding tunneling matrix element. If voltages V α are applied to the leads, their Hamiltonians take the form where H α and the number operators N e α are given by the proper generalizations of Eqs. (1) and (4) to two leads. In the absence of tunneling, the particle number is conserved separately for the two leads, [H α , N e α ] = 0. Tunneling allows the current to flow through a voltage-biased junction. This current is dissipative in the case of normal leads. There are many ways to evaluate it, below we present one of them.
It is convenient to change the representation, so that the excitation energies are measured from the respective electrochemical potential of each lead. For a generic, time-dependent unitary transformation U (t), the transformed HamiltonianH is given bỹ Here we take U in the form (the specific value of constants φ L,R here is inconsequential). Clearly, U commutes with the lead Hamiltonians, In the used representation, this is a periodic-in-time perturbation of the HamiltonianH L +H R , which results in the absorption of energy quanta Ω = eV by the system. To find the absorption power P in the weak-tunneling limit, we may use the Born approximation for the tunneling amplitudes and then apply Fermi's Golden Rule to evaluate the transition rate; lastly, we multiply it by the energy quantum eV to obtain: We replaced the averages of c † n c m over the Gibbs distributions with the Hamiltonians H L,R by the respective Fermi functions f F (ξ n L ,n R ). This required neglecting the electron-electron interaction in the leads, which is fine for the conventional tunnel junctions between normalstate metallic electrodes; the extra factor of 2 in Eq. (40) accounts for the summation over the spin variable. One may worry that we applied a formalism developed for isolated islands to two leads attached to a voltage source. Indeed, the charge transfer in the DC regime between two otherwise isolated islands is not sustainable over an arbitrarily long time. On a technical level, the difficulty arises in the derivation of Fermi's Golden Rule, which we used in Eq. (39), for the transitions between the discrete spectrum levels, as one would end up with interlevel Rabi oscillations instead. A formal way around this difficulty is to consider a slightly broadened tone Ω which would compensate for the discreteness of the energies ξ n L,R and allow one to use the standard derivation [11] of the Fermi's Golden Rule. (That derivation involves the consideration of the initial regime of linear growth of the perturbation [11]). After that, one takes the limit V L,R → ∞ keeping the products V L V R |t n L ,n R | 2 finite, and then returns to a fixed Ω = eV . To proceed with the derivation of DC conductance G N , we use the relation (we prefer the latter form of this relation for its compactness). In derivation of Eq.(40) we also assumed low temperature, T E F and denoted an average over the states n L , n R close to Fermi energy by |t n L ,n R | 2 . In an alternative standard derivation of Eq. (40), one uses momentum eigenstates for two clean infinite-size leads, see e.g. [12].
The microscopic properties of the junction are encoded in the matrix elements t n L ,n R of the tunneling Hamiltonian and abbreviated to a single parameter, conductance, by Eq. (40). The same value of G N can be achieved by enlarging the cross-sectional area Σ of the junction or by increasing the transmission coefficient |t B | 2 of the tunnel barrier. The conductance of a largearea junction can be estimated as G N ∼ (e 2 /h)(Σ/λ 2 F )|t B | 2 . The second factor here represents the cross-sectional area of the junction in units of the Fermi-wavelength-squared and can be viewed as the (large) number of independent electron modes fitting into the junction's area. Equation (40) will allow us to express various quantities of interest in the superconducting state in terms of the normal conductance G N .

Josephson energy and current
We define current operatorÎ as a time derivative of eN e R , At zero bias, the phase difference entering in Eq.
. We may introduce the "superconducting" phase difference ϕ = 2(φ L − φ R ) and, by inspecting Eqs. (37) and (41), establish the relationÎ = (2e/ )∂H T /∂ϕ between the current operator and tunneling Hamiltonian. Next we may account forH L and H R being independent of ϕ, in order to arrive at the exact relation,Î = (2e/ )∂H/∂ϕ, between the current operator and the full Hamiltonian (35). In equilibrium, i.e., with no bias applied, we may average the latter relation over the Gibbs distribution for the entire system and find for the current We removed the tilde sign in Eq. (42), as the trace there is gauge-invariant. A finite phase difference ϕ for the electron states across the barrier may be introduced by including the junction in a conducting loop and threading a magnetic flux through it. This causes a current running in the loop in equilibrium (known as persistent current) even if the entire system is in the normal state. This mesoscopic effect, however, vanishes in the limit δ → 0 and turns out quite hard to measure in normal-metal rings [13]. On the contrary, in the superconducting state, the Josephson current (42) remains finite at δ → 0, i.e., in the limit of macroscopic leads.
We are interested in the Josephson current I J (ϕ) and energy at temperatures T ∆, so we may replace the free energy F (ϕ) with the ϕ-dependent part δE G (ϕ) of the groundstate energy. We will evaluate δE G perturbatively to the lowest non-vanishing order (t 2 n L ,n R ) in tunneling and within the BCS mean field approximation. In the "tilde" basis the phase dependence is delegated toH T , so the defined in Eq. (8) mean-field order parameters of the leads are purely real, ∆ → |∆ L,R |. Alternatively, we may use the original basis, in which case the ϕ L,R dependencies are carried by ∆ L,R . For definiteness, we use the latter gauge.
We now write the tunneling Hamiltonian H T , Eq. (33), in terms of the quasiparticle operators; using the Bogoliubov transformation, Eq. (10), we find where for simplicity we assumed no bias (V = 0), and the presence of time-reversal symmetry resulting in real-valued tunneling amplitudes, t n L n R = t * n L n R ; we also adjusted the gauge, multiplying u n and v n by e −iϕ/2 compared to the definitions in Section 1.1. Here we have separated the terms accounting for single quasiparticle tunneling, H qp T , from that describing pair breaking and recombination, H p T . It is convenient to show explicitly the phase dependence of the Bogoliubov amplitudes: Obviously, the zeroth-order average H T 0 = 0; in the next order one finds where the sum is over all possible excited states |λ of energy E λ determined from the proper generalization of quasiparticle Hamiltonian Eq. (14) onto two leads. The only non-zero contribution comes from the first term in the RHS of Eq. (45). Then |λ = |n L n R σ ≡ γ † n L σ γ † n Rσ |GS and E λ = n L + n R . Evaluation of Eq. (47) yields: Trading the summations here for the integration over the energies over the corresponding states and dispensing with the dependence of the tunneling matrix elements on n L , n R we find The latter simplification made the integral in Eq. (49) ultraviolet-divergent; the dependence of t on n L and n R provides one with a model-dependent cut-off at some energies of the order of E F . The meaningful in the context of superconductivity part of E 0 (∆ L , ∆ R ), however, is model-independent. It can be evaluated with the help of the following regularization, which makes the integral convergent at L,R ∼ ∆ L,R and provides a way to express const in Eq. (27) in terms of ∆ L,R and G N . The integral in the Josephson energy (50) is convergent.
The applicability of the perturbation theory used in the derivation of Eq. (51) requires the smallness of the matrix elements t n L n R , which in turn means that the transmission coefficient |t B | 2 across the tunnel barrier must be small. As we discussed at the end of Section 2.2, the small factor |t B | 2 1 in the conductance G N can be compensated by a large number of electron modes Σ/λ 2 F fitting in the junction's cross-sectional area. A similar compensation happens also for the Josephson energy Eq. (51). In fact, it can be expressed in terms of G N and ∆. Comparing Eq. (51) with Eq. (40) and using Eq. (42) at T → 0 we conclude: The expression of E J in terms of ∆ and G N may be viewed as a version of the Ambegaokar-Baratoff formula. We may convert it to the familiar [1] form, eI c R N = (π/2)∆, by introducing the critical current of the junction I c = max ϕ {I J (ϕ)} and its resistance R N = 1/G N in the normal state.

Real part of the AC admittance of a junction
Josephson energy E J is one of the main parameters defining the energy spectrum of an ideal qubit, see Eq. (31). In the previous Section, we related E J to the normal-state conductance of the junction. Dissipation in the Josephson junction is one of the factors limiting the qubit coherence. To quantify the dissipation, this Section develops the theory of the AC admittance of a junction; we will focus on its real part. Similar to our discussion of the DC conductance in Section 2.2, the shortest way to get the dissipative part of AC admittance, Re Y (ω, ϕ), is to evaluate the absorption power P of bias V (t) = V cos(ωt) applied across the junction. In contrast to a fixed bias, the alternating one does not cause winding of the phase difference φ L (t)−φ R (t) in Eq. (37). Instead, this difference exhibits small oscillations around a finite value, The oscillatory term in the expansion drives the transitions in which energy quanta ω are absorbed. For finding P in the second order in V and second order in t n L n R we may apply Fermi's Golden Rule to the perturbation −[2eV sin(ωt)/ ω]∂ ϕHT (ϕ) and proceed similar to the derivation of Eq. (39). Next, we find the dissipative part of admittance by casting the result of calculation in the form This program works equally well for the tunnel junctions with normal or superconducting leads. In the former case, the result is independent of ϕ. We find Re Y (ω) = G N for a junction between two normal leads with energy-independent electron density of states, cf. Eq. (40). The limits ω → 0 and V → 0 of the absorbed power do not commute in general; this is exemplified by the evaluation of the DC dissipative conductance of an SNS or Josephson junction [14]. In the following we concentrate on the case of a finite ω.
Considering superconducting leads, we use Eqs. (43)- (45) in order to derive the appropriate form of the perturbation −[2eV sin(ωt)/ ω]∂ ϕHT (ϕ) in terms of the quasiparticles creation and annihilation operators. This is achieved by replacing (44); the gauge invariance of the observables allows us to assign all the phase dependence to lead L. The operator structure of the H p T and H qp T parts of the perturbation is quite different: the former one describes creation or annihilation of pairs of quasiparticles, while the latter term corresponds to the quasiparticle tunneling. Therefore, absorption processes originating in H p T are effective only at frequencies exceeding the threshold, ω > 2∆. Ultimately, we are interested in the interaction of quasiparticles with the qubit degrees of freedom evolving with frequencies well below ∆/ . Therefore, we focus on the terms in P stemming from H qp T , Power P then is found by multiplying the energy quantum by the transition rate: Double angular brackets denote averaging over the initial (not necessarily equilibrium) quasiparticles states |η with energy E η ; sum is over possible final quasiparticle (qp) states |λ with energy E λ . The two terms in square brackets account for energy absorbed and emitted by quasiparticles, respectively. Performing the averaging and using Eqs. (40) and (53) we find for the dissipative admittance of the Josephson junction: The quasiparticle distribution functions here do not have to be equilibrium ones; they merely represent occupation factors of various energy states. An assumption of L/R symmetry allows us to simplify Eq. (56): The ϕ-dependence here comes from the interference between two processes of charge-e transfer across the barrier: the first one consists of forwarding an electron as a quasiparticle across the barrier; the second process involves forwarding a Cooper pair accompanied by returning an electron. The involvement of the condensate makes the result of interference phase-dependent; the closer the quasiparticle energy to the gap, the stronger the relative effect of interference, cf. the numerator of the integrand in Eq. (57). Now assume that the microwave energy quantum ω and the characteristic value T eff of quasiparticle energy measured from the gap edge, − ∆, are small, ω, T eff ∆. Then, by changing the variable of integration in Eq. (57) via = ∆(1 + x) we find Under an additional assumption T eff ω, the quasiparticle admittance here is The dimensionless quasiparticle density is introduced here consistently with Eq. (25); it represents the density of quasiparticles n qp = n CP x qp normalized by the density of Cooper pairs n CP = 2ν 0 ∆ and assumes symmetric junction (ν 0 = ν L = ν R ).
The factor x qp in Eq. (59) accounts for the fact that all low-energy quasiparticles can absorb energy and hence contribute to dissipation. In equilibrium, x qp is controlled solely by the ratio T /∆, as the chemical potential µ of quasiparticles is pinned to zero. The simplest model of a non-equilibrium distribution allows for some µ = 0 and effective temperature T eff , and treats x qp and T eff as independent parameters. Leaving a more detailed discussion of the quasiparticle kinetics for Section 5, we mention here that there are reasons to expect T eff = T , as a single quasiparticle may relax its energy by emitting a phonon. On the contrary, in order to recombine that quasiparticle has to meet another existing quasiparticle. Therefore, the recombination rate is smaller by a factor x qp 1 than the rate of quasiparticle energy relaxation.
Assuming the quasiparticle distribution is described by a Fermi function with some µ and effective temperature T eff , the imaginary part of the quasiparticle admittance Y qp can be recovered using Kramers-Krönig transform. This would account for the effect of itinerant quasiparticles, but miss the main contribution to the imaginary part of the junction admittance Y J , which originates from the response of the condensate (or, equivalently, from the contribution of Andreev bound states). We refer to Ref. [10] for a discussion of these points. In the next section we will relate the transition rates in superconducting qubits to the admittance of the Josephson junction.

Qubit interaction with quasiparticles
For a Josephson junction shunted by an inductive loop, Fig. 2, the low-energy effective Hamiltonian can be written as The first term determines the dynamics of the phase degree of freedom in the absence of quasiparticles, see Eq. (31). The contribution from pair tunneling, Eq. (45), is taken into account by the E J term in Eq. (31).

The second term in Eq. (61) is the sum of the BCS Hamiltonians for quasiparticles in the leads
with H α qp of Eq. (14). The last term in Eq. (61) is the single quasiparticle tunneling Hamiltonian defined in Eq. (44) with a caveat: now the phase entering in H qp T is an operator, This way, H qp T becomes a Hamiltonian of the quasiparticles-qubit interaction. This way of accounting for the interaction is fine, as long as the dynamics of the condensate involves frequencies ω/2π much smaller that ∆/π . Fortunately, this is the case for the typical devices controlled by microwaves (with frequency 10 GHz, while ∆/π ∼ 100 GHz in Al).
Within the described model, we can calculate the transition rate Γ if between qubit states |i and |f (i.e., eigenstates of the Hamiltonian H ϕ ) associated with tunneling of a quasiparticle across the junction similarly to the calculation of admittance in Sec. 2.4. That is, we treat H qp T as a perturbation and evaluate the transition rates using Fermi's Golden Rule: where ω if = E i − E f is the difference between the energies of the two qubit states. After averaging over initial quasiparticle states |η and summing over final quasiparticle states |λ , and assuming |∆ L | = |∆ R | ≡ ∆, we find It is important to note that the transitions between the qubit states are accompanied by the charge-e transfer across the junction. In some cases, see Section 3.5, that helps one to single out the qubit transitions driven by quasiparticles. For generic states and flux bias, the matrix elements of sin(ϕ/2) and cos(ϕ/2) have similar orders of magnitude; then, at low effective temperature of quasiparticles (T eff ∆) the cos ϕ/2 contribution to Eq. (65) is suppressed by a small factor ∼ ω if /∆. Important exceptions are (quasi)elastic transitions (see Sec. 3.5) and transitions at special values of flux bias fine-tuned to suppress the sin ϕ/2 matrix element (see Sec. 3.4).

Qubit energy relaxation
Here we focus on the relaxation rate from the first excited, |i = |1 , to the ground, |f = |0 , state in a generic setting, therefore neglecting the cos ϕ/2 contribution, see the last line in Eq. (65).
With the assumptions discussed above, the qubit relaxation rate due to quasiparticle tunneling can be expressed as where is the quasiparticle current spectral density, see the first two lines in Eq. (65). The quasiparticle states occupation factors f ( ) typically are small for all allowed energies, whether quasiparticles are at equilibrium or not, f ( ) 1. This allows us to replace 1 − f ( + ω) → 1 in Eq. (67). Then, at low effective temperatures and frequencies, T eff , ω ∆, the quasiparticle current spectral density takes the form By comparing this formula to Eq. (59), we find the relation For a quasi-equilibrium distribution of quasiparticles with some x qp and T eff , it may be viewed as a particular case of more general fluctuation-dissipation relations [10]. That allows one to conclude that Γ 01 = Γ 10 · e − ω 10 /T eff indicating that Γ 10 > Γ 01 , as T eff > 0 (there are no reasons to expect an inversion in the energy distribution of quasiparticles). Relations (66) and (68) may allow one to link the magnitude of the junction dissipation (the real part of Y qp ) to the qubit relaxation rate. In the next sections we explore similarities and differences between the phase dependence of the admittance, see Eq. (58), and the variation of the relaxation rate with the external flux which controls the qubit and implicitly enters in Eq. (66).

Energy relaxation of a weakly anharmonic qubit
The last two terms (the "potential energy") in Eq. (31) in general possesses multiple minima, whose positions ϕ 0 are solutions of So long as the external flux is tuned away from half-integer multiples of the flux quantum and phase fluctuations are small, we can treat the potential energy in the harmonic approximation, The assumption of small phase fluctuations corresponds to the condition where is the qubit frequency in the harmonic approximation. For the transmon [6], E L = 0 and ϕ 0 = 0, the condition (73) corresponds to the requirement of a large ratio between Josephson and charging energy, E J /E C 1, which also enables us to neglect the dimensionless voltage n g of Eq. (31).
Within the harmonic approximation, it is straightforward to calculate the matrix element in Eq. (66) by expanding sin ϕ/2 to linear order around ϕ 0 . This way we find 0| sin Substituting this expression into Eq. (66), and using Eq. (69) and the definition of the charging energy, we arrive at Therefore, the qubit relaxation rate is given by the inverse of the classical RC time of the junction, where 1/R is identified with the real part of its flux-dependent admittance [ϕ 0 is identified with the phase difference ϕ in Eq. (58)]. This result seems to indicate that, as it is often the case, the behavior of the quantum harmonic oscillator is analogous to that of its classical counterpart. However, the analogy has its limitations, set by the form of the perturbation causing the relaxation: the selection rules for matrix elements of an operator sinφ/2 are less restrictive than forφ. That opens a possibility, e.g., of a direct decay from the second level to the ground state; the corresponding rate is [10] The dependence on phase/flux in this expression clearly differs from that in Eq. (58). Moreover, we remind that Eq. (76) is restricted to fluxes away from half-integer multiples of the flux quantum, so the relation to Eq. (58) does not necessarily hold at arbitrary flux -we explore this issue further in the next section.

The cos ϕ problem
The result for the real part of the junction admittance, Eq. (58), shows that as the phase difference approaches π, dissipation is suppressed. This is a manifestation of quantum mechanical interference: a quasiparticle is a coherent superposition of electron and hole-like excitations [cf. Eq. (10)], and these two components interfere during a tunneling event in a way that can preclude the absorption of energy. The exact cancellation at ϕ = π is expected only in the limit of small temperature; more generally, one can write where ε → 1 as T → 0, see [15] and Ch. 2.6 in Ref. [12]. In the latter reference, experimental attempts to determine ε in the 1970s are summarized, and the discrepancy between theory and experiments was termed "the cos ϕ problem". It is interesting to consider the behavior of the admittance when the junction is part of a loop, so that the flux Φ biases the junction, ϕ = 2πf with f = Φ/Φ 0 . Expanding Eq. (78) around f = 1/2 we find A fluxonium qubit consists of a small junction shunted by an inductor (the latter can be either an array of junctions or a superconducting nanowire). The qubit transition frequency ω 10 (f ) depends on the external flux; moreover for this qubit, assuming ε = 1, one can show that the relaxation rate near f = 1/2 takes the form [15] Γ 10 = F 2 4 where the dimensionless prefactor F can be calculated given the qubit parameters E C , E J , and E L . With regard to the flux dependence, this expression extends the validity of Eq. (76) to the region near half flux quantum. Being a consequence of fluctuation-dissipation relations [cf. Eq. (69)], Eq. (80) can be expected to hold also for ε = 1. Therefore, measuring the flux dependence of the relaxation rate makes it possible to estimate the value of ε. The result of the measurements together with theoretical curves for different values of ε are shown in Fig. 3; conservatively, one can estimate ε > 0.99.

Quasiparticles-driven e-jumps in a transmon
A qubit without an inductor is described by Hamiltonian (31) with E L = 0 acting in a space of 2π-periodic wave functions. If the ratio E J /E C is not too large (typically, less than about 25 for a transmon), the dependence of the energy levels on n g is well-resolved in experiments.
As we already discussed in Sec. 2.1, n g exhibits uncontrollable jumps by ±1/2 associated with the quasiparticle tunneling (e-jumps). Therefore, at a given n g a transmon is not a two-level system, but in fact is a four-level system: two states differing by the charge parity represent the logical ground state, and another pair forms the logical excited state. Quasiparticle tunneling results in the parity-changing transitions within the four levels. The relaxation rates between the logical states considered in the previous section are one example of the e-jumps, but transitions which do not change the qubit logical state while changing its parity are also possible, see Fig. 4. In our notations, see Eq. (65), these rates are Γ 00 and Γ 11 , for the transitions within the logical ground and excited states, respectively; they are also known as parity switching rates. Rates Γ 10 and Γ 01 contribute to the energy relaxation rate 1/T 1 of the qubit. The presence of finite rates Γ 00 and Γ 11 contribute to the qubit dephasing. The dephasing manifestation depends on the type of experiment. Suppose first the phase evolution of a qubit may be measured over time intervals shorter than 1/(Γ 00 + Γ 11 ), followed by averaging over many measurements. The energy levels of the qubit E 1 and E 0 would fluctuate from one measurement to another due to the e-jumps occurring between the measurements (cf. Section 2.1) and therefore the frequency ω 10 would fluctuate by some δω = (δE 1 + δE 0 )/ . The averaged over the measurements result then would yield the phase relaxation time T 2 ∼ 1/δω, in analogy with the inhomogeneous broadening [16] of magnetic resonance (assuming we identify δω with the spread of the magnetic resonance frequencies caused by static disorder in a solid). In a further analogy with the solid-state magnetic resonance, this relaxation mechanism is successfully countered by the echo technique [17], if the time 1/(Γ 00 + Γ 11 ) is long enough to allow for the echo pulse sequence. In the opposite case of high rates Γ 00 and Γ 11 , fast e-jumps would lead to the phase diffusion with the diffusion constant D ϕ ∼ (δE 2 1 /Γ 11 + δE 2 0 /Γ 00 )/ 2 . In this Section, we evaluate all four rates, Γ 10 , Γ 01 , Γ 00 , and Γ 11 for a transmon. We assume the quasiparticle distribution is described by an effective temperature T eff ∆ and the dimensionless density x qp 1, which may or may not correspond to zero chemical potential.
Aiming at the most realistic case, we take T eff ω 10 , which allows us to use Eq. (68) in Eq. (66). Furthermore, disregarding small anharmonicity we utilize Eqs. (74) and (75) with E L = 0 and ϕ 0 = 0, respectively, to find the qubit matrix element, entering Eq. (66). As a result, the latter equation leads to Generalization to arbitrary T eff / ω 10 amounts to multiplication of the right-hand side of Eq. (82) by 4 ω 10 /πT eff exp( ω 10 /2T eff )K 0 ( ω 10 /2T eff ), where K 0 (z) is a modified Bessel function, see [18]. Rate Γ 01 can be obtained from Eqs. (82) and (70). Within the harmonic approximation, the level-preserving transitions (0 → 0 and 1 → 1) occur between states which, up to exponentially small corrections [10], have the same parity of the wave function ψ(ϕ). Therefore, rates Γ 00 and Γ 11 come only from the i| cos(ϕ/2)|i ≈ 1 Figure 4: (adapted from Ref. [19]) Two lowest energy levels of the transmon qubit. The "even" and "odd" labels mark states of opposite charge parity. term in Eq. (65). A straightforward evaluation [18] yields In the derivation of these rates we accounted for the near-degeneracy between the states connected by the transitions. Indeed, the corresponding energy difference (divided by ) is at most a few MHz, so it satisfies the condition ω T eff , even is we assume that quasiparticles equilibrate at the fridge temperature (10 mK ≈ 200 MHz).
Comparing Eqs. (82) and (83), we find the ratio of the rates, Γ 11 /Γ 10 = (8E J /E C ) 1/2 · ( ω 10 T eff /π∆ 2 ) 1/2 . The first factor here is large, while the second one is small. For parameters of a typical transmon, the second factor wins the competition, so that Γ 11 /Γ 10 1; therefore it is unlikely for e-jumps to cause any phase diffusion.

Generation of quasiparticles by photons
So far we have considered transitions in a qubit due to quasiparticles, but neglected any effect of the external environment. In general, a qubit is coupled to the environment via a cavity or a waveguide resonator which support a number of modes; in other words, the qubit is coupled to photons whose frequencies ω ν depend on the geometry of the device. We can distinguish two kinds of photons: those with frequency lower than the pair-breaking energy, ω ν < 2∆/ , can be absorbed or emitted during a tunneling event of a quasiparticle already present, similarly to the absorption/emission of the qubit energy by the quasiparticles. Assuming that only a small number of such low-frequency photons are present in the cavity, this type of photon-assisted tunneling can only contribute a small correction (which we neglect) to the rates calculated in the previous section. In contrast, higher frequency photons, ω ν > 2∆/ , can always be absorbed, even in the absence of quasiparticles, as they have enough energy to break a Cooper pair and thus create two quasiparticles. The photon wave length, even at the photon energy somewhat exceeding 2∆, is still comparable to the size of the qubit. Therefore it is fair to assume that the alternating voltage generated by a photon is applied across the junction. Thus breaking of a Cooper pair generates two quasiparticles, one on each side of the junction. Here we consider the e-jump rates associated with such pair-breaking events.

Theory of photon-assisted e-jumps
The inclusion of the qubit-photon interaction in our model can be accomplished [20] by adding to Eq. (61) the Hamiltonian for the photon (we focus on one mode of a cavity here for simplicity) where b † ν and b ν are the creation and annihilation operators for the photon, and by replacing in Eqs. (46). Here φ ν is the amplitude of zero-point fluctuation of the phase due to the electric field E ν (0) at the junction position: with d ν being the effective dipole length that relates the electric field to the voltage drop U ν across the junction, U ν = d ν E ν (0). With these definitions, one can recognize that the replacement in Eq. (85) originates from the relation between phase and voltage in Eq. (36), see also the text preceding Eq. (53). For our purposes, it is sufficient to perform the replacement (85) in Eq. (45) and then consider the first term in the expansion over the small parameter φ ν 1. In this way, we obtain the following quasisparticle-qubit-photon interaction term: Using as before Fermi's golden rule, we can find the transition rate Γ ph if between the initial state with qubit in state |i , no quasiparticles, and one photon in the cavity, and the final state with qubit in |f , two quasiparticles, and no photon: can be related to the coupling strength between qubit and cavity [20], and with E and K the complete elliptic integrals of the second and first kind, respectively. These structure factors have the following properties: S ± (x) = 0 for x < 2, S ± (x) ≈ x for x 2, and Note the similarities between Eqs. (65) and (88): in both there are a prefactor that accounts for the coupling strength, and the squared matrix element of sin ϕ/2 and cos ϕ/2 multiplied by qubit-frequency-dependent structure factors. The latter additionally depend on the quasiparticle distribution function in Eq. (65) or on the photon frequency in Eq. (88); this difference relates to the different physical origin of the transitions: those with rates in Eq. (65) require quasiparticles to be present but no photons, while those in Eq. (88) require the presence of photons but not quasiparticles.
Let us consider again the case of a single junction transmon; then using Eq. (88) we find for the parity switching, relaxation, and excitation rates: The relaxation and excitation rates are generally close: for ω 10 < min{ω ν − 2∆/ , 2∆/ } we find with the lower bound saturated as ω ν → 2∆/ and the upper one for ω ν → ∞. This is in contrast to "cold" quasiparticles, in which case Γ qp 01 /Γ qp 10 1. Finally, the ratio between parity switching and relaxation rates can be large, since for ω ν 2∆/ . On the other hand, for photons near the pair-breaking threshold, ω ν −2∆ ∆, we find Γ ph and the large prefactor on the right hand side can be compensated by the small, final factor originating from S − (x).

Comparison with experiments
Various experiments have reported measurement of qubit rates for transitions that can be caused by quasiparticles. In the work introducing the 3D transmon architecture [21], the T 1 time was measured as function of temperature T . At low temperatures, this time was roughly independent of T , while it became quickly shorter at higher temperature, see Fig. 5.
A possible explanation of this behavior is that at low temperatures there are non-equilibrium, cold quasiparticles with the density x qp ≈ 3 · 10 −7 . This value is exceeded by the density x eq qp of equilibrium thermally-activated quasiparticles at T 120mK. Consistently with that, the relaxation time drops quickly upon the further increase of temperature. The quasiparticledriven relaxation rate scales linearly with the quasiparticle density, so we may separate the contributions of x qp and x eq qp to 1/T 1 , Alternatively, 1/T 0 1 could be dominated by a different, non-quasiparticle mechanism. One may attempt to distinguish between the mechanisms by measuring T 1 as a function of flux [15] or attempting to separate the relaxation processes involving e-jumps from those preserving the charge parity [22]. The fluxonium experiment indicates the presence of quasiparticle-driven relaxation at low temperatures and yields x qp ≈ (1−32)·10 −7 . In the transmon experiment [22] temperature-independent relaxation was dominated by mechanisms other than e-jumps.
A more recent experiment [19] in a similar setting (low E J /E C ratio) extracted all the six transition rates between the four qubit states. We show in Fig. 6 the parity-changing, e-jump rates. Again we see that they are independent of temperature at low temperatures, and quickly increase at higher temperatures. Moreover, Γ 10 /Γ 01 ∼ 1, giving a strong evidence for the photon-assisted e-jumps. In fact, except for Γ 11 , we can fit the data assuming that the rates are given by the sums of the contributions stemming from thermal quasiparticles and  Figure 6: (data points from Ref. [19]) Symbols: experimental data for the transition rates vs temperature. Dashed lines: theoretical rates calculated using Eq. (100) and the parameters given in the text.

Energy relaxation, recombination, and trapping
While multiple experiments indicate a low-temperature saturation of the quasiparticle density at a level x qp ∼ 10 −7 − 10 −6 , the source of such non-equilibrium population is not uniquely identified. Monitoring of the occupation of the fluxonium states [23] over a ∼ 10 min time span indicates that quasiparticles arrive in bunches, which leads to a non-Poissonian statistics of the quantum jumps between the qubit states. The qubit temperature (measured by the relative occupation of the states |0 and |1 ) remains low, favoring the assumption that the non-equilibrium quasiparticles are also cold, T eff ≈ 40 − 60 mK, see Figs. 7(a) and 7(b).
Quasiparticles may relax their energy in a superconductor by emitting phonons. Recombination is also accompanied by emission of a phonon that carries away the energy (∼ 2∆) released in the annihilation of two quasiparticles (in this discussion, we focus on the zerotemperature limit and low quasiparticle densities). Recombination requires a meeting of two quasiparticles, therefore the corresponding rate equation has the form dx qp /dt ∝ −x 2 qp . The recombination rate per quasiparticle scales as 1/τ r ∝ x qp . The relaxation rate, on the other hand, is independent of x qp but has a strong dependence on the energy E of the quasiparticle measured from the gap. The electron-phonon relaxation rate in metals is strongly affected by disorder. For thin films, the "dirty limit" in which the superconducting coherence length exceeds the electron elastic mean free path, is an adequate approximation. The theory of these rates is beyond the scope of these lectures; its summary can be found in Ref. [24]. The same work provides a detailed theory of the relaxation (1/τ E ) and recombination (1/τ r ) rates for quasiparticles in superconductors (with or without applied magnetic field). Using the results of [24] [specifically, Eqs. (57) at zero field, T eff ∆, and E ∆] we find Here 1/τ N (∆) is the rate of relaxation of an electron with energy ∆ (measured from the Fermi level) in the normal state. It depends, among other parameters on the electron mean free path. We extract the estimate 1/τ N (∆) = 4 · 10 8 s −1 for Al with electron diffusion constant D = 20 cm 2 /s from the data of Ref. [25]. In finding the rates, we considered, respectively, relaxation of a quasiparticle by a spontaneous phonon emission and recombination of a quasiparticle with a background quasiparticle distribution characterized by x qp and T eff . We also assumed that the thickness of the superconducting film exceeds the wavelength of a phonon with energy E; violation of that condition reduces [24] by 1 the exponent of the E/∆ factor in 1/τ E . Comparing the two rates of Eq. (101), we find that despite the precipitous drop of 1/τ E with energy, this rate still exceeds by an order of magnitude the recombination rate 1/τ r at E = T eff = 35 mK and x qp = 10 −6 . The recombination rate may be further reduced by the re-absorption of phonons in the superconductor, recreating pairs of quasiparticles [26]. That provides one with the grounds for assuming for quasiparticles a Gibbs distribution with a finite chemical potential and T eff = T . This assumption is further justified by the expectation that qubit manipulation with microwaves does not heat the quasiparticles [27].
The considered above processes may occur regardless the presence of macroscopic inhomogeneities in a superconductor. Inclusion of inhomogeneities, however, brings about two more elements of the quasiparticle dynamics, i.e., their diffusion and trapping. Experiments with qubits have opened new ways to investigate these processes, based on monitoring the relaxation rate of a qubit. Indeed, an excess density of quasiparticles created by a AC bias jolt applied to the junction affects the qubit relaxation rate, see Eqs. (66) and (68). The total, time-dependent relaxation rate Γ(t) of a transmon qubit can be written in the form where x qp (0, t) is the time-dependent excess quasiparticle density at the position of the junction, while Γ 0 accounts for all other relaxation mechanisms (including relaxation due to a possibly finite steady-state quasiparticle density). Using Eq. (75) with ϕ 0 = 0 for the matrix element of a transmon, we have γ 2∆ω 10 / /π for the constant in Eq. (102). In the next sections we consider first the effect of vortices, which hints at the possibility of affecting the dynamics and thus improve qubit performance. Motivated by these results, we then study the stronger effect of a normal-metal quasiparticle trap.

Single-vortex trapping power
When a thin superconducting film is cooled below its critical temperature in the presence of a perpendicular magnetic field, vortices are trapped in the film if the field is higher than a certain threshold. For a strip of width W , this threshold is of the order Φ 0 /W 2 , see Ref. [28] for a detailed discussion (further discussion of a ring and disk geometries can be found in Ref. [29]). The threshold field is usually small, amounting to a few milliGauss for a strip or disk with a width of few tens of microns. This implies that vortices can be avoided or permitted in certain regions of a superconducting circuit by properly choosing the dimensions of its features. For example, for the qubit design in Fig. 8, top panel, at sufficiently low field the vortices will only be trapped into the large square pads at the ends of the long and thin antenna wire.
In the presence of a vortex, the superconducting order parameter is suppressed over a core region of radius ∼ ξ, with ξ the coherence length. In this region, a quasiparticle can loose energy (e.g., by emitting a phonon), since there are states available with energy below the bulk gap, and thus get trapped in the vortex core. Experiment [30] revealed the effect of a single vortex on the quasiparticle dynamics and measured the relevant quantity, the trapping power P of a vortex, which is an intrinsic property of a vortex, independent of the device geometry.
At a phenomenological level, we expect the dynamics of the quasiparticle density to be governed by the following generalized diffusion equation: where D qp is the (temperature-dependent) quasiparticle diffusion constant, P is the "trapping power" of a single vortex, and the sum is over all the N vortices at positions R i . The trapping by vortices leads to an exponential decay of the density; we calculate this decay rate in a simplified model of a (quasi-)1D wire of length L and width W L attached to a square with the boundary condition ∇ ⊥ x qp = 0 at the boundaries (here ∇ ⊥ denotes the gradient in the direction perpendicular to the boundary). While Eq. (104) can be easily solved analytically, this is not possible for Eq. (105). However, so long as the diffusion rate D qp /S inside the pad is fast compared to the density decay rate, the density within the pad can be taken to be approximately uniform and we can derive from Eq. (105) a boundary condition for Eq. (104) by integrating the former over the pad area to obtain The solution to Eq. (104) can be written in the form: x qp (y, t) = e −st α cos (ky) where s = D qp k 2 is the density decay rate (i.e., the total trapping rate), and the boundary condition at the origin is satisfied. 1 Then, using the boundary condition Eq. (106), we have Figure 9: (adapted from Ref. [30]) Left: quasiparticle density decay rate (or trapping rate) vs cooling field for two devices . Right: trapping rate times area vs field for three devices; at low field, the stepwise increase of the rate is evident.
the following equation for k: with z = kL, A W = W L the wire area, and τ D = L 2 /D qp the diffusion time along the wire. There are two distinct limits for the solution of Eq. (108).
In the limit of small number of vortices of weak trapping power, N P τ D /A W 1, we have z 2 ≈ N P τ D /A, where A = S + A W is the total device area. In this regime, the density decay rate s is then which is proportional to the number of vortices. It is the vortices which provide the bottleneck for the quasiparticle evacuation from the vicinity of the Josephson junction. The measured in experiment rate s exhibited step-wise increase with the external field, see Fig. 9. Each step corresponds to entering of a vortex in a pad. The step height allowed one to extract [30] the trapping power of a single vortex, P ≈ 0.067 cm 2 /s. Comparing this finding with a theory remains to be a challenge. A crude estimate of the electron-electron interaction effect within the vortex core yields [30] trapping power smaller than the observed value by a factor ∼ 10 2 . The additional effect of the periphery of the vortex was considered in [24]. At the vortex periphery, the gap is suppressed compared to its nominal value; a propagating quasiparticle may emit a phonon and get trapped in that region. The additional rate associated with such process does not resolve discrepancy with the experiment, but indicates an interesting and yet unexplored temperature dependence of the trapping power.
As the number of vortices increases, the bottleneck shifts to the diffusion along the wires connecting the junction to the antenna pads. In terms of Eq. (108), it means the existence of an upper limit for the solution: z ≈ π/2 for N P τ D /A W 1. In this case, the density decay rate is determined by the diffusion rate: with diffusion coefficient D q ≈ 20 cm 2 /s. Along with the stepwise increase of the decay rate for small vortex number, an upper bound for the decay rate was also measured, see Fig. 9.

Normal-metal traps
The core of a vortex can be thought of as a small (size ξ 2 ) normal-state region inside a superconductor. Since quasiparticles can be trapped there, one can expect that an actual normal-metal island can also act as a quasiparticle trap. Here we consider the case of such an island in tunnel contact with a superconductor -that is, the normal and superconducting layers are separated by a thin insulating barrier, so that the contact has low transparency. In this situation, we have again a generalized diffusion equation for the quasiparticle density: In the last term, the function a(r) is unity in the normal-metal/superconductor contact region and zero elsewhere. The effective trapping rate Γ eff accounts for the competition of three effects (see Fig. 10): a quasiparticle in the superconductor can tunnel into the normal metal at rate Γ tr ; once in the normal metal, the excitation can relax to states with energy below the gap at rate Γ r , or it can escape back to the superconductor at rate Γ esc ( ). The latter is energy dependent because in calculating such a rate via Fermi's golden rule, the bare (energy-independent) tunneling-out rate Γ esc is enhanced by singularity of the (normalized) BCS density of states, Identifying the quasiparticle effective temperature with T , we can distinguish two limiting regimes (see Ref. [31] for details) for the effective trapping rate Γ eff : if relaxation is fast, Γ r Γ esc ∆/T , then the "bottleneck" process is tunneling into the normal metal and Γ eff ≈ Γ tr ; if relaxation is slow, Γ r Γ esc ∆/T , then relaxation is the bottleneck and For both the slow and fast relaxation regimes, according to Eq. (111) the dynamics of the density is determined by diffusion and effective trapping rate. Similarly to the case of vortices, we can gain a qualitative understanding of the dynamics by studying a simplified model, see Fig. 11. Let us consider a superconducting strip of length L + d and width W L, with the region −d < y < 0 in contact with normal metal and the region 0 < y < L free; then Eq. (111) simplifies to with the boundary condition ∂x qp /∂y = 0 at y = −d, L. So long as the trap is small, d λ tr , compared with the trapping length defined as λ tr = D qp /Γ eff , we can treat the trap in the Right: depiction of the processes determining the effective trapping rate: tunneling into the normal metal with rate Γ tr , relaxation to below the gap ∆ with rate Γ r , and escape back into the superconductor with rate Γ esc ( ). While the normal-metal density of states is featureless (top), the superconductor's one is peaked at the gap and zero below it (bottom).
same way as we treated the pad in the previous section, and from integrating over the trap area obtain an effective boundary condition at y = 0: Taking the solution in the region y > 0 to be of the form x qp = e −st α cos[k(y − L)] with s = D qp k 2 , and using this boundary condition, we arrive at with z = kL and l 0 = πD qp /2Γ eff L = πλ 2 tr /2L. Up to different parameters, this equation for z has the same form as Eq. (108); 2 therefore, we find again two regimes, one for small, weak traps with the density decay rate given by valid for d l 0 , and one for large/strong traps (d l 0 ) in which the decay rate is limited by diffusion, see eq. (110).
The cross-over between the weak and strong trapping regime can be studied by increasing the trap length d while keeping everything else equal. Such experiments were carried out with transmon qubits of design similar to that used in the vortex experiments, see Fig. 11. In Fig. 12 we show the (normalized) density decay rate vs (normalized) trap length: after 2 in the regime d λtr L, for the slow modes one can neglect the second term in Eq. (116), which then reduces to the expression found in Ref. [31]. an initial linear increase with length, the decay rate saturates. From the linear part, one can extract the effective trapping rate, which turns out to increase with increasing fridge temperature. This finding (as well as independent measurement of the relaxation rate) is in qualitative agreement with the expectation of slow relaxation in the normal metal being the bottleneck for trapping, see Eq. (113).
A more accurate modelling of the qubit geometry was considered in Ref. [32], where optimization in the number and position of traps was also studied, together with the other advantages of using traps (increase in the steady-state T 1 time and reduction of its fluctuations over long time scales). It should be noted that normal-metal traps can also introduce new dissipation mechanisms for the qubit. For example, the inverse proximity effect broadens and soften the superconducting gap, introducing subgap states into which the qubit can loose energy; this effect weakens exponentially with the ratio between junction-trap distance over coherence length and can be neglected [33]. Current within the normal metal and the tunneling current through the barrier between the normal metal and superconductor can also dissipate energy. While the contribution of the latter to qubit relaxation is negligible, the former one imposes some constraints on trap design which are not relevant to current qubits, but could be limiting for qubits with improved coherence [34]. Such limitations of normal-metal traps can be largely sidestepped by using instead gap-engineered traps, obtained by good contacts between two different superconductors [35].

Quasiparticle trapping in Andreev levels
At the end of Section 2.2 we mentioned that the normal-state conductance scales proportionally to the product of the junction's cross-sectional area Σ and the electron transmission coefficient |t B | 2 of the tunnel barrier, G N ∝ Σ · |t B | 2 . In discussing the Josephson effect in Section 2.3 and thereafter, we concentrated on the low-transmission, large-area tunnel junctions with G N ∼ e 2 / . One may ask, if any new effects appear in smaller-area junctions, where the same value of G N is achieved by increasing |t B | 2 . The answer is affirmative, due to the increasing prominence of the sub-gap Andreev levels associated with a junction. Phase biasing of a junction of any |t B | 2 leads to the Andreev levels appearance. We may illustrate it with a simple example of a point contact with |t B | 2 1. In terms of tunneling Hamiltonian (33), a contact of an area Σ λ 2 F is modeled by a matrix t n L n R having only one non-zero eigenvalue (i.e., only a single electron mode may go through the junction). Without loss of generality, we may take (ν 0 V)t n L n R = t B independent of n L , n R [here V and ν 0 are, respectively, the volumes of and the electron density of states in the two leads which we assume identical, cf. Eq. (40)]. At |t B | 1, we concentrate on a single-quasiparticle sector, keep only the term (44) and dispense with other terms (which do not conserve the quasiparticle number) in the tunneling Hamiltonian (43). Furthermore, expecting shallow bound states just below ∆ at small |t B |, we replace u n L u * n R − v n L v * n R in Eq. (46) by i sin(ϕ/2) and therefore simplify Eq. (44) to: Now we perform a canonical rotation into a new quasiparticle basis defined by the operators γ nσ± = (1/ √ 2)(γ n R σ ± iγ n L σ ). In new variables, the Hamiltonian for low energy ( n ≈ Figure 13: (from Ref. [37]) Left panel: a sketch of the Andreev levels dispersion with ϕ in a single-mode highly transparent junction. Its length is assumed to somewhat exceed the superconducting coherence length, allowing for two Andreev levels. The level degeneracy at ϕ = 0 is the manifestation of time-reversal symmetry; the degeneracy at ϕ = ±π is due to the symmetry with respect to the product of time reversal and spatial inversion transformations.
At other values of ϕ, Kramers doublets are split due to the combination of a finite Josephson current and present spin-orbit coupling. Various dashed arrows indicate transitions involving promotion of a quasiparticle from a lower to higher-energy Andreev state. Right panel: spectroscopic lines corresponding to the indicated transitions. The lines intersection at Φ/Φ 0 = 0 reflects the ϕ = 0 Kramers degeneracy. ϕ = 0. Regardless the junction's length, Zeeman effect associated with an applied magnetic field breaks time-reversal symmetry and lifts the spin degeneracy even at ϕ = 0. A phase bias ϕ = 0 across the junction leads to a Josephson current and is another source of time-reversal symmetry breaking. In the presence of spin-orbit coupling, a finite Josephson current may cause the spin splitting of an Andreev level. Such spin-split structure of Andreev levels in a finite-length, high-transmission junction is sketched in the left panel of Fig. 13. The Andreev levels lie below the edge of the quasiparticle continuum, and therefore are prone to trap quasiparticles. In terms of occupation factors n i A of the Andreev levels, the corresponding correction to the energy of the junction is δE(ϕ) = i n i A E i A (ϕ). At each given time, the set of factors {n i A } is drawn from integers 0 and 1 (we assign different superscripts to the components of a degenerate level). The set {n i A } changes from time to time, due to the inelastic relaxation of quasiparticles interacting with phonons. Based on Eqs. (101) and the discussion in the beginning of Section 5, we expect their rate to be slower than 1/τ N (∆). The rate is reduced further by low total average number of trapped quasiparticles, in which case we also may expect the energy relaxation to occur faster than the recombination.
At fixed {n i A }, the trapped quasiparticles affect the inductance of the junction. If the latter is a part of an LC-circuit, trapping shifts down its resonance frequency. This shift provides one with a measurable "fingerprint" of {n i A }. This kind of experiment was performed [38] with an Al nanobridge which may be considered as a short junction of a cross-section passing about 700 electron modes. Under applied phase bias, each mode with a particular value of |t i B | 2 gives rise to an Andreev level with energy E i A given by Eq. (122). At a given bias ϕ, resonance traces of the reflection amplitude were accumulated over time exceeding the evolution time of the {n i A } set. The averaged trace therefore corresponds to an average over Figure 14: (from Ref. [38]) Left panel: resonance traces of an LC circuit containing an Al nanobridge. The estimate number of quantum modes propagating through the bridge is ∼ 700. Traces are averaged over many measurements spanning time interval exceeding the time needed for the rearrangement of the set of occupation factors n i A . At low temperatures, two shoulders are clearly seen on the low-frequency side of the main resonance. The shoulders correspond to the shifts of the resonance frequency caused by quasiparticle poisoning, respectively, of one or two Andreev states. Right panel: the average number of trapped quasiparticles (left axis, black circles) and x qp (right axis, red squares) extracted from experiment. all configurations {n i A }, weighted by their probabilities. As flux bias grows, energies E i A (ϕ) drop, and so the likelihood of a quasiparticle trapping in an Andreev state grows. The left panel of Fig. 14 shows the averaged traces at ϕ = 0.464π for a set of temperatures. At the lowest temperature, multiple shoulders on the low-frequency side of the main resonance are resolvable, indicating multiple quasiparticle trapping numbers. The most prominent shoulder corresponds to a single quasiparticle trapped by an Andreev level associated with one of the modes. Its width comes from the range of the E i A (ϕ) values of ∼ 700 Andreev levels. It is quite remarkable that quasiparticle "poisoning" of one out of 700 quantum modes is traceable in the experiment. A weaker feature in the trace appearing further away from the main resonance peak corresponds to trapping of two quasiparticles. At higher temperatures, first the 2-quasiparticle and then the 1-quasiparticle shoulder shrink, leading to a Lorentzian resonance at T ∼ 150 mK. The extracted from experiment temperature dependence of the average number of trapped quasiparticles n trap is shown in Fig. 14. At T 170 mK, the number n trap grows as temperature is reduced; this is characteristic for a non-equilibrium population.
Microwave technique was recently also applied to studying Andreev levels in atomic point contacts [39,40] and in proximitized semiconductor wires [37,[41][42][43]. Such junctions carry only one or a few electron modes which allows one to perform spectroscopy of individual levels. There is a simple rule of thumb for assessing the odds of quasiparticle poisoning of an Andreev level at low temperature T and small, temperature-independent x qp . Assuming Boltzmann distribution of the quasiparticles in energy, we find their chemical potential (measured from the edge of the quasiparticle continuum), µ = (T /2) ln(x 2 qp ∆/2πT ). The justification for the use of Boltzmann distribution is the inequality between the relaxation and recombination rates, 1/τ E 1/τ r . This, in turn, requires small occupation factors of the quasiparticle states. Once the energy E A (ϕ) − ∆ of an Andreev level drops below µ, one may expect its high occupation. This condition was clearly satisfied in experiments [37], where all the detected transitions, see Fig. 13, had n i A = 1 in the initial state.

Conclusions
Constructing a quantum information device calls for finding elementary building blocks capable of maintaining quantum coherence over extended time periods. Superconductors provide one with a head-start in the race for a perfect device: the superconducting ground state locks together a macroscopic number of degrees of freedom, leaving a small number of collective variables to build a qubit from. The coherence distributed over many particles is inherently more robust than that of a single spin or atom. This robustness allows one to shorten the preparation and readout times for a superconducting qubit. However, some hazards come along with the macroscopic dimensions of a qubit. Many of them are defeated by now ubiquitous circuit QED architecture [44][45][46]. That sharpened the attention to the unwanted influences of superconducting quasiparticles in the "conventional" circuit QED devices [47] and in putative topological qubits [48]. It is clear by now that the observed low-temperature quasiparticle density by far exceeds its equilibrium values in a broad variety of devices. Their sources are still not fully identified, with photons [20], phonons [49,50], and even cosmic rays [47,51,52] being contenders. Meanwhile, superconducting qubits have provided one with an unrivaled technique for time-domain experiments. In many cases, it is by far the most sensitive tool for investigation of elementary processes in quasiparticle dynamics. It is this tool that allowed one to resolve such subtle effects as the tiny dissipative cos ϕ-component of the Josephson current and quasiparticle trapping rate by a core of a single vortex line. Improving the qubit performance goes hand-in-hand with the ever-increasing capability of the techniques they provide for physics research.