# High-Sensitivity ac-Charge Detection with a MHz-Frequency Fluxonium Qubit

Baldo-Luis Najera-Santos, Rémi Rousseau, Kyrylo Gerashchenko, Himanshu Patange, Angela Riva, Marius Villiers, Tristan Briant, Pierre-François Cohadon, Antoine Heidmann, José Palomo, et al.

► To cite this version:

Baldo-Luis Najera-Santos, Rémi Rousseau, Kyrylo Gerashchenko, Himanshu Patange, Angela Riva, et al.. High-Sensitivity ac-Charge Detection with a MHz-Frequency Fluxonium Qubit. Physical Review X, 2024, 14 (1), pp.011007. 10.1103/physrevx.14.011007. hal-04777888v2

HAL Id: hal-04777888

https://hal.science/hal-04777888v2

Submitted on 12 Nov 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
High-Sensitivity ac-Charge Detection with a MHz-Frequency Fluxonium Qubit

B.-L. Najera-Santos,¹ R. Rousseau,¹,² K. Gerashchenko,¹ H. Patange,¹ A. Riva,³,⁴ M. Villiers,³ T. Briant,¹

P.-F. Cohadon,¹ A. Heidmann,¹ J. Palomo,³ M. Rosticher,³ H. le Sueur,⁴ A. Sarlette,³ W. C. Smith,³,*

Z. Leghtas,³ E. Flurin,⁴ T. Jacqmin,¹ and S. Deléglise¹,†

¹Laboratoire Kastler Brossel, Sorbonne Université, ENS-Université PSL, CNRS,
Collège de France, 4 place Jussieu, Paris F-75252, France

²Alice & Bob, 53 Bd du Général Martial Valin, 75015 Paris, France

³Laboratoire de Physique de l'Ecole normale supérieure, ENS-PSL, CNRS, Inria,
Centre Automatique et Systèmes (CAS), Mines Paris, Université PSL, Sorbonne Université,
Université Paris Cité, Paris, France

⁴Quantronics group, Université Paris-Saclay, CEA, CNRS, SPEC, 91191 Gif-sur-Yvette Cedex, France

<figure><img src="image_2.png" /></figure>

(Received 28 July 2023; revised 3 November 2023; accepted 5 December 2023; published 24 January 2024)

Owing to their strong dipole moment and long coherence times, superconducting qubits have demonstrated remarkable success in hybrid quantum circuits. However, most qubit architectures are limited to the GHz frequency range, severely constraining the class of systems they can interact with. The fluxonium qubit, on the other hand, can be biased to very low frequency while being manipulated and read out with standard microwave techniques. Here, we design and operate a heavy fluxonium with an unprecedentedly low transition frequency of 1.8 MHz. We demonstrate resolved sideband cooling of the “hot” qubit transition with a final ground state population of 97.7%, corresponding to an effective temperature of 23 µK. We further demonstrate coherent manipulation with coherence times $T_1 = 34$ µs, $T_2^* = 39$ µs, and single-shot readout of the qubit state. Importantly, by directly addressing the qubit transition with a capacitively coupled waveguide, we showcase its high sensitivity to a radio-frequency field. Through cyclic qubit preparation and interrogation, we transform this low-frequency fluxonium qubit into a frequency-resolved charge sensor. This method results in a charge sensitivity of $33 \ \mu\text{e}/\sqrt{\text{Hz}}$, or an energy sensitivity (in joules per hertz) of $2.8\hbar$. This method rivals state-of-the-art transport-based devices, while maintaining inherent insensitivity to dc-charge noise. The high charge sensitivity combined with large capacitive shunt unlocks new avenues for exploring quantum phenomena in the 1–10 MHz range, such as the strong-coupling regime with a resonant macroscopic mechanical resonator.

DOI: 10.1103/PhysRevX.14.011007

Subject Areas: Quantum Physics

I. INTRODUCTION

Superconducting qubits consist of engineered quantum systems with lowest-level spacings designed to host a two-level system which can be manipulated and read out via its dipolar interaction with electromagnetic fields. Their strong dipole moment is also beneficial to interface them with other physical systems. For instance, fluorescence from individual electronic spins was successfully detected using a superconducting qubit-based microwave-photon detector [1] operating close to 7 GHz. Additionally, in the

realms of circuit quantum acoustodynamics (cQAD), the coupling between a qubit and a piezoelectric resonator is used to detect and manipulate the phononic state, typically within the 2–10 GHz range [2–5]. However, adapting these sensing schemes to lower frequencies, below the conventional operating frequency of superconducting qubits, introduces distinct challenges.

First, in a typical circuit quantum electrodynamics setup, the qubit state is read out thanks to the dispersive shift imparted to a nearby superconducting resonator. As the dispersive shift quickly drops for a cavity detuning exceeding the qubit anharmonicity, weakly anharmonic qubits, such as transmons, require nearly resonant resonators with dimensions scaling inversely with the frequency (as an illustration, a 1 MHz $\lambda/2$ coplanar cavity requires a 100-m-long waveguide). Second, systems with frequencies lower than $k_B T/h$, where T stands for the environmental temperature, are coupled to a hot thermal bath with which they exchange photons randomly, quickly turning pure quantum states into statistical mixtures.

<footnote>*Present address: Google Quantum AI, Santa Barbara, California, USA.</footnote>

<footnote>†samuel.deleglise@lkb.upmc.fr</footnote>

<footnote>Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI.</footnote>
In recent years, significant progress has been made in overcoming these challenges. Notable contributions include the development of a 14 MHz heavy-fluxonium qubit with a long coherence time and fast manipulation through fast-flux gates [6]. Furthermore, operation of a fluxonium qubit dispersively coupled to a 690 MHz piezoelectric mechanical system was demonstrated earlier this year [7].

In this work, we demonstrate a fluxonium qubit with a transition frequency as low as 1.8 MHz, achieving coherent operation and a charge sensitivity of $33 \mu e/\sqrt{Hz}$, reflecting its potential for coupling with other devices operating in the MHz range. We achieve single-shot readout and direct preparation in the qubit state basis using sideband cooling, attaining a preparation fidelity above 97%. Based on this fidelity, we calculate an effective temperature of approximately 23 µK. We also demonstrate direct resonant manipulation of the qubit state with a charge drive as low as $5 \times 10^{-3}$ Cooper pairs. This value corresponds to a single-shot charge sensitivity of $10^{-2}$ e. In order to compare the sensitivity of our qubit-based detection scheme to other charge sensors, we accumulate statistical data through a cyclic qubit preparation and interrogation sequence. The charge sensitivity demonstrated here, at $\delta q \approx 33 \mu e/\sqrt{Hz}$, rivals that of the most advanced transport-based devices [8–18], while maintaining intrinsic insensitivity to dc-charge noise. Furthermore, the capacitance $C \sim 50$ fF of the fluxonium electrode demonstrated in this work exceeds the typical gate capacitance of single-electron transistors (SETs) [9,12] by approximately 2 orders of magnitude, resulting in a record-low energy sensitivity $\delta q^2/2C \sim 2.8 \hbar$. This metric can be regarded as more fundamental than the bare charge sensitivity, which gets diluted when the self-capacitance of the probed system exceeds that of the sensor in a real-world charge-sensing scenario.

For instance, the demonstrated charge sensitivity, combined with large gate capacitance demonstrated here, are well suited to explore quantum phenomena with low-frequency mechanical systems. For one, the frequency and electrode capacitance demonstrated in our work align with those found in cutting-edge vacuum-gap dispersive electromechanical systems [19]. Additionally, the single-shot charge sensitivity demonstrated here is sufficient to detect the zero-point motion of such a system placed in a dc-biased vacuum gap capacitor [20]. Achieving the strong coupling between a low-frequency mechanical resonator and a superconducting qubit would enable us to test the superposition principle in a regime where general relativity and quantum mechanics interplay [21].

## II. CIRCUIT DESIGN

The heavy-fluxonium circuit is shown in Fig. 1. The qubit itself is composed of a small Josephson junction (energy $E_J = \phi_0^2/L_J$, $\phi_0$ denoting the reduced quantum of flux $\phi_0 = \Phi_0/2\pi$, with $\Phi_0 = h/2e$) shunted by a large

<figure><img src="image_3.png" /><figcaption>FIG. 1. Circuit diagram and implementation of the fluxonium qubit, controls, and readout. Optical micrograph (a) and circuit diagram (c) of the fluxonium qubit composed of a capacitor (green), an array of 360 Josephson junctions (blue), and a single junction (red). The qubit is capacitively coupled to the readout resonator (purple). The magnetic flux through the superconducting loop can be rapidly tuned via the current passing through the flux line (orange), and the qubit is capacitively coupled to a charge line (cyan). (b) Scanning electron micrograph of the fluxonium single junction (red) and four junctions of the array (blue). (d) Two-tone spectrum centered on the flux frustration point $\phi_{ext}/2\pi = 0.5$. The color scale in the image represents the phase of the reflected probe pulse. The fit to the data (dotted lines) yields the qubit parameters $E_J/h = 5.178$ GHz, $E_C/h = 0.4144$ GHz, and $E_L/h = 0.18$ GHz. Central inset: detailed view of the avoided crossing near the frustration point. The left- and right-hand insets represent the energy-level diagrams of the four lowest eigenstates, with the potential in gray and the wave functions in colors for an external flux of $\phi_{ext}/2\pi = 0.5$ and $\phi_{ext}/2\pi = 0.6$, respectively.</figcaption></figure>

capacitance to ground (capacitive energy $E_C = e^2/2C$) and a superinductance (inductive energy $E_L = \phi_0^2/L$) formed by 360 large Josephson junctions in series. We ensure that each junction of the array has a negligible phase-slip rate by taking $E_{J,A}/E_p \gtrsim 3$, where $E_{J,A}/h = 65$ GHz is the Josephson energy of each array junction, and $E_p/h = 17.9$ GHz is the junction plasma frequency [22]. In this regime the junction chain behaves as a linear inductor and the circuit Hamiltonian writes

$$ \hat{H}_Q = -E_J \cos(\hat{\phi} - \phi_{\text{ext}}) + 4E_C[\hat{n} - n_g(t)]^2 + \frac{E_L}{2} \hat{\phi}^2. \quad (1) $$
In this equation, $\hat{\phi}$ represents the superconducting phase across the junction, and $\hat{n}$ denotes its conjugate variable (the Cooper pair number), $\phi_{\text{ext}} = 2\pi\Phi_{\text{ext}}/\Phi_0$, where $\Phi_{\text{ext}}$ stands for the magnetic flux threading the superconducting loop, and $n_g(t)$ is the offset charge on the capacitor pad. $\phi_{\text{ext}}$ can be controlled by an on-chip flux line, and $n_g(t)$ can be controlled by a capacitively coupled coplanar waveguide. While the fluxonium spectrum is insensitive to a dc-charge offset [23], the main goal of this work is to evaluate the sensitivity of the qubit to a nearly resonant ac-charge modulation.

### III. QUBIT SPECTRUM

The circuit operates in the heavy-fluxonium regime, characterized by the two conditions $E_J \gg E_L$ and $E_J \gtrsim 10E_C$. The first condition ensures that the potential experienced by the positionlike variable $\hat{\phi}$ consists of multiple wells with distinct minima. The second condition ensures that the lowest energy eigenstates are well localized within each well, with a small tunneling probability between neighboring wells. The magnitude of the tunneling rate $E_S/\hbar$ is exponentially suppressed as a function of $\sqrt{8E_J/E_C}$, which relates the height of the potential barrier $2E_J$ to the zero-point energy $\frac{1}{2}\sqrt{8E_J E_C}$.

We denote $|g\rangle$ and $|e\rangle$ ($|f\rangle$ and $|h\rangle$) as the fundamental (first excited) states of the two lowest wells. Two families of transitions are observed in the two-tone spectroscopy of Fig. 1(d): intrawell (or plasmonic) transitions, $|g\rangle \to |f\rangle$ and $|e\rangle \to |h\rangle$, that are only weakly dependent on the external flux $\phi_{\text{ext}}$, and interwell transitions $|e\rangle \to |f\rangle$ and $|g\rangle \to |h\rangle$, that feature a linear dependence with $\phi_{\text{ext}}$.

Away from the two sets of flux points $\phi_{\text{ext}} \equiv 0[2\pi]$ and $\phi_{\text{ext}} \equiv \pi[2\pi]$ for which the potential is symmetric with respect to the position-like variable, the interwell transition $|g\rangle \to |e\rangle$ is exponentially suppressed, acting as a selection rule that can be used to protect a microwave qubit against relaxation [24] [right-hand inset of Fig. 1(d)].

At the flux-frustration point $\phi_{\text{ext}} = \pi$, the eigenstates undergo a transition, switching from localized modes around a single potential well to symmetric and antisymmetric superpositions of these well states. This transition results in a significant overlap of the flux wave functions, as evidenced by the magnitude of the flux matrix element $|\langle g|\hat{\phi}|e\rangle| \sim \pi$ [left-hand inset of Fig. 1(d)]. Importantly, at this point, the weakness of the tunneling element leads to a reduced qubit transition frequency $\omega_{ge} = E_S/\hbar$. The value of $\omega_{ge}$ can be tuned over several orders of magnitude by adjusting the circuit parameter $E_J/E_C$. In our specific case, we have chosen a transition frequency of 1.8 MHz, which closely matches the oscillation frequency of existing macroscopic mechanical systems based on suspended membranes [25,26]. Notably, this frequency is approximately one order of magnitude lower than the lowest frequency ever reported using superconducting qubits [6].

### IV. SIDEBAND COOLING

With this low frequency, the qubit has almost equal ground and excited state populations at thermal equilibrium. Inspired by experiments with trapped ions [27] and optomechanical systems [28], we initialize the qubit in a pure state by driving the readout cavity with a detuned tone. By sweeping the reset tone frequency $\omega_p$ in the vicinity of the cavity resonance $\omega_R/2\pi = 5.64$ GHz, we observe two distinct processes at the sideband frequencies $\omega_R \pm \omega_{ge}$, corresponding to the transitions $|g0\rangle \to |e1\rangle$ and $|e0\rangle \to |g1\rangle$ [see Fig. 2(a)]. More quantitatively, the qubit-resonator Hamiltonian can be linearized around the intracavity drive amplitude $\alpha$. For large drive amplitude and dropping all terms rotating at the drive frequency, we arrive at the Hamiltonian (see Appendix C)

$$ \hat{H} = \hat{H}_{\text{Q}} + \hbar\Delta_{\text{R}}\hat{a}^{\dagger}\hat{a} + \hbar g \cos(\hat{\phi} - \phi_{\text{ext}})(\alpha\hat{a}^{\dagger} + \alpha^{*}\hat{a}), \quad (2) $$

where $\hat{a}$ is the annihilation operator for photons in the readout cavity, $\Delta_R = \omega_p - \omega_R$ the pump-cavity detuning, and $\hbar g = E_J \phi_{zpf,R}^2$, with $\phi_{zpf,R}$ the zero-point fluctuations of the readout mode quantifying the energy-participation ratio of the cavity in the fluxonium junction. This Hamiltonian, expressed in a frame rotating at the drive frequency, describes the interaction between the fluxonium and an effective cavity mode of frequency $\Delta_R$. When $\Delta_R$ matches the frequency $\omega_{ge}$ ($-\omega_{ge}$), the interaction reduces to the Jaynes-Cummings (anti-Jaynes-Cummings) model between the cavity and the qubit. Furthermore, owing to the large cavity damping rate $\kappa/2\pi = 2.4$ MHz $\gg g|\alpha|$, the cavity field dynamics can be adiabatically eliminated, leading to the Purcell-like loss operators (see Appendix C),

$$ L^{\pm} = \frac{2g|\alpha|}{\sqrt{\kappa}} \langle e | \cos(\hat{\phi} - \phi_{\text{ext}}) | g \rangle \sigma^{\pm}, \quad (3) $$

where the $\pm$ sign depends on the sideband addressed by the drive pulse $\Delta_R \simeq \pm\omega_{ge}$.

At the flux-frustration point $\phi_{\text{ext}} = \pi$, the matrix element $\langle e | \cos(\hat{\phi} - \phi_{\text{ext}}) | g \rangle$ cancels due to the opposite parity of the qubit wave functions. Prior to the 10 µs reset pulse, we thus offset the flux by about $10^{-3} \times \Phi_0$ which corresponds to $\omega_{ge}/2\pi \simeq 10$ MHz, and we ramp it back to the frustration point afterward [see Fig. 2(b)]. In order to avoid undesired mixing of qubit states caused by nonadiabatic effects [6] while minimizing qubit decay, we have chosen a ramp duration of 2 µs.

The qubit population is then detected by standard dispersive readout. We were unable to directly resolve the qubit states $|g\rangle$ vs $|e\rangle$, due to a too small dispersive shift of the readout cavity. The qubit population was thus obtained by first mapping the population from $|g\rangle$ to $|h\rangle$, thanks to a $\pi$ Rabi pulse. The population in $|h\rangle$ is then measured by standard dispersive readout [see Fig. 2(b)].
<figure><img src="image_4.png" /><figcaption>FIG. 2. Sideband preparation of the fluxonium qubit. (a) Level diagram illustrating the sideband reset protocol: The qubit is reset to either $|g\rangle$ or $|e\rangle$ by driving one of the sideband transitions: $|e0\rangle \rightarrow |g1\rangle$ or $|g0\rangle \rightarrow |e1\rangle$. The subsequent rapid decay of the cavity photon ensures a directional transition toward the desired qubit state. (b) Pulse sequence for the reset protocol. The flux bias (blue) is set to the target value within 2 µs. A reset tone (purple) is applied to the readout resonator port for 10 µs. The flux is then reset to $\phi_{\text{ext}} = \pi$ within 2 µs. The qubit state is read out by transferring the population from $|g\rangle$ to $|h\rangle$ (orange $\pi$ pulse) followed by single-shot dispersive readout in the *eh* manifold (yellow pulse). (c) Qubit population (color scale) as a function of flux bias (x axis) and reset pulse detuning (y axis). The reset pulse power is adjusted to maintain a constant intracavity field at $\omega_R \pm \omega_{ge}(\phi_{\text{ext}})$. The horizontal dashed line indicates the readout frequency. The nominal working points for $|e\rangle$ and $|g\rangle$ preparation are denoted by the orange and blue dots, respectively. The orange and blue dot-dashed lines represent the predicted frequencies based on the Hamiltonian parameters in Fig. 1. (d) Final qubit population as a function of reset pulse duration [other parameters are the nominal parameters indicated in (c)]. The population extracted from single-shot readout distributions is shown on the left-hand axis, while the right-hand axis displays the population corrected for decay and mislabeling during readout (see Appendix D). Exponential fits to the data (yellow and blue lines) yield preparation times of 1.9 µs ($|e\rangle$) and 1.7 µs ($|g\rangle$) with final occupations of 97.7% and 97.7%, respectively.</figcaption></figure>

Figure 2(c) shows a map of the final qubit population as a function of the external flux and detuning of the reset tone. The expected location of the upper and lower sideband processes, as predicted by diagonalizing the Hamiltonian Eq. (1), are overlaid as dashed orange and blue lines, respectively. The dependence of the final qubit population as a function of the duration of the reset pulse [see Fig. 2(d)], for the preparation in $|e\rangle$ (preparation in $|g\rangle$), is well fitted by an exponential curve of time constant 1.9 µs (1.7 µs). The raw single-shot probability of detection is given by $P_g^{\text{prep } g} = 86.4\%$ for a qubit prepared in $|g\rangle$ ($P_e^{\text{prep } e} = 92.8\%$ for preparation in $|e\rangle$). By correcting for mislabeling and decay during readout (see Appendix D), we infer a state-preparation efficiency of 97.7% for the qubit preparation in $|g\rangle$ and 97.7% for the preparation in $|e\rangle$.

V. QUBIT COHERENCE

After having established this preparation process, we first investigate the energy relaxation of the $|e\rangle$ and $|g\rangle$ fluxonium states, toward a thermal state: Figure 3(a) displays the qubit population versus the delay time after preparation in either $|e\rangle$ or $|g\rangle$. When a qubit interacts with a thermal environment of occupation $n_{\text{th}}$, it experiences two relaxation channels, described by the operators $\sqrt{\Gamma_{\uparrow}}\hat{\sigma}^{+}$ and $\sqrt{\Gamma_{\downarrow}}\hat{\sigma}^{-}$, where $\Gamma_{\uparrow} = \gamma n_{\text{th}}$ and $\Gamma_{\downarrow} = \gamma(n_{\text{th}} + 1)$. Here, we have introduced the intrinsic damping of the circuit $\gamma = \omega_{ge}/Q_{\text{intr}}$, corresponding to the energy decay rate in an ideal zero-temperature environment, with $Q_{\text{intr}}$ the intrinsic circuit quality factor. In the case of low-frequency transitions, such as the $\{|e\rangle, |g\rangle\}$ manifold, where $\hbar\omega_{ge} \ll k_B T$, both emission and absorption are nearly equal $\Gamma_{\uparrow} \approx \Gamma_{\downarrow} \equiv \Gamma$ since $n_{\text{th}} \sim k_B T/\hbar\omega_{ge} \gg 1$. This leads to an exponential relaxation toward the statistical mixture $\rho^{\text{th}} = (|e\rangle\langle e| + |g\rangle\langle g|)/2$ at a rate $2\Gamma$ [see Eq. (F4) of Appendix F for details). By fitting exponential curves to the data of Fig. 3(a), we obtain $T_1 = 1/2\Gamma = 34$ µs.

As the qubit frequency explored in the current work extends well below the values reported in the literature so far [6], it is important to determine whether the qubit transition couples with a thermal environment or is primarily constrained by technical noises (e.g., $1/f$ charge or flux noise). To examine this, we conducted $T_1$ measurements similar to those shown in Fig. 3(a), while varying the cryostat base temperature (see Appendix E). To extract the effective temperature of the circuit, which is often observed to be higher than the cryostat temperature [29], we use the measured residual populations of the $\{|f\rangle, |h\rangle\}$ manifold as a local probe for the circuit temperature $T_{\text{circ}}$ (see Appendix E for details). Although the built-in temperature sensor of the cryostat indicates a minimal temperature of 7 mK, we have found that the circuit only thermalizes to $T_{\text{circ}} \approx 59$ mK. Furthermore, we observe a nearly linear
<figure><img src="image_5.png" /><figcaption>FIG. 3. Qubit coherence. (a) Energy relaxation time ($T_1$) measured at the flux point $\phi_{ext} = \pi$. The raw qubit population is plotted as a function of delay time after preparation in $|e\rangle$ (yellow dots) or $|g\rangle$ (blue dots) as the qubit relaxes toward a thermal mixture. A common exponential fit yields $T_1 = 34$ µs. The apparent imbalance in the equilibrium population (dashed line) is due to the residual decay of the intermediate state $|f\rangle$ used for readout. (b) Ramsey experiment after state preparation in $|g\rangle$. The exponential decay of the Ramsey fringes (black line) yields $T_2^* = 39.7$ µs and a transition frequency of 1.8 MHz. (c) $T_2^*$ (yellow curve, right-hand axis) and transition frequency (blue curve, left-hand axis) as a function of flux bias. The points highlighted in red correspond to the data presented in (b). The blue solid line is the predicted qubit frequency for the Hamiltonian parameters given in Fig. 1.</figcaption></figure>

relationship between $\Gamma$ and $T_{circ}$ in the range $59 \le T_{circ} \le 318$ mK, suggesting an intrinsic damping $\gamma = 12$ s$^{-1}$, or $Q_{intr} = 9.6 \times 10^5$ independent of temperature, with a marginal impact of technical noise approximately equal to the thermal contribution at $T_{circ} = 59$ mK (see Appendix E). This observation indicates that, despite its ultralow operational frequency, our qubit is marginally impacted by $1/f$ noise. This outcome stands in contrast with recent studies on frequency-tunable fluxonium [30] and may be attributed to the small superconducting loop area used in our circuit, limiting the influence of flux noise.

Finally, we probe the qubit dephasing time, denoted as $T_2^*$, as a function of external flux. To achieve this, we conducted Ramsey sequences on the $|g\rangle \rightarrow |e\rangle$ transition. As seen in Fig. 3(c), the coherence time reaches its maximal value of approximately 40 µs at the flux frustration point, $\phi_{ext} = \pi$. Indeed, as shown in Fig. 3(c), the qubit frequency is to first order insensitive to fluctuations in the external magnetic flux at this point. The Ramsey fringe measurement at $\phi_{ext} = \pi$ is depicted in Fig. 3(b). Notably, the measured coherence is not too far from the upper limit of $2T_1$, suggesting a pure dephasing rate of $\Gamma_{\phi} = 1/2T_1 - 1/T_2^* = (97$ µs)$^{-1}$.

## VI. AC-CHARGE SENSITIVITY OF THE FLUXONIUM QUBIT

In the following, we evaluate the sensitivity of the fluxonium to a nearly resonant ac-charge drive. We delve first into the theoretical advantages of the fluxonium qubit over other qubit implementations, before introducing a practical scheme for the experimental detection of weak charge modulation.

### A. Advantage of the heavy fluxonium over other capacitively shunted qubits

In this section, we aim to maximize the Rabi rate for a single-mode qubit subjected to a nearly resonant offset

charge of fixed oscillation amplitude $N_{drive}$ and frequency $\omega_d$. This thought experiment will provide a clearer understanding of why the heavy fluxonium holds an advantage over other capacitively shunted qubits.

Consider a single-mode qubit with a capacitive energy given by $4E_C[\hat{n} - n_g(t)]^2$, which interacts with a classical offset charge $n_g(t) = N_{drive} \cos(\omega_d t)$. For small charge modulations $N_{drive} \ll 1$, the Hamiltonian can be linearized. In a frame rotating at the drive frequency, it writes

$$ \hat{H}_{\text{int}} = -8N_{\text{drive}}E_C \langle e | \hat{n} | g \rangle \hat{\sigma}_x. \quad (4) $$

Using the relation between charge and flux matrix elements [22], namely, $8E_C|\langle e | \hat{n} | g \rangle| = \hbar\omega_{ge}|\langle e | \hat{\phi} | g \rangle|$, we derive the Rabi frequency:

$$ \Omega_r = 2N_{\text{drive}}\omega_{ge}|\langle e | \hat{\phi} | g \rangle|. \quad (5) $$

In a resonant coupling scenario, where the drive frequency $\omega_d$ is imposed by the resonance of an auxiliary system to probe, the qubit frequency needs to fulfill $\omega_{ge} = \omega_d$. In such a situation, maximizing the third factor $|\langle e | \hat{\phi} | g \rangle|$ is crucial. Indeed, only this term depends on the specifics of the qubit implementation, while the first two terms $N_{drive}$ and $\omega_{ge}$ are characteristics of the auxiliary system to be detected. For instance, in cQAD, the frequency $\omega_d$ is set by the mechanical resonance frequency, whereas the amplitude $N_{drive}$ depends on the details of the mechanical-electrical transduction. Consider the scenario of a silicon nitride membrane, which is a promising candidate for testing Penrose gravitational collapse due to its long coherence time and large zero-point fluctuations [21]. In this case, we expect an ac-charge modulation of $N_{drive} \sim 10^{-2}$ at a resonance frequency of $\omega_d/2\pi = \Omega_m/2\pi \approx 2$ MHz (see Appendix G).
While the matrix element $|\langle e|\hat{\phi}|g\rangle|$ is typically suppressed exponentially in the heavy-fluxonium regime, a radically different scenario emerges at the flux-frustration point. Here, the wave functions recover a large overlap $|\langle e|\hat{\phi}|g\rangle| \sim \pi$. This value compares favorably with weakly anharmonic qubits, where $|\langle e|\hat{\phi}|g\rangle| \sim (2E_C/E_J)^{1/4} \ll 1$, or even the Cooper pair box $|\langle e|\hat{\phi}|g\rangle| \sim 4E_C/E_J \sim 1$. In essence, the unique characteristics of fluxonium eigenstates at the flux-frustration point—manifesting as Schrödinger’s cat-like superpositions of persistent current states—endow it with a larger charge sensitivity compared to a transmon or Cooper pair box operating at the same transition frequency.

## B. Rabi oscillations of the qubit transition

In Fig. 4, we directly drive the qubit, biased at $\varphi_{\text{ext}} = \pi$, with a MHz pulse on the charge drive. We observe a Rabi oscillation pattern with maximum contrast for $\omega_d = \omega_{\text{ge}}$. The inset shows the Rabi frequency $\Omega_r/2\pi$ for a resonant drive at 1.8 MHz. As expected from Eq. (5), we observe a linear dependence of the Rabi oscillations with the drive amplitude, up to $\Omega_r/2\pi \sim 1$ MHz. For larger amplitude of the drive, the rotating-wave approximation breaks down as $\Omega_r$ approaches $\omega_{\text{ge}}$, leading to a deformed pattern with reduced contrast at the resonance drive condition. We use

<figure><img src="image_6.png" /><figcaption>FIG. 4. Direct Rabi manipulation of the radio-frequency qubit transition. After initial preparation in $|g\rangle$, the flux is reset to $\varphi_{\text{ext}} = \pi$, and the qubit is driven via the charge port with a MHz pulse of variable duration and frequency. The final qubit population is read out using the technique described in Fig. 2. The negative frequency part of the graph is here to highlight the validity range of the rotating-wave approximation. The inset shows the Rabi frequency for a resonant drive at 1.8 MHz, extracted from a sinusoidal fit, as a function of the drive voltage amplitude (upper horizontal axis). A linear fit of Eq. (5) to the data provides the lower horizontal axis calibration, where the drive amplitude is expressed in Cooper pairs on the fluxonium electrode. The red dot is obtained for the parameters of the main figure.</figcaption></figure>

Eq. (5) to relate the voltage amplitude on the digital-to-analog converter to the equivalent number of Cooper pairs $N_{\text{drive}}$ on the fluxonium electrode. We also deduce from this relation the minimum charge modulation $N_{\text{min}}$ required to observe coherent Rabi flipping, which is obtained for $\Omega_r \approx \pi/T_1$:

$$N_{\text{min}} = \frac{2\pi}{|\langle e|\hat{\phi}|g\rangle|\omega_{\text{ge}}T_1} \approx 5 \times 10^{-3}. \quad (6)$$

The ability to manipulate the qubit state with less than 1% of a Cooper pair shows the extreme sensitivity of the fluxonium to a resonant ac-charge modulation. For instance, this value would be sufficient to reach the strong-coupling regime with a dc-biased mechanical membrane in a resonant coupling scenario (see Appendix G).

The aforementioned value of $5 \times 10^{-3}$ Cooper pairs corresponds to a single-shot charge sensitivity of $10^{-2}$ e. However, through the implementation of quantum sensing protocols, like those routinely used in nitrogen-vacancy-center magnetometry [31] and similar methodologies [32], we are able to accrue substantial statistical data. This allows us to measure charge sensitivity within a 1 s integration period and subsequently compare these findings with other charge-sensing methods.

## C. Frequency-resolved ac-charge sensitivity

In a quantum sensing experiment, we can leverage the ability to swiftly prepare and read out the qubit state to detect a weak charge signal through repeated interaction with the two-level system. This involves preparing the qubit in $|g\rangle$, after which it interacts for an interrogation time $\tau_I$ with the weak continuous signal to be detected (referred to as the “calibration tone” henceforth), of frequency $\omega_{\text{cal}}$, applied to the charge port. For weak enough calibration tone, the Bloch vector undergoes a small rotation away from the south pole. We then probe this displacement by mapping the transverse component of the Bloch vector to the $\sigma_z$ basis with a $\pi/2$ pulse, before performing a single-shot readout of the qubit in the $|g\rangle, |e\rangle$ basis. In this scheme, the probability to detect the qubit in $|e\rangle$ slightly deviates from 1/2, by an amount that depends on the phase and amplitude of the calibration tone. Furthermore, the mismatch $\Delta$ between the calibration tone and qubit frequencies gives rise to a shot-to-shot rotation of the Bloch vector by an angle $\theta_k = k\Delta\tau$, where k is the repetition index and $\tau$ the repetition period of the experiment. Even though each measurement result $m_k \in \{0, 1\}$ only contains one bit of information, the complete measurement record $\{m_k\}_{0\le k < N_{\text{tot}}}$ can be used to reconstruct the spectrum of the charge modulation by the periodogram method [33].

Performing the $\pi/2$ rotation along a unique axis would lead to an ambiguity between positive and negative detuning $\Delta$. We thus perform the qubit rotations along an axis picked up sequentially in the set $(+X, +Y, -X, -Y)$.
<figure><img src="image_7.png" /><figcaption>FIG. 5. ac-charge sensing. (a) A weak monochromatic charge drive (also referred to as calibration tone) is detected thanks to a repeated pulse sequence: The qubit is prepared in $|g\rangle$ [black arrow in the Bloch sphere (b)]. After interacting for a time $\tau_I$ with the tone, a partial information on the qubit state is obtained by performing a $\pi/2$ pulse in one of the four directions $+X, +Y, -X, -Y$, followed by a qubit state readout in the eg basis. From the measurement samples $m_k \in \{0, 1\}$, a complex telegraphic signal $\sigma_k = i^k(m_k - 1/2)$ is constructed. The noise spectrum centered around the qubit frequency is estimated by the Bartlett's method, with periodograms of 1000 nonoverlapping consecutive samples. (c) The estimated noise spectrum presents a residual-bandwidth-limited peak at the calibration tone frequency $\nu_{cal} = 1.853$ MHz. Red inset: enlargement on the calibration peak and sinus-cardinal fit (solid line). Left- and right-hand insets: signal-to-noise-ratio (SNR) for the calibration peak as a function of interrogation time $\tau_I$ and calibration peak amplitude, respectively. The red dots in the insets correspond to the parameters used in the main graph of (c). The solid lines are the results of an analytic model taking into account the evolution of the qubit during the interrogation time. Signal cancellation occurs when the calibration tone amplitude is a multiple of that of a $\pi$ pulse. The spectrum in (c) is calibrated using the known variance of the calibration tone. RO stands for readout.</figcaption></figure>

This ensures a nonambiguous correspondence between discrete and continuous time frequencies over the interval $[-\Omega_{Ny}/2, +\Omega_{Ny}/2]$, where $\Omega_{Ny}$ is the Nyquist angular frequency $\Omega_{Ny} = \pi/\tau$ (see Appendix F). The charge-noise spectrum over this interval is then reconstructed by performing fast-Fourier transforms over adjacent windows of $N = 1000$ consecutive samples. Figure 5(c) shows an example of such an experimentally reconstructed spectrum. The calibration tone is visible as a sinus-cardinal-shaped peak, centered around $\omega_{cal}$ and of width $\Omega_{RBW} = 2\pi/N\tau$. This value is the residual bandwidth of our quantum spectrum analyzer, and it can be tuned by adjusting the window length $N$. The spectrum is normalized in units of elementary charge $e/\sqrt{Hz}$ using the known amplitude of the calibration tone, as determined from the linear fit of Fig. 4.

The calibration peak sits on a flat noise background, which is attributable to the sampling noise of the quantum sensor [34]. An analytic model for the signal-to-noise ratio (SNR) as a function of the experimental parameters has been derived (see Appendix F) and shows good agreement with the measured data [see Fig. 5(c)]. Qualitatively, the SNR increases linearly for $\tau_I \ll T_1$, as the initial Bloch vector accumulates a transverse component $2|\langle\sigma\rangle| = \Omega_r\tau_I$. On the other hand, due to the interaction with the thermal bath, the Bloch vector relaxes eventually toward the origin of the Bloch sphere such that the SNR vanishes for $\tau_I \gg T_1$. In practice, around the optimal value $\tau_I \sim 20$ µs, the detector achieves a noise level as low as 33 µe/√Hz. This value approaches that of the most sensitive electrometers such as the radio-frequency quantum point contact [13,35] or the radio-frequency
SET [11,12]. Yet, these transport-based sensors are very different in nature from the current qubit-based quantum protocol. The shunt capacitor on which the charge is detected in our system is typically 2 orders of magnitude larger than the superconducting islands employed in those systems [9,12,13]. This is of utmost practical importance when it comes to connecting the sensor to an auxiliary quantum system. As an example, when trying to detect the charge modulation of an electromechanical system such as Ref. [19], the 50 fF capacitor of the vacuum-gap system would perfectly match the value employed in this work, whereas traditional sensors would suffer a large dilution of the signal. The challenge of detecting extremely small charge signals while maintaining a large island capacitance is more directly captured by the energy sensitivity [12] $δq²/2C ≈ 2.8ħ$ which is below the sensitivity of any other charge detectors operating at MHz frequencies. Contrary to rf SETs, whose sensitivity is fundamentally bounded by the shot-noise limit $δq²/2C ≥ ħ$ [9], the quantum sensor introduced in this work could in principle operate below this bound, for instance, by improving the qubit coherence time as demonstrated in Appendix F. Furthermore, in stark contrast with transport-based measurements, featuring a flat frequency response from dc to several tens of MHz, our resonant detector features a narrow frequency response around the qubit frequency, the full bandwidth being given by $Ω_{full} = 2π/τ_I ∼ 50$ kHz (see Appendix F). This peculiar frequency response is highly advantageous when coupling the fluxonium to a nearly resonant system, as it guarantees perfect immunity to low-frequency environmental charge noise while maximizing charge sensitivity at the MHz region of interest.

## VII. CONCLUSION

In conclusion, we have demonstrated high-fidelity preparation, manipulation, and single-shot readout of a heavy-fluxonium qubit with a transition frequency as low as 1.8 MHz. To the best of our knowledge, this is the lowest frequency reported so far for a superconducting qubit. As demonstrated in earlier work [6], this circuit represents a realistic alternative to the transmon in a quantum computing architecture. Our work furthermore demonstrates the potential of this circuit in sensing experiments. This can be routed from the peculiar frequency response of the circuit which filters efficiently the environmental noise at audio frequency while being maximally sensitive at the resonant qubit frequency in the MHz range. The high charge sensitivity combined with the large capacitive shunt demonstrated in this work opens up avenues in hybrid circuits, where the fluxonium can be used as a resonant probe to manipulate other physical systems. As an example, we show (Appendix G) that the coherence time and electric dipole achieved in the current work are sufficient to attain the strong-coupling regime in a hybrid electromechanical system involving a dc-biased nanomechanical resonator.

## ACKNOWLEDGMENTS

We thank Natalia Pankratova for her valuable insights and advice regarding the fabrication of Josephson junctions. The authors acknowledge support from ANR project MecaFlux (ANR-21-CE47-0011), the CryoParis project from the Région Ile-de-France DIM Sirteq, and the HyQuTech project from the Emergence Sorbonne Université program. This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (Grant Agreement No. 851740). K. G. is funded by the Quantum Information Center Sorbonne (QICS) Ph.D. program, H. P. is funded via the CNRS/UArizona joint Ph.D. program. E. F. acknowledges funding from the European Research Council under Grant No. 101042315 (INGENIOUS). A. S. acknowledges support from ANR project Hamroqs and the Plan France 2030 through the PEPR NISQ2LSQ Project No. (ANR-22-PETQ-0006). B.-L. N.-S. acknowledges support from ANR QFilters (ANR-18-JSTQ-0002).

## APPENDIX A: MICROFABRICATION

The large circuit parts, i.e., the coplanar waveguide resonator, the flux line, the charge-drive electrode, and the fluxonium coplanar capacitor, were fabricated with standard UV laser lithography: starting with a 280-µm-thick silicon (100) wafer with a resistivity of 20 kΩ cm, and coated with a 150-nm-thick layer of niobium (Nb), we spin coat S1805 positive resist, and bake it at 115 °C for 1 min. The resist is then exposed to UV light with a dose of 100 mJ/cm², and developed using MF-319. Nb is etched using reactive ion etching with a $SF_6$ plasma. Any remaining resist is finally removed with acetone in an ultrasound bath at 50 °C for 15 min, rinsed in isopropanol (IPA), and dried.

The small Josephson junction and the numerous large junctions of the superinductor were fabricated using the Dolan-bridge technique. First, the sample is spin coated with MMA EL13 at 4000 rpm, and subsequently baked for 1 min at 195 °C. The sample is then spin coated with PMMA A3 at 5000 rpm and baked for 30 min at 195 °C. Electron-beam lithography with a dose of 280 µC cm⁻² is used to create the freestanding bridges on the MMA-PMMA bilayer. The development is performed with a mixture of IPA and deionized water at a temperature of 6 °C for 90 s. The Al/AlOₓ/Al junctions are fabricated by evaporation of a 31-nm-thick aluminum (Al) layer with a -22° angle, followed by oxidization with a mixture 9:1 of argon and O₂, at 200 mbar and for 12 min. Finally, a 100-nm-thick Al layer is evaporated at +22° angle. A lift-off is then performed in an N-methyl-2-pyrrolidone bath at 80 °C for 20 min. The sample is then dried after rinsing with acetone and IPA.
<figure><img src="image_8.png" /><figcaption>FIG. 6. Room-temperature rf and dc circuitry. IR mixer: Image rejection mixer; IF: intermediate frequency; RT amplifier: room temperature amplifier; and TWPA: traveling wave parametric amplifier.</figcaption></figure>

## APPENDIX B: EXPERIMENT SCHEMATIC

The room-temperature and cryogenic rf and dc connections are depicted in Figs. 6 and 7, respectively. A fast data acquisition system Operator X (OPX) from quantum machine is used in combination with a microwave source (Anapico APMS40G) to generate the rf and microwave pulses. The magnetic flux control is obtained by combining a stabilized voltage source (Yokogawa 7651) with a fast analog output of the OPX. This setup allows us to scan the flux bias over more than $\Phi_0$, while enabling fast control over a range of approximately $3 \times 10^{-2} \cdot \Phi_0$ with a sub-microsecond time resolution.

## APPENDIX C: SIDEBAND COOLING

The starting point of our analysis is the Hamiltonian of the fluxonium, including the driven readout cavity. The Hamiltonian is expressed in the normal mode basis, as obtained by diagonalizing the classical equations of motion for $E_J \to 0$ [36]. The normal modes are labeled “R” for readoutlike and “Q” for qubitlike. The Hamiltonian writes

<figure><img src="image_9.png" /><figcaption>FIG. 7. Cryogenic rf and dc circuitry. The color code matches that of Fig. 1. VLFX 500: 500 MHz low-pass filter from Minicircuits. K&L: 12 GHz low-pass filter (model 6L250-12000) from K&L. SLP 1.9+: 1.9 MHz low-pass filter from Minicircuits. JTWPA: Josephson traveling wave parametric amplifier.</figcaption></figure>

$$ \hat{H} = \hbar\omega_Q\hat{c}^\dagger\hat{c} + \hbar\omega_R\hat{a}^\dagger\hat{a} + \hbar(\hat{a}\epsilon_d^*e^{i\omega_p t} + \hat{a}^\dagger\epsilon_d e^{-i\omega_p t}) \quad (C1) $$

$$ -E_J \cos(\hat{\phi}_Q + \hat{\phi}_R - \varphi_{\text{ext}}), \quad (C2) $$

where $\hat{\phi}_R = \varphi_{zpf,R}(\hat{a} + \hat{a}^\dagger)$ and $\hat{\phi}_Q = \varphi_{zpf,Q}(\hat{c} + \hat{c}^\dagger)$ represent the normal mode positionlike operators in the absence of the Josephson term ($E_J = 0$). The bosonic annihilation operators for the qubit and cavity modes are denoted by $\hat{c}$ and $\hat{a}$, with respective frequencies $\omega_Q$ and $\omega_R$. The zero-point fluctuations of the readout and qubit modes, as seen by the Josephson junction, are given by $\varphi_{zpf,R}$ and $\varphi_{zpf,Q}$. Note that the participation ratio of the resonator mode in the Josephson junction is very small [36], such that $\varphi_{zpf,R} \ll 1$.

To account for photon loss in the readout resonator, we model the dynamics with a master equation,

$$ \frac{d}{dt}\hat{\rho} = -i[\hat{H}, \hat{\rho}] + \kappa\mathcal{D}_{\hat{a}}(\hat{\rho}), \quad (C3) $$
where the Lindbladian writes $D_{\hat{a}}(\hat{\rho}) = \hat{a}\hat{\rho}\hat{a}^{\dagger} - \frac{1}{2}[\hat{a}^{\dagger}\hat{a}\hat{\rho} + \hat{\rho}\hat{a}^{\dagger}\hat{a}]$.

In the following, we demonstrate how Eq. (C3) simplifies to the effective qubit dissipative dynamics, as represented by the loss operator in Eq. (3). We proceed by going in a frame rotating at the drive frequency $\omega_p$, and displaced around the mean amplitude of the cavity field $\alpha$, by doing the substitution $\hat{a} \rightarrow (\alpha + \hat{a})e^{-i\omega_p t}$. The steady state value $\alpha$ is chosen such that it cancels the three drift terms:

$$ \hat{H}_1 = (\hbar\kappa/2i)(\alpha^*\hat{a} - \alpha\hat{a}^\dagger), \quad (C4) $$

$$ \hat{H}_2 = \hbar\Delta(\alpha\hat{a}^\dagger + \alpha^*\hat{a}), \quad (C5) $$

$$ \hat{H}_3 = \hbar(\epsilon_d^* \hat{a} + \epsilon_d \hat{a}^\dagger), \quad (C6) $$

where $\hat{H}_1$ is an effective Hamiltonian dynamics stemming from the expression of the Lindbladian in the displaced frame, $\hat{H}_2$ comes from the linearization of the term $\hbar\Delta_R \hat{a}^\dagger \hat{a}$ in the rotating frame Hamiltonian, with $\Delta_R = \omega_p - \omega_R$ being the drive detuning, and $\hat{H}_3$ is the drive term. We thus obtain the value of $\alpha$:

$$ \alpha = \frac{-\epsilon_d}{\Delta_R + i\kappa/2}. \quad (C7) $$

The Hamiltonian in the displaced frame becomes

$$ \hat{H} = \hbar\omega_Q \hat{c}^\dagger \hat{c} + \hbar\Delta_R \hat{a}^\dagger \hat{a} - E_J \cos(\hat{\phi}_Q + \tilde{\phi}_R - \varphi_{\text{ext}}), \quad (C8) $$

where $\tilde{\phi}_R = \varphi_{zpf,R}(\alpha + \hat{a})e^{-i\omega_p t} + \varphi_{zpf,R}(\alpha^* + \hat{a}^\dagger)e^{i\omega_p t}$ is the resonator coordinate in the new frame. Since $\varphi_{zpf,R}|\alpha| \ll 1$, we Taylor expand this expression to second order with respect to $\tilde{\phi}_R$:

$$ \begin{aligned} \hat{H} = \hbar\omega_Q \hat{c}^\dagger \hat{c} + \hbar\Delta_R \hat{a}^\dagger \hat{a} - E_J \cos(\hat{\phi}_Q - \varphi_{\text{ext}}) & \\ + E_J \sin(\hat{\phi}_Q - \varphi_{\text{ext}}) \times \tilde{\phi}_R & \\ + \frac{E_J}{2} \cos(\hat{\phi}_Q - \varphi_{\text{ext}}) \times \tilde{\phi}_R^2 & \end{aligned} \quad (C9) $$

The first line of Eq. (C9) corresponds to the resonator and unperturbed qubit Hamiltonian, the second line, which corresponds to the first order Taylor expansion, can be safely neglected as it only consists of terms rotating at $\pm\omega_p$. On the other hand, the third line, corresponding to the second order Taylor expansion, reduces to $\cos(\hat{\phi}_Q - \varphi_{\text{ext}}) \times \varphi_{zpf,R}^2(\alpha\hat{a}^\dagger + \alpha^*\hat{a})$, once fast rotating terms have been neglected, and linearizing for $\alpha \gg \hat{a}$, $\hat{a}^\dagger$. We thus obtain Eq. (2), with

$$ \begin{aligned} \hat{H}_Q = \hbar\omega_Q \hat{c}^\dagger \hat{c} & \\ - E_J \left[ 1 - \left( \frac{2\pi\Phi_{zpf,R}}{\Phi_0} \right)^2 \left( |\alpha|^2 + \frac{1}{2} \right) \right] \cos(\hat{\phi}_Q - \varphi_{\text{ext}}) & \\ \approx \hbar\omega_Q \hat{c}^\dagger \hat{c} - E_J \cos(\hat{\phi} - \varphi_{\text{ext}}). & \end{aligned} \quad (C10) $$

In the last expression, we have used $\varphi_{zpf,R}|\alpha| \ll 1$, and $\hat{\phi}_Q \approx \hat{\phi}$ since the participation ratio of the readout resonator in the junction is small.

Let us now project the Hamiltonian on the qubit subspace, with the projector $\hat{\Pi}_{eg} = |e\rangle\langle e| + |g\rangle\langle g|$:

$$ \begin{aligned} \hat{H}_m = \hat{\Pi}_{eg}\hat{H}\hat{\Pi}_{eg} = \frac{\hbar\omega_{ge}}{2}\hat{\sigma}_z + \hbar\Delta_R \hat{a}^\dagger \hat{a} & \\ + \hbar g(c_0 + c_x\hat{\sigma}_x + c_z\hat{\sigma}_z) \times (\alpha\hat{a}^\dagger + \alpha^*\hat{a}) & \end{aligned} \quad (C11) $$

with

$$ c_0 = (\langle g|\hat{\beta}|g\rangle + \langle e|\hat{\beta}|e\rangle)/2, $$

$$ c_x = \langle g|\hat{\beta}|e\rangle, $$

$$ c_z = (\langle e|\hat{\beta}|e\rangle - \langle g|\hat{\beta}|g\rangle)/2, $$

where $\hat{\beta} \equiv \cos(\hat{\phi}_Q - \varphi_{\text{ext}})$. The interesting processes occur when the cavity drive is nearly resonant with one of the two sidebands, $\Delta_R \sim \pm\omega_{ge}$. We treat separately the two cases by going to the interaction picture with respect to $\hat{H}_0^\pm = \hbar\omega_{ge}\hat{\sigma}_z/2 \pm \hbar\omega_{ge}\hat{a}^\dagger \hat{a}$:

$$ \begin{aligned} \hat{H}_m^\pm = \hbar\Delta_R^\pm \hat{a}^\dagger \hat{a} & \\ + \hbar g(c_0 + c_x[\hat{\sigma}^- e^{\mp i\omega_{ge} t} + \hat{\sigma}^+ e^{\pm i\omega_{ge} t}) + c_z\hat{\sigma}_z] & \\ \times (\alpha\hat{a}^\dagger e^{\pm i\omega_{ge} t} + \alpha^*\hat{a} e^{\mp i\omega_{ge} t}), & \end{aligned} \quad (C12) $$

where $\hat{H}_m^\pm (\hat{H}_m^-)$ is the Hamiltonian in the rotating frame $\hat{H}_0^\pm (\hat{H}_0^-)$, and $\Delta_R^\pm = \Delta_R \pm \omega_{ge}$ is the drive detuning with respect to the upper or lower sideband. Since the qubit-cavity system operates deep in the resolved sideband regime (at the bias point chosen for sideband preparation, $2\omega_{ge}/\kappa \approx 10$), we can safely neglect fast rotating terms, which yields:

$$ \hat{H}_m^\pm \approx \hbar\Delta_R^\pm \hat{a}^\dagger \hat{a} + \hbar g c_x (\alpha \hat{\sigma}^\pm \hat{a}^\dagger + \alpha^* \sigma^\mp \hat{a}). \quad (C13) $$

We proceed with the adiabatic elimination of the cavity field [37], since the cavity dissipation $\kappa D_{\hat{a}}(\hat{\rho})$ dominates over the coupling $g c_x |\alpha|$. We define the parameter $\epsilon \equiv g c_x |\alpha|/\kappa$, with respect to which we can expand the density matrix $\hat{\rho}$, with $\hat{\rho}_{mn} = \langle m|\hat{\rho}|n\rangle$ acting on the qubit's subspace:

$$ \begin{aligned} \frac{1}{\kappa} \frac{d\hat{\rho}}{dt} = & -\frac{i}{\hbar\kappa} [\hat{H}_m^\pm, \hat{\rho}] + D_{\hat{a}}(\hat{\rho}), \\ \hat{\rho} = & \hat{\rho}_{00}|0\rangle\langle 0| + \epsilon(\hat{\rho}_{10}|1\rangle\langle 0| + \rho_{01}|0\rangle\langle 1|) \\ & + \epsilon^2(\hat{\rho}_{11}|1\rangle\langle 1| + \hat{\rho}_{02}|0\rangle\langle 2| + \hat{\rho}_{20}|2\rangle\langle 0|) + O(\epsilon^3). \end{aligned} \quad (C14) $$
The goal of the adiabatic elimination procedure is to obtain the reduced dynamics of the qubit alone,

$$ \frac{d\hat{\rho}_Q}{dt} = \mathrm{Tr}_R \left[ \frac{d\hat{\rho}}{dt} \right] = \frac{d\hat{\rho}_{00}}{dt} + e^2 \frac{d\hat{\rho}_{11}}{dt} + O(e^3), \quad (C15) $$

which can be done by projecting the Lindblad evolution of Eq. (C14) on the resonator's elements $|0\rangle\langle0|$, $|0\rangle\langle1|$, $|1\rangle\langle1|$, so that

$$ \frac{1}{\kappa} \frac{d\hat{\rho}_{00}}{dt} = -i e^2 (\hat{\sigma}^{\pm} \hat{\rho}_{10} - \hat{\rho}_{01} \hat{\sigma}^{\mp}) + e^2 \hat{\rho}_{11} + O(e^3), \quad (C16) $$

$$ \frac{1}{\kappa} \frac{d\hat{\rho}_{10}}{dt} = -i\hat{\sigma}^{\mp}\hat{\rho}_{00} - \left(i\frac{\Delta_R^{\pm}}{\kappa} - \frac{1}{2}\right)\hat{\rho}_{10} + O(\epsilon), \quad (C17) $$

$$ \frac{1}{\kappa} \frac{d\hat{\rho}_{11}}{dt} = -i(\hat{\sigma}^{\mp}\hat{\rho}_{01} - \hat{\rho}_{10}\hat{\sigma}^{\pm}) - \hat{\rho}_{11} + O(\epsilon). \quad (C18) $$

Equation (C16) shows that $\hat{\rho}_{00}$ is slowly varying, since its derivative is of order $O(e^2)$. In the rhs of Eq. (C17), the first term is a source term, and the second one a damping term. As the source term is slowly varying, we can assume that $\hat{\rho}_{10}$ is always in its stationary state: $(1/\kappa)(d\hat{\rho}_{10}/dt) \approx 0$. The same argument applies to Eq. (C18), such that $(1/\kappa)(d\hat{\rho}_{11}/dt) \approx 0$. In the end, we obtain

$$ \hat{\rho}_{10} = -\frac{\kappa}{\Delta_R^{\pm} + i\frac{\kappa}{2}} \hat{\sigma}^{\mp} \hat{\rho}_{00} + O(\epsilon), \quad (C19) $$

$$ \hat{\rho}_{11} = \frac{\kappa^2}{\Delta_R^{\pm 2} + \frac{\kappa^2}{4}} \hat{\sigma}^{\mp} \hat{\rho}_{00} \hat{\sigma}^{\pm} + O(\epsilon). \quad (C20) $$

By inserting these expressions in Eq. (C16), we recognize the Lindblad evolution associated to the following effective loss operators on the qubit:

$$ \hat{L}^{\pm} = gc_x |\alpha| \sqrt{\frac{\kappa}{\Delta_R^{\mp 2} + \frac{\kappa^2}{4}}} \hat{\sigma}^{\pm}. \quad (C21) $$

When the drive is set at resonance with one of the sidebands ($\Delta_R^{\pm} = 0$), we retrieve the expressions Eq. (3).

Following a more rigorous treatment of the adiabatic elimination using a systematic pertubative approach, as described in Ref. [38], leads to the same results.

## APPENDIX D: PREPARATION FIDELITY ESTIMATION

To evaluate the preparation fidelity, we prepare the qubit in either $|g\rangle$, $|e\rangle$, or a thermal state. The density matrix reduced to the $ge$ manifold reads

$$ \rho_0^{\mathrm{prep} j} = P_g^{\mathrm{prep} j} |g\rangle\langle g| + (1 - P_g^{\mathrm{prep} j}) |e\rangle\langle e|, $$

where $P_g^{\mathrm{prep} j}$ is the probability to be in the $|g\rangle$ state after preparation in state $j \in \{g, e, h\}$. Since the dispersive shift of the readout cavity is too small to directly distinguish $|g\rangle$ from $|e\rangle$, we apply a 64 ns pulse, resonant with the $|g\rangle \to |h\rangle$ transition. The initial $|g\rangle$ population undergoes Rabi oscillations of angle $\theta$, while the $|e\rangle$ population remains unaffected:

$$ \begin{aligned} \rho^{\mathrm{prep} j}(\theta) = & P_g^{\mathrm{prep} j} \frac{1 + \cos \theta}{2} |g\rangle\langle g| + P_g^{\mathrm{prep} j} \frac{1 - \cos \theta}{2} |h\rangle\langle h| \\ & + (1 - P_g^{\mathrm{prep} j}) |e\rangle\langle e| \\ & + P_g^{\mathrm{prep} j} \frac{|\sin \theta|}{2} (|g\rangle\langle h| + |h\rangle\langle g|). \end{aligned} \quad (D1) $$

Decoherence has been neglected in this process as the pulse duration is short compared to the decoherence rates of the $|g\rangle-|h\rangle$ transition. We then read out the state of the qubit through the dispersive shift of the readout resonator. Histograms of the real ($I$) and imaginary part ($Q$) of the reflection coefficient, measured with a 600 ns pulse are plotted in Fig. 8(a). The continuous variable $I$ is then compared to a threshold $I_{\mathrm{threshold}} = 0$ to yield a Boolean detection result left $I < I_{\mathrm{threshold}}$. By averaging a large number of repetitions, we measured the probability $P^{\mathrm{prep} j}[\mathrm{left}](\theta)$ to obtain $I < I_{\mathrm{threshold}}$ after a Rabi pulse of angle $\theta$ [see Fig. 8(b)], for each preparation protocol. This probability is given by

$$ P^{\mathrm{prep} j}[\mathrm{left}](\theta) = \sum_{x \in \{g, e, h\}} P[\mathrm{left}|x] \langle x | \rho^{\mathrm{prep} j}(\theta) | x \rangle, \quad (D2) $$

where $P[\mathrm{left}|g]$, $P[\mathrm{left}|e]$, and $P[\mathrm{left}|h]$ are the conditional probabilities of measuring $I < I_{\mathrm{threshold}}$ knowing that the qubit was in the state $|g\rangle$, $|e\rangle$, or $|h\rangle$, respectively. In a perfect detection scenario, $P[\mathrm{left}|g] = P[\mathrm{left}|e] = 1$ and $P[\mathrm{left}|h] = 0$. Combining Eqs. (D1) and (D3), we arrive at

$$ \begin{aligned} P^{\mathrm{prep} j}[\mathrm{left}](\theta) = & P_g^{\mathrm{prep} j} \frac{P[\mathrm{left}|g] + P[\mathrm{left}|h]}{2} \\ & + (1 - P_g^{\mathrm{prep} j}) P[\mathrm{left}|e] \\ & + P_g^{\mathrm{prep} j} \cos(\theta) \frac{P[\mathrm{left}|g] - P[\mathrm{left}|h]}{2}. \end{aligned} \quad (D3) $$

Considering the large occupation of the thermal bath, we assume equal populations in $|g\rangle$ and $|e\rangle$: $P_g^{\mathrm{prep} \mathrm{th}} = 0.5$. We proceed by fitting the three curves of Fig. 8 with the free parameters $P[\mathrm{left}|g]$, $P[\mathrm{left}|e]$, $P[\mathrm{left}|h]$, $P_g^{\mathrm{prep} g}$ and $P_g^{\mathrm{prep} e}$. We obtain the conditional readout probabilities $P[\mathrm{left}|g] = 94.04 \pm 0.04\%$, $P[\mathrm{left}|e] = 95.87 \pm 0.03\%$, and $P[\mathrm{left}|h] = 10.99 \pm 0.05\%$. The values of $1 - P[\mathrm{left}|g] = 6\%$, $1 - P[\mathrm{left}|e] = 4\%$ correspond to the mislabeling, due to the overlap of the Gaussian distributions. The larger value of $P[\mathrm{left}|h] = 11\%$ is due to the decay of the $|h\rangle$ state (of lifetime 7 µs) during the 600 ns readout pulse. The
<figure><img src="image_10.png" /><figcaption>FIG. 8. (a) Histograms of the real (I) and imaginary (Q) parts of the cavity reflection coefficient, as obtained by demodulating and integrating a 600 ns pulse. The x and y axes have been rescaled by the standard deviation of the Gaussian envelope. Top plot is the histogram obtained when the system is initially prepared in its thermal state. Even though $|g\rangle$ and $|e\rangle$ are almost equally populated, the dispersive shift is insufficient to separate the two distributions. Bottom plot is the histogram obtained after a $|g\rangle \rightarrow |h\rangle$ $\pi$ pulse. The left-hand blob corresponds to the $|e\rangle$ population, unaffected by the $\pi$ pulse, while the right-hand blob is due to the population transferred in $|h\rangle$. Setting a threshold at $I = 0$ (white dashed line) implements a single-shot readout. (b) Calibration of the state-preparation fidelity. The qubit is prepared through sideband cooling in $|e\rangle$ (yellow), $|g\rangle$ (blue), or a thermal state (red). We apply a 64 ns pulse at the $|g\rangle-|h\rangle$ transition frequency with a varying amplitude. The dynamics is then fitted using Eq. (D3) to extract the preparation fidelity (see main text for details).</figcaption></figure>

extracted preparation fidelity for the $|g\rangle$ and $|e\rangle$ states are $P_g^{\text{prep} g} = 97.67 \pm 0.05\%$ and $1 - P_g^{\text{prep} e} = 97.69 \pm 0.08\%$. The quoted error intervals are obtained by a bootstrap technique: The fit is repeated on a subset of the data obtained by random sampling with replacement of the data point, from which the mean value and standard deviation of each parameter is extracted. The effective temperature of the qubit after the preparation is

$$T_g^{\text{prep} g} = \frac{\hbar\omega_{ge}}{k_B[\log(P_g^{\text{prep} g}) - \log(1 - P_g^{\text{prep} g})]} = 23 \text{ µK.}$$

## APPENDIX E: NOISE TEMPERATURE OF THE QUBIT ENVIRONMENT

In this appendix, we determine how the relaxation rate varies with temperature. To achieve this, we heat the mixing chamber of the cryostat with a resistor. The temperature $T_{RuO_2}$, as measured by a ruthenium oxide probe built in with the cryostat (model Bluefors BF-LD250), is stabilized thanks to a feedback loop to various setpoints ranging from 7 to 100 mK. For each point, we measure the decay rate $2\Gamma$ of the states $|e\rangle$ and $|g\rangle$, akin to the measurement presented in the main text [refer to Fig. 3(a)].

To precisely correlate our data with the circuit's temperature, we harness the residual thermal population of the higher qubit excited states, $\{|f\rangle, |h\rangle\}$. This functions as a local probe, sensitive to the noise temperature which directly reflects the noise spectrum experienced

by the circuit near its second transition at a frequency of 3.7 GHz [39].

In practice, we let the circuit thermalize with its environment, and then record a histogram of the real (I) and imaginary part (Q) of the readout cavity reflection coefficient, as visible in Fig. 9(b). Three peaks are visible on the histogram, corresponding to the population in the manifold $\{|g\rangle, |e\rangle\}$, the state $|f\rangle$, and the state $|h\rangle$, respectively. We assume a Boltzmann distribution for the population in the various qubit states: $p_k \propto e^{-E_k/k_B T}$, where $E_k$ is the energy of state k ($k \in \{|g\rangle, |e\rangle, |f\rangle, |h\rangle\}$. Furthermore, by neglecting the small transition frequencies $\omega_{ge} / \sim 2\pi \times 1.8$ MHz and $\omega_{fh} / \sim 2\pi \times 50$ MHz, compared to $\omega_{ge-fh} \equiv (\omega_{ef} + \omega_{gh})/2 \sim 2\pi \times 3.7$ GHz, we get $p_{|g\rangle} = p_{|e\rangle} = p_{g,e}/2$ and $p_{|f\rangle} = p_{|h\rangle} = p_{f,h}/2$. We extract the populations $p_{g,e}$ and $p_{f,h}$ by a triple Gaussian fit to the readout histogram, where the Gaussian peaks corresponding to $|f\rangle$ and $|h\rangle$ are constrained to the same area. From the values $p_{g,e}$ and $p_{f,h}$, we determine the effective circuit temperature:

$$T_{\text{circ}} \equiv \frac{\hbar\omega_{ge-fh}}{k_B \log\left(\frac{p_{g,e}}{p_{f,h}}\right)} \quad (\text{E1})$$

Figure 9(c) presents the evolution of $T_{\text{circ}}$ with respect to $T_{RuO_2}$. Throughout the entire temperature range, there is a notable disparity between the two temperatures. Specifically, the circuit temperature remains relatively stable at $T_{\text{circ}} \sim 59$ mK for $7 \le T_{RuO_2} \le 40$ mK, before rising to
<figure><img src="image_11.png" /><figcaption>FIG. 9. Temperature dependence of the decoherence rate. (a) Two-dimensional histogram of the I, Q quadratures of the readout reflection coefficients for a qubit at thermal equilibrium with the environment. (b) Histogram of the I quadrature. The population in the manifolds {|g⟩, |e⟩⟩, |f⟩, |h⟩⟩} are determined by fitting the various peaks with Gaussian functions (see text for details). Qubit effective temperature $T_{\text{circ}}$, as determined from Eq. (E1), as a function of cryostat temperature $T_{\text{RuO}_2}$. (d) Energy decay rate $2\Gamma$ in the {|g⟩, |e⟩⟩} manifold measured via $T_1$ relaxometry [see Fig. 3(a)], as a function of the effective temperature $T_{\text{circ}}$. The dashed line is a linear fit to the data yielding $2\Gamma = \xi \cdot T_{\text{circ}} + \zeta$, with $\xi = 0.27 \text{ mK}^{-1} \text{ms}^{-1}$ and $\zeta = 16 \text{ ms}^{-1}$. The filled point in (c) and (d) is the one extracted from the histograms in (a) and (b).</figcaption></figure>

$T_{\text{circ}} = 318 \text{ mK}$ when $T_{\text{RuO}_2} = 100 \text{ mK}$. Such discrepancies could arise from inadequate sample thermalization, temperature gradients within the mixing chamber plate, or potential miscalibration of the RuO₂ temperature sensor. This underscores the value of our independent temperature assessment using the higher qubit excited states.

We then plot the decay rate $2\Gamma$ as a function of $T_{\text{circ}}$ in Fig. 9(d). The data are well fitted by an affine relation $2\Gamma = \xi T_{\text{circ}} + \zeta$, with $\xi = 0.27 \text{ mK}^{-1} \text{ms}^{-1}$ and $\zeta = 16 \text{ ms}^{-1}$. The first term, proportional to $T_{\text{circ}}$, can be interpreted as the thermally induced relaxation, occurring at a rate $\gamma(2n_{\text{th}} + 1) \sim 2\gamma n_{\text{th}}$, where the second equality results from $n_{\text{th}} \gg 1$ for the 1.8 MHz qubit transition. The second term, which can be traced back to a technical noise independent of $T_{\text{circ}}$, has a relatively small contribution over the studied temperature range—$\zeta \sim \xi T_{\text{circ}}$ for the lowest temperature $T_{\text{circ}} \sim 59 \text{ mK}$. This points toward the fact that technical noise is only marginally affecting the qubit relaxation, in spite of the very low transition frequency $\omega_{ge}/2\pi = 1.8 \text{ MHz}$. Finally, from the linear dependence, we can extract the intrinsic circuit damping $\gamma = \xi\hbar\omega_{ge}/2k_B \sim 12 \text{ s}^{-1}$, or equivalently, $Q_{\text{intr}} = \omega_{ge}/\gamma = 9.6 \times 10^5$.

## APPENDIX F: CHARGE SPECTRUM ANALYZER

In this appendix, we develop a theoretical model for the expected signal-to-noise ratio in the frequency-resolved charge detection experiment.

### 1. Qubit evolution during the interrogation time

We first model the evolution of the qubit during the interrogation time, by taking into account the interaction with the calibration tone (Rabi frequency $\Omega_r$, finite detuning of the calibration tone $\Delta = \omega_{ge} - \omega_{cal}$). Since the qubit is coupled to a thermal bath with a large occupation, we choose an equal rate $\Gamma$ for the loss and gain of qubit excitations. From the empirical finding $T_1 \approx T_2$ (see Fig. 3), we also assume a dephasing rate $\Gamma_{\phi} \approx \Gamma/2$. The full evolution of the qubit's density matrix $\rho$ is thus, in a frame rotating at the drive frequency,

$$ \frac{d\hat{\rho}}{dt} = -\frac{i}{\hbar}[\hat{H}, \hat{\rho}] + \Gamma \left[ \hat{\sigma}\rho\hat{\sigma}^{\dagger} - \frac{1}{2}(\hat{\sigma}^{\dagger}\hat{\sigma}\hat{\rho} + \hat{\rho}\hat{\sigma}^{\dagger}\hat{\sigma}) \right] \\ + \Gamma \left[ \hat{\sigma}^{\dagger}\rho\hat{\sigma} - \frac{1}{2}(\hat{\sigma}\hat{\sigma}^{\dagger}\hat{\rho} + \hat{\rho}\hat{\sigma}\hat{\sigma}^{\dagger}) \right] \\ + \frac{\Gamma}{2} \left[ \hat{\sigma}_{z}^{\dagger}\rho\hat{\sigma}_{z} - \frac{1}{2}(\hat{\sigma}_{z}\hat{\sigma}_{z}^{\dagger}\hat{\rho} + \hat{\rho}\hat{\sigma}_{z}\hat{\sigma}_{z}^{\dagger}) \right], \quad (\text{F1}) $$

with $\hat{H} = \hbar\Delta\hat{\sigma}_z/2 + \hbar\Omega_r\hat{\sigma}_x/2$. We proceed by calculating the Bloch equations for the three components of the qubit pseudo-spin:

$$ \frac{d\langle\hat{\sigma}_x\rangle}{dt} = -\Delta\langle\hat{\sigma}_y\rangle - 2\Gamma\langle\hat{\sigma}_x\rangle, \quad (\text{F2}) $$
$$ \frac{d\langle\hat{\sigma}_y\rangle}{dt} = -\Omega_r\langle\hat{\sigma}_z\rangle + \Delta\langle\hat{\sigma}_x\rangle - 2\Gamma\langle\hat{\sigma}_y\rangle, \quad (F3) $$

$$ \frac{d\langle\hat{\sigma}_z\rangle}{dt} = \Omega_r\langle\hat{\sigma}_y\rangle - 2\Gamma\langle\hat{\sigma}_z\rangle. \quad (F4) $$

These equations describe a rotation around an axis $\Delta\mathbf{e}_z + \Omega_r\mathbf{e}_x$ combined with an isotropic relaxation toward the origin of the Bloch sphere at a rate $\Gamma$ due to the various relaxation channels. We solve for a qubit initially prepared in $|g\rangle$ ($\langle\hat{\sigma}_z\rangle = -1$, $\langle\hat{\sigma}_x\rangle = \langle\hat{\sigma}_y\rangle = 0$), and obtain

$$ \langle\hat{\sigma}_x\rangle = e^{-2\Gamma t} \left(\cos\left(\sqrt{\Omega_r^2 + \Delta^2} t\right) - 1\right) \frac{\Delta\Omega_r}{\Omega_r^2 + \Delta^2}, $$

$$ \langle\hat{\sigma}_y\rangle = e^{-2\Gamma t} \sin\left(\sqrt{\Omega_r^2 + \Delta^2} t\right) \frac{\Omega_r}{\sqrt{\Omega_r^2 + \Delta^2}}, $$

$$ \langle\hat{\sigma}_z\rangle = -e^{-2\Gamma t} \left(\Delta^2 + \Omega_r^2 \cos\left(\sqrt{\Omega_r^2 + \Delta^2} t\right)\right) \frac{1}{\Omega_r^2 + \Delta^2}. $$

At the end of the interrogation time, we can thus obtain the magnitude of the pseudo-spin projection in the x, y plane:

$$ |2\langle\hat{\sigma}_0\rangle| = |\langle\hat{\sigma}_x\rangle + i\langle\hat{\sigma}_y\rangle| = \Omega_r\tau_I e^{-2\Gamma\tau_I} f(\Delta), \quad (F5) $$

where $f(\Delta)$ is the frequency response function of the detector, given by

$$ f(\Delta) = \sqrt{\frac{\Delta^2 \mathrm{sinc}^2\left(\sqrt{\Omega_r^2 + \Delta^2 \frac{\tau_I}{2\pi}}\right) + \Omega_r^2 \mathrm{sinc}^2\left(\sqrt{\Omega_r^2 + \Delta^2 \frac{\tau_I}{\pi}}\right)}{\Omega_r^2 + \Delta^2}} \\ \approx \mathrm{sinc}\left(\frac{\Delta}{\Omega_{\mathrm{full}}}\right), \quad (F6) $$

where $\Omega_{\mathrm{full}} = (2\pi/\tau_I)$ and the convention $\mathrm{sinc}(x) = \sin(\pi x)/\pi x$ has been used. The second equation is valid in the limit $\Omega_r\tau_I \ll 1$.

## 2. Signal processing

At the end of the interrogation time, a projective measurement of one of the transverse components of the pseudo-spin is performed in the qubit frame:

$$ \langle m_k \rangle = 1/2 + \mathrm{Re}[\langle\sigma\rangle_0 (-i)^k e^{i\Delta k \tau}]. \quad (F7) $$

The term $(-i)^k$ encodes for the alternating measurement basis $\{X, Y, -X, -Y\}$. The term $e^{i\Delta k\tau}$ describes the phase difference between the frames of the qubit and calibration tone. Without loss of generality, we can ignore the phase of $\langle\sigma\rangle_0$ and assume $\langle\sigma\rangle_0 \in \mathbb{R}^+$, such that

$$ \langle m_k \rangle = 1/2 + \langle\sigma\rangle_0 \mathrm{Re}[-(-i)^k e^{i\Delta k \tau}]. \quad (F8) $$

Finally, samples undergo the transformation $\sigma_k = i^k(m_k - 1/2)$. We thus get

$$ \langle\sigma_k\rangle = \begin{cases} \langle\sigma\rangle_0 \cos(\Delta k \tau) & k \text{ even} \\ i\langle\sigma\rangle_0 \sin(\Delta k \tau) & k \text{ odd.} \end{cases} \quad (F9) $$

Hence, the real and imaginary parts of the complex values $\langle\sigma\rangle_0 e^{i\Delta t}$ are encoded pairwise on the successive samples $\sigma_k$. The records are then grouped by windows of $N = 1000$ consecutive samples, and Fourier transformed to yield periodograms. In order to reduce the spacing between adjacent frequency bins, we perform the Fourier transform on a zero-padded version of the samples $\{z_k\}_{0\le k<N_pN}$, with

$$ z_k = \begin{cases} \sigma_k & 0 \le k \le N-1 \\ 0 & N \le k < N_pN. \end{cases} \quad (F10) $$

The padding factor $N_p$ represents the number of frequency bins in each measurement bandwidths. We typically use $N_p = 5$ in our data analysis. We denote $\{Z_k\}_{0\le k<N_pN}$ the Fourier transform of the samples $\{z_n\}_{0\le n<N_pN}$:

$$ Z_n = \sum_{k=0}^{N_p N-1} z_k e^{-2i\pi kn/N_p N}. \quad (F11) $$

Following Bartlett's method, the spectrum is then estimated by taking the mean value $S_n = \langle|Z_n|^2\rangle$ over a large number of periodograms.

## 3. Response to the calibration tone and frequency aliasing

Because of the calibration tone, the samples $z_k$ have a nonzero expectation value [see Eq. (F9)]. We now estimate the line shape $S_n^{\mathrm{signal}} = |\langle Z_n \rangle|^2$ resulting from this signal. By combining Eq. (F9) with Eq. (F11), and separating the contribution of even and odd index $k$ in the sum, we get

$$ \langle Z_n \rangle = \langle \sigma \rangle_0 \sum_{k=0}^{N/2-1} \cos \left( 2\pi \frac{2k\Delta}{2\Omega_{Ny}} \right) e^{-i2\pi(2k\Delta_n/2\Omega_{Ny})} \\ + i \sin \left( 2\pi \frac{(2k+1)\Delta}{2\Omega_{Ny}} \right) e^{-i2\pi[(2k+1)\Delta_n/2\Omega_{Ny}]}, $$

where the $n$th frequency bin is given by $\Delta_n = (2\pi n/\tau N_p N)$ and the Nyquist frequency $\Omega_{Ny} = \pi/\tau$. After elementary arithmetic manipulations, we arrive at

$$ \langle Z_n \rangle = \frac{\langle \sigma \rangle_0}{2} \left[ e^{i2\pi[(\Delta - \Delta_n)/\Omega_{Ny}][(N-1)/4]} \frac{\sin(2\pi \frac{\Delta - \Delta_n}{\Omega_{Ny}} \frac{N}{4})}{\sin(2\pi \frac{\Delta - \Delta_n}{4\Omega_{Ny}})} \right. \\ \left. + i e^{-i2\pi[(\Delta + \Delta_n)/\Omega_{Ny}][(N-1)/4]} \frac{\sin(2\pi \frac{\Delta + \Delta_n}{\Omega_{Ny}} \frac{N}{4})}{\cos(2\pi \frac{\Delta + \Delta_n}{4\Omega_{Ny}})} \right]. $$
For large values of $N$, we have

$$S_n^{\text{signal}} \approx \frac{\langle \sigma \rangle_0^2}{4} \left[ \left( \frac{\sin(2\pi \frac{\Delta - \Delta_n}{\Omega_{\text{Ny}}})}{\sin(2\pi \frac{\Delta - \Delta_n}{4\Omega_{\text{Ny}}})} \right)^2 + \left( \frac{\sin(2\pi \frac{\Delta + \Delta_n}{\Omega_{\text{Ny}}})}{\cos(2\pi \frac{\Delta + \Delta_n}{4\Omega_{\text{Ny}}})} \right)^2 \right].$$

This expression is peaked around the values $\Delta \equiv \Delta_n (\text{mod } 2\Omega_{\text{Ny}})$ and $\Delta \equiv \Omega_{\text{Ny}} - \Delta_n (\text{mod } 2\Omega_{\text{Ny}})$. Figure 10 shows the two families of peaks in the $(\Delta, \Delta_n)$ plane. In the aliasing-free region $-\Omega_{\text{Ny}}/2 \le \Delta \le \Omega_{\text{Ny}}/2$ highlighted by the gray square, the signal is given in a good approximation by

$$S_n^{\text{signal}} \approx \left[ \langle \sigma \rangle_0 \frac{N}{2} \text{sinc} \left( \frac{\Delta - \Delta_n}{\Omega_{\text{RBW}}} \right) \right]^2 . \quad (\text{F12})$$

In this expression, we use the definition $\text{sinc}(x) = \sin(\pi x)/\pi x$, and the residual bandwidth of the measurement is given by

$$\Omega_{\text{RBW}} = 2\pi/\tau N. \quad (\text{F13})$$

## 4. Signal-to-noise ratio

Owing to the quantum nature of our sensor, the measurement records are in essence discrete, such that a fundamental sampling noise, of spectral shape $S_n^{\text{sampling}}$, affects our measurement. Indeed, the spectrum estimator can be decomposed according to

$$S_n = S_n^{\text{signal}} + S_n^{\text{sampling}}, \quad (\text{F14})$$

with $S_n^{\text{signal}} = |\langle Z_n \rangle|^2$ and $S_n^{\text{sampling}} = \langle |Z_n - \langle Z_n \rangle|^2 \rangle$. To calculate the sampling noise, we can consider the situation where no calibration tone is applied, such that the samples $\{z_k\}_{0 \le k < N}$ are independent, with $\langle z_k \rangle = 0$ and $\langle z_k z_k^* \rangle = \delta_{k,k'}/4$. Combined with the relation (F11), we get

$$S_n^{\text{sampling}} = \langle |Z_k|^2 \rangle = \sum_{k=0}^{N-1} \langle z_k z_k' \rangle = N/4. \quad (\text{F15})$$

By combining the relations (F12) and (F15), we get the signal-to-noise ratio:

$$\text{SNR} = \sqrt{\max(S_n^{\text{signal}})/S_n^{\text{sampling}}} = \sqrt{N\sigma_0}. \quad (\text{F16})$$

The blue curve in the left inset of Fig. 5(c) is calculated using Eqs. (F16) and (F5), with $2\Gamma = (24~\mu\text{s})^{-1}$, in qualitative agreement with the value obtained with a more direct measurement (see main text), and an 84% scaling factor to account for finite readout efficiency.

<figure><img src="image_12.png" /><figcaption>FIG. 10. Aliasing diagram of the charge spectrum analyzer. (a) The blue lines indicate the position of the peaks in discrete frequency space $\Delta_n$ as a function of the continuous frequency $\Delta$ of the applied tone, $\Delta \equiv \Delta_n (\text{mod } 2\Omega_{\text{Ny}})$ and $\Delta \equiv \Omega_{\text{Ny}} - \Delta_n (\text{mod } 2\Omega_{\text{Ny}})$ (see main text for details). (b) Experimental spectrogram with an applied tone of fixed frequency. The detuning $\Delta$ is swept by adjusting the qubit local oscillator frequency. In this particular instance, the repetition time of the experiment was set to $\tau = 14.4~\mu\text{s}$, corresponding to a Nyquist frequency $\Omega_{\text{Ny}}/2\pi = 35~\text{kHz}$. The frequency axis in the Fig. 5(c) has been cropped to only display frequencies from $-\Omega_{\text{Ny}}/2 \le \Delta_n \le \Omega_{\text{Ny}}/2$.</figcaption></figure>

## 5. Approximate expression for the optimal charge sensitivity

The noise spectrum in units of $e^2/\text{Hz}$ is calibrated such that the area under the calibration peak matches the known modulation amplitude:

$$\int_0^\infty S_{ee} d\omega / 2\pi = (2N_{\text{drive}})^2, \quad (\text{F17})$$

where the factor 2 accounts for the number of elementary charges in each Cooper pair. The left-hand side of Eq. (F17) is approximately given by $S_{ee}[\omega_d] \cdot \Omega_{\text{RBW}}/2\pi$, such that the peak of the noise spectrum is given by

$$\max(S_{ee}^{\text{signal}}) = (2N_{\text{drive}})^2 \frac{2\pi}{\Omega_{\text{RBW}}}. \quad (\text{F18})$$
We can now use the definition of the signal-to-noise ratio (in conjunction with the linear relationship between $S_{ee}$ and $S_n$):

$$S_{ee}^{\text{sampling}} = \max(S_{ee}^{\text{signal}})/(\text{SNR})^2. \quad (\text{F19})$$

Additionally, by combining Eq. (F5) with Eq. (F16), we obtain the approximate expression of the signal-to-noise ratio for a calibration tone well within the detector bandwidth ($\Delta \ll \Omega_{\text{full}}$):

$$(\text{SNR})^2 = N\left(\frac{\exp(-\tau_I/T_1)\Omega_R\tau_I}{2}\right)^2. \quad (\text{F20})$$

Finally, by inserting Eq. (F20) into Eq. (F19), and using the expressions (5) for $\Omega_r$ and (F13) for $\Omega_{\text{RBW}}$, we derive

$$S_{ee}^{\text{sampling}} = \delta q^2 = \frac{4\tau}{\tau_I^2 \omega_{ge}^2 \pi^2 \exp(-2\tau_I/T_1)}. \quad (\text{F21})$$

To minimize $\delta q$, it is beneficial to maximize the duty cycle $\tau_I/\tau$. In our experiment, the total preparation and readout time is approximately amounts to $\tau_{\text{prep}} \sim 13$ µs, rendering the total cycle time as $\tau = \tau_I + \tau_{\text{prep}}$. Figure 11 illustrates the evolution of $S_{ee}^{\text{sampling}}$ as a function of the interrogation time $\tau_I$ for various $T_1$ values under two distinct scenarios.

In the first scenario (dotted lines), we have considered the ideal case $\tau_{\text{prep}} = 0$ µs. In this ideal case, the optimal sensitivity is obtained for $\tau_I = T_1/2$, reaching a value

$$\delta q_{\min}^2 = \frac{8 \exp(1)}{T_1 \omega_{ge}^2 \pi^2}. \quad (\text{F22})$$

Remarkably, $\delta q_{\min}$ only depends on the qubit frequency $\omega_{ge}$ and coherence time $T_1$. This stems from the observation that at the flux-frustration point, the Rabi frequency depends only on the product $\omega_{ge}N_{\text{drive}}$ [see Eq. (5)], and not on the specific qubit parameters, as long as the system operates in the heavy-fluxonium regime.

In the second scenario (full lines in Fig. 11), we consider a realistic preparation and readout time $\tau_{\text{prep}} = 13$ µs. As evident from Eq. (F21), for a given interrogation time $\tau_I$, the sensitivity is degraded by a factor $\sqrt{1/\eta}$, where $\eta = \tau_I/\tau$ denotes the duty cycle of the experiment, in comparison to the ideal case. However, the optimal sensitivity, defined as

$$\delta q_{\min} = \min_{\tau_I}(\delta q), \quad (\text{F23})$$

remains close to the ideal one [given in formula (F22)] as long as $\tau_{\text{prep}} \ll T_1$, as visible in Fig. 11(b).

<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /><figcaption>FIG. 11. Sensitivity of the heavy-fluxonium charge spectrum analyzer. (a) Charge sensitivity as a function of interrogation time $\tau_I$, for various qubit lifetimes $T_1$ (see legend). The dashed line corresponds to an ideal scenario where the duty cycle $\tau_I/\tau = 1$ [see Eq. (F22)]. The full line corresponds to a fixed preparation or readout time $\tau_{\text{prep}} = 13$ µs. (b) Minimal sensitivity (obtained for an optimal value of $\tau_I$) as a function of $T_1$. The two scenarios considered in (a) are still represented by dashed and full lines, respectively.</figcaption></figure>

## APPENDIX G: ESTIMATE OF THE CHARGE MODULATION BY A dc-BIASED MEMBRANE

In this appendix, we assess the possibility for the heavy fluxonium to reach the strong-coupling regime with state-of-the-art macroscopic electromechanical systems. For this, we estimate the magnitude of the charge modulation induced by the zero-point fluctuations of a dc-biased vacuum-gap capacitor. In this scenario, we consider that the out-of-plane vibrations of a silicon nitride membrane modulate the capacitance between two parallel electrodes subjected to a dc bias voltage $V_g$.

Table I summarizes the main geometric parameters of the membrane. The membrane lateral dimensions are chosen such that the fundamental mechanical mode matches the qubit frequency $\omega_{ge}$ [40]. The area of the electrodes is chosen to obtain a capacitance $C = 50$ fF matching the value reported in our fluxonium implementation. We assume an electrode separation $h = 500$ nm, which is a conservative estimate based on flip-chip assemblies already
<table><caption>TABLE I. Estimated parameters for a macroscopic electro-mechanical system.</caption><tr><td>Membrane side</td><td><i>l</i></td><td>150 µm</td></tr><tr><td>Membrane stress</td><td><i>σ</i></td><td>1 GPa</td></tr><tr><td>Mechanical mode frequency</td><td><i>Ω</i><sub>m</sub>/2π</td><td>1.8 MHz</td></tr><tr><td>Motional mass</td><td><i>m</i></td><td>3 ng</td></tr><tr><td>Zero-point fluctuations</td><td><i>x</i><sub>zpf</sub></td><td>7 fm</td></tr><tr><td>Silicon nitride density</td><td><i>ρ</i></td><td>3200 kg m<sup>-3</sup></td></tr><tr><td>Capacitor electrodes distance</td><td><i>h</i></td><td>500 nm</td></tr><tr><td>Electrode surface</td><td><i>S</i></td><td>(90 µm)<sup>2</sup></td></tr></table>

reported in the literature [19]. The mechanical resonator undergoes the sum of the restoring force and the electro-static force:

$$F = -m\Omega_m^2(z - h) - V_g^2\epsilon_0S/2z^2. \quad (G1)$$

Mechanical stability requires $V_g < \sqrt{m\Omega_m^2h^3/\epsilon_0S} \approx 50$ V. If we assume a conservative bias voltage $V_g = 5$ V, we obtain

$$N_{\text{drive}} = \frac{V_g}{2e} x_{\text{zpf}} \frac{dC}{dx} \approx 0.01, \quad (G2)$$

where $x_{zpf} = \sqrt{\hbar/2m\Omega_m}$ is the mechanical mode zero-point motion amplitude, and $(dC/dx) \approx C/h$.

[1] Z. Wang, L. Balembois, M. Rančić, E. Billaud, M. L. Dantec, A. Ferrier, P. Goldner, S. Bertaina, T. Chanelière, D. Estève, D. Vion, P. Bertet, and E. Flurin, *Single-electron spin resonance detection by microwave photon counting*, Nature (London) **619**, 276 (2023).

[2] A. D. O'Connell, M. Hofheinz, M. Ansmann, R. C. Bialczak, M. Lenander, E. Lucero, M. Neeley, D. Sank, H. Wang, M. Weides, J. Wenner, J. M. Martinis, and A. N. Cleland, *Quantum ground state and single-phonon control of a mechanical resonator*, Nature (London) **464**, 697 (2010).

[3] Y. Chu, P. Kharel, T. Yoon, L. Frunzio, P. T. Rakich, and R. J. Schoelkopf, *Creation and control of multi-phonon Fock states in a bulk acoustic-wave resonator*, Nature (London) **563**, 666 (2018).

[4] K. J. Satzinger, Y. P. Zhong, H.-S. Chang, G. A. Peairs, A. Bienfait, M.-H. Chou, A. Y. Cleland, C. R. Conner, E. Dumur, J. Grebel, I. Gutierrez, B. H. November, R. G. Povey, S. J. Whiteley, D. D. Awschalom, D. I. Schuster, and A. N. Cleland, *Quantum control of surface acoustic-wave phonons*, Nature (London) **563**, 661 (2018).

[5] P. Arrangoiz-Arriola, E. A. Wollack, Z. Wang, M. Pechal, W. Jiang, T. P. McKenna, J. D. Witmer, R. Van Laer, and A. H. Safavi-Naeini, *Resolving the energy levels of a nanomechanical oscillator*, Nature (London) **571**, 537 (2019).

[6] H. Zhang, S. Chakram, T. Roy, N. Earnest, Y. Lu, Z. Huang, D. K. Weiss, J. Koch, and D. I. Schuster, *Universal fast-flux control of a coherent, low-frequency qubit*, Phys. Rev. X **11**, 011010 (2021).

[7] N. R. A. Lee, Y. Guo, A. Y. Cleland, E. A. Wollack, R. G. Gruenke, T. Makihara, Z. Wang, T. Rajabzadeh, W. Jiang, F. M. Mayor, P. Arrangoiz-Arriola, C. J. Sarabalis, and A. H. Safavi-Naeini, *Strong dispersive coupling between a mechanical resonator and a fluxonium superconducting qubit*, arXiv:2304.13589.

[8] A. N. Korotkov and M. A. Paalanen, *Charge sensitivity of radio frequency single-electron transistor*, Appl. Phys. Lett. **74**, 4052 (1999).

[9] A. Aassime, D. Gunnarsson, K. Bladh, P. Delsing, and R. Schoelkopf, *Radio-frequency single-electron transistor: Toward the shot-noise limit*, Appl. Phys. Lett. **79**, 4031 (2001).

[10] S. J. Angus, A. J. Ferguson, A. S. Dzurak, and R. G. Clark, *A silicon radio-frequency single electron transistor*, Appl. Phys. Lett. **92**, 112103 (2008).

[11] W. Lu, Z. Ji, L. Pfeiffer, K. West, and A. Rimberg, *Real-time detection of electron tunnelling in a quantum dot*, Nature (London) **423**, 422 (2003).

[12] R. Schoelkopf, P. Wahlgren, A. Kozhevnikov, P. Delsing, and D. Prober, *The radio-frequency single-electron transistor (RF-SET): A fast and ultrasensitive electrometer*, Science **280**, 1238 (1998).

[13] M. C. Cassidy, A. S. Dzurak, R. G. Clark, K. D. Petersson, I. Farrer, D. A. Ritchie, and C. G. Smith, *Single shot charge detection using a radio-frequency quantum point contact*, Appl. Phys. Lett. **91**, 222104 (2007).

[14] C. Volk, A. Chatterjee, F. Ansaloni, C. M. Marcus, and F. Kuemmeth, *Fast charge sensing of Si/SiGe quantum dots via a high-frequency accumulation gate*, Nano Lett. **19**, 5628 (2019).

[15] M. Gonzalez-Zalba, S. Barraud, A. Ferguson, and A. Betz, *Probing the limits of gate-based charge sensing*, Nat. Commun. **6**, 6084 (2015).

[16] J. J. Viennot, M. R. Delbecq, M. C. Dartiailh, A. Cottet, and T. Kontos, *Out-of-equilibrium charge dynamics in a hybrid circuit quantum electrodynamics architecture*, Phys. Rev. B **89**, 165404 (2014).

[17] H. Brenning, S. Kafanov, T. Duty, S. Kubatkin, and P. Delsing, *An ultrasensitive radio-frequency single-electron transistor working up to 4.2 K*, J. Appl. Phys. **100**, 114321 (2006).

[18] M. P. Blencowe and M. N. Wybourne, *Sensitivity of a micromechanical displacement detector based on the radio-frequency single-electron transistor*, Appl. Phys. Lett. **77**, 3845 (2000).

[19] Y. Seis, T. Capelle, E. Langman, S. Saarinen, E. Planz, and A. Schliesser, *Ground state cooling of an ultracoherent electromechanical system*, Nat. Commun. **13**, 1507 (2022).

[20] J. J. Viennot, X. Ma, and K. W. Lehnert, *Phonon-number-sensitive electromechanics*, Phys. Rev. Lett. **121**, 183601 (2018).

[21] M. F. Gely and G. A. Steele, *Superconducting electromechanics to test Diósi-Penrose effects of general relativity in massive superpositions*, AVS Quantum Sci. **3**, 035601 (2021).
[22] V. E. Manucharyan, *Superinductances*, Ph.D. thesis, Yale University, 2012.

[23] V. E. Manucharyan, J. Koch, L. I. Glazman, and M. H. Devoret, *Fluxonium: Single Cooper-pair circuit free of charge offsets*, *Science* **326**, 113 (2009).

[24] Y.-H. Lin, L. B. Nguyen, N. Grabon, J. San Miguel, N. Pankratova, and V. E. Manucharyan, *Demonstration of protection of a superconducting qubit from energy decay*, *Phys. Rev. Lett.* **120**, 150503 (2018).

[25] Y. Tsaturyan, A. Barg, E. S. Polzik, and A. Schliesser, *Ultracoherent nanomechanical resonators via soft clamp- ing and dissipation dilution*, *Nat. Nanotechnol.* **12**, 776 (2017).

[26] E. Ivanov, T. Capelle, M. Rosticher, J. Palomo, T. Briant, P.-F. Cohadon, A. Heidmann, T. Jacqmin, and S. Deléglise, *Edge mode engineering for optimal ultracoherent silicon nitride membranes*, *Appl. Phys. Lett.* **117**, 254102 (2020).

[27] F. Diedrich, J. C. Bergquist, W. M. Itano, and D. J. Wineland, *Laser cooling to the zero-point energy of motion*, *Phys. Rev. Lett.* **62**, 403 (1989).

[28] J. D. Teufel, T. Donner, D. Li, J. W. Harlow, M. S. Allman, K. Cicak, A. J. Sirois, J. D. Whittaker, K. W. Lehnert, and R. W. Simmonds, *Sideband cooling of micromechanical motion to the quantum ground state*, *Nature (London)* **475**, 359 (2011).

[29] X. Y. Jin, A. Kamal, A. P. Sears, T. Gudmundsen, D. Hover, J. Miloshi, R. Slattery, F. Yan, J. Yoder, T. P. Orlando, S. Gustavsson, and W. D. Oliver, *Thermal and residual excited-state population in a 3D transmon qubit*, *Phys. Rev. Lett.* **114**, 240501 (2015).

[30] H. Sun, F. Wu, H.-S. Ku, X. Ma, J. Qin, Z. Song, T. Wang, G. Zhang, J. Zhou, Y. Shi *et al.*, *Characterization of loss mechanisms in a fluxonium qubit*, *Phys. Rev. Appl.* **20**, 034016 (2023).

[31] C. Bonato, M. S. Blok, H. T. Dinani, D. W. Berry, M. L. Markham, D. J. Twitchen, and R. Hanson, *Optimized quantum sensing with a single electron spin using real-time adaptive measurements*, *Nat. Nanotechnol.* **11**, 247 (2016).

[32] E. Polino, M. Valeri, N. Spagnolo, and F. Sciarrino, *Photonic quantum metrology*, *AVS Quantum Sci.* **2**, 024703 (2020).

[33] Alan V. Oppenheim, Ronald W. Schafter, and John R. Buck, *Discrete-Time Signal Processing* (Prentice Hall, Upper Saddle River, NJ, 1999), Chap. 2.

[34] C. L. Degen, F. Reinhard, and P. Cappellaro, *Quantum sensing*, *Rev. Mod. Phys.* **89**, 035002 (2017).

[35] D. Reilly, C. Marcus, M. Hanson, and A. Gossard, *Fast single-charge sensing with a rf quantum point contact*, *Appl. Phys. Lett.* **91**, 162101 (2007).

[36] W. C. Smith, A. Kou, U. Vool, I. M. Pop, L. Frunzio, R. J. Schoelkopf, and M. H. Devoret, *Quantization of inductively shunted superconducting circuits*, *Phys. Rev. B* **94**, 144507 (2016).

[37] Z. Leghtas, S. Touzard, I. M. Pop, A. Kou, B. Vlastakis, A. Petrenko, K. M. Sliwa, A. Narla, S. Shankar, M. J. Hatridge, M. Reagor, L. Frunzio, R. J. Schoelkopf, M. Mirrahimi, and M. H. Devoret, *Confining the state of light to a quantum manifold by engineered two-photon loss*, *Science* **347**, 853 (2015).

[38] R. Azouit, *Adiabatic Elimination for Open Quantum Systems*, Ph.D. thesis, MINES ParisTech, 2017.

[39] A. A. Clerk, M. H. Devoret, S. M. Girvin, F. Marquardt, and R. J. Schoelkopf, *Introduction to quantum noise, measurement, and amplification*, *Rev. Mod. Phys.* **82**, 1155 (2010).

[40] P.-L. Yu, T. P. Purdy, and C. A. Regal, *Control of material damping in high-Q membrane microresonators*, *Phys. Rev. Lett.* **108**, 083603 (2012).