Circuit-QED with phase-biased Josephson weak links

By coupling a superconducting weak link to a microwave resonator, recent experiments probed the spectrum and achieved the quantum manipulation of Andreev states in various systems. However, the quantitative understanding of the response of the resonator to changes in the occupancy of the Andreev levels, which are of fermionic nature, is missing. Here, using Bogoliubov-de Gennes formalism to describe the weak link and a general formulation of the coupling to the resonator, we calculate the shift of the resonator frequency as a function of the levels occupancy and describe how transitions are induced by phase or electric field microwave drives. We apply this formalism to analyze recent experimental results obtained using circuit-QED techniques on superconducting atomic contacts and semiconducting nanowire Josephson junctions.


I. INTRODUCTION
Circuit-QED (cQED) was born from the transcription of the concepts of Cavity Quantum Electrodynamics [1,2], which describes the strong coupling of atoms to the photons of a cavity, to superconducting circuits in which the discrete spectrum of collective electromagnetic modes mimics the energy levels of an atom [3,4]. These techniques are nowadays at the core of quantum information processing with superconducting circuits [5,6]. Furthermore, the coupling of superconducting cavities with a variety of systems is now intensively used to probe and manipulate other degrees of freedom at the single quantum level [7][8][9]. These can be either collective modes, like the modes of a nanomechanical oscillator [10], or microscopic ones like the spin of an electron in semiconducting quantum dots or in defects in an insulator [11]. Most superconducting circuits probed so far using cQED correspond to the former case, as the relevant modes are surface plasmons in nonlinear oscillators, rendered dissipationless by the existence of the superconducting gap. A very interesting exception arises in phased-biased Josephson weak links or in quantum dots with superconducting electrodes, where localized discrete subgap states (Andreev bound states (ABS) [12][13][14][15] or Yu-Shiba-Rusinov states [16,17]) develop, giving rise to a truly atomic-like spectrum of levels for fermionic quasiparticles. Two types of transitions, both conserving fermion parity, can be driven in an Andreev system: pair transitions, in which two quasiparticles are excited at once, and single quasiparticle transitions that correspond to atomic-like transitions of a quasiparticle from one Andreev level to another. Pair transitions, first observed in atomic contacts [18], give rise to the Andreev qubit [19][20][21]. Single quasiparticle transitions, first observed in InAs weak links [22], offer a route, alternative to quantum dots, to couple a single fermionic spin with a microwave resonator and develop an Andreev spin qubit [23][24][25][26][27]. Most of these experiments were performed using cQED techniques [19,22,23,27,28]. The aim of the present work is to develop a theory of the coupling of a microwave resonator to a multi-level many-body fermionic system of Andreev levels, able to describe the main features of the spectra measured in those recent experiments, in particular the intensity of the transition lines, and analyze the possible existence of selection rules associated to the spin.
The article is organised as follows. In Sec. II we introduce the theoretical model and discuss the resonator frequency shift in the absence of driving fields. Starting from the microscopic Bogoliubov-de Gennes (BdG) equations for the weak link, we express the shift in the resonator frequency as a function of the occupancy of individual Andreev levels. In Sec. III we introduce a model to describe the driving through either an ac flux or an ac gate voltage and find under which conditions spin-nonconserving transitions can occur. In Sec. IV we derive the resonator frequency shifts in the presence of driving fields. In Sec. V, we compare the predictions of the theory to experimental data on superconducting atomic contacts and nanowire weak links.

II. RESONATOR-WEAK LINK COUPLING
The system we consider comprises a weak link of length L embedded in a superconducting loop that is inductively coupled to a microwave resonator (Fig. 1). The microwave cavity is represented as a lumped-element LC resonator with bare resonance frequency f r = (2π √ L r C r ) −1 . Introducing the photon annihilation (creation) operators a (a † ), it can be described by the Hamil-tonianĤ r = hf r a † a.
The dc phase difference δ = 2πΦ/Φ 0 across the weak link, with Φ 0 = h/2e the flux quantum, is imposed by the magnetic flux Φ threading the superconducting loop (we assume in this section that the loop inductance is negligible). An electrostatic gate fixes with its dc value V g the electrochemical potential in the weak link. We consider applying a microwave drive either on the gate voltage V ac g or on the flux Φ ac . The Hamiltonian for the weak link can be written in the formĤ T is the Nambu bispinor field operator and x is the position along the weak link. We denote by |Φ iσ the eigenstates of the Bogoliubov-de Gennes (BdG) equation H 0 (δ)|Φ iσ = E iσ |Φ iσ , which correspond to Andreev states when |E iσ | < ∆, where ∆ is the superconducting gap in the leads. In this notation the subscript iσ refers to the level i with spin σ. Levels labeled with positive i are above the Fermi level. Due to the electron-hole symmetry implicit in the BdG formalism, each state iσ is associated to a state with opposite spin at opposite energy −iσ. Notice that when spin-orbit interaction is at play, as can be the case in a nanowire weak link, spin is no longer a good quantum number and σ has to be understood as a pseudospin index noted "u" or "d" in the following. Otherwise σ corresponds to spin up (↑) or spin down (↓). In Fig. 2 we show the energy of the Andreev states arising in a weak link with a single occupied channel, for three cases of increasing complexity. (a) For a zero-length junction, there is only one pair of subgap spin-degenerate levels i = ±1. Each level can be occupied by 0, 1 or 2 quasiparticles. This case describes well atomic contacts weak links with length L ξ, where ξ is the superconducting coherence length. For ballistic conduction channels, ξ = v F /∆, where v F is the Fermi velocity at the weak link. (b) For finite length weak links, without spin-orbit coupling, the parameter λ = L/ξ determines how many spin-degenerate Andreev pairs appear in the gap: depending on δ and channel transmission τ , this number is 1 + 2λ/π or 2 + 2λ/π ( x is the integer part of x). For the parameters of Fig. 2b (λ = 1.7, τ = 0.97), i = ±1, ±2. Compared to short junctions, a new type of excitation arises: a quasiparticle in state 1 can absorb a photon and be excited to state 2. (c) The spin character of these excitations becomes relevant when spin-orbit interaction is present in the weak link, which leads to a lifting of the levels spin-degeneracy when δ = 0, π. This lifting results from a spin-dependent Fermi velocity, leading to different values of λ for the two spin textures [22,26]. In Fig. 2(c), the gray level of the lines encodes the state pseudospin. Among all possible transitions, some conserve the pseudospin, others do not. This regime describes well InAs nanowire weak links [22,23].
The coupling between resonator and weak link occurs through current fluctuations in the resonator, assumed to be in its ground state, which induce phase fluctuations across the weak link, so that δ → δ +δ r , whereδ r = δ zp (a + a † ) with δ zp the amplitude of zero-point phase fluctuations. In accordance with experiments, we assume δ zp 1. In Ref. [29] we derived a general expression to compute the frequency shift of the resonator when a single level |Φ iσ is occupied (we write here explicitly the spin indices) by expanding the HamiltonianĤ 0 up to second order inδ r : , leftmost arrow in (b) and in (c)). Single particle transitions are possible when quasiparticle excitations are already present in the system. They correspond to atomic-like transitions between two levels both at either positive or negative energies (rightmost arrow in (b), rightmost arrows in (c)).
where we have introduced the transition energies E iσ,jσ = E jσ − E iσ , the curvature E iσ = ∂ 2 E iσ /∂δ 2 and M iσ,jσ = | Φ i |H 0 |Φ j | the modulus of the matrix element of the current operator H 0 = ∂H 0 /∂δ between states iσ and jσ . The coupling strength g iσ,jσ is related to M iσ,jσ by g iσ,jσ = δ zp M iσ,jσ . The term associated to virtual transitions from iσ to jσ is noted V iσ,jσ . We note that the terms in Eq. (2) arise from the H 0 = ∂ 2 H 0 /∂δ 2 term in the expansion of Hamiltonian (1), whereas the remaining terms that result from the H 0 term can also be derived from the Jaynes-Cummings Hamiltonian [29]. According to Ref. [29] two different regimes can be identified: the adiabatic one, which occurs when hf r |E iσ,jσ | for all iσ, jσ , leading to hδf (iσ) r ∝ E iσ ; and the dispersive one when hf r ∼ |E iσ,jσ | for a set of iσ, jσ , in which case the terms which involve exchange of virtual photons dominate. We shall analyze the occurrence of the two regimes for the examples presented below.
It can be seen from Eq. (2) that all transitions which couple a given state iσ with other states jσ via H 0 are relevant to calculate the shift of the resonator. Note that jσ cannot be taken equal to −iσ since electron-hole symmetry anticommutes with H 0 .
Fermionic statistics: The actual resonator frequency shift is determined not by a single but by all Andreev levels which are populated in a given many-body state of the weak link. We consider first the ground state |g , in which all negative energy levels are occupied. The frequency shift is then where the factor 1/2 compensates for the redundancy of the BdG description. Note how we differentiate in the notation the shift δf (iσ) r (with parentheses) associated to the occupancy of a single level iσ and the shift δf |Ψ r (with a ket) associated to a many-body state |Ψ . When combining Eq. (2) and Eq. (3), and taking into account that V iσ,jσ = −V jσ ,iσ , only virtual transitions to positive energy levels contribute where E |g = (1/2) i<0,σ E iσ is the energy of the ground state. Further simplification occurs in the absence of a magnetic field and in the presence of a mirror symmetry, where the operator H 0 does not mix opposite pseudospins (σ andσ) [26], so that V iσ,jσ = 0. Any many-body state |Ψ can be built by creating electron-like γ † iσ |g (i > 0) or hole-like γ iσ |g (i < 0) excitations from the ground state. Here γ † iσ is the Bogoliubov quasiparticle creation operator. Notice also that γ † −iσ = −sγ iσ due to double counting in the semiconducting picture that we are using, where s = 1(−1) for σ = u(d). The frequency shift in |Ψ is where n iσ = 0, 1 is the occupancy of the state iσ. More generally, n iσ has to be understood as the average occupancy of this state. The number of fermionic quasiparticle excitations in the weak link given by N |Ψ = i>0,σ [n iσ + 1 − n −iσ ] can be even or odd, but states with different parity are not coupled by photons.
We now illustrate these ideas in a simple case. Zero-length junctions: In the limit L → 0, accessed experimentally in superconducting atomic contacts, the BdG equation can be solved analytically. For a single conduction channel of transmission τ , the resulting pair of Andreev states within the gap, represented in Fig. 2(a), [13][14][15]. In Fig. 3 the frequency shift of a resonator at f r = 0.2∆/h as given by Eq. (5) is shown for τ = 0.8 and τ = 0.999, and for three many-body states: the ground state |g , the odd parity state |o obtained by creation of one quasiparticle |o = |1σ = γ † 1σ |g , and the lowestin-energy excited state with even parity |e = |1↑1↓ = γ † 1↑ γ † 1↓ |g . In each state, the resonator frequency shift (dashed red line in Fig. 3) results from the sum of four contributions. The first one corresponds to the contribution of the curvature E |Ψ of the many-body state |Ψ (green lines in Fig. 3). The second one (blue lines) is associated with virtual transitions between the Andreev levels -1 and +1 (see arrow in Fig. 2(a)), coupled by the matrix element of H 0 [19,30] The corresponding term reads with f A = 2E A /h the Andreev transition frequency. The third type of contribution (orange lines), is associated with virtual transitions between an Andreev level and states in the continuum C − at energies E < −∆ or C + at energies E > ∆. Using the expressions for the matrix elements of H 0 given in Refs. [30,31], one finds that the associated shift grows positive from δ = 0, presents a maximum, and exhibits a negative dip close to the condition where ∆ − E A matches hf r . This is characteristic of a threshold behavior associated with the continuum, also discussed in Ref. [32].
The last contribution, which results from virtual transitions from states in C − to states in C + , is negligible.
In the ground state ( Fig. 3(a) and (b)), the level -1 and all levels of C − are doubly occupied, so that the factor 1 2 in Eq. (4) cancels out with a factor 2 for the spin, and hδf |g r where we dropped the spin indices since H 0 conserves the spin for zero-length junctions. In the zero-length limit, the energy of the states in the continuum does not depend on phase [33], and E |g = −E A . The second term reads i<0 j>0 In practice, because of the large energy E −1,j > E A + ∆ for a transition to C + , the last term in Eq. 10 can always be neglected (11) When f A f r , the three terms from V −1,1 compensate and the frequency shift is entirely due to E A , as shown in Fig. 3(a) and far from δ = π in Fig. 3 hf A (green and blue lines in Fig. 3(b,b')), i.e. the contribution due to H 0 vanishes [29], and the frequency shift is essentially the one that can be derived from the Jaynes-Cummings Hamiltonian [34,35] δf |g ,JC with g(δ) = Mδ zp / . At the scale of Fig. 3(b'), δf |g ,JC r and the exact δf |g r coincide within the linewidth. The rotating-wave approximation (RWA), which consists in neglecting the second term: overestimates the little bump of δf |g r at δ = π by a factor ∼ 2.
The odd parity state |o = γ † 1σ |g has energy E g +E A = 0. The shift of the resonator in this case is Using and Although this state does not disperse with δ (E |o = 0), one obtains a finite shift associated to transitions from Andreev states to continuum states, which becomes sizable when E A approaches ∆ (see Fig. 3(c)). The excited state |e = γ † 1↑ γ † 1↓ |g has energy E g + 2E 1σ = E A . The shift in this state (shown in Fig. 3(d,e)) is and one gets When E A ∆ the continuum contributions can be neglected and δf |e r ≈ −δf |g r . This is no longer the case when E A approaches ∆, a situation in which both 2 j∈C + V 1,j and E A contribute to the shift, as shown in Fig. 3(d,e).
In spectroscopy experiments [19], transitions |g → |e are observed. The frequency shift that governs the measured signal is ∆f r = δp |e δf |e r − δf |g r , with δp |e the population change in the excited state due to the microwave excitation. In Fig. 3(f,g), δf |e r −δf |g r is encoded in the color of the line showing 2E A (δ). For τ = 0.8, it is dominated by the curvature term 2E A , except when E A approaches ∆ and virtual transitions enter in δf |e r . For τ = 0.999, close to δ = π, the terms associated to virtual transitions between ABS causes a change of sign of ∆f r when 2E A crosses hf r , as expressed by the dispersive approximation (Eq. (12)).
Finite-length junctions, with spin-resolved levels: The situation becomes richer when there are several Andreev levels within the gap as in Fig. 2(b,c). Furthermore, in the presence of spin-orbit ( Fig. 2(c)) the subgap states are spin split, which gives rise to a larger number of possible transitions. The terms V −iσ,jσ depend on the matrix elements of H 0 which do not have analytic forms , and (f,g) transition energy 2EA with color-coded frequency change when driving the system from |g to |e , all as a function of the phase δ. The bare resonator frequency, shown as a black line in (f,g), was taken at fr = 0.2∆/h. In the left panels, τ = 0.8, and the transition frequency is always larger than hfr; in the right ones, τ = 0.999, and 2EA crosses hfr. Total frequency shifts δf |g ,|o ,|e r in state |g (a,b,b'), |o (c) and |e (d,e) are shown with dashed red lines. They are decomposed into three contributions: states' curvature E" |g ,|o ,|e (green lines), virtual transitions among Andreev levels (blue lines) and virtual transitions from Andreev levels to continuum levels (orange lines). In each panel, the contributions that have a negligible contribution are not shown, for clarity. in this case. They can be obtained by solving numerically the BdG equation, for which we use two complementary approaches: the scattering model of Ref. [22] and a discretized tight-binding model of the nanowire (see Appendices B and C for technical details). As these methods rely on different approximations one cannot expect a one-to-one correspondence of their results. For instance, the scattering method is based on a linearization of the electrons and holes dispersion relations around the Fermi level (Andreev approximation) which is not as-   sumed in the tight-binding model. On the other hand the tight-binding model only includes two sites to describe the nanowire cross-section. We have checked, however, that the methods yield qualitatively similar results for the limits where their approximations are both valid.
In Fig. 4, we show the predictions for the frequency shifts in the case of a weak link with three spin-split manifolds of Andreev levels (same parameters as Fig. 2(c)), at zero Zeeman field. Two values of the bare resonator frequency are considered: f r = 0.07∆/h (panels (a-e)) and f r = 0.18∆/h (panels (f-j)). The frequency shift in the ground state |g is first evaluated using Eq. (4). All matrix elements are computed with a tight-binding model, as described in Appendix C. We assume that scattering takes place only in the longitudinal direction, and hence does not mix the subbands. Thus, in absence of magnetic field the matrix elements of H 0 are then zero for all pseudospin-non-conserving transitions [26]. Frequency shifts in the other states are found from Eq. (5). Transitions from |g create pairs of excitations (pair transitions), leading for example (leftmost arrow in Fig. 2 Because of the redundancy between negative-and positiveenergy states, we use here only labels corresponding to positive energies ( Fig. 4(a,f)). The states accessible from |g and involving only the two lowest subgap levels are therefore those shown in Fig. 4(b,g). We also consider the closest-in-energy states that can be reached from the single-particle state |1u ( Fig. 4(c,h)). On the one hand, states with a single quasiparticle are accessible through single quasiparticle transitions, like |2u = γ † 2u γ 1u |1u , or |2d . On the other hand, the same fermion parity is also maintained with pair transitions that lead to states with three quasiparticles: |1u1d2u = −γ † 1d γ −2d |1u , |1u1d2d , and |1u2u2d . For all these possible states, frequency shifts are given by Eq. (5), which simplifies to Figure 4(b,c,g,h) shows the transition energies from |g and |1u , with line colors encoding the resonator frequency shift for each transition (color scale in (d) or (e)). The phase asymmetry of the transition energies shown in panels (c) and (h) comes from the fact that we consider an initial state with a given pseudospin (|u ). The mirrored spectra about δ = π would be obtained when considering transitions from |1d . The situation is the simplest when the resonator photons energy hf r is smaller than the energy of all the virtual transitions entering in the calculation of δf The results look more complicated in Fig. 4(g,h), with many sign inversions of the frequency shift when sweeping δ. Sign inversions occur when the energy of one of the virtual transitions entering in the calculation of the frequency shift in the initial or the final state coincides with hf r . These coincidences are marked in panel (f) with small vertical bars linking levels with same spin distant by hf r . For example, there is one of them at and δf (2u) r present abrupt changes of sign at this phase, which is seen in all the lines involving 1u or 2u in Fig. 4(g,h). Similarly, there is another such coincidence at δ/π ≈ 1.27, where E 3d − E 2d = hf r , leading to color changes in the transition lines |g → |...2d . Detailed plots of the frequency shift for pair and single particle transitions are shown in Fig. 4(i,j), with divergences when energy differences match hf r .

III. GATE AND FLUX DRIVING
As illustrated in Fig. 1 transitions between Andreev states can be driven either with an ac electric field, using a gate [22], or with a magnetic flux, by means of an ac current in a conductor placed nearby the loop [28]. The magnetic flux can also be modulated with an excitation applied through the resonator coupled to the superconducting loop [19]. The driving can be modeled by the following term added to the system Hamiltonian where ω d = 2πf d is the driving frequency. In the case of a flux driving, which acts on the phase δ, A iσ,jσ ∝ Φ iσ |H 0 |Φ jσ . In the absence of magnetic field and for a ballistic model which preserves the transverse spatial symmetry [26], the H 0 operator does not mix the transverse channels of the weak link and thus only pseudospinconserving transitions are allowed. Notice, however, that whenever the drive or the scattering breaks the transverse spatial symmetry spin-flip transitions can take place [27].
In the case of the gate driving the ac signal induces a displacement δV ( r) in the electrostatic potential experienced by the electrons in the junction region. The corresponding matrix elements in the driving Hamiltonian are thus A iσ,jσ = Φ iσ |δV ( r)τ z |Φ jσ , where τ z is a Pauli matrix in electron-hole space.
In connection to recent experiments [22,23,28], we analyse the case of spin-split Andreev states in semiconducting nanowires. As described in Ref. [26], the pseudospin of the Andreev states comes from nanowire's transverse modes with different spins hybridized by Rashba spin orbit coupling. A perturbation δV uniform in the transverse direction does not couple different transverse modes and therefore pseudospin flip transitions are not allowed, i.e. A iu,jd = 0. Only a non-uniform perturbation couples transverse modes and allows pseudospin flip transitions.
The fact that all possible transitions between two Andreev manifolds have been observed in the experiments of Ref. [22] indicates that the non-uniform component of the induced potential δV by the gate electrode was significant. More insight into the possibility to engineer the selection rules using gate driving can be obtained by considering the model of Ref. [26] for the nanowire's transverse channels. Within this model the nanowire confining potential is assumed to have cylindrical symmetry. Thus, the modes in the lowest subband have zero angular momentum along the nanowire axis (l = 0) and they have l = 1 on the first excited subband. A lateral gate would impose a perturbation δV ( r) which typically breaks the rotational symmetry and therefore would couple states on different subbands, naturally leading to both pseudospin flip transitions and pseudospin conserving transitions. One could think, however, of a more general gate configuration like the one in Fig. 5(a), where two lateral gates can be set such that δV (y) = −δV (−y) (or δV (y) = δV (−y)) in an anti-symmetric (or symmetric) configuration as indicated in the right panel. In the anti-symmetric case, the δV matrix elements vanish for states on the same subband, but are finite for states in different subbands. As a consequence we would have A iu,ju = A id,jd = 0. The allowed transitions between spin split Andreev states are indicated in Fig. 5(c) with arrows of different colors for symmetric (blue), antisymmetric (green) or an intermediate (magenta) configuration. The phase-dependent matrix elements for each case are shown in Fig. 5(d) for the parameters that give the spectrum in (c) which correspond to a fit of the data in Fig. 11 discussed below.

IV. RESONATOR FREQUENCY SHIFTS FOR THE DRIVEN WEAK LINK
While a driven two-level system can be described using Bloch equations [36], we develop here the theory for the general multi-level case in the weak drive regime. In order to obtain the frequency shift for the resonator coupled to the driven weak link we analyze the resonator spectral We take first the interaction picture in which the time evolution of the resonator and the weak link are provided by solving master equations including dissipation, and then we treat the resonator-weak link coupling and the drive as perturbations. We calculate the perturbation terms up to second order in both δ zp and A iσ,jσ (see Appendix A for details).
The frequency shift for a single quasiparticle transition where , ω d is the driving field frequency and Γ a(b) are phenomenological parameters to account for the finite linewidths in the transition spectrum, which in our present theory are associated to the states relaxation. For a pair transition from the ground state |g , one obtains where {jσ, kσ } means a set of indices jσ and kσ corresponding to positive energy levels ordered in energy, and does not contain a permutation of the indices. jσ and −jσ are for a pair of particle-hole symmetric Andreev levels.

A. Atomic contacts
We first focus on the simplest experiments, reported partly in Ref. [19], dealing with atomic contacts hosting a small number of transport channels. The superconducting loop containing the atomic contact is coupled to a microwave resonator at f r = 10.1 GHz measured in reflection. In the experiment, atomic contacts with various channel transmissions are formed and probed with the same sample.
We first discuss single-tone continuous-wave (CW) spectroscopy data taken on seven different atomic contacts, as shown in Fig. 6. To acquire these data the microwave response of the resonator is probed as a function of flux ϕ = 2πΦ/Φ 0 , with no drive applied on the weak link. Over a small flux range around ϕ = π, the amplitude of the reflection coefficient |S 11 | displays up to three distinct local minima (in dark) as a function of frequency, as shown in Fig. 6(a-g). The positions f 1,3 of these minima were extracted by fitting |S 11 |(f ) with the linear combination 11 |(f, f i ) corresponding to a single resonance centered at f i . The extracted f i are shown with symbols in Fig. 6(a'-g'). These data can essentially be understood by considering the contribution of just one dominant channel of transmission τ 1 such that the corresponding Andreev frequency f A1 comes very close to the resonator frequency f r (a,b) or crosses it (c-g). The data are ordered with increasing τ 1 from (a) to (g). The three resonances are attributed to partial occupancy of ground, odd and excited states for the corresponding channel. As shown in Eqs. (10,18), and since, according to the analysis of Section II, contributions of the continuum can safely be neglected for phases close to π, the frequency shifts associated with ground and excited state are opposite ( Fig. 6(a,b)), and the frequency shift associated with the odd state is zero ( Fig. 6(a,e)). Finally, in order to explain a small residual global shift, it is necessary to consider the contribution of one or two additional channels with smaller transmissions, as explained below.
The resonator frequency shift is obtained from Eqs. (10,16,18), adding the contributions of all channels. The bare resonator frequency f r was determined from measurements of the resonator with open contacts. The phase δ across the contact differs slightly from ϕ = 2πΦ/Φ 0 due to the phase drop across the loop inductance (see Appendix D). In a first step, the value of δ zp = 0.0042, common to all contacts, was determined by fitting Fig. 6(d), taken on a contact for which a fit of the two-tone spectroscopy yielded τ 1 = 0.992 ( [19] and Fig. 7). For the other contacts, the fitting parameters are the transmissions τ i of two or three channels (fit parameters given in Appendix E). The transmission τ 1 of the most-transmitted channel determines the overall shape of the spectra. In the small phase interval considered here, the effect of the other channels is simply a constant overall shift of the order of a MHz. The predictions for δf r assuming the most-transmitted channel to be in |g , |o or |e are shown with blue, green and red lines.
The δ− and τ − dependence of the coupling constant g is an essential ingredient to obtain a consistent fit of all the data at once. The difference between the full theory and the JC contribution is very small at the scale of Fig. 6. In Fig. 6(a-c), since |f A1 (π)−f r | f A1 (π)+f r , the counter-rotating term ∝ 1/(f A1 + f r ) can also be neglected, and the RWA is sufficient. When τ 1 is closer to 1 (Fig. 6(d-g)), f A1 becomes significantly smaller than f r near δ = π and the counter-rotating term must be taken into account, as illustrated in Fig. 6(g") where the RWA prediction shown with a dashed line departs clearly from the data.
We present in Fig. 7 the two-tone spectroscopy measured for the atomic contact of Fig. 6(d) [19]. A single transition appears in the measurement window, corresponding to a channel with transmission τ ≈ 0.992. To acquire this data a strong microwave pulse drives the twolevel system, during a time exceeding its relaxation and dephasing times, immediately followed by a microwave tone probing the resonator. The two quadratures of the signal reflected by the resonator are measured by homodyne detection. In Fig. 7, the color scale represents the 7. Two-tone spectrum measured on the same atomic contact as for Fig. 6(d). Left: experimental data. Color scale corresponds to the change of amplitude of a quadrature of the reflected signal. Right: calculated spectrum, line color codes change in S11, scaled in order to fit at best the data.
change of amplitude δA of one of these quadratures, relatively to its value in the absence of excitation. It depends on the steady state occupancy of the states after the excitation pulse and on the frequency shift of the resonator in each state. For this data taken on a very high-transmission contact and around δ = π, the resonator frequency shifts are dominated by the dispersive shifts, with negligible effects of the states in the continuum, and the rotation wave approximation (Eq. (13)) applies: δf (6), when E A ∆). When the resonator frequency is f r + δf r , the complex reflection coefficient of a measurement tone at frequency f r + δf m on the resonator is [19] S 11 (δf r ) = 1 − 1 + e iθ /(1 + Q ext /Q int ) with In this expression, Q int , Q ext and Q t are the internal, external and total quality factors of the resonator. When driving the system, if it is in an even state (probability 1 − p o ), the occupancies of the ground and excited states change by δp g and δp e = −δp g , resulting in a change of S 11 by δS 11 = δp g (S 11 (δf |g r ) − S 11 (δf |e r )). Putting everything together, one gets, for δf m = 0, with u = Q t δf |g r /f r . To compute δp g , one uses the result of the Bloch equations [36] adapted for the presence of the odd state: with T 1 and T 2 the life time and coherence time, δω = 2π(f 1 − f A ) the detuning between the drive frequency f 1 and the qubit frequency f A , and ω R the Rabi frequency. For simplicity, we assume here that T 1 and T 2 are constant: T 1 = 4 µs and T 2 = 38 ns (values measured at δ = π). Since the excitation acts on the phase across the contact, the Rabi frequency depends on δ: ω R ∝ M 2 (see Eq. 6). The drive tone being sent through the resonator, its amplitude is filtered: In the fitting of the data, ω R is set to 2π × 4.2 MHz at δ = π in order to reproduce the measured line width.
Using Q t = 2200, Q int = 4800, f r = 10.13 GHz, and g(π)/2π = 72 MHz, one obtains the fit shown on the right hand side of Fig. 7, with the color scale of δQ adapted to match the data. Not only does the change in the Q quadrature reproduce the changes δA in the measured quadrature on the resonance line f A (δ), but one also predicts a signal at f 1 = f r , which has its origin in the very large Rabi frequency when the drive signal is not filtered by the resonator. At this precise frequency, the strong detuning δω is compensated by the large ω R in Eq. (26), and δp g is non-zero even if the drive is far from the resonant frequency f A . This feature is clearly visible in the data, although not as strong as in the calculation for f A > f r perhaps due to an effect of the resonator nonlinearity not included in the model. In Fig. 8, we show the two-tone spectroscopy and the single-tone CW data of one particular contact obtained with the same sample and showing a double avoided crossing. In the two-tone spectroscopy ( Fig. 8(a)), one observes three Andreev transitions (labeled f A1 , f A2 and f A3 ) corresponding to channels with transmissions τ 1 =0.998, τ 2 =0.992 and τ 3 =0.980, with minimum transition frequencies (at δ = π) of 4.1, 7.7, and 12.6 GHz.
Two transition lines f A1 and f A2 cross the resonator at  Fig. 8(a)). Depending on the state of the two first channels, the change of Q takes very different values, the largest change corresponding to |oog → |ooe (in magenta). In Fig. 8(a), δQ for the four transitions is further scaled by the probability of each initial state, e.g. po1po2(1 − po3) for |oog → |ooe . Among the four transition lines, only this one is clearly visible in Fig. 8(a).
10.1 GHz. The experimental data (left half of Fig. 8(a)) show split transition lines at δ = π. This is accounted for by the relation between the phase δ across the contact and the reduced flux ϕ, which involves the phase drop across the loop inductance, proportional to the current through the atomic contact (see Appendix D). In the presence of several channels which can be either in the ground or the odd state (we neglect the occupancy of the excited states), this current can take several values. As a consequence, the phase across the contact, common to all the channels, depends on the states occupancy, which gives rise to an effective coupling between the channels. To calculate the spectrum from theory, we use Eq. (23) (with δf m = −0.4 MHz) for each state |Ψ , and compute the weighted average S 11 = |Ψ p |Ψ S 11 (|Ψ ) using the probabilities p |Ψ for each state. The probabilities to find each channel in the odd state were taken constant, p o1 = 0.55, p o2 = 0.5, p o3 = 0.4, determined by fitting the phase dependence of S 11 (φ) at an excitation frequency where no transition is observed (f 1 = 16 GHz). In presence of the excitation pulse, the change in each p |Ψ was obtained from Eq. (25). For simplicity, we assumed T 1 = 4 µs and T 2 = 38 ns at each phase and for each state. Depending on the state of the other channels, a transition from |g to |e in one channel leads to different changes in S 11 , as illustrated in Fig. 9. The resulting calculated spectrum is shown in the right half of Fig. 8(a), in which we represent the changes δQ in the Q quadrature. Fainter lines in the data are multi-photon transitions (at f Ai /n), or transitions involving two channels (f A1 + f A2 and (f A2 + f A3 )/2), not included in the theoretical plot. Despite the strong simplifications in the analysis, the changes δA in the amplitude of the reflected signal in Fig. 8(a) are well reproduced by the calculated changes in the Q quadrature of S 11 . Similar agreement is found on the other quadrature. In Fig. 8(b) we show the single-tone spectroscopy of the resonator. The resonator frequency shift calculated under the assumption that the three channels are in the ground state is shown as a black solid line. Features associated with configurations in which one (red and green lines) or both (blue line) of the two most-transmitted channels are in the odd state are also observed. Horizontal shifts of the red and green lines with respect to the black ones result from the phase drop across the loop inductance, different in each configuration. The analysis of the data of Fig. 8 illustrates how, in a multi-channel weak link, the frequency shift associated with a transition in one channel depends strongly on the occupancy of the Andreev states in the others.

B. Nanowire junctions
The two-tone spectra shown in the left-half of Figs. 10 and 11 were measured on the same InAs nanowire device at two different gate-voltages [22]. In this experiment, the coupling between the weak link and the resonator was two orders of magnitude smaller than in the atomic contacts case, and consequently the resonator frequency shift was always smaller than its linewidth. In this limit, ∆S 11 depends linearly on the sum of the frequency shifts associated to each level, and the complications of the well-coupled system described in the previous section can be ignored. Hence, the color scales in Figs. 10 and 11 directly corresponds to the resonator frequency shift. Both spectra correspond, like in Fig. 4(a-e), to a situation in which the resonator frequency f r = 3.26 GHz is very low as compared to most of the observed transition lines.
A pair transition and several single-particle transitions are clearly recognized in the spectrum of Fig. 10. The signal at each point being integrated over 100 ms, a time much longer than the one for parity changes due to quasiparticle poisoning [28], transitions from states of different parities appear in the spectrum. From the analysis illustrated by Fig. 4(a-e), one understands that the frequency shifts corresponding to transitions above the resonator frequency are essentially given by the curvature of the transition lines: with the color scale of Fig. 10, lines are red where they have positive curvature, and blue when negative. For a more quantitative comparison with theory, we first fit the position of the bundle of single-particle transitions appearing in the range 3-10 GHz, using the scattering model presented in Ref. [22] (see parameters in Appendix F). The calculated single particle transition lines shown in the right-half of Fig. 10 reproduce well the observed single particle transition energies, but the pair transition predicted from the same Andreev levels disperse less than in the experimental data. Using these parameters, we calculate the matrix elements for H 0 (see Appendix B). Since the two-tone spectroscopy data was 10. Two-tone spectroscopy of an InAs nanowire weak link (from Ref. [22]). Left: experimental data. Colorscale represents the resonator frequency shift (sign corrected compared to Ref. [22]). A pair transition (PT) and a bundle of single particle transitions (SQPT) are pointed at. Right: calculation for a single occupied channel (see text). Colorscale is the difference in frequency shift between initial and final state. Solid line at 3.26 GHz indicates the resonator frequency. In the calculation, it was assumed that δzp = 1.2 × 10 −5 and the dissipation rate is Γ1σ + Γ 2σ = 0.62 GHz. The sign of the frequency shift in the experiment has been corrected compared to Ref. [22].
taken at very small power, the theory of weak-driving is fully applicable. We evaluate the matrix elements for the weak-driving through the gate and compute the resonator shift shown in the figure, using Eqs. (21,22). Globally, the shifts calculated for the four single particle transitions reproduce quite well the observed ones. However, some details differ, notably for the highest single quasiparticle transition, with shifts near ϕ = 0 larger in the data than in the calculation. The shift for the pair transition is reproduced only at a qualitative level. A similar procedure was used to fit the data in Fig. 11, taken at another value of the gate voltage [22], but again in a situation where the effect of states' curvature dominates. The fit parameters are given in Appendix F. In this case, the calculated pair transition lies outside the frequency range of the graph and only the bundle of four single-particle transitions is clearly recognisable (transitions in the range 13-23 GHz). Theory captures most of the features of the experiment for this set of single particle transition lines.
In the experimental results reported in Ref. [22], the resonators shift was remarkably low (tens of Hz) as com- 11. Comparison between experimental results on Andreev transitions in a nanowire junction [22] and a full calculation of the shift taking into account that the transitions were induced by microwaves applied on the gate. The parameters for the calculation are those that allowed fitting the spectrum [22] and the matrix elements for the microwaves shown in Fig.  5(c) with magenta are used. The same values of δzp and the dissipation rate are used as in Fig. 10.
pared to that observed for atomic contacts (tens of MHz) [19]. There are two reasons for this. On the one hand, the geometry of the circuit, which determines the phase fluctuations the resonator induces in the loop. It can be optimized with the circuit design. On the other hand, and more fundamental, the reduction of the matrix element of H 0 in the long-junction limit. As a rough approximation (more accurate close to δ = π), the matrix element expression reduces to that for a short junction with an effective gap ∆ eff = ∆ 1+λ ; this contributes to a (1 + λ) 2 reduction in the coupling.
In recent experiments [41] we coupled the nanowire through a shared inductance, increased the resonator frequency and diminished its impedance, leading to an enhanced δ ZP by three orders of magnitude. The data from Figs. 12, 13 and 14 have been taken using this new setup. The sample comprises a L ∼ 550 nm-long weak link from a full shell InAs-Al nanowire (same batch as the one in Ref. [22]), coupled to a resonator with f r = 6.6 GHz. The coupling being much larger than in Figs. 10 and 11, the measured signal cannot be easily converted in a resonator shift, and we show the change in one quadrature, δQ.
In the range between 5 and 10 GHz of the otherwise very busy spectrum of Fig. 12, one recognizes three pair transitions and a bundle of single particle transi- 12. (a) Two-tone spectrum measured on an InAs nanowire weak link of length L ∼ 550 nm, and using a resonator at fr = 6.6 GHz. The color-coded quadrature of the measured signal shows many sign changes along the transition lines, qualitatively in agreement with the behavior illustrated in Fig. 4(g,h): the sign changes are attributed to situations where the energy of some virtual transitions match hfr. (b) Associated single-tone spectrum.
tions. The color-coded quadrature shows many abrupt sign changes along the transition lines, like the behavior illustrated in Fig. 4(g,h): the sign changes are attributed to situations where the energy of some virtual transition matches the resonator frequency (as indicated by the black dashed lines). This spectrum has other remarkable features, like the occurrence of sets of pair transitions very close in energy. This is something we have observed many times in these nanowire weak links (see for example Fig. 13 between 11-20 GHz). Although we do not have a closed explanation yet, it could originate from subband orbital quasidegeneracies due to the approximate rotational symmetry of the weak link, or to multiple subbands occupancy. Moreover, the shape of several other lines in the spectrum does not correspond to what our simplistic model predicts, and levels anticrossings might be required to explain them.
In Fig. 13 we show data measured in a different cooldown. In this case, changes of sign of the displayed quadrature occur when the lowest transition of a singleparticle transition bundle (drawn in black based on the signal of both quadratures) crosses the resonators line. The resonator shift measured in the CW single-tone spectroscopy (Fig. 13(c)) shows a dominant contribution from the pair transition that lies at 12 GHz at δ = π, but the For each phase value, the resonator frequency is extracted from the raw data ; the shift from its bare value is shown with gray disks and compared to the calculated shift due to a single channel (dashed red), whose pair transition towards the lowest Andreev manifold fits the transition line fA shown in (a). Blue dashed line: contribution of an effective second channel (parameters in Appendix F; because of its low transmission, the corresponding pair transition does not fall in the frequency range of (a)). Black line: total shift due to both channels.
anti-crossings expected at the position of the dotted lines in Fig. 13(a) are not visible, indicating the very small occupancy of the initial state for the corresponding transition. This is in agreement with the difference in intensity between the transition lines in the two-tone spectroscopy. A remarkable feature in this spectrum is the presence of two very low-frequency lines (below 2 GHz), better seen in Fig. 13(b). By analyzing their position in energy, they can be identified as spin-flip intra-manifold transitions. The lines labelled a,b,c,d correspond to transitions between the first and the second manifold, at energies E 2d,1d , E 2u,1d , E 2d,1u , E 2u,1u . Their differences, labelled a-c, b-d, a-b and c-d, coincide two by two. They are shown with red lines and perfectly match the observed low-frequency transitions. The two first ones correspond to a transition energy E 1u,1d , the two last ones to E 2d,2u . The lines are dimmer at low frequency because the matrix elements go to zero at phases 0 and π (see Fig. 5(d)), and because the difference in occupancy of the two spin states diminishes when their energy difference is comparable to temperature: k B T /h ≈ 0.8 GHz. Note that such transitions have been recently driven indirectly through Raman processes [27]. In general, a complete fit of the spectra found in nanowire weak links is not possible with a simple modelling of the weak link. However, in the absence of a drive, the frequency shift of the resonator is often dominated by the contribution from a single channel, which allows for a simpler description. In Fig. 14(a) and (b), we show the two-tone spectroscopy and single-tone measurement of the same nanowire weak link as in Figs. 12, again from another cooldown. Among the observed transition lines, there is a high-contrast pair transition that crosses the resonator at ϕ = π(1 ± 0.12). Within the scattering model (Appendix B), it can be fitted as a pair transition towards the lowest of three Andreev manifolds arising from a high-transmission channel (τ = 0.996). In Fig. 14(a) we indicate this fit with a red line on the right hand side (parameters in Appendix F), as well as two replicas shifted by ±f r , also visible in the data and associated to a strong measurement tone. The corresponding shift of the resonator, fitted with Eq. (2) and (3) with δ zp = 0.012 and using f r = 6.60762 GHz (bare frequency measured when the nanowire is fully depleted), is shown in dashed red in the right hand side of Fig. 14(b). Although it contains contributions from both the continuum and the three Andreev manifolds, it is mainly dominated by the H 0 contribution at energy E −1,1 associated to the transition to the lowest manifold. Therefore it can be well approximated within a simplified Jaynes-Cummings description, taking Eq. (12) with a renormalized gap ∆ eff /h = 15.4 GHz and f A = 2E −1,1 (not shown in the figure for clarity, because it coincides almost exactly with the full theory for the channel shown in dashed red). Although it does not fit perfectly, it offers a simple analytical form that captures well the main features of the data around the anticrossing. The small discrepancies with the experimental data are attributed to the other channels. Many other transitions are indeed visible in the two-tone spectroscopy, which we model with an effective second channel (parameters in Appendix F). Its contribution, shown with a dashed blue line in Fig. 14(b), produces the smooth background that, added to the shift from the main transition, accounts precisely for the data for all phases (black line).

VI. CONCLUSIONS
We have developed a general formalism to describe the readout of the states of a phase-biased superconducting weak link coupled to a microwave resonator in a circuit-QED setup, based on Ref. [29]. The fermionic character of the weak link excitations gives rise to a rich phenomenology in the presence of a driving field, with different responses depending on the states' parity. We show how both spin-conserving or non-conserving transitions can occur, depending on how the system is driven.
We have applied this theory to analyze experimental results obtained on atomic contacts and semiconducting nanowire junctions. For the former, only one pair of spin-degenerate Andreev levels exists per channel, which greatly simplifies the theory, leading to analytical results. Quantitative agreement with the data, both for the ground state properties and the spectroscopy, is reached. For the nanowire case, the junctions host many Andreev states in each channel and spin-degeneracy is broken by spin-orbit coupling. This requires a more involved analysis and the comparison with experiments is less straightforward. However, with a simple modeling of these junctions we were able to interpret several features observed in the experiment, such as the effect on the resonator shift of the levels' curvature and the crossing of virtual transitions with the resonator frequency. In addition, we reported data on longer nanowire weak links with a stronger coupling to a resonator, in which direct intra-manifold spin-flip transitions of a single quasiparticle in a superconducting weak link have been observed in the absence of a Zeeman field. The spectra also suggest that Andreev levels in different channels are coupled, and display groups of similar transitions that require further modelling.
Appendix A: Theory for the frequency shift of a resonator coupled to a driven weak link We provide a detailed derivation of Eqs. (21) and (22) in Sec. IV.

Dissipation in the resonator and the weak link
For more realistic description, we assume that external baths are coupled to the system, leading to photon losses in the resonator and relaxation in the weak link. For the resonator, the time evolution of an operator O including the photon loss can be described by the master equation where E g = (1/2) jσ <0 E jσ is the ground state energy.
In this representation, in order to consider only states with positive energy we used the relation γ † iσ = −sγ −iσ , where s = 1(−1) for σ = u(d). The master equation taking into account the relaxation of a quasiparticle from one energy level to the next lower one is given by where ω iσ = E iσ / , Γ jσ is the relaxation rate and L jσ = γ † jσ γ j−1σ . If we truncate the equations at first order in Γ jσ , we obtain approximate solutions, and n iσ (0) is the occupancy of the iσ Andreev level at t = 0. For simplicity, we assume that n i−1σ (0) and n i+1σ (0) are zero, that is,Γ iσ = Γ iσ .

Spectral function of the resonator
The shift of the resonator coupled to the driven weak link in a many-body state can be obtained by using a retarded Green's function of the resonator, In the absence of the coupling to the weak link, the Green's function is given by D R (ω) = 1/(ω − ω r + iκ/2) where we used the form of a(t) in Eq. (A2). Let us now consider the resonator-weak link coupling Hamiltonian [29] including the driving field, whereĤ ( ) 0 (δ) is given bŷ and the driving HamiltonianÂ(t) iŝ In the interaction picture, the expectation value in Eq. (A7) is written as where the subscript 'I' denotes the interaction picture where creation and annhilation operators evolve in time according to Eqs. (A2) and (A5), and ρ = ρ R ⊗ρ WL is the density matrix at time t = 0 with ρ R for the resonator and ρ WL for the weak link, where Here, N is the mean photon number of the resonator and the many-body state of the weak link at t = 0 is either the ground (|g ) or the excited odd (|i 0 σ 0 = γ † i0σ0 |g ) state. The time evolution operatorÛ I in Eq. (A12) iŝ where T denotes the time-ordered product. We approximate the evolution operator to obtain δ 2 zp -order term (the shift without the driving field) and δ 2 zp |A iσ,jσ | 2 -order term (the shift with the driving field), requiring to expand up to fourth order inP I (t ). We assume that the coupling coefficient δ zp , the strength of the driving field and the dissipation in each system are sufficiently weak such that a leading order correction to the expectation value of (A12) can be written as, Below, we calculate this shifted frequency δω r , and provide the results given in the main text. Hereafter, all operators will be in the interaction picture (described by Eqs. (A2) and (A5)) and we omit the subscript 'I' for simplicity unless stated otherwise.
The perturbation expansion ofÛ † (t)a(t)Û (t) and keeping terms that contribute to the frequency shfit lead tô where C H (t) and C H 2 (t) are δ 2 zp -order terms, where θ ij = θ(t i − t j ) is the Heaviside step-function and The field-dependent terms C H A 2 (t) and C H 2 A 2 (t) are given by where the subscript numbers n = 1, ...4 indicate the time dependence, Substituting Eq. (A17) into Eq. (A16) and performing the integrations with respect to the given many-body state, we find the shifted frequency δω r , where δω X , with X ∈ H , H 2 , H A 2 , H 2 A 2 , are the frequency assciated with C X (t), For ρ WL = |g g|, they are given by Here, E jσ ,iσ = E iσ − E jσ is the energy difference between the levels iσ and jσ , andẼ r is the resonant energy of the resonator including the effect of the dissipation, E r = hf r − i(Γ iσ + Γ jσ + κ)/2. The frequency shift for ρ WL = |i 0 σ 0 i 0 σ 0 | is obtained as First, note that in absence of drive, the coupling to the weak link induces a shift in the resonator's frequency which is which is a generalized version of Eq. (4) taking into account the intrinsic resonator's linewidth and the inherited contribution due to the finite life-time of Andreev states. To see this, note that the last term can be written as −|h iσ,jσ | 2 1 E iσ,jσ − hf r + i(Γ iσ + Γ jσ + κ)/2 and has the form of Jaynes-Cummings with dissipation. The contribution from the drive is which for example, if the drive frequency is close to |g → |1u1d transition, the term with the denominator which corresponds to the shift in |1u1d multiplied by a lorenztian coefficient associated with the probability of the transition, here the matrix element is the Rabi frequency. So the result is very intuitive: once the shifts δω |Ψ r dr off for the different many-body configurations are determined, the effect of the drive is to take into account selection rules for transitions to occur.
The continuum states Φ iσ at |E| ≥ ∆ are degenerate with i = (i p , i d ) where i p = e, h denote electron-like or hole-like state and i d = l, r are the left or right source. They are normalized by the local density at the source region, where the group velocity v g (E) at energy E is given by v 1 √ E 2 − ∆ 2 /|E| for (σ, i p , i d ) = (u, e, l), (u, h, r), (d, e, r), (d, h, l), and v 2 √ E 2 − ∆ 2 /|E| for (σ, i p , i d ) = (u, e, r), (u, h, l), (d, e, l), (d, h, r).
The current operator H 0 = ∂H 0 /∂δ of the linearized band is given by where the factor 1/2 reflects the particle hole redundancy andτ 0 is the identity matrix in the Nambu space. We find the matrix element of the current operator between states iσ and jσ . For pseudospin-up states, σ, σ = u, we obtain and for pseudospin-down sates, σ, σ = d,
In runs I and III, the bare resonator frequency fr was measured when the contact was open, whereas for run II it was a fit parameter. The largest transmission τ1 is the essential fit parameter, determining the overall shape. In the last column, we indicate the minimal value of the Andreev frequency associated to this transmission f min A1 = fA1(π) = 2∆ √ 1 − τ1. When f min A1 < fr (c-g), one observes an avoided crossing. A second channel with transmission τ2, and for Fig. 6(a',b') a third one with transmission τ3 taken equal to τ2, essentially accounts for an overall shift of the resonator frequency. From the estimated loop inductance ≈ 0.1 nH [40], we obtain the screening parameter β = 0.03 used to fit the data.

Appendix F: Spectra parameters
Most of the spectra (illustrative of fit) were obtained from Eq. (B12). We give in Table II the parameters corresponding to the various figures.

ACKNOWLEDGMENTS
Technical support from P. Sénat is gratefully acknowledged. We thank P. Orfila and S. Delprat for nanofabrication support, and our colleagues from the Quantronics group and C. Strunk for useful discussions. We