Maxime Boissonneault

A C Doherty

F R Ong

P Bertet

D Vion

D Esteve

A Blais

Superconducting qubit as a probe of quantum fluctuations in a nonlinear resonator

In addition to their central role in quantum information processing, qubits have proven to be useful tools in a range of other applications such as enhanced quantum sensing and as spectrometers of quantum noise. Here we show that a superconducting qubit strongly coupled to a nonlinear resonator can act as a probe of quantum fluctuations of the intra-resonator field. Building on previous work [M. Boissoneault et al. Phys. Rev. A 85, 022305 (2012)], we derive an effective master equation for the qubit which takes into account squeezing of the resonator field. We show how sidebands in the qubit excitation spectrum that are predicted by this model can reveal information about squeezing and quantum heating. The main results of this paper have already been successfully compared to experimental data [F. R. Ong et al. Phys. Rev. Lett. 110, 047001 (2013)] and we present here the details of the derivations.

I. INTRODUCTION

Nonlinearity in oscillators was first observed by Huygens who discovered that large oscillations in pendulum clocks introduced inaccuracies because of the resulting change in natural oscillation frequency [START_REF] Kovacic | The Duffing Equation: Nonlinear Oscillators and their Behaviour[END_REF]. It is however most famously Duffing, in the context of combustion engines, who tackled the problem of nonlinear oscillators in a systematic way [START_REF] Duffing | Erzwungene Schwingungen bei Veränderlicher Eigenfrequenz[END_REF]. Although they have a long history in physics, nonlinear oscillators still manage to surprise and are the focus of intense research [START_REF] Boyd | Nonlinear Optics[END_REF][START_REF] Dykman | Fluctuating Nonlinear Oscillators[END_REF]. This is particularly true in optics where optical nonlinearity can be realized by taking advantage of the change in index of refraction of certain media with light intensity. This nonlinearity can lead to frequency down-and up-conversion, and parametric oscillations and amplification [START_REF] Boyd | Nonlinear Optics[END_REF]. Nonlinearities produced in this context are however rather weak, and nonlinear phenomena at optical frequencies are therefore revealed mostly under high pumping intensities.

The situation is quite different with superconducting circuits where very strong nonlinearities at microwave frequencies can be achieved [START_REF] Bourassa | [END_REF], revealing nonlinear behaviour even at the single-photon level [6][7][8]. These circuits are based on Josephson junctions embedded in otherwise linear circuit elements to create superconducting qubits and nonlinear microwave resonators. These resonators can take various forms, ranging from LC-circuits where the inductance is replaced by a Josephson junction [9], stripline resonators with an embedded Josephson junction [START_REF] Boaknin | Dispersive microwave bifurcation of a superconducting resonator cav[END_REF], and to metamaterial resonators where the central resonator conductor is replaced by an array of Josephson junctions [START_REF] Castellanos-Beltran | [END_REF].

These superconducting nonlinear resonators have proven themselves to be valuable tools the context of circuit quantum electrodynamics (cQED) where one couples * maxime.boissonneault@usherbrooke.ca microwave resonators to superconducting qubits [12,13]. In this context, nonlinear resonators have for example been used has parametric [14][15][16] or bifurcation amplifiers [START_REF] Siddiqi | Quantum Computing in Solid State Systems[END_REF][START_REF] Lupaşcu | [END_REF][19] in qubit state measurement. There, one is interested in the information contained about the state of the qubit in the field at the output of the resonator. In the present paper we take the converse point of view: we show how the qubit can be used as a probe of quantum fluctuations of the field inside the resonator.

The theory presented below was developed in parallel to, and already tested against, the experimental results of Ref. [20]. The goal of the present paper is thus to give the details of the derivation of the model whose main results can be found in Ref. [20]. Moreover, the present work is based on the same experimental setup as studied in Ref. [21] and builds on previous calculations presented in Ref. [22] -which we will refer to as Paper I from now on.

In Paper I, we have developed a model describing the measurement backaction of a driven nonlinear resonator on a qubit strongly coupled to the resonator. This model went beyond many approximations that are standard in the literature. First, we considered a many-level instead of a two-level Hilbert space for the superconducting qubit. Second, we took into account the fact that the ac-Stark shift on the qubit caused by a strong pump on the resonator depends on the detuning of the qubit to the pump, and not the qubit-resonator detuning as is usually assumed [12]. Finally, our model went beyond the standard linear response theory for the qubit-state dependance of the resonator state. This model was compared with the experimental results presented in Ref. [21] and was found to be in quantitative agreement with the measured qubit's ac-Stark shift before and after bifurcation of the resonator. We also found excellent agreement with the nontrivial dependence of the qubit's measurementinduced dephasing on the pump power. From these results, we have concluded that the system is close to the quantum limit for measurement in a parameter range.

The model developed in Paper I however made two main approximations that we relax here: small separation between the two pointer states of the resonator corresponding to the two qubit states, and weak squeezing of the resonator field. We still assume that the pointer state separation and the squeezing are both relatively small, but take into account first order corrections to these approximations. As we will show, relaxing the first approximation leads to combined qubit-resonator transitions, i.e. red and blue sideband transitions. Relaxing the second approximation allows us to take into account squeezing of the resonator field and, as predicted by the theory of quantum heating [23][24][25][26], this leads to an effective temperature of the resonator field. Using sideband spectroscopy, a standard tool in ion-trapping experiments [27,28], we then discuss how the qubit can act as an absolute thermometer of this effective temperature. In practice and as discussed in more details below, the main difference between the experimental results of Ref. [21] and those of Ref. [20] is a drive of increased amplitude on the qubit.

The production of squeezed light by nonlinear microwave resonators has of course already been studied before [29,30]. While these studies focussed on the light at the output of the resonator, as mentioned above here we are focussing on the light inside the resonator. Moreover, our work adds to an already quite extensive literature concerning nonlinear resonators -see for example Refs [START_REF] Boyd | Nonlinear Optics[END_REF][START_REF] Dykman | Fluctuating Nonlinear Oscillators[END_REF]23,[31][32][START_REF] Dykman | Applications of Nonlinear Dynamics: Model and Design of Complex Systems[END_REF][START_REF] Yurke | [END_REF][35][36]. Here however the usual assumptions, such as very small nonlinearities, small qubitresonator dispersive coupling, or strictly two-level qubits cannot be made when comparing to the experimental results of Ref. [20], and they are avoided in this paper.

The paper is organized as follow. In section II, we present the system's bare Hamiltonian and master equation, and introduce the notation for the nonlinear resonator and the qubit. In section III, we summarize the calculation presented in Paper I and highlight the main approximations and results that were obtained. In section IV, we relax the small distinguishability and small squeezing approximations and obtain first order corrections. In section V, we compare our model to experimental results first presented in Ref. [20]. Concluding remarks are made in section VI while details of some of the calculations can be found in Appendices A and B.

II. PRESENTATION OF THE SYSTEM

As discussed in the Introduction, we consider a nonlinear resonator strongly coupled to a superconducting qubit. An example of such a system is illustrated in Fig. 1, where a transmon qubit [37] is coupled to a coplanar transmission-line resonator rendered nonlinear by Josephson junction embedded in the resonators' center conductor. We introduce in section II A the notation used for the nonlinear resonator and in section II B we FIG. 1. (Color online) Schematic representation of a possible implementation of the system considered in this paper and realized in Ref. [20]. This represents a coplanar resonator (blue) made nonlinear using an embedded Josephson junction (dark green) and capacitively coupled to a transmon qubit [37] placed between the central conductor and the ground planes. The model described in this paper however applies to other nonlinear resonators and qubits (see text).

focus on the qubit and its coupling to the resonator.

A. Nonlinear resonator

Following the notation of Yurke and Buks [START_REF] Yurke | [END_REF] and also used in Paper I, we define the Hamiltonian of the nonlinear resonator as ( = 1)

H r = ω r a † a + K 2 a † a † aa + K ′ 3 a † 3 a 3 , (2.1)
where ω r , K and K ′ are respectively the resonator's bare resonance frequency, its Kerr coefficient and a higher order nonlinearity Kerr coefficient. The operator a ( †) annihilates (creates) an excitation in the resonator. An important aspect of the experiment described in Ref. [20] is the presence of multiple drives on the resonator. We will denote with the subscript d any drive, of amplitude ǫ d and frequency ω d , far detuned from the qubit transition frequency. We will allow for many such drives in our description. In addition, we will denote with the subscript s a spectroscopy drive of amplitude ǫ s and frequency ω s close to the qubit's transition frequency. The presence of these drives can be represented by the usual Hamiltonians [START_REF] Gardiner | Quantum Noise : A Handbook of Markovian and Non-Markovian Qu 3rd[END_REF],

H d = d ǫ d e -iω d t a † + ǫ * d e iω d t a , (2.2a
)

H s = ǫ s e -iωst a † + ǫ * s e iωst a.
(2.2b)

The drives d, far from the qubit resonance, are used to populate the resonator and will not drive transitions of the qubit. They will result in dispersive shifts of the qubit frequency. On the other hand, the spectroscopy drive is aimed specifically at driving the qubit. Because of their different influence, we treat these various drives very differently below. We finally introduce photon loss in the resonator at the rate κ. Together with the above Hamiltonian, this is In the region H (above the full line), the resonator is in a highamplitude state. In the region L, (below dashed line), it is in a low-amplitude state. Between the dashed and full lines, the resonator is bistable. The dashed vertical line corresponds to the detuning at which the data presented in this paper were taken.

captured by the Lindblad form master equation

ρr ≡ L r ρ r = -i [H r + H d + H s , ρ r ] + κD[a]ρ r , (2.3)
where we have introduced the usual dissipation superoperator

D[A]ρ ≡ 1 2 (2AρA † -A † Aρ -ρA † A). (2.4)
Nonlinear resonators in circuit QED typically have a negative Kerr constant K [START_REF] Bourassa | [END_REF]. As one drives the resonator, the nonlinearity therefore causes a back-bending of the resonator's response as illustrated in Fig. 2 (a). For large drive amplitudes ǫ d , the resonator becomes bistable, with the two stable solutions denoted L and H, and of respectively low and high amplitude of oscillation. To simplify the description, it is useful to introduce the reduced detuning Ω = 2(ω r -ω d )/κ [21]. As a function of this reduced detuning and of the drive amplitude ǫ d , the resonator's response is captured by the stablity diagram illustrated in Fig. 2 (b). If Ω is smaller than the critical detuning Ω C = √ 3 the resonator can be in the L state at low power and in the H state at high power, with an intermediate region where both states are stable. Since both states go from stable to unstable at different powers, the resonator in this region is hysteretic, and can be used as a sample-and-hold detector [START_REF] Lupaşcu | [END_REF]19,[START_REF] Siddiqi | [END_REF].

B. Qubit and qubit-resonator coupling

Most superconducting qubits require a larger Hilbert space than their logical subspace {|0 , |1 } for an accurate description. This is the case for example for the transmon [37], the capacitively-shunted flux qubit [40], the fluxonium [41], the tunable-coupling qubit [42] and the phase qubits [43]. Here, we consider up to M levels and write the free qubit Hamiltonian as

H q = M-1 i=0 ω i Π i,i ≡ Π ω , (2.5)
where ω i is the frequency of the qubit eigenstate |i , Π i,j ≡ |i j|, and where we have introduced the shorthanded notation

Π x ≡ M-1 i=0 x i Π i,i . (2.6)
In practice, the values of the eigen-frequencies ω i should either be extracted from experiments, or computed by diagonalizing the full qubit Hamiltonian. For the results presented in this paper, we considered a transmon qubit and carefully calibrated the qubit frequencies from experimental spectroscopic data. More details about the sample and its parameters can be found in Refs. [20,21] and in Paper I.

As usual, we consider a dipolar qubit-resonator coupling

H I = M-2 i=0 g i (a † + a)(Π i,i+1 + Π i+1,i ), (2.7)
where each qubit transition i ↔ j is coupled to the resonator if and only if i = j ± 1. This restrictive condition is often made true either by selection rules or because the other transitions are too far detuned from the resonator frequency to have an impact [37]. We stress that we use here the full Rabi Hamiltonian rather than its Jaynes-Cummings counter-part since, as will be seen below, the counter-rotating terms in Eq. (2.7) will play a predominant role in the sideband transitions.

Taking into account qubit damping and pure dephasing, we finally write the master equation describing the coupled system as

ρ = -i [H, ρ] + κD[a]ρ + γ M-2 i=0 g i g 0 2 D [Π i,i+1 ] ρ + 2γ ϕ D [Π ε ] ρ, (2.8)
where the total Hamiltonian is

H = H r + H q + H I + H d + H s .
In this master equation, γ is the qubit |1 → |0 decay rate and γ ϕ is the qubit pure dephasing rate for the same states. With the above form, we have assumed that qubit decay between two consecutive states scales as the coupling of this transition to the resonator. This assumption is not essential to this work but is convenient and realistic. Moreover, to describe pure dephasing of this multi-level system, we have defined Π ε with

ε i ≡ ∂(ωi-ω0) ∂X × ∂(ω1-ω0) ∂X -1
the X-dispersion and with ε 0 = 0, ε 1 = 1 by definition.. Here, X represents some control parameter (for example flux or charge) whose fluctuations cause dephasing. The above master equation is the same as the one used in Paper I, with the exception that we have set here the resonator's rate of two-photon loss [START_REF] Yurke | [END_REF] to zero for simplicity [44].

Our goal in the next two sections is to obtain a reduced qubit model that includes squeezing of the resonator field and captures qubit-resonator sideband transitions in the high-power regime of the spectroscopy drive.

In section III, we first summarize the results obtained in Paper I. In section IV, we then build on these results and consider first order corrections to two main approximations that were used in Paper I and mentioned in the Introduction. We will show that in the presence of strong spectroscopy drive ǫ s , these corrections will yield a qubit spectrum displaying red and blue sidebands in addition to the main qubit line, and that the amplitude of these sidebands reveals information about squeezing of the intra-resonator field.

III. SUMMARY OF PREVIOUS RESULTS

In Paper I, we have performed a series of unitary transformations on the master equation (2.8) and have obtained an effective master equation for the qubit only. The first step is to transform the master equation using a polaron transformation [45]

P = M-1 i=0 Π i,i D(α i ), (3.1)
where D(α) = exp[αa † -α * a] is the displacement operator [START_REF] Gardiner | Quantum Noise : A Handbook of Markovian and Non-Markovian Qu 3rd[END_REF]. This transformation displaces the resonator field in a qubit-state-dependent manner, such that a → a+Π α . If the pointer states α i are chosen properly, the intraresonator field in this transformed frame is in -or close to -the vacuum. In this situation, it is simple to trace over the resonator states to obtain an effective equation for the qubit only. This can be done exactly within the linear dispersive approximation [45], but unfortunately not when taking into account the full Jaynes-Cummings coupling such as in Paper I. In this situation, the additional complexity arises from transforming the operator Π i,i+1 in Eq. (2.7) which yields

P † Π i,i+1 P = Π i,i+1 D(α i+1 -α i )e -iIm[α * i+1 αi] . (3.2)
This transformed operator is problematic since it contains all powers of the ladder operators a ( †) throught the displacement operator and these will not leave the resonator field in its vacuum state in the transformed frame.

To simplify the situation we assumed in Paper I that the distinguishability |α i+1 -α i | is very small and took

P † Π i,i+1 P ≈ Π i,i+1
D C = exp M-2 i=0 ξ * i Π i,i+1 -ξ i Π i+1,i , (3.3)
where ξ i is a classical (scalar) analog of the operator

λ i a † = [g i /(ω i+1,i -ω r )
]a † found in the usual dispersive transformation of the Jaynes-Cummings Hamiltonian [46]. After these two steps, the result is a transformed master equation containing the ac-Stark shift of the qubit frequency, dressed-dephasing of the qubit [47] and measurement-induced dephasing [48]. These various quantities are related by the nonlinear equations for the pointer states

0 = (ω r -ω d -i κ 2 )α i,d + K|α i | 2 α i,d + K ′ |α i | 4 α i,d + ǫ d + S d i + 1 3! K d i |α i | 2 α i,d , (3.4)
where the expressions for S d i and K d i are given below. With this formulation of Eq. (3.4) we have assumed that α can be written as α = i α i = d,i α i,d e -iω d t + α i,s e -iωst . This form assumes that the multiple drives are spread out enough in frequency such that one drive does not contribute significantly to the field oscillating at another drive's frequency.

In the third step, we apply one last transformation, the quantum dispersive transformation, which takes here the form

D = exp M-2 i=0 λ i a † Π i,i+1 -λ * i aΠ i+1,i , (3.5)
with λ i = g i /(ω i+1,i -ω r ). Since the polaron transformation moves the system to a frame where the photon population is small, this transformation can safely be performed to lowest order. The resulting master equation now contains the Lamb shift of the qubit frequency as well as Purcell decay. After these three transformations and projecting the qubit into its logical subspace {|0 , |1 }, the effective Hamiltonian takes the form H ′′′ = H ′′′ 0 +H ′′′ 2 (each prime indicating a transformation) where

H ′′′ 0 = ω ′′′ 1,0 2 σ z + g 0 α 0,s e -iωst σ + + h.c. , (3.6a)
H ′′′ 2 = ω ′ r (α) + Π S(α) a † a + Υa † 2 + Υ * a 2 . (3.6b)
In these expressions, the ac-Stark and Lamb shifted qubit frequencies are given by

ω ′′′ i (α) ≡ ω ′′ i (α) + L i (α), ω ′′ i (α) ≡ ω i + d S d i |α d | 2 + 1 4 d K d i |α d | 4 (3.7)
where

S d i ≡ -(X d i -X d i-1 ), K d i ≡ -4S d i (|Λ d i | 2 + |Λ d i-1 | 2 ) -(3X d i+1 |Λ d i | 2 -X d i |Λ d i+1 | 2 ) + (3X d i-2 |Λ d i-1 | 2 -X d i-1 |Λ d i-2 | 2 ), (3.8)
are the classical Stark shift coefficients, with

Λ d i ≡ -g i /(ω i+1,i -ω d ) and X d i ≡ -g i Λ d i ,
and where

L i (α) ≡ χ i-1 (α), (3.9a) S i (α) ≡ -(χ i (α) -χ i-1 (α)), (3.9b)
are the Lamb shift and the cavity-pull, with χ i (α) ≡ -g i λ i (α) and

λ i (α) ≡ -g i /[ω ′′ i+1,i (α) -ω ′ r (α)].
We have also defined

Υ ≡ K 2 + K ′ |Π α | 2 Π 2 α . (3.10)
Because of its a ( †)2 dependence, the term proportional to Υ in the transformed Hamiltonian leads to squeezing of the resonator field. This contribution was assumed to be small in Paper I and dropped. Putting all of this together, the resulting transformed qubit-resonator master equation is then given by ρ′′′

= -i [H ′′′ 0 + H ′′′ 2 , ρ ′′′ ] + κD[a]ρ ′′′ + γ ′′′ ↓ D[σ -]ρ ′′′ + γ ′′′ ϕ 2 D[σ z ]ρ ′′′ , (3.11)
where

γ ′′′ ↓ = γ + λ 2 0 (α)κ, (3.12a)
γ ′′′ ϕ = γ ϕ + Γ ϕm , (3.12b)
Γ ϕm = κ|α 1 -α 0 | 2 2 , (3.12c)
are the modified rates having neglected dressed dephasing as well as two-photon losses, both of which were included in Paper I. Using this master equation, we have showed in Paper I that the measurement-induced dephasing rate given by Eq. (3.12c), with the pointer states given by Eq. (3.4), is in quantitative agreement with measured qubit spectroscopic linewidth within the limits of the approximations that were made. There are few key points behind this good agreement between theory and experiments. First, contrary to what is usually used in circuit QED, our expression for the ac-Stark shift depends on the pump drive frequency rather than the resonator frequency. Our treatment moreover goes beyond linear response theory for the resonator state. Indeed, in Paper I we show that for typical circuit QED parameters, describing the measurementinduced dephasing quantitatively with a nonlinear resonator always requires going beyond a linear response. That is, whenever the gain of the resonator is large, the susceptibility of the resonator response to a shift in its resonance frequency is large. Because of this large susceptibility, the qubit cannot be treated as a simple perturbation causing a small shift of the resonator frequency.

IV. SQUEEZING AND SIDEBANDS

In this section, we relax the two main approximations that are discussed above (small distinguishability and negligible squeezing) and include first order corrections. First, instead of approximating Eq. (3.2) as P † Π i,i+1 P ≈ Π i,i+1 , we now take

P † Π i,i+1 P ≈ Π i,i+1 1 + β i a † -β * i a , (4.1)
where β i ≡ α i+1 -α i is assumed to be small. Taking these terms into account when transforming H I yields a term in the Hamiltonian that was neglected in Paper I and that is given by

H ′ SB ≡ M-2 i=0 g i (Π * α + Π α ) β i a † -β * i a Π i,i+1 + h.c. . (4.
2) This Hamiltonian generates multi-photon qubitresonator sideband transitions that will appear in the qubit spectrum.

With the proper choice of polaron frame (i.e. of pointer states α i ), the average field a is small. Assuming that β i is also small, we consider H ′ SB to be itself a correction to the transformed system Hamiltonian. We therefore omit to apply the classical and quantum dispersive transformations on H ′ SB since this would only yield even smaller corrections.

Taking into account this correction, the master equation (3.11) describing the system in the three-times transformed frame now reads

ρ′′′ = -i [H ′′′ 0 + H ′′′ 2 + H ′′′ SB , ρ ′′′ ] + κD[a]ρ ′′′ + γ ′′′ ↓ D[σ -]ρ ′′′ + γ ′′′ ϕ 2 D[σ z ]ρ ′′′ , (4.3)
where H ′′′ SB = H ′ SB .

A. Bogoliubov transformation

In Paper I, we had only H ′′′ 0 + H ′′′ 2 as the Hamiltonian in the transformed frame and assumed that H ′′′ 2 was a perturbation small enough that the state of the resonator was the vacuum and could be readily traced out. In practice however the terms proportional to a ( †)2 in H ′′′ 2 lead to squeezing of the stationary state of the resonator. Here, we do not drop these terms and take care of them using a Bogoliubov transformation before adiabatically eliminating the resonator.

Indeed, assuming Hamiltonian H ′′′ 2 of Eq. (3.6b) depends only weakly on the qubit state, it can be diagonalized using a Bogoliubov transformation, which takes the form

S = e 1 2 ξ * a 2 - 1 2 ξa † 2 (4.4)
and whose action on the field operator a is [START_REF] Gardiner | Quantum Noise : A Handbook of Markovian and Non-Markovian Qu 3rd[END_REF] S † aS = cosh(r)a -e i2θ sinh(r)a † , (

where ξ = re 2iθ is the squeezing parameter.

To keep the analytical calculations tractable, we will keep in this transformation only the time-dependance of θ with θ(t) ≡ -ω p t + Θ. This implies that the field has reached a steady state in a frame rotating at ω p . We show in Appendix A that under transformation by S the qubit-resonator master equation now takes the form ρ(4) = -i H (4) s , ρ (4)

+ γ ′′′ ↓ D [σ -] ρ (4) + γ ′′′ ϕ 2 D[σ z ]ρ (4) + κ sinh 2 (r) + 1 D[a]ρ (4) + κ sinh 2 (r)D[a † ]ρ (4) , (4.6)
where

H (4) s = H (4) 0 + H (4) 2 + H (4)
SB , with

H (4) 0 = g 0 α s,0 σ + e iδt + h.c. , H (4) 2 = ∆r (α) a † a + 1 2 , H (4)
SB = F (4) σ + e iδt + h.c., (4.7) and

F (4) = g 0 α s,0 (ca † -c * a), c ≡ β cosh(r) + β * e i2θ sinh(r), ∆r (α) ≡ [ω ′ r (α) + Π S(α) -ω p ]/ cosh(2r), δ = ω ′′′ 1,0 -ω s . (4.8)
In obtaining this master equation, we have moved to a frame rotating at ω ′′′ 1,0 for the qubit and at ω p for the resonator. We have also neglected rapidly oscillating terms in H (4) SB using the rotating wave approximation. Finally and as presented in more details in Appendix A, we have assumed that the photon population in the transformed frame is small. As can be seen from the term sinh 2 (r)D[a † ]ρ responsible for heating in Eq. (4.6), this assumption will only be true in the limit of small squeezing. The squeezing coefficient r is given by the solution of the equations cos

[arg(Υ p ) -2θ] = ω ′ r (α) + Π S(α) -ω p 2|Υ p | tanh 2r, (4.9a) sin [arg(Υ p ) -2θ] = κ 4|Υ p | sinh(2r). (4.9b)
While complicated and not reproduced here, the solution of these equations is analytical, and yields a maximum squeezing coefficient

r max = 1 2 arcsinh 4|Υ p | κ . (4.10)
Maximal squeezing is reached for ω p = ω ′ r (α)+Π S(α) . For the parameters of interest here and in Ref. [20], r max ∼ 0.5 for all qubit states corresponding to sinh 2 r max ∼ 0.3. Although we would have prefered to have sinh 2 r ≪ 1 to justify our approximation, we will see below that this model nevertheless semi-quantitatively compares with experimental results.

B. Adiabatic elimination of the resonator

We now adiabatically eliminate the resonator to obtain a master equation for the qubit only. As described in Appendix B, we use the projector formalism [START_REF] Gardiner | Quantum Noise : A Handbook of Markovian and Non-Markovian Qu 3rd[END_REF]49] to obtain the following reduced master equation:

ρq = -i H, ρ q + γ↓ D[σ -]ρ q + γ↑ D[σ + ]ρ q + γϕ 2 D[σ z ]ρ q .
(4.11) In this expression, we have defined the rates

γ↓ = γ ′′′ ↓ + |g 0 α s,0 c| 2 [L(-δ) + L(δ)] sinh 2 r + L(-δ) , (4.12a) γ↑ = |g 0 α s,0 c| 2 [L(-δ) + L(δ)] sinh 2 r + L(δ) , (4.12b) γϕ = γ ′′′ ϕ , (4.12c)
and the Hamiltonian

H = δ 2 σ z + g 0 (α s,0 σ + + α * s,0 σ -), (4.13)
where

L(δ) = Re[f (δ)] is a Lorentzian with f (ω) = κ/2 + i[ ∆r (α) + ω] κ 2 /4 + [ ∆r (α) + ω] 2 . (4.14) Moreover, δ = ω ′′′ 1,0 (α) -ω s + Im[S ↓ (δ) -S ↑ (-δ)]
where S ↓ (δ) and S ↑ (-δ) are defined in Appendix B. The above master equation is one of the central results of this paper and the significance of the various terms is discussed below.

C. Steady state qubit population

We are now almost ready to compare the above model to the experiments presented in Ref. [20]. It is however useful to give some more details of the experiment so as to compute the appropriate qubit observable for comparison with theory. In the experiment presented in Ref. [20], a transmon qubit [37] is coupled to a nonlinear transmission resonator. The nonlinearity is provided by a Josephson junction embedded in the central conductor of the resonator. The resonator is first driven with a pump drive of frequency ω p and power P p pushing the resonator field out of its ground state. The drive is chosen to be long enough for the coupled resonator-qubit system to reach its steady state. Under this steady state resonator field, the qubit is ac-Stark shifted as described in Eq. (3.7). This frequency shift is then measured in order to reveal information about the internal resonator field.

For this, a second (spectroscopy) drive at frequency ω s close to the qubit frequency is turned on while the pump field is still present. This spectroscopy drive excites the qubit only if ω s matches the shifted qubit transition frequency. Both drives are then turned off for a time longer than the resonator decay time, but shorter than the qubit's relaxation time. The qubit state is then measured using the standard bifurcation readout procedure [19]. This process is then repeated multiple times for different spectroscopy frequency ω s and pump power P p to extract the probability P (|1 ) of the qubit to be in its excited state.

To compare our model to experimental results, we therefore compute from Eq. (4.11) the steady state probability of the qubit to be in its excited state. We find: The above expression for P (|1 ) can be understood by focussing on three different contributions, leading to three peaks in P (|1 ) versus pump power P p and spectroscopy frequency ω s . The first peak is obtained for γ↓ ≫ γ↑ such that Π 1,1 eq ∼ 0. In this regime, P (|1 ) reduces to

P (|1 ) = Π 1,1 eq γ2 2 + δ2 + 2γ 2 |g 0 α s,0 | 2 /(γ ↑ + γ↓ ) γ2 2 + 4γ 2 |g 0 α s,0 | 2 /(γ ↑ + γ↓ ) + δ2 , ( 4
P (|1 ) ≈ 2γ 2 |g 0 α s,0 | 2 /γ ↓ γ2 2 + 4γ 2 |g 0 α s,0 | 2 /γ ↓ + δ2 . (4.17)
That is, we find a Lorentzian centered at δ = 0 with width γ2 in the absence of power broadening. This Lorentzian is the main qubit line and is power-broadened by the spectroscopy field α s,0 .

The two other contributions are found when γ↑ is large. From the definition of γ↑ in Eq. (4.12b), this requires that ± δ ∼ ∆r (α) which corresponds to L(±δ) taking its maximal value. If |δ| is sufficiently large and if the undriven decay rate γ ′′′ ↓ is negligible compared to the effective rates arising from squeezing [∝ |g 0 α s,0 c| 2 /κ in Eqs. (4.12a) and (4.12b)], then the dominant contribution is Π 1,1 eq

P (|1 ) ≈ Π 1,1 eq ≈ L(-δ) sinh 2 r + L(δ)(sinh 2 r + 1) [L(-δ) + L(δ)](2 sinh 2 r + 1
) .

(4.18) In this situation, corresponding to the resolved sideband limit, P (|1 ) takes the form of red and blue sidebands on either side of the main qubit line. Interestingly these two sidebands depend on the squeezing parameter r. In the next section, we compare this three-peak excitation spectra to experimental data and analyze the amplitude, width and position of the peaks to extract information about the internal resonator state and in particular about the squeezing parameter.

V. COMPARISON TO EXPERIMENTS

As discussed in Paper I and in Ref. [21], only the centre qubit line is observed under low spectroscopy power P s . This can be understood from the expressions (4.12a) and (4.12b) for γ↑ and γ↓ . Indeed, for negligible spectroscopy power |α s,0 | 2 → 0 and γ↓ ≫ γ↑ such that as discussed above only the centre line is apparent. The position and width of this peak were analyzed in details in Paper I and in Ref. [21]. Here we focus on the situation where the spectroscopy power P s is important such that γ↑ cannot be neglected with respect to γ↓ . The experimental results in this situation for P (|1 ) versus pump power and spectroscopy frequency -first presented in Ref. [20] -are reproduced in Fig. 3a). On the right side of Fig. 3, we show the result of P (|1 ) as given by the fully analytical expression of Eq. (4.15). In both cases, we see a sudden jump in the qubit transition frequency. This corresponds to switching of the resonator state from its L to its H state and was discussed in details in Ref. [21]. At larger pump power, sidebands are clearly resolved and the agreement between experiments and theory is excellent. In producing this figure, all parameters except one have been extracted independently and are given in Ref. [20] and Paper I. The adjustable parameter is an ad-hoc multiplicative coefficient to the effective sideband driving amplitude F (4) defined in Eq. (4.8). This parameter accounts for the large dependence of the coefficient c in F (4) on the squeezing coefficient ξ = re 2iθ . For our model to reproduce the experimental sideband amplitudes, we have multiplied c by 2. Considering the number of approximations that are done in obtaining the model, we consider such a factor to be an acceptable correction. It is important to stress that this correction does not affect the position (frequency) of any of the linesmadeMoreover, while it changes the absolute amplitude of the sidebands relative to the main line, it does not change their width nor does it change the ratio of their amplitude since both sidebands are affected in the same way by this correction.

For a more quantitative comparison, we show in Fig. 4 four linecuts of Fig. 3 (a) and (b), for increasing pump power P p . In these plots, the dots correspond to the experimental data, while the lines are the analytical predictions. With the above single correction, the agreement is almost quantitative at all powers. We note however that at P p = -5 dBm (top left panel) the analytical results predict two peaks while a single one is observed experimentally [START_REF]This sideband was observed experimentally with another sample[END_REF]. In fact, in all cases the sidebands are better resolved in the analytical model than is observed Sideband splitting fifc (a), width wi (b) and ratio of amplitude A b /Ar (c) as a function of the pump amplitude ǫp, in logarithmic scale. In (a) and (b), red circles (full red lines), blue crosses (dotted blue lines) and black squares (dashed black lines) are experimental data (theoretical predictions) for the red sideband (i = r), blue sideband (i = b) and main line (i = c), respectively. In (c), black squares (orange double-dotted line) are extracted from fits to the experimental (analytical) spectrum, while the dotted grey line and dashed-dotted green line correspond to Eq. (5.1a) and Eq. (5.1b), respectively. experimentally, especially close to the bifurcation threshold. We attribute this discrepancy to the rotating wave approximation made in Appendix B where we neglected terms oscillating at the sideband-detuning frequency. As is discussed there and observed in Fig. 4, making this approximation corresponds to considering the resolvedsideband limit.

While the analytical expression for P (|1 ) in terms of ω s is not simply that of three superposed Lorentzians, it is useful to compare fits of both the analytical expressions and the experimental data to such a simplified model. The position, amplitude and widths of the three peaks were therefore extracted by fitting the sum of three Lorentzian curves both to the experimental and analytical spectra. We plot the results of these fits in Fig. 5, where we show the sideband detuning from the main line [panel (a)], the width of the three peaks [panel (b)], and the ratio of amplitude of the blue to the red sidebands [panel (c)]. Due to the breakdown of the resolvedsideband approximation and as mentioned above, we see in (a) that the sidebands are slightly more separated in the analytical model (lines) than in the experimental data (dots). We see from panel (b) that the width of the red sideband (full line and circles) and the centre peak (dashed line and squares) are in quantitative agreement. The model however shows a blue sideband that is much wider than obtained experimentally. We attribute this discrepancy to the small signal-to-noise ratio of the blue sideband (i.e. it is not very distinguishable from the noise). In fact, the statistical error on the fit parameters of the blue sidebands are about 30%.

In panel (c), we show the ratio of the amplitude of the blue sideband to that of the red sideband. In addition to the data from the fit to the experimental data (black squares) and to the analytical expression (orange double-dotted line), we have extracted this ratio directly from the analytical expression for P (|1 ), i.e. without assuming a Lorentzian profile. Since the blue (red) sideband is at a frequency corresponding to δ = ∆r (δ = -∆r ), the ratio of amplitudes is given by A b /A r = P (|1 )| δ=-∆r /P (|1 )| δ= ∆r . Using the simplified expression Eq. (4.18) for the amplitude of the sidebands, we find

A b A r ≈ [L( ∆r ) + L(-∆r )] sinh 2 r + L(-∆r ) [L(-∆r ) + L( ∆r )] sinh 2 r + L( ∆r ) (5.1a) ≈ sinh 2 r sinh 2 r + 1 , (5.1b)
where in the last approximation, we have assumed that L(-∆r ) ≫ L( ∆r ), which is valid in the limit ∆r ≫ κ, i.e. in the resolved-sideband limit. Eq. (5.1a) above corresponds to the dotted grey line in panel (c) and Eq. (5.1b) to the dashed-dotted green line. Because of space constrains, only the simpler expression Eq. (5.1b) was used for comparison to the experimental data in Ref. [20].

We emphasize once more that the correction factor applied to the theoretical results does not affect any of the results showed in Fig. 5. As a result, given the quantitative agreement displayed in panel c, it is possible to accurately determine the squeezing coefficient r of the intra-resonator field. As already mentioned above, we find a maximal value of sinh 2 r max ∼ 0.3 corresponding to r max ∼ 0.5.

Alternatively, the ratio of the two sidebands also allows us to extract the effective temperature of the oscillator as described in the quadruply transformed frame by the master equation (4.6). Indeed, the last line of this master equation takes the standard form

κ(n th + 1)D[a]ρ (4) + κn th D[a † ]ρ (4) , (5.2)
where we have identified sinh 2 (r) with an effective thermal occupation number of the oscillator, n th = sinh 2 (r). Using the Bose-Einstein distribution for n th then allows us to define an effective temperature T eff for the system due to squeezing of the resonator field. This corresponds to the essential prediction of the quantum heating theory [23][24][25][26]. As discussed in more details in Ref. [20], T eff corresponds to a few tenths of quanta, much larger than the temperature expected from the base temperature of the dilution refrigerator and filtering of the lines. Moreover, the non-monotonic dependence with pump power, especially the decrease of the effective temperature with increasing pump power above bifurcation exclude a classical heating effect (due to the contacts for example) from being the cause of this effective temperature.

VI. CONCLUSION

We have developed a theoretical model for a multilevel qubit coupled to a pumped nonlinear resonator. The model holds within the dispersive regime of circuit QED and for pumping powers well above the bifurcation threshold of the resonator. The reduced qubit master equation that we have obtained contains Purcell decay, measurement-induced dephasing, dressed dephasing, quadratic ac-Stark shift, Lamb shift, as well as the first order correction of the resonator's squeezing on the qubit. This allows us to obtain quantitative agreement with experimental data in a wide range of parameters, without adjustable parameters. In this way, we show how the qubit can be used as a probe of squeezing of the intraresonator field. By comparing the ratio of the red and blue sidebands in the qubit excitation spectrum, we have extracted the squeezing coefficient of the field or equivalently the effective temperature of the so-called quasioscillator, providing a direct demonstration of quantum heating.

Interesting extensions to this model include a finer treatment of the intra-resonator squeezing, especially the squeezing angle, as well as a qubit-dependent squeezing. This latter aspect could provide insights on how to improve dispersive qubit measurement using squeezing, or on how measurement-induced dephasing is changed by squeezing of the intra-resonator field.

Note added.-After this work was completed, we became aware of related work [START_REF] Peano | [END_REF].

No. EQuS CE110001013). We thank Calcul Québec and Compute Canada for computational resources. and the complementary projector R ≡ I -Q, where I is the identity. There are a number of useful identities that one can prove with these projectors. Among these are L r Q = QL r = 0, L q Q = QL q , RL r = L r R and RL q = L q R. The first property arises from the adiabatic approximation and the trace-preserving nature of Lindbladians, while the three others are consequences of the definition of Q and R [START_REF] Gardiner | Quantum Noise : A Handbook of Markovian and Non-Markovian Qu 3rd[END_REF].

Using the above properties, and assuming -without loss of generality -that QL c (t)Q = 0, we then search for the evolution equations of v(t) = Qρ(t) and w(t) = Rρ(t). We obtain v(t) = QL c (t)w(t) + L q v(t), ẇ(t) = (L r + RL c (t) + L q )w(t) + RL c (t)v(t).

(B5)

Assuming that L r contains the dominant dynamic and that the integration time is long compared to 1/κ yet short compared to 1/γ, 1/γ ϕ , we obtain the approximate solution for w(t)

w(t) ≈ ∞ 0 d t exp [L r t] RL c (t -t)v(t). ( B6
)
We can then replace this solution into the equation for v(t) and, using the definition of Q, obtain an equation for ρq ρq = L q ρ q (t) + ∞ 0 dt ′ Tr r L c (t)e Lrt ′ L c (t -t ′ )ρ(t) .

(B7) Given the form of L c (t), one can then show that the qubit reduced master equation can be written as ρq = L q ρ q -i [δH, ρ q ] + Re[S ↑ (-δ)]D[σ + ]ρ q + Re[S ↓ (δ)]D[σ -]ρ q , (B8)

where we have defined the spectra

S ↑ (ω) = ∞ 0
dt ′ e iωt ′ F (4) † (t ′ )F (4) (0) Importantly, in obtaining the above master equation, we have assumed that we could neglect terms oscillating at a frequency 2δ. As discussed in section V, this approximation corresponds to the well resolved sideband limit. Finally, we use the quantum regression theorem [START_REF] Gardiner | Quantum Noise : A Handbook of Markovian and Non-Markovian Qu 3rd[END_REF] to obtain S ↑ (ω) = |g 0 α s,0 c| 2 f (ω) sinh 2 r + f * (-ω)(1 + sinh 2 r) , (B11) where we have defined the complex function

f (ω) = κ/2 + i[ ∆r (α) + ω] κ 2 /4 + [ ∆r (α) + ω] 2 , (B12)
and S ↓ (ω) = S ↑ (ω). With this spectrum, we easily get the reduced master equation (4.11) for the qubit only.

FIG. 2 .

 2 FIG. 2. (Color online) (a) Amplitude of the resonator internal field (arbitrary units) in response to a drive of reduced frequency Ω = 2(ωrω d )/κ (ΩC = √ 3) for increasing drive amplitudes ǫ d . (b) Stability diagram of the resonator.In the region H (above the full line), the resonator is in a highamplitude state. In the region L, (below dashed line), it is in a low-amplitude state. Between the dashed and full lines, the resonator is bistable. The dashed vertical line corresponds to the detuning at which the data presented in this paper were taken.

FIG. 3 .

 3 FIG. 3.Qubit spectrum for a strong spectroscopy amplitude ǫs as a function of the pump amplitude ǫp, in logarithmic scale. Experimental results (left) are compared to the analytical spectrum (right). Results are for pump frequency ωp/2π = 6439 MHz, corresponding to Ω/ΩC = 1.74 and spectroscopy power Ps = -4 dBm, corresponding to ǫs/2π ≈ 25 MHz. All other system parameters can be found in Ref.[20].

FIG. 4 .

 4 FIG. 4. Line cuts of Fig. 3 for four pump powers. Black lines are theoretical predictions. Dots are experimental data.

  FIG. 5.Sideband splitting fifc (a), width wi (b) and ratio of amplitude A b /Ar (c) as a function of the pump amplitude ǫp, in logarithmic scale. In (a) and (b), red circles (full red lines), blue crosses (dotted blue lines) and black squares (dashed black lines) are experimental data (theoretical predictions) for the red sideband (i = r), blue sideband (i = b) and main line (i = c), respectively. In (c), black squares (orange double-dotted line) are extracted from fits to the experimental (analytical) spectrum, while the dotted grey line and dashed-dotted green line correspond to Eq. (5.1a) and Eq. (5.1b), respectively.

  iωt ′ F (4) (t ′ )F (4) † (0) with • s = Tr r (•ρ s r) and whereδH = Im[S ↑ (-δ)]σ -σ + + Im[S ↓ (δ)]σ + σ -.(B10)

  . With this approximation, the interaction Hamiltonian H I transforms into a detuned drive acting directly on the qubit.The second step in Paper I is to remove this effective detuned qubit driving using what we called a classical dispersive transformation

ACKNOWLEDGMENTS

We acknowledge discussions with J. Gambetta, M. Dykman and within the Quantronics group. We acknowledge support from NSERC, FQRNT, the Alfred P. Sloan Foundation, CIFAR, the European project Scaleqit and the ITN network CCQED and the Australian Research Council Centre of Excellence Scheme (Grant

Appendix A: Squeezon transformation

In this Appendix, we diagonalise the sideband Hamiltonian by applying the transformation S of Eq. (4.4) on the master equation (4.3). Given the action of S on a [see Eq. (4.5)], we find for the transformed dissipator

As we can see in the first line, what was pure damping in the original frame now sees heating in the transformed frame. The second line of the above expression takes the form of a commutator and can be added to the Hamiltonian part of the master equation. Finally, we will neglect the last line because it vanishes if the resonator is in its ground state. Because of the presence of the heating term, this is an approximation that restricts this theory to low squeezing (i.e. low effective temperature).

Taking into account the above contribution from the dissipation, we can then transform the system Hamiltonian

where F (4) is

) We see that choosing r and θ such that

yields vanishing squeezing terms, leaving only a renormalized harmonic oscillator and a driven qubit. The solution to this equation can be expressed in the simpler form of Eq. (4.9) when considering the real and imaginary parts separately.

Assuming that the squeezing coefficient does not depend on the qubit state, and moving to a frame rotating at ω p for the resonator and ω ′′′ 1,0 for the qubit, we then obtain the master equation (4.6).

Appendix B: Adiabatic elimination through the projector formalism

In this Appendix we use the projector formalism [START_REF] Gardiner | Quantum Noise : A Handbook of Markovian and Non-Markovian Qu 3rd[END_REF]49] to adiabatically eliminate the resonator's degrees of freedom. We first note that the master equation (4.6) can be expressed as

where to simplify the notation used in this Appendix we have dropped the index on ρ and defined the resonator, qubit and coupling Lindbladians as

respectively. We then assume that the resonator relaxes much faster than the qubit (i.e. κ ≫ γ, γ ϕ ). This allows us to assume a stationary state for the resonator given by ρs r = L r ρ s r = 0. In the same way, we also assume that the qubit-resonator state can be written as ρ(t) = ρ s r ⊗ ρ q (t), (B3)

where ρ q (t) is the reduced qubit density matrix. In order to find the evolution equation for ρ q (t), we define the projector on the qubit subspace Q: Qρ ≡ ρ s r ⊗ Tr r (ρ) , (B4)


Sébastian Minjeaud

Richard Pasquetti
email: richard.pasquetti@unice.fr

High order approximation of a tokamak edge plasma transport minimal model with Bohm boundary conditions

Keywords: Magnetic confined fusion, Edge plasma, Bohm boundary conditions, Hyperbolic system, Tokamaks

In the recent paper [START_REF] Isoardi | Penalization modeling of a limiter in the tokamak edge plasma[END_REF], a penalty method is proposed to address the so-called Bohm boundary conditions which are generally imposed to model the limiters, i.e. the walls that intercept in a tokamak the magnetic field lines. This approach is validated by considering a simple one-dimensional hyperbolic system, that constitutes a minimal transport model for ionic density and momentum. In [START_REF] Ph | An optimal penalty method for a hyperbolic system modeling the edge plasma transport in a tokamak[END_REF] this penalty model is discussed and an improved formulation of the penalty method is proposed. More recently, in [START_REF] Paredes | A penalization technique to model plasma facing components in a tokamak with temperature variations[END_REF], a similar model is considered but completed with evolution equations for the ionic and electronic temperatures. Again, a penalty method is proposed to enforce the boundary conditions at the walls.

In the present note, we consider the system studied in [START_REF] Ph | An optimal penalty method for a hyperbolic system modeling the edge plasma transport in a tokamak[END_REF][START_REF] Isoardi | Penalization modeling of a limiter in the tokamak edge plasma[END_REF] and show that as soon as the solution is smooth, a spectrally accurate approximation can be obtained. To this end, a stabilized spectral element method (SEM) is prefered to the Godunov schemes used in [START_REF] Ph | An optimal penalty method for a hyperbolic system modeling the edge plasma transport in a tokamak[END_REF][START_REF] Isoardi | Penalization modeling of a limiter in the tokamak edge plasma[END_REF][START_REF] Paredes | A penalization technique to model plasma facing components in a tokamak with temperature variations[END_REF]. Moreover, instead of the penalty approach, which in the multidimensional context can only offer a first order algebraic convergence, we use a direct imposition of the Bohm boundary condition in an explicit time marching. The system studied in [START_REF] Ph | An optimal penalty method for a hyperbolic system modeling the edge plasma transport in a tokamak[END_REF][START_REF] Isoardi | Penalization modeling of a limiter in the tokamak edge plasma[END_REF] writes:

∂ t N + ∂ x Γ = S N (1)
∂ t Γ + ∂ x Γ 2 N + N = S Γ (2)
with t ∈ R + and x ∈ Ω = (x min , x max ) for the time and space variables and where N and Γ stand for the dimensionless ion density and ion momentum, respectively. In this model, the ion and electron temperatures are assumed to be constant, so that when using the perfect gas law, in dimensionless form the pressure is equal to the density. We then recognize the usual conservation equations of mass and momentum, where some source terms S N and S Γ are considered. Moreover, the dimensionless velocity equals the Mach number, i.e. M = Γ/N . It is easy to check that the system (1)-( 2) is strictly hyperbolic. In convective form it writes:

∂ t N Γ + 0 1 -M 2 + 1 2M ∂ x N Γ = S N S Γ (3)
and so the wave speeds, defined as the eigenvalues of the matrix in (3), equal λ ± = M ± 1. As well known, the number of boundary conditions to be considered depends on the number of ingoing waves. Thus, at x max two boundary conditions are needed if M < -1, one boundary is needed if -1 < M < 1 and no boundary conditions are needed if M > 1.

In order to model the limiter by a particles sink, the Bohm boundary condition states that M (x max ) ≥ 1 [START_REF] Ph | Transition to supersonic flows in the edge plasma[END_REF]. Note that M (x max ) = 1 corresponds to the case where the boundary x = x max is characteristic, so that it is not clear what should be done. Thus, in [START_REF] Ph | An optimal penalty method for a hyperbolic system modeling the edge plasma transport in a tokamak[END_REF] it is suggested to slightly modify the boundary condition by imposing M (x max ) = 1 -η and M (x min ) = 1 + η, with a fixed (small) η > 0 in order to consider a wellposed hyperbolic problem. From our point of view, here this problem is mainly theoretical and probably not relevant as soon as numerical time and space discretizations are involved. Hereafter, we enforce the following Bohm conditions:

M (x min ) ≤ -1 and M (x max ) ≥ 1 . (4)
To check this approach in the frame of high order methods, we use in space the SEM, i.e. a high order nodal finite element method such that the approximation space, say E h , contains all C 0 functions whose restriction in each element is a polynomial of degree p, see e.g. [START_REF] Karniadakis | Spectral hp element methods for CFD[END_REF]. In each element the basis functions are Lagrange polynomials based on the Gauss-Lobatto-Legendre (GLL) points and these GLL interpolation points are also used as quadrature points (quadratures are then exact for polynomials of degree 2p -1). Moreover, since the problem is hyperbolic, shocks may develop and a stabilization technique is required. We use the Spectral Vanishing Viscosity (SVV) technique introduced in [START_REF] Maday | Legendre pseudo-spectral viscosity method for nonlinear conservation laws[END_REF][START_REF] Tadmor | Convergence of spectral methods for nonlinear conservation laws[END_REF] (see also [START_REF] Xu | Stabilized spectral element computations of high Reynolds number incompressible flows[END_REF] in the multidimensional context). With (N h , Γ h ) for the numerical approximations of (N, Γ), we thus consider the following variational formulation:

xmax xmin ∂ t N h v h dx + xmax xmin ∂ x Γ h v h dx + V N = xmax xmin S N h v h dx , ∀v h ∈ E h , (5)
xmax xmin ∂ t Γ h v h dx + xmax xmin ∂ x ( Γ 2 h N h + N h )v h dx + V Γ = xmax xmin S Γh v h dx , ∀v h ∈ E h (6)
where the SVV stabilization terms write:

V N = ǫ hp xmax xmin Q p (∂ x N h )∂ x v h dx , V Γ = ǫ hp xmax xmin Q p (∂ x Γ h )∂ x v h dx . (7)
In these expressions, ǫ hp controls the amplitude of the SVV term and Q p is the so-called SVV operator whose kernel contains all low frequency Legendre components, say for k ≤ m p < p of the spectral approximation, see e.g. [START_REF] Maday | Legendre pseudo-spectral viscosity method for nonlinear conservation laws[END_REF] for details. In the computations, we have systematically used ǫ hp = h/p, where h is the grid-size, and m p = [ √ p] ([.], for nearest integer).

For the integration in time, we use a standard fourth order Runge Kutta scheme (RK4 scheme). Then, to enforce the Bohm conditions we impose at the end of each RK4 step:

N h (x min ) := -Γ h (x min ) if Γ h N h (x min ) > -1 ; N h (x max ) := Γ h (x max ) if Γ h N h (x max ) < 1 . (8)
On the basis of this SEM approximation, we have considered three different test cases.

Test case 1: This test-case is the one introduced in [START_REF] Isoardi | Penalization modeling of a limiter in the tokamak edge plasma[END_REF]. The length of the domain is unitary, Ω = (0, 1). At the initial time t = 0 the fluid is at rest, M 0 = Γ 0 = 0, and the density is constant, N 0 = 1. For the source terms one has S N = 2 and S Γ = 0. In [START_REF] Isoardi | Penalization modeling of a limiter in the tokamak edge plasma[END_REF] the asymptotic solution, i.e. obtained at t ≫ 1, is compared to the analytical one, say (N ∞ , Γ ∞ ), such that:

N ∞ = S N (0.5 + x(1 -x)) , Γ ∞ = S N (x -0.5) . (9)
The computation has been carried out with K = 20 elements and with the polynomial degree p = 6. In Fig. 1, we show the solution (N h , Γ h /N h ) at different times. At the final time of the computation, the numerical solution clearly agrees with the analytical one. Despite the use of the SVV technique, one can discern that the solution is not perfectly smooth. This is mainly associated to the fact that for this test case the initial profiles N 0 and Γ 0 do not verify the Bohm conditions, since M (0) = M (1) = 0 at t = 0.

Test case 2: This test-case is similar to test-case 1, i.e. S N = 2 and S Γ = 0, but to check our algorithm with a supersonic -subsonic transition at the boundary, the initial condition for the momentum is now Γ 0 = 2(x -0.5) -2 sin(πx). Then, at the initial time, the flow is sonic at the boundaries and partially subsonic and supersonic inside the domain. This situation may be met in tokamaks [START_REF] Ph | Transition to supersonic flows in the edge plasma[END_REF]. Using again K = 20 and p = 6, the solution thus obtained is shown at different times in Fig. 2, where we focus on the shortest times, when occurs the transition. Later on, the evolution is similar to the test-case 1. It should be mentionned that:

-for the present computation the SVV stabilization is not mandatory. This was not the case previously, with p = 6, probably because at t = 0, |M (±1)| < 1 for test 1.

-we have also obtained satisfactory results by enforcing the Bohm conditions at the end of each RK4 step, using:

Γ h (x min ) := -N h (x min ) if Γ h N h (x min ) > -1 ; Γ h (x max ) := N h (x max ) if Γ h N h (x max ) < 1 , (10)
which points out that in the 1D case, there is no predominant role of one of the two variables. With p = 6 this was however not true for the test 1. Thus, for robustness reasons we prefer to use [START_REF] Xu | Stabilized spectral element computations of high Reynolds number incompressible flows[END_REF] in the RK scheme.

-Using a more complex initial condition, e.g. (N 0 = 1, Γ 0 = -cos(3πx)), one can observe shocks. Such a non physical situation is well supported by the SVV stabilization, even if as expected one observes the Gibbs phenomenon.

Test case 3: This test-case is essentially the one proposed in [START_REF] Ph | An optimal penalty method for a hyperbolic system modeling the edge plasma transport in a tokamak[END_REF]; For x ∈ Ω = (-0.4, 0.4)

N (x, t) = exp -x 2 0.16(t + 1) , Γ(x, t) = c sin πx 0.8 -x 2 0.16(t + 1) (11)
with c = 1 in [START_REF] Ph | An optimal penalty method for a hyperbolic system modeling the edge plasma transport in a tokamak[END_REF]. The source terms S N and S Γ are those obtained by plugging these expressions in the governing equations. The initial profiles N 0 and Γ 0 are those obtained for t = 0. Clearly, M = Γ/N is constant in time and shows a sinusoidal behaviour in space such that M (±0.4) = ±1 if c = 1. Such a very smooth solution is of interest to carry out a convergence study. We have introduced the parameter c in order to check our SEM approximation in the case c > 1, i.e. for a supersonic outgoing flow. The computations were carried out with K = 20 elements and various polynomial degrees, 1 ≤ p ≤ 8. The time-step is chosen small (∆t = 10 -4 ) so that the time approximation errors are negligible. Fig. 3 shows the variation of the error, in the max norm (i.e. in the l ∞ norm), at t = 1, with respect to the polynomial degree, for the sonic (c = 1) and supersonic flows (c = 2), with and without using the SVV stabilization technique. Even if such convergence curves are not really smooth, one clearly discerns the expected SEM exponential rate of convergence.

Figure 1 :

 1 Figure 1: Density and Mach number at different times and exact asymptotic solution. Computation done with K = 20 elements and a polynomial degree p = 6. SVV stabilization: mp = 2 and ǫ hp = 1/120.

Figure 2 :

 2 Figure 2: Density and Mach number at different times. Computation done with K = 20 elements and a polynomial degree p = 6. SVV stabilization: mp = 2 and ǫ hp = 1/120.

Figure 3 :

 3 Figure 3: Error versus the polynomial degree with and without SVV stabilization, for the sonic (at left) and supersonic (at right) flows. The space discretization is based on K = 20 elements. If SVV is used, then ǫ hp = 0.8/(20p) and mp = [ √ p].

Acknowledgments: This study was supported by the FR-FCM, project "Fluid numerical modelling of edge plasma of tokamaks. Application to ITER" and by the EUROFUSION, project "Synergetic numericalexperimental approach to fundamental aspects of turbulent transport in the tokamak edge". We also thank all partners of the ANR project "ESPOIR" for fruitful discussions.


M Stern

G Catelani

Y Kubo

C Grezes

A Bienfait

D Vion

D Esteve

P Bertet

Flux Qubits with Long Coherence Times for Hybrid Quantum Circuits

We present measurements of superconducting flux qubits embedded in a three dimensional copper cavity. The qubits are fabricated on a sapphire substrate and are measured by coupling them inductively to an on-chip superconducting resonator located in the middle of the cavity. At their flux-insensitive point, all measured qubits reach an intrinsic energy relaxation time in the 6-20 µs range and a pure dephasing time comprised between 3 and 10 µs. This significant improvement over previous works opens the way to the coherent coupling of a flux-qubit to individual spins.

Electronic spins in semiconductors such as NV centers in diamond or phosphorus donors in silicon can reach coherence times up to seconds [1][2][3] and are therefore promising candidates for quantum information processing. However, the main obstacle to an operational spinqubit quantum processor is the difficulty of coupling distant spins. To solve this issue, it has been proposed to couple each spin to a superconducting circuit which acts as a quantum bus and mediates the spin-spin interaction [4,5]. This approach requires to reach the strong coupling regime where the coupling strength g between the spin and the circuit is larger than their respective decoherence rates. Among all superconducting circuits, the largest coupling constants could be obtained with flux qubits (FQ) [6][7][8][9][10][11][12][13] due to their large magnetic dipole. They can reach up to g/2π ∼ 100 kHz for realistic parameters, which is much larger than the best reported spin decoherence rates. This brings a strong motivation for developing FQs with a coherence time T 2 > 2/g ∼ 4 µs.

FQ coherence times reported up to now are limited to T 2 2µs, with a sizeable irreproducibility [10,13]. The reasons for these relatively short coherence times are numerous but stem in part from the poor control of the electromagnetic environnement of the qubit in previously used dc-SQUID readout setups [9][10][11]. A better control of the environnement was recently demonstrated in the case of another superconducting qubit, the transmon, by using a three dimensional (3D) cavity that allows reading out the qubit and protecting it from spontaneous emission [9]. A natural question is therefore whether it is also possible to increase the coherence times of FQs and their reproducibility by using such a setup.

In this work, we present the first measurements of FQs in a 3D cavity. The six qubits measured reach reproducible coherence times T 2 between 2 and 8 µs, which would be already sufficient to reach the strong coupling regime with a single spin. In addition, our results shed light on decoherence of FQs, giving evidence that charge noise is the dominant decoherence mechanism at their flux-insensitive point.

A scheme of the three dimensional cavity used in our experiment is shown in Fig. 1a. The cavity is made of copper to enable the application of an external magnetic field B to the FQs. Its dimensions are chosen for its first mode to be at 5.6 GHz. The sample inserted in the cavity is a sapphire chip with an LC resonator inductively coupled to six FQs, with a coupling constant ∼ 50 MHz. The LC resonator acts as an intermediate coupler [12] between the FQs and the first cavity mode. It appears as a resonance peak at ω LC /2π = 4.643 GHz in the transmission spectrum (Fig. 1b) with a quality factor Q LC ∼ 1.5 × 10 4 determined by the length of the input and output antennas inside the cavity (see Fig. 1a).

Figure 1d presents a colorized SEM micrograph of one of the FQs. It consists of a superconducting aluminum loop of area A intersected by four Josephson junctions. Three of the junctions are identical with a Josephson energy E J /h = 250 GHz and a single electron charging energy E C /h = 3.6 GHz; the fourth junction area is smaller than others by a factor α (see Table I). When the flux threading the loop Φ = B A is close to half a flux quantum Φ 0 /2, the two states characterized by a persistent current ±I P in the loop become degenerate, hybridise and give rise to an energy splitting ∆ called the fluxqubit gap. This circuit behaves therefore as a two-level system and its transition frequency is

ω 01 = √ ∆ 2 + ε 2 with ε = 2 I P Φ b / and Φ b = Φ -Φ 0 /2 [6, 7].
The qubits are fabricated by double angle-evaporation of Al-AlOx-Al on sapphire . We use a tri-layer PMMA-Ge-MAA process [18], which provides a good precision and reproducibility of the junction size and a rigid germanium mask, robust to the O 2 ashing and ion milling cleaning steps, which evacuates efficiently the charges during e-beam lithography. The measurements are performed in a cryogen-free dilution refrigerator at a temperature of 33 mK. The device is magnetically shielded with 2 Cryoperm boxes surrounding a superconducting enclosure. The cavity is also closed using Ecosorb corks and seals, in order to protect the sample from electromagnetic radiation that could generate quasiparticles [19]. The input line is attenuated at low temperature to minimize thermal noise and filtered with impedance-matched radiation-absorbing filters. The readout output line includes several filters, two isolators and a cryogenic HEMT amplifier. Qubit state manipulations are performed by injecting in the input line microwave pulses at ω d ∼ ω 01 , followed by a readout pulse at ω m ∼ ω LC whose amplitude and phase yield the qubit excited state probability [20].

An advantage of the cavity readout is the possibility to measure several qubits in a single run, by fabricating them with different loop areas A i so that the field B i = Φ 0 /2A i at which the flux reaches Φ 0 /2 is different for each qubit. Figure 1e presents the amplitude of the transmitted signal at frequency ω LC as a function of B, showing a dip in the amplitude of the transmitted signal when the frequency of any of the six FQs comes close to ω LC .

Each qubit is characterized by its spectroscopic parameters ∆ and I P , extracted from the dependence of its resonance frequency on the applied flux. These values, given in Table 1, are in good agreement with the predictions of the model described in [21] using both the measured values of α and of the tunnel resistance of the junctions. The coherence properties of each qubit are measured with the appropriate microwave pulse sequence [22]: the energy-relaxation time T 1 , the Ramsey coherence time T 2R from which one gets the Ramsey pure dephasing time (T ϕR ) -1 = (T 2R ) -1 -(2T 1 ) -1 , and the echo decay time T 2E yielding the echo pure dephasing time

(T ϕE ) -1 = (T 2E ) -1 -(2T 1 ) -1 .
We now present detailed measurements on the qubits Q3 and Q4 having the longest coherence times. The flux dependence of their frequency, shown in Fig. 2a (b), yields ∆/2π = 8.47 (9.01) GHz and I P = 169 (160) nA for Q3 (Q4). Since both qubits were designed to have the same parameters, this shows good control of our e-beam lithography and oxidation parameters. We now turn to the coherence times at the so-called optimal point Φ b = 0, where the qubit frequency ω 01 = ∆ is insensitive to first order to flux-noise [10,11]. Energy relaxation (see Fig. 2cd) is exponential with T 1 = 8 µs for Q3 and 18 µs for Q4. Ramsey fringes also show an exponential decay for Q3 with T ϕR = 8.5 µs, and an exponentially decaying beating pattern for Q4 with T ϕR = 10 µs. These features are consistent with the qubit spectra measured after an excitation pulse of ∼ 20µs with a power well below saturation (see Fig. 2a): Q3 line is Lorentzian with a full-width-halfmaximum (FWHM) of 85 kHz, while Q4 line consists of a doublet of two Lorentzians separated by 150 kHz and having a width of 40 kHz, whose origin is discussed further below.

The amplitude of the spin-echo signal decays exponentially (see Fig. 2 c,d) with finite pure dephasing times T ϕE = 17 µs for Q3 and T ϕE = 16 µs for Q4. This moderate improvement compared to the Ramsey pure dephasing times points out to the presence of high-frequency noise in our circuit in contrast to previous reports [11,13]. We attribute this effect to fluctuations in the photon number (photon noise) of one or several cavity modes inducing fluctuations of the qubit frequency due to the dispersive shift [10,23,24]. This noise cannot be compensated by the echo protocol because its correlation time (∼ 100 ns -1µs) given by the mode energy damping is shorter than the echo sequence duration. By plunging the antennas deeper in the cavity for reducing the cavity damping time, the observation of a lower T 2E confirms this explanation [23,24]. Interestingly, removing the photon-noise contribution from the Ramsey pure dephasing time yields a "low-frequency Ramsey dephasing time" ( TϕR ) -1 = (T ϕR ) -1 -(T ϕE ) -1 with TϕR = 16 µs for Q3 and TϕR = 30 µs for Q4. This one order of magnitude improvement compared to previous flux-qubit experiments that reported TϕR in the 0.2 -2.5µs range at the optimal point [10,11,13] is discussed later.

Away from the optimal point, decoherence is governed mainly by flux-noise in agreement with previous works [11,13]. The Ramsey and spin-echo damping become Gaussian as |Φ b | increases, which is characteristic of 1/f noise [25]. Fitting the Ramsey (or echo) envelope as f R,E (t) = e -t/(2T1) e -(Γ ϕR,E t) 2 , we observe a linear dependence of Γ ϕR,E = (T ϕR,E ) -1 on |Φ b |, with Γ ϕR ∼ 4.5 Γ ϕE (see Fig. 3), consistent with dephasing caused by flux-noise. Indeed, assuming a flux-noise power spectral density S Φ (ω) = A Φ /ω, one can show [11,25] that ln(1/ω IR t) ∼ 3.7 predicting Γ ϕR ∼ 4.5 Γ ϕE in agreement with the measured value. We find an amplitude A Φ = (2.5 µΦ 0 ) 2 comparable to previously reported values [11,13].

Γ ϕE = √ A Φ ln 2 |∂ω 01 /∂Φ b | and Γ ϕR = A Φ ln(1/ω IR t) |∂ω 01 /∂Φ b |, with ω IR an infrared cut-
All six qubits were characterized in this way, over several cooldowns (see Table I). Energy relaxation times were found to change from cooldown to cooldown, and also to occasionally change abruptly in the course of one cooldown. Several mechanisms contribute to relaxation; among them, spontaneous emission of a photon by the qubit in the measurement lines (the so-called Purcell regime [26]). Because this spontaneous emission rate T -1 P is also the rate at which a photon coming from the measurement line is absorbed by the qubit, it can be quantitatively determined by measuring the qubit Rabi frequency Ω R for a given microwave power P in at the cavity input. For a qubit coupled symmetrically to the input and output lines, one gets Comparing these estimates with the measured T 1 times (see Table I), we find that Q2 and Q6 are almost Purcelllimited. The intrinsic energy relaxation time 1/(T -1

T P = 2 Ω 2 R P in ω 01 . (1)
1 - T -1
P ) of all six measured qubits at their optimal point is thus in the 6 -20 µs range. This significant improvement over previous reports of T 1 in the 0.5 -4 µs range [9][10][11]27] (except one sample for which T 1 = 12 µs [13]) is probably due to a combination of several factors: good control of the electromagnetic environment in the 3D cavity [9], careful filtering and shielding against infrared radiation [19], low-loss sapphire substrate and different fabrication process.

The frequency dependence of the relaxation rate Γ 1 = (T 1 ) -1 of Q4 in the vicinity of the optimal point is shown in Fig. 4. Large variations are observed, with in particular a reproducible increase of the relaxation rate by a factor 2 over 1 MHz, as also recently observed for a transmon qubit [28]. No corresponding anomaly in the Rabi frequency was observed at this point, which excludes spontaneous emission into the measurement lines (see Eq. ( 1)). We attribute therefore this peak to one resonant microscopic two-level system (TLS) weakly coupled to the qubit [28]. The remaining constant background ∼ (20 µs) -1 could be due to dielectric losses, vortex motion, or out-of-equilibrium quasiparticles.

To estimate the quasiparticle contribution to relaxation, the same measurements were performed at 150 mK, a temperature at which the quasiparticle density is expected to be close to its thermal equilibrium value. The relaxation rate increases to Γ (150 mK) 1

(5 µs) -1 due to quasiparticles with a similar frequency dependence as at 33 mK, although less prononced. Assuming that quasiparticles are mainly generated in the pads of the LC resonator and diffuse into the galvanically coupled qubits, we estimate the density of quasiparticles in the vicinity of the qubit n qp (150 mK) = 1 µm -3 [29]. This density yields a theoretical relaxation rate Γ (qp) 1 = (14 µs) -1 [21], lower than the measured value by a factor 3, a discrepancy which we attribute to the crudeness of the modelling of quasiparticle diffusion. Since Γ (qp) 1 is proportional to the quasiparticle density [29], we conclude that an out-of-equilibrium quasi-particle density of ∼ 0.3 µm -3 would be sufficient to explain the measured energy relaxation times at 33mK, which seems a plausible value in view of earlier reports in other superconducting qubit circuits [9][10][11]. However, the dielectric loss contribution is also important. Taking into account reported values of dielectric loss tangents ∼ 2 × 10 -5 [21], we find Γ (dielectric) 1

∼ (25 µs)

-1 , which is comparable to the measured values. Along these lines, we note that flux qubits fabricated on a high resistivity silicon chip and measured with the same setup at 33 mK showed a five-fold increase in relaxation rate.

Another interesting aspect of our experiment is the long Ramsey pure dephasing time T ϕR measured for Q3 and Q4 at their flux optimal point. Although not quite as long, all measured qubits have T ϕR in excess of 3 µs (see Table I). The origin of decoherence at the optimal point for FQs has so far not been identified. One striking feature is the large variability of reported values of T ϕR at the optimal point for rather similar FQ samples, ranging from 0.2 µs [10] up to 10 µs in this work, whereas T ϕR = 2.5 µs in [13]. A doublet structure in the qubit line was frequently observed at the optimal point [10,11,31], with a greatly varying splitting (20 MHz in [10] and 150 kHz in this work, as seen in Fig. 2b) which was also found to vary in time.

All these features are consistent with charge noise being the dominant noise source limiting T ϕR at the optimal point. The sensitivity to this noise is exponentially dependent on the ratio E J /E C [7] which is twice smaller in [10] compared to the present work. Using the model described in [21], we estimate a charge modulation amplitude of ∼ 50 kHz for Q3/Q4 compared to ∼ 120 MHz for [10] (∼ 300 kHz for [13]), yielding a three orders of magnitude lower charge noise sensitivity which explains qualitatively the difference in dephasing time. The doublet lineshape of Q4 can be attributed to slow fluctuations of the electron number parity on one of the qubit islands [31] as observed for transmon qubits [11,32].

In conclusion, we have characterized the coherence properties of six FQs in a three-dimensional microwave cavity. We consistently find intrinsic energy relaxation times T 1 ranging between 6 and 20 µs, a significant improvement over previous FQ measurements that we attribute to good control of the electromagnetic environment provided by the 3D cavity, low-loss substrate, and careful filtering. We identify weakly coupled twolevel systems, quasiparticles and dielectric losses as likely sources of energy damping. At the optimal point, long Ramsey pure dephasing times up to 10 µs are measured, limited by a combination of photon noise and charge noise with roughly equal contribution. We argue that charge noise is the dominant microscopic dephasing mechanism for FQs at the optimal point, and that its effect can be greatly minimized by chosing proper qubit parameters. Our results prove that FQs can reliably reach long coherence times, which opens new perspectives for the field of hybrid quantum circuits, in particular for the coherent coupling of single spins to superconducting circuits. Based on the geometry of Q4, we estimate that at this temperature there are on average just 0.1 thermal quasiparticles in island I 1 . This means that the with high probability qubit decay can only be induced by a quasiparticle tunneling from the LC inductor gavanically connected to the qubit into one of the small island. Therefore, indicating with x qp = n qp /2ν 0 ∆ 0 the normalized quasiparticle density in the inductor (ν 0 is the aluminum density of states at the Fermi energy and ∆ 0 = 200 µeV), we have

Γ (qp) 1 4 π E J x qp 2∆ 0 ω 01 1| sin ϕ 1 2 |0 2 + 1| sin ϕ 2 2 |0 2 (B1)
The matrix elements can be estimated numerically after having obtained the approximate eigenstates of the Hamiltonian in Eq. (A1) and we find

1| sin ϕ 1 2 |0 1| sin ϕ 2 2 |0 0.21 (B2)
Since in thermal equilibrium we have x qp = 2πT /∆ 0 e -∆0/T , using the values E J = 250 GHz and ω 01 = 9.01 GHz, from Eq. (B1) we arrive at the value Γ (qp) 1

= 0.07 µs -1 reported in the main text. This calculation neglects completely the role of quasiparticles having tunneled into the other islands of the qubit and therefore underestimates Γ

Γ Ci 1 = 16 tan δ i E Ci 3 j,k=1 A ji A ki N jk (B3)
where tan δ i is the inverse quality factor of the material causing the loss, the Hermitian matrix N jk is defined, using the matrix elements of the number operators, as

N jk = 1|n j |0 0|n k |1 (B4)
and the dimensionless matrix A ij accounts for the different coupling strengths between voltage fluctuations V i in each capacitors and the number operators; specifically, this coupling adds the term 2e i,j

n j A ji V i (B5)
to the Hamiltonian in Eq. (A1). For a simple order-of-magnitude estimate, we take the entries of matrix N ik to be approximately identical

N ik ≈ 0.09 (B6)
with the numerical value calculated numerically using the (approximate) eigenstates of the Hamiltonian in Eq. ( 1). Then summing over indices j and k and over all four junctions, we find for the junction capacitors

A ij A ik E C 4 c r [1 + 3α + c r (1 + α)] 2 (B8)
From the measured decay rate at 33 mK, we can put upper bounds on the loss tangents for junction and ground capacitors; using Eqs. (B3) and (B7)-(B8). Assuming that the relaxation of the qubit would be due mainly to dielectric losses, we find tan δ J ∼ 2.1 × 10 -6 , (B9) and tan δ G ∼ 2.4 × 10 -5 .

(B10)

The value in Eq.( 10) is an order of magnitude larger than recent estimates [12], where tan δ J ≤ 3 × 10 -7 excluding the junction dielectric loss as the likely source of relaxation in our qubits. On the contrary, the number in Eq.( 11) is close the the bound tan δ G ≤ 2 × 10 -5 of Ref. [13], indicating that dielectric losses in the substrate or interface oxides can be one of the main sources of low-temperature relaxation. We note that similar experiments, realized in the same conditions but on a high-resistivity silicon chip, yielded an increased relaxation rate by a factor of ∼ 5.

Figure 1 :

 1 Figure 1: (a) Cut-away representation of the 3D cavity, with the LC circuit (in blue) on its sapphire chip. The green arrow represents the applied magnetic field B. The red arrow represents the ac electric field E(t) of the first mode of the cavity. (b) Transmission spectrum of the cavity coupled to the LC resonator. The first peak at frequency ωLC /2π = 4.643 GHz corresponds to the resonance of the LC resonator while the two other peaks correspond to the first modes of the cavity. (c-d) Circuit diagram and colorized SEM micrograph showing the FQ (in red) inductively coupled to the LC resonator (in blue). (e) Amplitude of the transmitted signal at frequency ωLC as a function of B, showing the signal from the six qubits.

Figure 2 :

 2 Figure 2: Characterisation of Q3 (top panels) and Q4 (bottom panels): (a-b) (left panels) Measured qubit frequency ω01(Φ b ) (blue dots) and fit (red curve) yielding the qubit parameters ∆ and IP . (Right panels) Spectroscopy data at Φ b = 0 (blue dots). Q3 spectrum is fitted with a single Lorentzian peak ; Q4 spectrum is fitted with a sum of two Lorentzian peaks separated by 150kHz. (c-d) (left panels) Qubit energy relaxation and spin-echo measurements. The excited state probability is plotted as a function of the delay between the π pulse and readout pulse (blue dots) or between the two π/2 pulses of the echo sequence (purple dots). Red (orange) solid line is an exponential fit to the energy relaxation (spin-echo) data. (Right panels) Measured Ramsey fringes (blue solid line), with fit (red solid line) to exponentially damped single (top) and double (bottom) sine curves.

Figure 3 :

 3 Figure 3: Pure dephasing rates of Q4 as a function of Φ b : Experimental (dots) and fitted (line, see text) echo (a) and Ramsey (b) dephasing rates.

Figure 4 :

 4 Figure 4: Frequency dependence of the relaxation rate Γ1 of Q4 in the vicinity of its optimal point at T = 33mK (blue squares) and T = 150mK (red circles), showing an increased relaxation rate caused by quasi-particles.

Figure 5 :

 5 Figure 5: (a) Circuit diagram of the flux qubit including the main contributions of the geometric capacitance. (b) SEM micrograph of one of the flux qubits showing the different islands. The island I3 (in blue) is connected galvanically to the pads of the LC resonator. (c) Persistent current IP vs its gap ∆. For the blue solide line, the parameters of the model were EJ /h = 250 GHz, EC /h = 3.6 GHz, cr = 0.1102. The upper (lower) blue dotted line corresponds to EJ /h = 265 GHz and EC /h = 4 GHz (EJ /h = 235 GHz and EC /h = 3.2 GHz). The green circles are the experimental values. (d) Persistent current IP vs α for the same parameters as (c). The horizontal error bars correspond to the confidence interval in the measure of α by SEM observation.

1 .

 1 Turning now our attention to the dielectric loss relaxation mechanism, we use Fermi's golden rule and the quantum fluctuation-dissipation relation to write each capacitor's C i contribution to the decay rate in the form[2]

A 2 [ 1 + 2 (

 212 ij A ik E C (1 + 2c r ) 2 + 2 + α(3 + c r ) 3α + c r (1 + α)] B7)Here we neglect for simplicity the renormalization of parameters, discussed above, due to the geometrical capacitances not included in the diagram of Fig. a. Similarly, summing over the geometric capacitors we get i∈G

Table I :

 I Parameters of the different measured FQs. Here ∆ is the FQ gap, IP is the persistent current flowing in the loop

		∆/2π (GHz) IP (nA) α T1(µs) TP (µs) TϕR(µs) TϕE(µs)
	Q1	2.70	245 0.61 6-10 1.1 × 10 5	2	7
	Q2	4.91	207 0.55	2	3	-	-
	Q3	8.47	169 0.49 6.5-8	30	8	17
	Q4	9.01	160 0.49 13-18	270	10	16
	Q5	9.71	171 0.43 5.5-12	90	5	>100
	Q6	15.15	140	0.4	4	12	-	-

of the qubit, α corresponds to the measured ratio between the small and big junctions, T1 is the relaxation time, TP the Purcell limit time due to the coupling of the qubit with the cavity, TϕR and TϕE the coherence times obtained by Ramsey and Echo sequences respectively.

We would like to acknowledge fruitful discussions with M. Devoret, A. Lupascu and within the Quantronics group. We thank P. Sénat, P.-F. Orfila, D. Duet, J.-C. Tack, P. Pari, D. Bouville, P. Forget and M. de Combarieu for their technical support. This work was supported by the ANR CHIST-ERA project QINVC, the C'Nano IdF project QUANTROCRYO, the ERC project CIRQUSS, the JSPS and by the EU under REA Grant Agreement No. CIG-618258 (GC).

Supplementary Materials

Appendix A: Hamiltonian of the 4-junction flux qubit

In this section, we present a model of a 4-junction flux qubit taking into account the main geometric capacitance terms and compare the calculation results with our experimental data.

Fig. 5a presents the circuit diagram of our model. The loop of the qubit is formed by four islands I i=1,4 shown on Fig. 5b and connected by 4 Josephson junctions. Three junctions are identical while the area of the junction between islands I 1 and I 4 is smaller than the others by a factor α.

The full (4×4) geometric capacitance matrix of the circuit was calculated using a three-dimensional electromagnetic simulator [1]. We include in the diagram three capacitors which account for the main contribution of the geometric capacitance. Capacitor C c represents the capacitance between the largest qubit island I 1 and I 3 ; the latter is connected galvanically to the pads of the LC resonator. Capacitors C g1 and C g3 are respectively the capacitances to ground of I 1 and I 3 .

For qubit Q4, the values of these capacitances were found to be C g1 = 0.04 fF, C g3 = 477 fF and C C = 0.57 fF.

In the limit C g3 C g1 relevant to our experiment, the circuit Hamiltonian can be written as:

where ϕ i are the phase differences across the three large junctions of Josephson energy E J , the operators n i are conjugated to ϕ i and count the number of Cooper pairs tunneling across the junctions, C is the capacitance of the large junctions, E C is given by

, Φ is the magnetic flux through the qubit loop, and Φ 0 is the magnetic flux quantum.

A complete description of the system includes all the inter-island capacitances, the capacitances to ground of the two small qubit islands and the inductance of the qubit loop. Assuming symmetry between the small islands, it can be shown [2] that the main effect of the additional capacitances is to slightly renormalize the parameters E C , α, and c r entering the charging energy part of Eq. (A1). In our numerical calculations we take this renormalization into account, using for the geometric capacitances the values extracted from the full capacitance matrix of our system but neglected the effect of the loop inductance.

We approximately solve the Schrödinger equation in the charge basis by truncating the Hilbert space to a large but finite number of charge states (∼ 25 3 ) and compare in Figs. c andd the results of the calculation to the measured values of the flux qubit gap ∆ and persistent current I P .

The Josephson energy E J was estimated from Ambegaokar-Baratoff relation E J = h∆ 0 /(8e 2 R N ) where ∆ 0 is the superconducting gap of a 25-nm-thin aluminium film and R N the tunneling resistance of the junction at low temperature and in the normal state. Assuming ∆ 0 = 200 µeV [3] and measuring the room temperature tunneling resistance of the junctions, we find E J /h = 250 ± 15 GHz.

The charging energy of the junctions E C was estimated by measuring their size and using standard values given in the literature for the capacitance per unit area. Assuming c = 90 ± 10 fF/µm 2 [4-6], we find C = 5.4 ± 0.6 fF and E C /h = 3.6 ± 0.4 GHz.

The experimental data match well with the calculation and the major influence of adding the geometric capacitance terms is to reduce the flux qubit gap ∆ by around 1 GHz.

Appendix B: Relaxation Rate Calculations

In this section, we discuss the calculations of the relaxation rates due to quasiparticles and dielectric losses; more details will be presented elsewhere [2].

Our estimate of the contribution Γ (qp) 1

of quasiparticle tunneling to the relaxation rate of Q4 is based on the general theory of quasiparticle effects developed in Refs. [7,8]. It was shown there that in a multi-junction qubit, Γ (qp) 1 is obtained by summing the contributions of each individual junction, and that at low temperatures each term in the sum is proportional to the density of quasiparticles. A number of experiments [9][10][11] have shown that in aluminum devices at temperatures 150 mK, quasiparticles are in thermal equilibrium.