Broken selection rule in the quantum Rabi model

Understanding the interaction between light and matter is very relevant for fundamental studies of quantum electrodynamics and for the development of quantum technologies. The quantum Rabi model captures the physics of a single atom interacting with a single photon at all regimes of coupling strength. We report the spectroscopic observation of a resonant transition that breaks a selection rule in the quantum Rabi model, implemented using an LC resonator and an artificial atom, a superconducting qubit. The eigenstates of the system consist of a superposition of bare qubit-resonator states with a relative sign. When the qubit-resonator coupling strength is negligible compared to their own frequencies, the matrix element between excited eigenstates of different sign is very small in presence of a resonator drive, establishing a sign-preserving selection rule. Here, our qubit-resonator system operates in the ultrastrong coupling regime, where the coupling strength is 10% of the resonator frequency, allowing sign-changing transitions to be activated and, therefore, detected. This work shows that sign-changing transitions are an unambiguous, distinctive signature of systems operating in the ultrastrong coupling regime of the quantum Rabi model. These results pave the way to further studies of sign-preserving selection rules in multiqubit and multiphoton models.

. The strength of the interaction g between qubit and oscillator defines different regimes of coupling, reflecting the dominating energy scales in the dynamics of the system. The limit of coupling strength g much smaller than qubit/resonator frequencies ω ω  g { , } r q is the well-known Jaynes-Cummings (JC) model 3 , where the rotating-wave approximation (RWA) holds. The ultrastrong coupling (USC) regime is defined by g/ω r ≥ 0.1. Such condition implies that counter-rotating terms in the interaction Hamiltonian cannot be neglected. The USC regime presents interesting features such as a non-trivial ground state with photonic excitations 4 , the generation of qubit-oscillator nonclassical states 5 and ultrafast gates for quantum computation 6 . The intricate ground-state structure has motivated various proposals for quantum simulation 7,8 . Despite its fundamental interest and applications, the experimental study of the USC regime is challenging. In this sense, superconducting quantum circuits 9 are suitable systems to achieve large coupling strengths, despite some intrinsic limitations for most superconducting qubit types 10 . In fact, recent ideas have emerged [11][12][13] on how to reach the USC regime using superconducting flux qubits. The latter can attain a large galvanic coupling to resonators 14 and have a huge anharmonicity, typically larger than ω π  20GHz /2 q , both of which are crucial to reach g/ω r > 0.1. USC in an open system has recently been reported using a superconducting flux qubit 15 .
In the JC model ω  g/ 1 r , the energy-level structure consists of a manifold of dressed-state doublets each containing n excitations |n, ± 〉 3 . Each doublet |n〉 is labeled with a quantum number ± corresponding to the relative sign in the superposition of uncoupled qubit-resonator states {|e, n − 1〉 , |g, n〉 }: c os( /2) , 1 sin( /2) , , The mixing angle is defined as θ ω ω = − g n tan 2 /( ) n q r . In the case of resonator driving represented by the operator H d ~ (a + a † ), transitions between dressed states of different sign, ± ↔ ±  n n , 1 , JC JC , are much less favored than those preserving it, |n, ± 〉 JC ↔ |n ± 1, ± 〉 JC . For USC rates, g/ω r ~ 0.1, one can still effectively describe the spectrum of the system using dressed-state doublets ± ∼ n, with correspondingly larger transition matrix elements and corrections of order g/(ω r + ω q ) to manifolds with total number of excitations n ± 2 [16][17][18][19] . The considerations on sign-changing transitions in the USC regime do not apply. Thus for coupling strength much smaller than qubit/resonator frequencies, one can establish a selection rule for the case of driven resonator H d ~ (a + a † ). In the USC regime, due to the different nature of the system eigenstates, the selection rule may be broken. For qubit driving H d ~ σ x , the selection rules are slightly different than the resonator driving case, as detailed in the results section.
Observations of sign-preserving transitions have been reported for superconducting transmon qubits coupled to transmission line resonators with g/ω r ~ 0.02 20,21 . To date, only a few experimental systems have reached the USC regime. In superconducting circuits, artificial atoms have been ultrastrongly coupled to resonators 14,[22][23][24] and transmission lines 15 , where new physics beyond the RWA were reported. Polaritons in quantum wells have been ultrastrongly coupled to microcavities, causing large frequency shifts due to enhanced nonlinearities 25 . Finally, electrons in a two-dimensional electron gas undergoing cyclotron transitions ω c were coupled to metamaterial cavities where the coupling strength achieved was g/ω c ~ 0.58 26 . None of the above cited experiments reported observations or was able to provide evidence of sign-changing transitions.
Here, we report the spectroscopy of a qubit-resonator system in the USC regime of the quantum Rabi model. We use the fact that our system is at finite temperature to excite several transitions between dressed states. One of the transitions is shown to break the sign-preserving selection rule, demonstrating another paradigmatic example of the unique physics occurring in the USC regime. Using a superconducting flux qubit coupled to a lumped-element LC resonator we are able to reproduce the quantum Rabi model of a qubit coupled to a single oscillator 1 . In our experimental system, the qubit and the resonator have no higher energy levels within the relevant frequency range, which is advantageous for observing the breakdown of the sign-selection rule.

Results
Theoretical model. Our experiment operates in the regime of g/ω r ~ 0.1, which will be referred to as the quantum Bloch-Siegert (QBS) regime. The eigenstates of the system in this regime are analogous to the doublet structure of dressed states of the JC model (see the methods section for more details): With respect to the JC model, a correction of order λ = g/(ω q + ω r ) appears connecting to states containing n ± 2 number of excitations due to the effect of the counter-rotating terms. The mixing angle is here defined as n n n which is different from the JC in a photon-number dependent detuning δ n ≡ δ + 2nω BS and a photon-number dependent coupling strength ω ω ω ≡ − + g n g n g /( ) n r q 3/2 BS . δ ≡ ω q − ω r is the detuning in the JC model, ω BS ≡ g 2 /(ω r + ω q ) is the Bloch-Siegert shift that originates from the counter-rotating terms. The existence of a selection rule for transitions between dressed states of different sign depends on the relative magnitude of matrix elements for sign-changing/preserving transitions.  cos( /2) cos( /2) 1 (1 ( 2)) sin( /2) cos( /2),  Clearly, the sign-preserving transitions are always possible while the sign-changing transitions are less likely to be excited. For vanishing qubit-resonator coupling g and δ > 0 (δ < 0), ϕ n → 0 (ϕ n → π), so sin(ϕ n /2) → 0 (sin(ϕ n /2) → 1), cos(ϕ n /2) → 1 (cos(ϕ n /2) → 0). This implies that sign-changing transitions disappear when the coupling strength is much smaller than qubit/resonator frequencies, establishing a selection rule for this type of transitions.
For single qubit driving H d / = A q σ x cos(ω d t) the transition matrix elements σ ↔ T i j x between dressed states read For δ < 0, as is the case in our experiment, ϕ → as g → 0, implying a sign-preserving rule. However, the other transitions do not follow the same rule since Clearly for single-qubit driving, the selection rule considerations of driven resonator do not apply. In the rest of this work we only discuss the transitions |1, − 〉 ↔ |2, − 〉 and |1, − 〉 ↔ |2, + 〉 that follow the sign-preserving rule both for resonator-driven as well as qubit-driven transitions. The case of single-qubit σ z driving is described in the Methods. Therefore all our conclusions on broken selection rules still hold. Throughout the rest of the text we will drop the tilde to refer to the dressed states in the USC regime simply as |n, ± 〉 .
Analyzing spectra in the USC regime requires accounting for the out-of-equilibrium quantum dynamics featuring the coloured nature of dissipative baths. One possibility is to project the master equation in the dressed-state basis |n, ± 〉 19 . Another possibility is the second-order time-convolutionless projection operator method (TCPOM) 34 , which has been proven useful to correctly describe the open system dynamics in the USC regime 18,35 . Here, we follow the latter approach where the master equation reads r a x z r r r r r r r r 0 d rive , , We consider thermal baths at temperature T for all dissipative channels, where the relaxation coefficients depend on the spectral density of the baths d k (ω) and the system-bath coupling strength α k (ω). Moreover, the system Hamiltonian includes the qubit-resonator interaction as well as the presence of classical microwave fields, that is The qubit mixing angle θ q is defined as tan θ q = Δ /ε, with Δ being the tunnel coupling between the qubit eigenstates and ε = 2I p (Φ − Φ 0 /2) the magnetic energy of the qubit (see Methods). The out-of-equilibrium dynamics assumes a perturbative treatment of driving mechanisms which is valid for amplitudes 36 , as is the case throughout this work. In all numerical simulations we use up to 6 Fock states in the resonator basis which is sufficient for g/ω r ~ 0.1.
Scientific RepoRts | 6:26720 | DOI: 10.1038/srep26720 Sign-preserving transition. In order to explore dressed-state transitions, we use a superconducting flux qubit galvanically coupled to an LC resonator. The results presented here are obtained with the same device where the Bloch-Siegert was first reported 14 . All relevant experimental details can be found there. The resonator has frequency ω r /2π = 8.13 GHz, while the qubit has persistent current I p = 500 nA and a minimal splitting of Δ /2π = 4.2 GHz, so the detuning is δ = Δ − ω r = − 2π × 3.93 GHz at the symmetry point of the qubit. The coupling strength is g/2π = 0.82 GHz. The latter implies ω .  g/ 01 r , so the system operates in the USC regime. The qubit and the resonator are thermally anchored to the mixing chamber stage of a dilution refrigerator at 20 mK. The qubit state is read out using a DC-SQUID magnetometer. An external superconducting coil is used to set the magnetic flux in the qubit loop near Φ = Φ 0 /2. At this flux bias point the effects of the counter-rotating terms of the Hamiltonian are strongest (equation (16) with sin θ q = 1).
A numerical diagonalization of the Hamiltonian in equation (16) leads to the spectrum shown in Fig. 1(a). The solid lines correspond to transitions from the ground state, while the dashed lines are transitions from the first excited state. The intermediate dashed transition |1, − 〉 ↔ |2, − 〉 (see Fig. 1(b)) around 8 GHz is the one studied in detail in this section and corresponds to a sign-preserving transition. The highest-energy, dashed transition |1, − 〉 ↔ |2, + 〉 is studied in the next section and is the first transition in ascending energy to break the sign-preserving selection rule, given our choice of parameters. At the qubit symmetry point Φ = Φ 0 /2 only transitions connecting states of different parity are allowed 37 . Therefore, the red sideband |1, − 〉 ↔ |1, + 〉 and blue sideband |0, g〉 ↔ |2, − 〉 are forbidden at that point. Outside the symmetry point all transitions are permitted due to the asymmetry in the flux qubit potential 38,39 . At the qubit symmetry point Φ = Φ 0 /2, the system eigenstates correspond to those in equations (3) and (4) and the selection rules introduced in the previous section are therefore expected. Figure 2(a) shows the spectrum of the system at low power near the resonator frequency ω r /2π = 8.13 GHz. The solid and dashed curves are a fit of equation (16). The dotted line is the JC model with the same fitted parameters. The difference between the JC and Rabi models is the Bloch-Siegert shift ω BS 14 , with a maximum of 55 MHz at Φ = Φ 0 /2. A negative (positive) shift of the frequency corresponds to a resonator (qubit)-like dressed state transition. Figure 2(b) shows a spectroscopy trace at Φ = Φ 0 /2. The measurements are performed following a similar protocol as the one developed previously for a tunable-gap flux qubit 40 (see also Methods). Following from Fig. 1, we find that the higher frequency resonance in Fig. 2(b) corresponds to the transition |0, g〉 ↔ |1, + 〉 , which given the detuning δ/2π = − 3.93 GHz is mostly a resonator-like excited dressed state. The lower frequency resonance corresponds to |1, − 〉 ↔ |2, − 〉 , therefore it is a sign-preserving transition. The finite temperature in our system enables the excitation of this transition. The signal from this resonance is weak given the small amount of population of qubit-like state |1, − 〉 in thermal equilibrium. The fact that we can observe this transition at low driving strength indicates a large matrix element as predicted by equations (6)(7)(8)(9)(10)(11)(12)(13). The low signal-to-noise ratio leads to spurious peaks and dips which do not represent additional transitions, as verified by the spectrum in Fig. 1. The geometry of our system is such that qubit and resonator driving amplitudes have a comparable value, , therefore both systems are simultaneously driven in all measurements presented. The resonances in this spectrum look broad due to the low relaxation time of this particular device, mostly due parasitic capacitance to the SQUID readout circuitry. Other flux qubit experiments have shown much better performance, even with SQUID detectors 41,42 . There is a clear asymmetry between the two resonances in Fig. 2(b). This difference allows us to approximately calibrate the electronic temperature of the system. Solving the master equation iteratively to reach steady-state is a lengthy process, making it impractical to run a fit of the observed spectrum. Instead, the unknown parameters are swept in order to find the optimal set that leads to a calculated spectrum closest to the measurements. We plot the expectation value of the qubit current operator ρσ ∼ˆÎ Tr( ) z circ , where ρ is the steady-state solution of the TCPOM master equation, equation (14) using the eigenstates from equation (16). More details can be found in the Methods. The optimal set of parameters is Scientific RepoRts | 6:26720 | DOI: 10.1038/srep26720 . For the temperature, we plot T = 90 mK (blue asterisks) and T = 150 mK (red circles) to indicate the possible range of temperatures. Other experiments using superconducting qubits have observed comparable finite electronic temperatures 43 .
The results of the simulation clearly indicate an asymmetry between the two resonances near 8 GHz. A global scaling factor and an offset have been applied to the numerical solution to match the amplitude of the experimental resonances. The scaling factor gives the transduction between average persistent current in the qubit to SQUID switching probability. In the experiment, the offset value is chosen to be 50% probability of the DC-SQUID to switch to the running state where a finite voltage is generated 42 . The differences between simulation and experiment for the resonance at 8.02 GHz is most likely due to the low signal-to-noise, leading to an imprecise determination of the temperature of the system. Broken sign-preserving selection rule. In this section, we perform spectroscopy of the qubit-resonator system with driving amplitude approximately 5 times larger than the one used in Fig. 2. The enhanced driving strength enables additional transitions to occur, both one-and two-photon transitions. The high-power spectrum of the system is shown in Fig. 3(a), showing a rich spectral structure. For frequencies below 8 GHz, multi-photon processes can be identified. The frequency of each multi-photon transition follows the proper scaling ω ω θ = − − n g / s in / 2 n r q , as was already reported for a transmon qubit coupled to a resonator 44 . θ q is the qubit mixing angle defined below equation (16). The resonance at 8.25 GHz from Fig. 2 has been power-broadened, overlapping with the sign-preserving transition |1, − 〉 ↔ |2, − 〉 at 8.02 GHz. The large signal generated by this resonance may be related to a buildup of the population in the resonator due to the strong drive 45 .
Above 8 GHz, there are two distinctive resonances near Φ Φ  /2 0 . The transition increasing in frequency away from the symmetry point corresponds to the blue sideband, |0, g〉 ↔ |2, − 〉 . The reason it vanishes at Φ = Φ 0 /2 is due to parity selection rules 37 . The other weaker resonance corresponds to the transition |1, − 〉 ↔ |2, + 〉 and is, therefore, a sign-changing transition that breaks the sign-preserving selection rule explained in the introductory section. To our knowledge, this is the first observation of this kind of transitions between excited states of the quantum Rabi model. The small signal from this resonance evidences the difference in magnitude of the transition matrix elements as compared to the one for |1, − 〉 ↔ |2, − 〉 given the difference in driving amplitudes used to acquire the data in Figs 2 and 3.
The spectroscopic sign of any resonance detected by the SQUID magnetometer relates to the magnetic field generated by the qubit in steady-state when an external drive is present. The transition |1, − 〉 ↔ |2, + 〉 increases the switching probability above the reference at 50%, compared to all other resonances where the switching probability decreases below 50%. This implies that in steady-state the qubit is more polarized to the ground state than in thermal equilibrium, contrary to the other resonances where the qubit ends up with more population in the excited state. The resonance |1, − 〉 ↔ |2, + 〉 therefore cools the qubit down by transferring its excess thermal energy to the resonator, since the state |2, + 〉 is a resonator-like excited dressed state. A spectroscopy trace at the symmetry point can be observed in Fig. 3(b). The switching probability at 12.3 GHz increases with respect to the background set at P sw = 0% (referenced to 50%).

Discussion
In Fig. 4(a,b)   . (b) Spectrum at Φ = Φ 0 /2. The switching probability of the SQUID detector is referenced to 50%. A numerical simulation of master equation equation (14) is shown for T = 90 mK (T = 150 mK) with blue asterisks (red circles). The low amplitude of the sign-preserving transition at 8.02 GHz is related to the qubit thermal population in equilibrium. equations (7), (8), (11) and (12) and the red-dotted line is the Jaynes-Cummings model calculated from the states in equations (1) and (2). In Fig. 4(c,d) we also plot the matrix elements T |1−〉↔|2−〉 , T |1−〉↔|2+〉 for the case of single-qubit driving H d ~ σ x . The difference in matrix elements for sign changing/preserving transitions is very clear for both types of driving. Figure 4(b,d) show that in our experiment the broken-sign transition is mostly driven by the direct qubit drive. As shown in Fig. 4(e), the relative weight T |1−〉↔|2−〉 /T |1−〉↔|2+〉 decreases with increasing g/ω r for both types of driving since we consider finite detuning δ/2π = − 3.91 GHz. If we had the resonance condition ω r = ω q = 2π × 8.13 GHz, the relative weight of the matrix elements for resonator (qubit) driving would reach T |1,−〉↔|2,−〉 /T |1,−〉↔|2,+〉 ≈ 7(1) at g/ω r = 0.1, instead of 234(4.8) as shown in Fig. 4(e). The dashed vertical line in Fig. 4(a-e) marks the operating point of our device. The observation of the sign-changing transition |1, − 〉 ↔ |2, + 〉 in this experiment can therefore be attributed to the system operating in the USC regime.
In Fig. 5 we compare the measured spectrum near 12 GHz (Fig. 5(a)) with a simulation (Fig. 5(b)) with the TCPOM master equation (14) using larger drive amplitude, A qb /2π = A r /2π = 50 MHz, keeping the rest of terms equal as those found in the previous section, with T = 100 mK. The increased driving amplitude used here is approximately a factor of 5 larger compared to the amplitude used for the simulations in Fig. 2(b), closely matching the ratio of driving amplitudes used in the measurements of Fig. 3 relative to those in Fig. 2. The results are plotted in Fig. 5(b) for the probability to find the system excited, where ρ is the density matrix of the coupled system calculated with the master equation in steady-state. The position of the main two resonances appearing in the experiment in Fig. 5(a) are clearly reproduced up to ~− 1 mΦ 0 , particularly at the symmetry point Φ = Φ 0 /2 where only the sign-changing transition |1, − 〉 ↔ |2, + 〉 is excited. Beyond − 1 mΦ 0 , the experimental signal in Fig. 5(a) disappears due to an increase of the dephasing rate of the qubit away from the symmetry point 41,46,47 which is not taken into account in the numerical simulations.

Conclusions
We report the observation of transitions between excited states of a superconducting flux qubit-resonator system in the ultrastrong coupling regime. The strength of the coupling combined with initial thermal population of the qubit permits the observation of a transition between excited states that was never detected before. Such transition changes the relative sign in the superposition of bare qubit-resonator states in the dressed-state level structure. We developed a theoretical model based on the time convolutionless projection operator method that reproduces all transitions observed spectroscopically. Our work, therefore, verifies the existence of sign-changing transitions in the ultrastrong coupling regime of the quantum Rabi model, despite their weak strength. Until now, evidence for superconducting qubit-resonator systems to operate in the USC regime was found as spectral deviations from the Jaynes-Cummings model 14,22,23 . We instead put forward the observation of sign-changing transitions in the quantum Rabi model as a direct signature of any physical system in the USC regime for resonator-type driving. Our work can be extended to multi-qubit 33 or multi-photon 48 quantum Rabi models where forbidden transitions play a key role in understanding and manipulating the energy-level structure of the system.

Methods
Quantum Bloch-Siegert Hamiltonian and eigenenergies. For a qubit coupled to a resonator the Hamiltonian that describes the full system dynamics is the quantum Rabi model 1 . When g/ω r , ω  g/ 1 q the rotating-wave approximation (RWA) holds and the Jaynes-Cummings (JC) model is obtained 3 . In the regime   (14). The z-scale corresponds to probability to find the system excited, where ρ is the steady-state solution of the master equation. The transition |1, − 〉 ↔ |2, + 〉 breaking the sign-preserving selection rule is clearly visible at all fluxes while the blue sideband generates much more signal away from the symmetry point. For this simulation T = 100 mK, A qb /2π = A r /2π = 50 MHz.
Scientific RepoRts | 6:26720 | DOI: 10.1038/srep26720 g/ω r ~ g/ω q ~ 0.1 the RWA is no longer valid but a perturbative treatment still permits deriving the Hamiltonian of the system analytically 16,17 , as well as a model on dissipation 18,19,36 .
The Bloch-Siegert Hamiltonian derives from the quantum Rabi Hamiltonian via perturbative treatment of the parameter λ ω ω = +  g/( ) 1 r q 16,17,19 . In this perturbative regime, the counter-rotating terms are treated as off-resonant fields and the quantum Rabi model can be transformed in an effective dispersive picture. The resulting Hamiltonian is: H a a a a gn a a / 2 where the Bloch-Siegert shift ω BS is defined as with photon-dependent coupling strength ĝ n ( ) Here ω ε = ∆ + q 2 2 is the flux qubit splitting with ε = 2I p (Φ − Φ 0 /2) being the magnetic energy of the qubit in its truncated two-state Hilbert space. The qubit mixing angle θ q is defined as tan θ q = Δ /ε.
Projecting the Hamiltonian on the product state |Ψ 〉 = |i, n〉 , i = {e, g} has a box-diagonal representation that simplifies the full diagonalization: For box n, the following eigenvalues are found: g ,0 BS where n ≥ 1 represents the number of total excitations in the uncoupled basis. All levels are shifted by ω BS as compared to the JC model. In addition, both the effective detuning between levels δ + 2nω BS as well as the effective coupling constant (second term in the square root of equation (22)) depend on the number of excitations. The latter is the second-order correction in this dispersive picture. Using the following definitions: r n r , B S the eigenenergies take the form  (28) and (29) show that the energy eigenstates of the Bloch-Siegert Hamiltonian correspond to doublets of superpositions exactly as with the Jaynes-Cummings model, with coefficients having a different dependence on the number of photons n and corrections of order λ = g/(ω r + ω q ) to other states with different number of excitations. Therefore the selection rule considerations between transitions that change the sign quantum number explained in the introduction break down in the USC regime.

Dissipation dynamics.
Here, the parameters of the qubit-resonator system described in the main text are derived using numerical simulations. Besides the temperature of the bath T, there are a set of unknown parameters that need to be determined: the qubit relaxation and decoherence rates Γ 1 , Γ 2 , the resonator decay rate Γ r and the amplitudes of the external fluxes A qb and A r driving the qubit and the resonator, respectively. Estimates of the mutual inductance between qubit/resonator and microwave line permit setting A qb = A r by geometry (see Fig. 6). We also assume that qubit decoherence is governed by relaxation alone Γ 2 = Γ 1 /2 since this device shows very short T 1 ~ 10 ns times near the symmetry point Φ = Φ 0 /2.
A technique to characterize the temperature of a transmon qubit by driving Rabi oscillations was already developed 43 , but due to the short relaxation time of our qubit this technique is not suitable. We would like to emphasize that the coherence of the flux qubit used in this experiment is limited by spontaneous emission to the readout circuitry of the SQUID detector 42 . Using proper filtering 49 and a more symmetric circuit design 40 would allow better quantitative study of the spectroscopic resonances. Other experiments have shown much longer values of qubit lifetime at the symmetry point 41,46,47 , exceeding 10 μs.
We run a numerical simulation of the master equation presented in the theory section and compare the steady-state solution to the spectroscopy measurements. The numerical simulation of the master equation (Equation (14) in the main article) is performed by a Runge-Kutta method considering a Fock space of up to 6 Fock states in the resonator, sufficient for the relative coupling strength g/ω r = 0.1. The operators Û and Ŝ are numerically built by means of projecting them into the eigenstates of the qubit-resonator Hamiltonian, Notice that taking into account Equation (16) from the main article in the above expression will lead to integrals of the form Scientific RepoRts | 6:26720 | DOI: 10.1038/srep26720 The second term in the above equation stands for the Cauchy principal value. The latter leads to small Lamb shifts that we have neglected in our numerical simulation.
An example of the dynamical evolution governed by the master equation can be seen in Fig. 7. The plots show the probability to populate the lowest four energy levels from Fig. 1 Fig. 7(a,b). The value of the driving frequency ω D is resonant with the transition |0, g〉 ↔ |1, − 〉 . Figure 7(a) features small oscillations due to a driving amplitude stronger than the qubit decay rate. Figure 7(b) shows a smooth behavior given that the dissipative mechanisms surpass the driving microwave fields.
In order to compare the numerical simulations with the actual data we need to take into account that the DC-SQUID detector is sensitive to the flux generated by the qubit, Φ qb = M SQ−qb 〈 I circ 〉 , where M SQ−qb is the SQUID-qubit mutual inductance and 〈 I circ 〉 the expectation value of the circulating current operator of the qubit. The current operator has the form α ϕ = αÎ I sin( ) C circ with ϕ α the phase operator across the α junction. Using the two-level approximation of the flux qubit, σ =Î I p z circ calculated in the energy basis of the qubit, with I p being the persistent current. In this case it is easy to show that for a pure state 〈 σ z 〉 = ± cos θ q , with θ q defined in the previous section. The sign depends on whether the qubit is in the ground or the excited state, respectively. For a mixed state ρ = | | + | | p p 0 0

Measurement protocol.
To perform measurements at the symmetry point we use a protocol developed in previous experiments 40 . The qubit is initiated in an external bias flux Φ i = − 2.4 mΦ 0 where it produces a net magnetic field. Then a flux square pulse is applied to the qubit via the microwave line, followed by a long microwave burst. The length of the burst is sufficient to bring the system in steady-state for all resonances. At the end of the burst the flux pulse is ramped down until the qubit is brought back to the initial point. The rise time of the pulse (~2 ns) is slow enough that the Landau-Zener tunneling probability is below 1% when the qubit is brought across the avoided-level crossing with the resonator, so the qubit population is transferred adiabatically from the symmetry point to the readout point.

Selection rule for longitudinal driving.
A longitudinal single-qubit driving characterized by a σ z operator will not break the parity symmetry of the quantum Rabi model (QRM) at the qubit symmetry point. Correspondingly, the driving may induce transitions between states of the same parity subspace. For instance, transitions will be allowed between states − ∼ n, to + ∼ n, or to ± ± ∼ n 2, . Such transitions are activated by a proper choice of the resonance condition, with the signal amplitude being much smaller than the energy differences between higher-level states in the QRM. This assures the validity of the rotating-wave approximation. Away from the symmetry point, the QRM breaks its parity symmetry and consequently its selection rule associated with cavity and/or qubit driving.
More explicitly, let us consider the QRM with a flux qubit at its symmetry point driven by a microwave field along the longitudinal axis σ z . The total Hamiltonian reads:  Fig. 2(b) in the main article to obtain the spectroscopy traces to reproduce the experimental data.