# Probing the quantum motion of a macroscopic mechanical oscillator with a radio-frequency superconducting qubit

Kyrylo Gerashchenko, Rémi Rousseau, Léo Balembois, Himanshu Patange,

Paul Manset, W. Clarke Smith, Zaki Leghtas, Emmanuel Flurin, Thibaut

Jacqmin, Samuel Deléglise

► To cite this version:

Kyrylo Gerashchenko, Rémi Rousseau, Léo Balembois, Himanshu Patange, Paul Manset, et al. Prob-
ing the quantum motion of a macroscopic mechanical oscillator with a radio-frequency superconducting
qubit. 2025. [hal-05088558]

HAL Id: hal-05088558

https://hal.science/hal-05088558v1

Preprint submitted on 28 May 2025

**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
Probing the quantum motion of a macroscopic mechanical oscillator with a
radio-frequency superconducting qubit

K. Gerashchenko,¹, * R. Rousseau,¹, ², * L. Balembois,¹, * H. Patange,¹ P. Manset,¹

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

¹Laboratoire Kastler Brossel, Sorbonne Université, CNRS, ENS-Université
PSL, Collège de France, 4 place Jussieu, 75005 Paris, 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, Sorbonne Université, Université Paris Cité, Centre Automatique et Systèmes, Mines Paris, Université PSL, Inria, Paris, France

⁴Google Quantum AI, Santa Barbara, CA

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

(Dated: May 28, 2025)

Long-lived mechanical resonators like drums oscillating at MHz frequencies and operating in the quantum regime offer a powerful platform for quantum technologies and tests of fundamental physics [1]. Yet, quantum control of such systems remains challenging, particularly owing to their low energy scale and the difficulty of achieving efficient coupling to other well-controlled quantum devices. Here, we demonstrate repeated, and high-fidelity interactions between a 4 MHz suspended silicon nitride membrane and a resonant superconducting heavy-fluxonium qubit [2–4]. The qubit is initialized at an effective temperature of 27 µK and read out in a single-shot with 77 % fidelity. During the membrane’s 6 ms lifetime the two systems swap excitations more than 300 times. After each interaction, a state-selective detection is performed, implementing a stroboscopic series of weak measurements that provide information about the mechanical state. The accumulated records reconstruct the membrane’s position noise-spectrum, revealing both its thermal occupation $n_{th} \approx 47$ at 10 mK and the qubit-induced back-action. By preparing the qubit either in its ground or excited state before each interaction, we observe an imbalance between the emission and absorption spectra, proportional to $n_{th}$ and $n_{th} + 1$ respectively, a hallmark of the non-commutation of phonon creation and annihilation operators. Since the predicted Diósi–Penrose gravitational collapse time [5, 6] is comparable to the measured mechanical decoherence time, our architecture enters a regime where gravity-induced decoherence could be tested directly [7].

Mechanical resonators, such as vibrating membranes or beams, are essential components of quantum technologies, serving as high-precision sensors [8–11] and trans-

ducers between distinct quantum systems [12–14]. They also provide a platform for probing fundamental physics at the intersection of quantum mechanics and gravity [7]. For instance, owing to their long lifetime and large zero-point fluctuations, state-of-the-art drums oscillating at MHz frequencies have been identified as good candidates to test the existence of a gravitational collapse phenomenon postulated by Diosi [5] and Penrose [6]. Measuring and controlling the quantum state of these systems remains challenging because it requires coupling them to a resonant two-level system [15]. In circuit Quantum Acoustodynamics (cQAD) [16–18], this is typically achieved with superconducting qubits that operate in the GHz range. Significant progress has been made toward bridging the gap with macroscopic mechanical resonators. For instance, mechanical squeezing was demonstrated in a 22 MHz aluminum drum dispersively coupled to a 4 GHz Cooper-pair box [19]. However, coherent interaction and efficient qubit readout are hindered by the short coherence times of this qubit architecture. More recently, techniques from cQAD have been extended to observe coherent exchange dynamics between a 700 MHz piezoelectric resonator and a fluxonium qubit [20]. However, owing to the resonator’s extremely small mass (only a few picograms), the gravitational decoherence time predicted by the Diosi–Penrose model is roughly six orders of magnitude longer than the environmental decoherence time.

In this paper, we report an experimental implementation of an idealized mechanical-cavity quantum electrodynamics system, where a mechanical resonator is repeatedly interrogated by a resonant two-level system. The mechanical resonator is a suspended silicon nitride membrane oscillating at $\omega_m/2\pi \sim 4.4$ MHz. It is coupled to a heavy-fluxonium circuit, which features simultaneously a low frequency qubit transition [3, 4], resonant with the MHz mechanical mode, and GHz transitions to higher excited states that can be used to readout and prepare the qubit state deterministically (see Fig. 1a). We rely exclusively on the binary outcomes of qubit measurements to infer the mechanical state. Owing to the long mechanical lifetime ($T_1^m \sim 6$ ms), the qubit interacts with the membrane over 300 times during that in-

<footnote>* These three authors contributed equally</footnote>

<footnote>† These authors jointly supervised the work samuel.deleglise@lkb.upmc.fr</footnote>
<figure><img src="image_2.png" /><figcaption>FIG. 1. **Experimental principle and design.** (a) Energy-level diagram of the mechanical resonator (left) and fluxonium circuit (right). The membrane mode at frequency $\omega_m$ is modeled as a quantum harmonic oscillator. The heavy fluxonium's potential (black line) at $\varphi_{ext} = \pi$ is shown with its four lowest eigenstates. The qubit transition (blue/red) is tunable into resonance with the mechanical oscillator. (b) 3D schematic of the mechanical-fluxonium device. A superconducting loop composed of a small Josephson junction (red) and a superinductor (purple), is threaded by an external flux $\varphi_{ext}$ (brown). The small junction is shunted by a capacitor (blue). One capacitor electrode couples to a readout resonator (light pink) and the other to a control waveguide (cyan). The suspended SiN membrane, metallized with Al (yellow), and the two fluxonium electrodes form a vacuum-gap capacitor. A charge gate (orange) provides DC bias and drives to the membrane. (c) Spatial profile of the 4.4 MHz mechanical mode of the membrane, plotted at an amplitude equal to its zero-point fluctuations and overlaid on the fluxonium capacitor electrodes layout (dashed grey lines). (d) Colored optical and electronic micrograph of a twin fluxonium sample. (e) Lumped element circuit diagram of the sample.</figcaption></figure>

terval, extracting information bit-by-bit about its state. This process realizes a periodic weak measurement of the mechanical mode, from which we reconstruct the noise spectrum $S_{xx}(\omega)$ of the dimensionless membrane position $x$ under the influence of the thermal environment of occupation $n_{th}$. By letting the membrane interact with the qubit either prepared in its ground, or excited state, we probe separately the emission and absorption noise spectra around $\pm\omega_m$. This measurement reveals the fundamental quantum asymmetry between emission and absorption processes, proportional to $n_{th}$ and $n_{th}+1$ respectively. While similar signatures have already been observed through sideband asymmetries under dispersive coupling to microwave [21] or optical [22–24] cavities, this is the first time, to our knowledge, the asymmetric mechanical noise spectrum is measured via coupling to a two-level system, as initially proposed theoretically [25–28]. Most importantly, given the mass of our system ($M \sim 5.3$ ng), the postulated Diósi-Penrose gravitational collapse time [29] $\tau_G \sim 0.5$ ms is comparable to the mechanical decoherence time $\tau_{th} = 2T_1^m/n_{th} \sim 0.3$ ms, positioning our experiment in a regime where gravitational effects could become experimentally relevant (see App. B).

Our system, shown in Fig. 1, consists of two distinct chips. The fluxonium qubit [2] is fabricated on an intrinsic silicon substrate using standard microfabrication tech-

niques (App. C). It consists of two coplanar electrodes shunted by a nonlinear Josephson junction ($E_J/h = 4.82(6)$ GHz) and a large inductor formed by 360 Josephson junctions ($E_L/h = 0.128(7)$ GHz). The mechanical resonator is a clamped 110 × 140 µm², 90 nm-thick rectangular silicon nitride film under $\gtrsim 800$ MPa in-plane tensile stress, coated with a superconducting aluminum pad [13, 30–32] (see App. C for details). The qubit primarily couples to the membrane mode at $\omega_m/2\pi = 4.4$ MHz, which exhibits an anti-symmetric displacement pattern across the Y-axis (see Fig. 2b). The coupling is achieved via a flip-chip technique, forming two vacuum gap capacitors $C_m^{\pm} = C_m \pm zdC_m^{\pm}/dx$ between the membrane and fluxonium electrodes, yielding a total fluxonium capacitive energy $E_C/h = 0.418(9)$ GHz. Static charges are forced onto the membrane via a large capacitive electrode, maintaining the membrane at a fixed potential $V_b$, while the fluxonium island is referenced to ground by capacitors $C_g$. This scheme enables the transmission of the membrane motion into a charge modulation on the fluxonium electrodes.

The external magnetic flux $\varphi_{ext}$ through the inductive loop can be controlled via an on-chip flux line. When $\varphi_{ext} = \pi$, the flux degree of freedom experiences a double-well potential, illustrated in Fig. 1a. The transition between the two lowest energy eigenstates, $|g\rangle$ and $|e\rangle$, is governed by the tunneling rate across the central barrier,
<figure><img src="image_3.png" /><figcaption>FIG. 2. Qubit-mechanical coupling. (a) Ramsey spectroscopy under mechanical drive. Top: Pulse sequence—each cycle begins with a mechanical displacement of the membrane, immediately followed by a fixed-delay Ramsey sequence, repeated to sample the qubit frequency shift $\delta\omega = \chi \langle \hat{a}^\dagger \hat{a} \rangle$. Bottom: Extracted qubit frequency shift versus bias voltage and drive detuning (color map). Dark squares mark the Lorentzian response peak at each $V_b$, with a quadratic fit (dashed line). The drive amplitude is adjusted at each voltage to keep the peak response level approximately constant (except for $V_b = V_{\text{offset}}$ where the driving force vanishes). (b) Mechanical lifetime. Top: Pulse sequence—after resonant membrane displacement, a series of $N = 121$ Ramsey measurements are performed, sampling the instantaneous qubit frequency at discrete times $k \times T$. Bottom: Averaged qubit frequency shift $\delta\omega$ plotted versus delay $k \times T$ (for $k = 1, \dots, N$), where data at each delay are ensemble-averaged. An exponential fit yields a mechanical lifetime $T_1^m = 5.9$ ms ($Q = 1.62 \times 10^5$). The right-hand axis translates $\delta\omega$ into the equivalent mechanical occupation in units of phonon number. (c) Coherent Rabi exchange. Top: Pulse sequence—after initializing the membrane in a coherent state $|\alpha\rangle$, with $|\alpha|^2 \simeq 1.6 \cdot 10^5$, the qubit is Stark-shifted to $\omega_q^{\text{st}}$ for an interaction time $\tau$, and measured along $\sigma_z$. Bottom: Ground-state population as a function of interaction time and Stark-shifted qubit frequency. The spheres on the top right of each figure show the qubit Bloch vectors immediately after the interaction for the data points highlighted by blue circles.</figcaption></figure>

which decreases exponentially with $\sqrt{E_J/E_C}$.

In the heavy-fluxonium regime $[3, 4] (E_J \gg E_C > E_L)$, the circuit is highly anharmonic, with a qubit transition frequency $\omega_q/2\pi = 2.35$ MHz, while the transition to the next excited state occurs at 3.45 GHz. This permits treating $\{|g\rangle, |e\rangle\}$ as an almost ideal two-level system, with higher levels used for fast reset and dynamic tuning. A two-photon preparation protocol yields fidelities $\eta_g = 98.5(2.9)\%$ and $\eta_e = 98.3(1.0)\%$ (App. E 5), while single-shot dispersive readout achieves measurement fidelities $F_g = 84.8(1.8)\%$ and $F_e = 69.1(1.2)\%$ (App. E 6). By driving the circuit near the higher excited state transitions at the flux sweet spot $\varphi_{\text{ext}} = \pi$, we dynamically tune $\omega_q$ via the AC-Stark effect to bring it into resonance with the 4.4 MHz membrane (App. E 2). Coherence times are $T_1 = 22.4$ µs and $T_2 = 17.1$ µs, decreasing to $T_1 = 7.4$ µs and $T_2 = 4.2$ µs when tuned to $\omega_m$ via the AC-Stark shift (App. E 3).

A membrane mode of frequency $\omega_m$ can be modeled as a harmonic oscillator of Hamiltonian $H_m = \hbar\omega_m(a^\dagger a + 1/2)$, where $a$ and $a^\dagger$ are the ladder operators (Fig. 1a). Projecting the system dynamics in the $\{|g\rangle, |e\rangle\}$ manifold and applying the rotating wave approximation yields the Jaynes-Cummings Hamiltonian

$$H_{JC}/\hbar = \omega_m a^\dagger a + \frac{\omega_q}{2} \sigma_z + i \frac{\Omega}{2} (a^\dagger \sigma - a \sigma^\dagger), \quad (1)$$

where $\sigma$ represents the lowering operator for the qubit

and the coupling strength is given by

$$\Omega = \omega_q \langle g | \varphi | e \rangle \left| \frac{\partial C_m^+}{\partial x} \beta \frac{V_b - V_{\text{offset}}}{2e} \right| \quad (2)$$

Here, $V_b$ denotes the membrane bias voltage, $\beta \sim C_g/(C_g + C_m) \sim 1/2$ is a dilution factor resulting from the capacitive biasing scheme, and $V_{\text{offset}} = Q_2/2C_g$ is an offset voltage resulting from the static charge $Q_2$ confined on the floating fluxonium island. With the maximum effective voltage of $|V_b - V_{\text{offset}}| \sim 11.5$ V practically achievable in our setup, we measure $\Omega/2\pi \sim 1.50(8)$ kHz (App. A).

In absence of AC-Stark effect, the qubit and the membrane modes are out-of-resonance: the large detuning $|\Delta| \gg \Omega$ suppresses direct energy exchange between the two systems. Instead, the qubit undergoes a frequency shift proportional to the membrane's phonon population $\chi a^\dagger a \sigma_z$, where $\chi = \Omega^2/4\Delta$ [33, 34]. This effect is detected via a Ramsey spectroscopy experiment (see Fig. 2a). The qubit is initialized in the $|+\rangle = (|g\rangle + |e\rangle)/\sqrt{2}$ state and evolves for a duration $\tau$ with the membrane in a coherent state, during which it accumulates a phase proportional to the phonon number. The experiment alternates between measuring the $\sigma_x$ and $\sigma_y$ components of the superposition over multiple realizations. Fig. 2a shows the qubit frequency shift as a function of the membrane drive frequency and applied DC
<figure><img src="image_4.png" /></figure>

FIG. 3. Spectrum analyzer experiment in the semi-classical limit. (a) Phase-space distribution of the mechanical mode amplitude, with an overlaid classical trajectory (black curve), and the qubit Bloch sphere. After initialization in $|g\rangle$ (blue) or $|e\rangle$ (red), the membrane mode acts as a coherent drive on the qubit. The complex mechanical mode amplitude is mapped onto the $(\sigma_x, \sigma_y)$ plane of the Bloch sphere. The average energy transferred from membrane to qubit during the interaction $\Delta E_{g/e}$ is displayed on the $\langle\sigma_z\rangle$ histogram on the right-hand inset. (b) Classical position $x(t)$ of the mechanical resonator, shown in a frame rotating at the qubit-mode detuning $\Delta$. The signal is an oscillation at frequency $\Delta$, with amplitude and phase undergoing random fluctuations at rate $\kappa_m = 1/T_1^m$. Blue bars represent qubit $\sigma_x$ measurement outcomes after interaction with the mechanical mode. (c) Experimental pulse sequence. The qubit, prepared in $|g\rangle$ or $|e\rangle$, is Stark-shifted to a detuning $\Delta$ from the membrane frequency for a duration $\tau$, followed by a $\sigma_x$ measurement. (d) Single-shot outcomes $m_k \in \{-1, 1\}$ (for initial state $|e\rangle$) are recorded, grouped into batches of length $N$, and each batch is discrete-Fourier-transformed. (e) Fourier transform of the first 10 batches represented as a waterfall plot. (f) Position noise spectra obtained for the qubit prepared in $|g\rangle$ (blue) and $|e\rangle$ (red), by averaging the squared magnitudes of the Fourier transforms shown in (e). The two Lorentzian peaks on the left are centered at the membrane frequency and correspond to the membrane thermal noise. Additionally, a weak calibration tone produces a sinc-shaped peak offset by 500 Hz towards higher frequencies. (g) Zoom on the Lorentzian peaks, with each spectrum normalized to the calibration-tone amplitude and flat background subtracted. The broader, lower-amplitude blue curve ($|g\rangle$) indicates qubit-induced cooling, while the narrower, higher-amplitude red curve ($|e\rangle$) indicates heating of the mechanical mode. Data in (a) and (b) are simulated for illustration purpose, whereas real data are displayed in (d), (e), (f), (g).

bias voltage $V_b$. The data reveal a quadratic decrease in the mechanical frequency with $V_b - V_{\text{offset}}$, explained by the electrostatic spring-softening effect (App. D 4) [35].

The mechanical resonator's relaxation time $T_1^m$, is extracted by preparing the membrane in a coherent state and measuring the exponential decay of its phonon population. This decay manifests as a time variation of the qubit frequency, tracked through repeated Ramsey measurements. An exponential fit to the data, shown in Fig. 2b, yields $T_1^m = 5.9$ ms. The corresponding quality factor, $Q_m = \omega_m T_1^m = 1.62 \cdot 10^5$, is consistent with similar studies on silicon nitride membranes [36].

We now focus on the resonant interaction regime. In these experiments, the qubit is prepared and measured at its bare resonance frequency $\omega_q$, while the resonant interaction, is activated for a duration $\tau$ by tuning the

qubit frequency $\omega_q^{\text{st}} = \omega_m + \Delta$ close to the membrane one, governed by the evolution operator $U_I = e^{-iH_{\text{C}}\cdot\tau/\hbar}$. When the membrane is initialized in a large enough coherent state $\alpha$, characteristic Rabi oscillations emerge at a rate $\Omega|\alpha| \gg 1/T_1$ (see Fig. 2c). In this limit, $U_I$ describes a two-level system undergoing a rotation by an angle $|\alpha|\theta_1$ about the equatorial axis defined by $i\alpha^*$ under the influence of a classical drive.

When the membrane is in a thermal state, the resonant interaction induces a qubit precession of typical angle $\theta_1 \sqrt{n_{\text{th}}}$, where $\theta_1 = \Omega\tau \sim 47$ mrad is the angle induced by a single phonon. Though small, this precession is detectable through sequential qubit measurements performed within a single mechanical lifetime [4, 37]. Starting from the prepared qubit state (either $|g\rangle$ or $|e\rangle$), we turn on the interaction for a duration $\tau = 5$ µs. We then
project the qubit onto $\sigma_x$ by applying a $\sqrt{X}$ rotation, followed by a $\sigma_z$ readout. This sequence is periodically repeated with the qubit reinitialized in the same state (Fig. 3c).

To build physical intuition, we first treat the membrane's degrees of freedom—its position $\mathbf{x}$ and momentum $\mathbf{p}$—classically. In the qubit drive frame, rotating at $\omega_q^{st}$, the complex mechanical amplitude $\alpha(t) = \tilde{\alpha}(t)e^{i\Delta t}$ appears as a signal oscillating at frequency $\Delta$, modulated by a slowly varying envelope $\tilde{\alpha}(t)$ that diffuses in phase space at a characteristic rate $\kappa_m$, as shown in Fig. 3a. This amplitude is sampled at discrete times $t_k = k \times T$ by the sequential qubit measurements of $\sigma_{x,k}$, where $T = 15.5$ µs. In the small angle limit $|\alpha|\theta_1 \ll 1$, the rotation $U_I$ maps the phase space coordinates (Re($\alpha(t_k)$), Im($\alpha(t_k)$)) of the mechanical oscillator onto the qubit's transverse components ($\sigma_{x,k}, \sigma_{y,k}$) (Fig. 3b). Therefore,

$$ \langle \sigma_{x,k} \rangle = \theta_1 x(t_k), \quad (3) $$

where $x(t_k) = \text{Re}(\alpha(t_k))$ is the dimensionless membrane position. By operating the qubit in a frame tuned by $\Delta > \kappa_m$, the free evolution between interaction pulses continually rotates the measurement axis. Consequently, successive $\sigma_{x,k}$ measurements sweep through all mechanical quadrature angles—probing both Re($\tilde{\alpha}(t_k)$) and Im($\tilde{\alpha}(t_k)$)—in direct analogy with a heterodyne detection.

Our focus being on characterizing the steady-state statistical properties, we directly compute the noise spectral density $S_{mm}^{g/e}$ of the telegraphic signal, where the superscript $g/e$ denotes the initial qubit state used in the sequence. To do so, we first group the measurement outcomes $\{m_k\}$ into batches of $N = 10000$ samples (see Fig. 3d). For each batch, we compute the discrete Fourier transform, $\mathcal{F}[m_k]_r$, where the frequency bins are defined as $\omega_r = 2\pi r/NT$ for $0 \le r < N$. We then average the squared magnitudes over all batches to obtain the noise spectral density $S_{mm}(\omega_r) = |\mathcal{F}[m_k]_r|^2 T/N$, with the overline indicating an average over batches corresponding to distinct classical trajectories $\{x(t_k)\}_{0 \le k < N}$ (see Fig. 3e). At our sampling rate, the aliasing-free frequency span is $\omega_{N/2} \approx 2\pi \times 32$ kHz. We therefore choose a small qubit-membrane detuning $\Delta/2\pi = 5$ kHz to ensure the oscillations of $x(t)$ remain entirely within the Nyquist zone.

In a semiclassical approximation that neglects qubit backaction on the membrane, the trajectory $x(t)$ acts as an independent external drive, so that the qubit's response at each sampling time depends only on the instantaneous $x(t_k)$ and not on prior measurements. Consequently, the two-point correlator factorizes:

$$ C_k \equiv \overline{\langle \sigma_{x,k} \sigma_{x,0} \rangle} = \begin{cases} \overline{\langle \sigma_{x,k} \rangle} \langle \sigma_{x,0} \rangle, & k \neq 0, \\ 1, & k = 0. \end{cases} \quad (4) $$

Recalling Eq. (3), and Wiener-Khinchin theorem, linking the measured spectrum to the Fourier-transform of $C_k$, we have

$$ S_{mm}^{g/e}(\omega) = \frac{T}{N} \mathcal{F}[C_k] = \frac{\theta_1^2}{4} S_{xx}(\omega) + T. \quad (5) $$

where $S_{xx}$ is the position noise spectrum and $T$ is the constant background arising from the sampling noise due to the binary nature of the qubit measurements. In our experimental data (see Fig. 3f), the spectra $S_{mm}^{g/e}$ display a Lorentzian peak centered at $\omega_m$ that reflects the thermal fluctuations of the membrane, with the sampling noise background $T$ indicated by a dashed line. Although Eq. (5) predicts identical spectra regardless of the initial qubit state, we observe significant differences in both the area and the linewidth of the thermal peaks.

These discrepancies arise because our simple model neglects the energy exchange between the membrane and the qubit during their interaction. In particular, when the qubit is prepared in $|g\rangle$ it absorbs energy from the membrane, whereas when prepared in $|e\rangle$, it emits energy to the membrane (see Fig. 3a). On average, the qubit energy is lowered or raised on the Bloch sphere by $\Delta E_{g,e}/\hbar\omega_q^{st} = \theta_1^2 n_{th}/4$ due to the interaction. Because of energy conservation, the membrane must experience a qubit-induced relaxation or amplification at a rate $\kappa = \theta_1^2/4T$. The mechanical linewidth is thus modified (broadened or narrowed) according to $\kappa'_{g/e} = \kappa_m \pm \kappa$. Furthermore the qubit being reinitialized in a pure state before each interaction, it provides a 0-occupation bath for the membrane, leading to an occupation $n'_{th,i} = \kappa_m n_{th}/\kappa'_i$. This qubit-induced cooling (or heating) mechanism is analogous to the dynamical backaction effects observed in cavity optomechanics [38, 39].

Lorentzian fits to the thermal peaks yield linewidths of $\kappa'_g = 2\pi \times (28.5 \pm 1.4)$ Hz for $|g\rangle$ preparations and $\kappa'_e = 2\pi \times (23.7 \pm 0.7)$ Hz for $|e\rangle$ preparations. The inferred linewidth $\kappa_m = (\kappa'_e + \kappa'_g)/2 = 2\pi \times (26.1 \pm 0.8)$ Hz is in agreement with the previously measured $1/T_1^m = 2\pi \times (27.16 \pm 2.6)$ Hz, confirming that the membrane exhibits no additional dephasing beyond its Fourier-limited linewidth. Moreover, the dynamical back-action rate $\kappa = 2\pi \times (2.4 \pm 0.8)$ Hz agrees with a full experimental model that includes finite qubit preparation and readout fidelities. To account for these imperfections, the spectrum is normalized using a calibration peak generated by a weak continuous tone applied on the qubit charge drive port (see Fig. 3f). The normalized, background-subtracted spectra are displayed in Fig. 3g. The experimentally obtained area ratio $A_e/A_g = 1.42 \pm 0.11$ is in a reasonable agreement with the expected value $n'_{th,e}/n'_{th,g} = \kappa'_g/\kappa'_e = 1.22 \pm 0.07$. Here the dynamical backaction is moderate due to the relatively low cooperativity $\kappa/\kappa_m \sim 1/10$. Nevertheless, we expect that ground-state cooling will be within reach in future experiments—either by further increasing the mechanical lifetime [40, 41] or by boosting the qubit-membrane coupling strength via a larger DC bias.

In a full quantum treatment, each interaction cycle constitutes a weak measurement of the membrane state
<figure><img src="image_5.png" /><figcaption>FIG. 4. Quantum position-spectrum asymmetry. (a) Position-noise spectrum of a harmonic oscillator at thermal equilibrium, showing absorption (blue peak at $-\omega_m$, area $\propto \bar{n}_{th}$) and emission (red peak at $+\omega_m$, area $\propto \bar{n}_{th} + 1$). (b) Experimental pulse sequence: to eliminate the average dynamical backaction, we alternately initialize the qubit in the ground state $|g\rangle$ and the excited state $|e\rangle$. From the resulting measurement records, we then reconstruct the emission spectrum $S_{xx}^g$ (when starting in $|g\rangle$) and the absorption spectrum $S_{xx}^e$ (when starting in $|e\rangle$). (c) Reconstructed spectra for $|g\rangle$ (blue) and $|e\rangle$ (red) preparations at varying cryostat temperatures at 10, 15, 20, 25, and 30 mK, as illustrated by the thermometer. (d) Lorentzian-peak areas at $\omega_m$ versus cryostat temperature: blue and red markers indicate $S_{xx}^g$ and $S_{xx}^e$ areas for each qubit preparation. A linear fit to the mean-area of the 5 temperature series is shown (dashed line), and the right-hand axis converts the fitted peak area into mean phonon number using this calibration. (e) Measured phonon-number difference between $|e\rangle$ and $|g\rangle$ spectra versus cryostat temperature. The horizontal dashed line at 1.37 quanta is the theoretical prediction, which differs from 1 due to qubit preparation infidelity and finite lifetime. All points include error bars representing the standard deviation from bootstrap fits to random FFT-batch sub-ensembles.</figcaption></figure>

$\rho$, described by the operators $M_{\pm}^i = \langle \pm | U_I | i \rangle$ where the superscript $i \in \{g, e\}$ denotes the initial qubit state, and the subscript $\pm$ the measurement outcome, occurring with probability $p_{\pm}^i = \text{Tr}[M_{\pm}^i \rho M_{\pm}^{i\dagger}]$. The mechanical state update, conditioned on the measurement outcome $\pm$ reads [42, 43]

$$ \rho \rightarrow \frac{M_{\pm}^{i} \rho M_{\pm}^{i\dagger}}{p_{\pm}^{i}}. \qquad (6) $$

The dynamical backaction effects described above (qubit induced cooling or amplification) follow naturally from writing the mean evolution, where qubit measurement outcomes are discarded:

$$ \rho_i \rightarrow \sum_{f=\pm} M_f^i \rho M_f^{i\dagger}. \qquad (7) $$

Expanding this expression to second order in $\theta_1$ yields a dissipation channel characterized by a loss operator that depends on the qubit state: $L_g = \sqrt{\kappa}a$ (i.e. phonon annihilation) when the qubit is prepared in $|g\rangle$ and $L_e = \sqrt{\kappa}a^\dagger$ (i.e. phonon creation) when prepared in $|e\rangle$ (see App. F 4).

Not only does each qubit-membrane interaction transfer energy on average, but it also entangles the two systems. A subsequent projective measurement of the qubit then steers the membrane into a conditional state that depends on the measurement outcome, imprinting information which persists over the mechanical lifetime $T_1^m$. These measurement-induced quantum correlations yield a two-point qubit correlator that departs from the semiclassical expression of Eq. (4) (see App. F 4), for $k \neq 0$

$$ \begin{aligned} C_k^g &= \theta_1^2 \text{Re}[\langle a^\dagger(t_k) a(0) \rangle] \\ C_k^e &= \theta_1^2 \text{Re}[\langle a(t_k) a^\dagger(0) \rangle]. \end{aligned} \qquad (8) $$

In these expressions $a(t)$ is the annihilation operator expressed in the Heisenberg picture. Upon Fourier transform, these correlation functions yield two spectra peaked around $\omega_m$, proportional to $\langle a^\dagger a \rangle$ and $\langle a^\dagger a \rangle + 1$ respectively, corresponding to the emission and absorption spectra (Fig. 4a).

In Fig. 3, the difference between $S_{xx}^g[\omega_m]$ and $S_{xx}^e[\omega_m]$ is largely due to classical backaction-cooling for $|g\rangle$ and amplification for $|e\rangle$. In other words, the equilibrium occupation $\langle a^\dagger a \rangle$ is different in the two experimental sequences. To mitigate this effect, and isolate the small quantum asymmetry, we implemented an experimental sequence where the qubit is alternatively prepared in $|g\rangle$ and $|e\rangle$ (see Fig. 4(b)). In this manner, the average energy exchanged with the qubit cancels out in a coarse-grained approximation. The spectra $S_{mm}^{g/e}$ are then computed from the corresponding subsequences $\{m_k^{g/e}\}$ (see Fig. 4b).

The experimental spectra $S_{mm}^{g/e}$, represented in blue and red in Fig. 4c, exhibit a small but statistically significant difference in their areas. To confirm that this offset
has a quantum origin, we vary the cryostat temperature, $T_{\text{cryo}}$, within the interval [10 mK, 30 mK]. The spectra are normalized using a calibration peak and corrected by subtracting the sampling noise background. The mechanical mode's thermalization with the cryostat is evidenced by the linear dependence of the mean peak area. From this linear dependence, we infer the mechanical mode temperature as $T_m = T_{\text{cryo}} + 2.2$ mK, allowing us to recalibrate the mode area in units of quanta (right axis of Fig. 4(d)). The constant offset of 2.2 mK observed can be attributed either to a miscalibration of the generic Bluefors temperature sensor, additional acoustic noise transduced into the mechanical mode, or a slight temperature gradient between the sensor location near the mixing chamber and the sample holder. Most notably, the area difference between $S_{mm}^g$ and $S_{mm}^e$ remains subtle yet statistically significant across the entire temperature range, with an average of 1.25(25) phonons. This experimental offset—slightly above the ideal quantum asymmetry of exactly 1 phonon—agrees well with the 1.37-phonon prediction of our full theoretical model, which incorporates finite state-preparation fidelity and qubit-lifetime asymmetry (see App. F 5).

To conclude, we have demonstrated, for the first time, direct resonant coupling between a superconducting qubit and a macroscopic mechanical membrane oscillating at 4 MHz—over two orders of magnitude below previous realizations [20]. By employing a heavy-fluxonium architecture, we bridge the gap between the GHz-scale of superconducting circuits and the low-frequency mechanical domain, while preserving coherence. Exploiting the membrane’s long lifetime and periodic weak measurements, we reconstruct its noise spectrum and reveal a clear asymmetry between emission and absorption—a direct signature of the non-commutativity of phonon creation and annihilation operators.

Importantly, operating a 5.3 ng membrane at 4 MHz brings us into the regime where a macroscopic spatial superposition would decohere via the Diósi–Penrose mechanism on timescales comparable to the membrane’s intrinsic coherence [7]. Looking ahead, this platform naturally lends itself to quantum and adaptive feedback protocols—using membrane drives or conditional qubit preparations informed by real-time readout—to achieve active stabilization, state steering, and generation of non-classical mechanical states [44]. Several optimizations could improve the device performance. Phononic-crystal engineering should boost the electromechanical cooperativity by three orders of magnitude. Increasing the DC bias above 11.5 V—by mitigating leakage currents—reducing the membrane-qubit separation to sub-µm scales, and applying standard surface passivation to extend the qubit coherence time $T_2$ will together maximize the coupling–coherence product $\Omega T_2$. These measures should close the current order-of-magnitude gap and enable operation in the strong-coupling regime, thereby enabling experimental tests of gravitational collapse theories.

# ACKNOWLEDGMENTS

We thank Emilio Rui and Philippe Campagne-Ibarcq for insightful discussions. We thank Jean-Michel Raimond, Antoine Heidmann, Pierre-François Cohadon, and Tristan Briant for carefully reading the manuscript. This work was supported by the Agence Nationale de la Recherche under projects MECAFLUX (ANR-21-CE47-0011), FERBO (ANR-23-CE47-0004), and RobustSuperQ (ANR-22-PETQ-0003). This work has been supported by Region Île-de-France in the framework of DIM SIRTEQ (project CryoParis) and DIM QuanTiP (project COCONUT), and by Sorbonne Université through the HYQUITECH “Emergence” program. K. G. acknowledges support from the Quantum Information Center Sorbonne (QICS doctoral fellowship), H. P. is funded by the CNRS–University of Arizona joint Ph.D. program.

# AUTHOR CONTRIBUTIONS

K.G., L.B. and R.R. jointly designed and simulated the qubit sample, L.B., K.G. and P.M. micro-fabricated the sample. K.G., S.D. and R.R. performed measurements, K.G. and R.R. performed data analysis, R.R. and E.F. developed the theoretical framework for sideband asymmetry. H.P. and T.J. developed the flip-chip process and fabricated the membrane. E.F, Z.L. and W.C.S. provided essential guidance on the experimental implementation, contributed key conceptual ideas, and critically read and revised the manuscript. K.G., R.R., T.J. and S.D. contributed equally to writing the article, with contributions from all authors. S.D and T.J. supervised the project.

# COMPETING INTERESTS

The authors declare no competing interests.

# DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding authors upon reasonable request.
# SUPPLEMENTARY MATERIALS

## CONTENTS

A. Device Parameters 8

B. Expected Diósi–Penrose collapse time 8

C. Device and Fabrication 10

1. Qubit Design 10

2. Qubit Chip Fabrication 10

3. Membrane Design 10

4. Membrane Fabrication 10

5. Flip-Chip Assembly 12

D. Hamiltonian derivation 12

1. Effective mass 12

2. Circuit quantization 13

3. Qubit–membrane coupling 14

4. Electrostatic spring softening 14

5. Offset voltage discussion 15

6. Qubit-membrane distance estimation 15

E. Characterization of the qubit 16

1. Qubit spectroscopy 16

2. Qubit AC-Stark shift 16

3. Decoherence and relaxation times 17

4. Qubit reset 18

5. Qubit Preparation Fidelity 18

6. Qubit readout fidelity 19

7. Flux stabilization 20

F. Spectrum analyzer theory 20

1. Quantum harmonic oscillator spectrum 20

2. System definition 21

3. Semi-classical treatment 21

4. Quantum treatment 22

5. Effect of the qubit imperfections 24

G. Experimental Setup 26

References 26

## Appendix A: Device Parameters

Table I summarizes the key device parameters; the fourth column indicates the characterization method—often citing the specific section of this document where each measurement is described. During the interval between the measurements for Fig. 3 and Fig. 4, an unexpected temperature rise in the dilution refrigerator reduced the mechanical lifetime from $T_1^m = 5.8$ ms to $T_1^m = 1.47$ ms. Concurrently, the interrogation time was shortened to 4 µs, yielding a single phonon rotation angle of $\theta_1 = 38(2)$ mrad.

## Appendix B: Expected Diósi–Penrose collapse time

In the Diósi–Penrose (DP) framework [5, 6], a spatial superposition of a massive object acquires a finite lifetime due to the difference in its gravitational self-energy. For a resonator of mass *m* prepared in the "Schrödinger-cat" superposition

$$|\Psi\rangle = \frac{|\alpha\rangle + |-\alpha\rangle}{\sqrt{2}}, \qquad (B1)$$

one defines the gravitational self-energy difference

$$\begin{gathered} \Delta E = 2E_{12} - E_1 - E_2, \\ E_{ij} = -G \iint \frac{\rho_i(\mathbf{r})\rho_j(\mathbf{r}')}{|\mathbf{r} - \mathbf{r}'|} d^3\mathbf{r} d^3\mathbf{r}', \end{gathered} \qquad (B2)$$

where $\rho_i$ and $\rho_j$ are the mass-density distributions of the two components of the superposition and *G* the gravitational constant. The collapse time is then

$$\tau_G = \frac{\hbar}{\Delta E}. \qquad (B3)$$

Assuming the mass is concentrated in an ensemble of spherical nuclei of radius *a* (for silicon, $a \approx 2.7$ fm), an essential experimental prerequisite for testing DP-induced collapse is that the two displaced wavepackets be separated by $\Delta X > 2a \sim 5$ fm, with the spatial separation:

$$\begin{gathered} \Delta X = 2|\alpha| X_{zpf}, \\ X_{zpf} = \sqrt{\frac{\hbar}{2m\omega_m}} \sim 1 \text{ fm.} \end{gathered} \qquad (B4)$$

Importantly, owing to its low frequency, our mechanical system exhibits comparatively large zero-point fluctuations ($X_{zpf} \sim 1$ fm), so that only a modest coherent amplitude of $|\alpha| \approx \sqrt{6}$ ($\sim 6$ phonons) meets this criterion. Building on established circuit-QED cat-state protocols [45], our tunable resonant coupling between a superconducting qubit and the 4 MHz mechanical mode opens the path to preparation and measurement of such few-phonon Schrödinger-cat states.

Once the separation criterion (B4) is satisfied, following the reasoning of ref. [7], the DP-self energy difference can be calculated explicitly, yielding

$$\tau_G = \frac{1}{m} \frac{5\hbar a}{48\pi G m_a}, \qquad (B5)$$

with $m_a$ the atomic mass. Using the parameters of Silicon ($m_a = A m_u$, with $A = 28$ and $m_u = 1.7 \times 10^{-27}$, the atomic mass unit), we obtain $\tau_G \sim 0.5$ ms.

In a thermal environment, with $\bar{n}_{th} \approx 50$ (at 10 mK, 4 MHz) and mechanical decay rate $\kappa_m = 1/T_1^m \approx 170$ s$^{-1}$, the thermal decoherence time of the cat is

$$\tau_{\text{cat}} = \frac{\tau_{\text{th}}}{8|\alpha|^2} = \frac{T_1^m}{4\bar{n}_{\text{th}}|\alpha|^2}, \qquad (B6)$$
<table><caption>TABLE I. Key device parameters, grouped by Membrane, Qubit, and Coupling.</caption><thead><tr><th>Parameter</th><th>Description</th><th>Value</th><th>Origin</th></tr></thead><tbody><tr><td colspan="4" style="text-align: center; font-weight: bold;">Membrane</td></tr><tr><td>ω<sub>m</sub>/2π</td><td>Membrane frequency</td><td>4.4 MHz</td><td>Membrane spectroscopy Fig. 2(a)</td></tr><tr><td>M</td><td>Physical mass</td><td>5.3 ng</td><td>Appendix D 1</td></tr><tr><td>m</td><td>Effective mass</td><td>2.3 ng</td><td>Appendix D 1</td></tr><tr><td>X<sub>zpf</sub></td><td>Zero-point fluctuations</td><td>0.9 fm</td><td>Appendix D 1</td></tr><tr><td>T<sub>1</sub><sup>m</sup></td><td>Mechanical lifetime</td><td>5.9 ms</td><td>Ringdown experiment Fig. 2(b)</td></tr><tr><td>Q</td><td>Mechanical quality factor</td><td>1.62 × 10<sup>5</sup></td><td>Ringdown experiment Fig. 2(b)</td></tr><tr><td>τ<sub>G</sub></td><td>Expected gravitational collapse time</td><td>0.5 ms</td><td>Appendix B</td></tr><tr><td>τ<sub>th</sub></td><td>Thermal decoherence time</td><td>0.3 ms</td><td>Appendix B</td></tr><tr><td colspan="4" style="text-align: center; font-weight: bold;">Qubit</td></tr><tr><td>E<sub>J</sub>/h</td><td>Josephson energy</td><td>4.82(6) GHz</td><td>Spectroscopy fit Sec. E 1</td></tr><tr><td>E<sub>L</sub>/h</td><td>Superinductance energy</td><td>0.128(7) GHz</td><td>Spectroscopy fit Sec. E 1</td></tr><tr><td>E<sub>C</sub>/h</td><td>Charging energy</td><td>0.418(9) GHz</td><td>Spectroscopy fit Sec. E 1</td></tr><tr><td>ω<sub>q</sub>/2π</td><td>Qubit |g⟩ ↔ |e⟩ frequency</td><td>2.35 MHz</td><td>Ramsey Fig. 11</td></tr><tr><td>T<sub>1</sub></td><td>Qubit T<sub>1</sub> (bare)</td><td>22.4(8) μs</td><td>Relaxation Fig. 11</td></tr><tr><td>T<sub>2</sub></td><td>Qubit T<sub>2</sub> (bare)</td><td>17.1(1) μs</td><td>Ramsey Fig. 11</td></tr><tr><td>T<sub>1</sub><sup>st</sup></td><td>Qubit T<sub>1</sub> (Stark-shifted)</td><td>7.4(2) μs</td><td>Relaxation Fig. 11</td></tr><tr><td>T<sub>2</sub><sup>st</sup></td><td>Qubit T<sub>2</sub> (Stark-shifted)</td><td>4.2(0) μs</td><td>Ramsey Fig. 11</td></tr><tr><td>η<sub>g</sub></td><td>Preparation fidelity |g⟩</td><td>98.5(2.9)%</td><td>Appendix App. E 5</td></tr><tr><td>η<sub>e</sub></td><td>Preparation fidelity |e⟩</td><td>98.3(1.0)%</td><td>Appendix App. E 5</td></tr><tr><td>F<sub>g</sub></td><td>Readout fidelity |g⟩</td><td>84.8(1.8)%</td><td>Appendix App. E 6</td></tr><tr><td>F<sub>e</sub></td><td>Readout fidelity |e⟩</td><td>69.1(1.2)%</td><td>Appendix App. E 6</td></tr><tr><td colspan="4" style="text-align: center; font-weight: bold;">Coupling</td></tr><tr><td>d</td><td>Membrane–qubit gap distance</td><td>2.5(3) μm</td><td>See Sec. D 6</td></tr><tr><td>C<sub>m</sub>(d)</td><td>Membrane–qubit capacitance</td><td>13.9 pF</td><td>Ansys simulation</td></tr><tr><td>Ω/2π</td><td>Qubit–membrane vacuum rabi rate</td><td>1.50(8) kHz</td><td>Sec. F 5</td></tr><tr><td>θ<sub>1</sub></td><td>Single-phonon rotation angle</td><td>47(3) mrad</td><td>From Ω</td></tr><tr><td>C<sub>q</sub></td><td>Qubit shunt capacitance</td><td>5.7 pF</td><td>Ansys simulation</td></tr><tr><td>C<sub>g</sub></td><td>Qubit–ground capacitance</td><td>44.3 pF</td><td>Ansys simulation</td></tr></tbody></table>
giving $\tau_{cat} \sim 5$ µs for $|\alpha|^2 \approx 6$. Consequently, gravitationally induced decoherence would contribute a few percent of the total decoherence, making it experimentally resolvable.

Employing the technique introduced in ref. [40], our group has demonstrated soft-clamped silicon nitride membranes—with defect sizes, fundamental-mode frequencies (~ 2.5 MHz), and aluminum metallization comparable to our current resonator—achieving $T_1^m \approx 2$ s [46]. This close correspondence suggests that replacing the current membrane with such a device would directly extend the mechanical lifetime to the seconds scale. Integrating such low-loss oscillators into our platform would yield $\tau_{cat} \sim 3\tau_G$, fully satisfying both $\Delta x > 5$ fm and $\tau_{cat} > \tau_G$ conditions for a decisive test.

# Appendix C: Device and Fabrication

## 1. Qubit Design

Photographs of the device are given in Fig. 1(e). The fluxonium qubit is realized as a superconducting loop containing a small Josephson junction in series with a large superinductance, with two coplanar capacitor pads providing the shunt capacitance across the junction-inductor branch. The qubit is manipulated via a direct waveguide that is capacitively coupled to one of the electrodes. It is measured via a readout line that is capacitively coupled to the other electrode. The readout line consists of a $\lambda/2$ waveguide resonator that is interrupted in the middle, forming an inductive coupling. This configuration creates a readout cavity at 5.59 GHz and Purcell filter at 5.92 GHz that protect the qubit from purcell losses via the readout line. Slow changes in the external flux are applied using an external coil placed on the chip holder, whereas small flux adjustments, which are important for locking onto $\phi_{ext} = \pi$, are made via an on-chip flux line that is inductively coupled to the superinductance loop. The membrane is biased via a waveguide located on the qubit chip.

## 2. Qubit Chip Fabrication

The qubit chip (shown in Fig. 1(e)) is fabricated on a 280-µm-thick silicon substrate. All superconducting elements, except for the junctions, are formed from a 200 nm niobium layer deposited by sputtering. The film was then patterned at the wafer scale by UV-laser lithography using S1805 photoresist and finally etched to define all superconducting elements. The gap between the qubit chip and the membrane is set by 600-nm-high aluminum pillars, which are not connected to the qubit ground. These pillars are fabricated using a bilayer of LOR20B and AZ5214 E photoresists, followed by UV lithography and aluminum evaporation.

All Josephson junctions are fabricated using the Dolan bridge technique. A bilayer of MMA EL13 and PMMA A3 is deposited and patterned by electron-beam lithography. After development, aluminum is deposited using a Plassys evaporator with double-angle evaporation (±22°). Unwanted aluminum is removed via a lift-off process in an NMP solution.

The wafer is then diced into 9.8×10.8 mm chips, which are then mounted on a printed circuit board and wirebonded.

## 3. Membrane Design

The membrane design is shown in Fig. 6(a). It consists of a free-standing 90-nm-thick silicon nitride rectangular membrane clamped at its edges to a silicon substrate. To bias the membrane, a square bias pad is positioned on the substrate below the membrane. The only metalized areas are the bias pad, the membrane, and the connection between them.

The mechanical modes of this device are measured at room temperature with a Michelson interferometer used as a position-meter. The sensitivity of the interferometer is sufficient to monitor the membrane out-of-plane motion due to thermal noise. The fundamental mode was observed at 2.54 MHz, the first asymmetric mode at 3.75 MHz, the second asymmetric mode, used in this work was observed at 4.30 MHz in this room-temperature measurement. Optical ringdown measurements—where the membrane is driven by modulating the radiation-pressure force at its resonance frequency—were not possible after evaporating the aluminum electrodes. Laser absorption in the metal during the driving phase induced rapid frequency drifts that prevented stable resonant excitation and hence ringdown acquisition. Instead, we characterized the quality factor $Q$ before metallization, obtaining $Q \approx 60\,000$ for the mechanical mode studied.

## 4. Membrane Fabrication

The membranes are fabricated from a 350-µm-thick (100) silicon wafer with a 100 nm silicon nitride layer deposited on both sides under a uniform in-plane stress $\ge 800$ MPa. First, the back side is masked with AZ5214E photoresist. Then UV lithography and reactive ion etching (RIE) of the SiN layer is performed. To release the membranes, we etch through the whole substrate from the backside opening: the wafer is placed in a PEEK holder with a seal to protect the front side and immersed in a 30 % KOH solution at 85 °C. At the end of the process, we obtain a free-standing SiN membrane. A view of the backside, after etching, is shown in Fig. 6(c). A similar lithography/RIE process is applied to the front side, etching the SiN layer except in the membrane and bias pad regions. The sample is then etched in KOH to remove 25 µm of silicon (recess seen on the front side
<figure><img src="image_6.png" /><figcaption>FIG. 5. **Schematic of the setup.** A detailed wiring schematic of the experimental setup, showing both the RF and DC lines inside the cryostat and the control electronics at room temperature.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>FIG. 6. (a) Top: suspended SiN rectangular membrane in green (only the non metalized edge is visible). Aluminum layer in gold. Silicon substrate in blue. Bottom: square biasing electrode lying on the substrate. (b) Optical microscope image of the front-side of the final device. The suspended membrane and the biasing electrode sit on a silicon mesa, the typical rough etch bottom, arising from KOH etch can be seen at the edges of the picture. (c) Optical image of the same sample, seen from the back side. The released metalized SiN membrane can be seen at the top. Mirror images on the four edges of the membrane come from the smooth crystalline planes obtained after etching with KOH trough the substrate.</figcaption></figure>

in Fig. 6(b)). This serves to increase the gap between chips, thereby facilitating the flip-chip process. Finally, the wafer is cleaned in a 10% HF solution for 5 minutes, cleaning and thinning the SiN membrane to a 90 nm thickness. A final UV-laser lithography step, followed by evaporation of 30 nm of aluminum and liftoff, defines the biasing electrodes.

## 5. Flip-Chip Assembly

The flip-chip assembly is performed using an MJB4 mask aligner for micropositioning and a custom-built holder. The chips are aligned with an estimated accuracy of 5 µm and brought into contact. Bonding is achieved using drops of Dymax OP-67-LC UV-curing glue at the corners of the membrane chip. The aluminum spacers on the qubit chip nominally set the inter-chip gap to 600 nm. While room-temperature Fizeau interferometry on unmetalized test samples confirmed this spacing, cryogenic qubit spectroscopy consistently reveals a gap of 1–3 µm, which we attribute to differential bending under thermo-mechanical stress.

# Appendix D: Hamiltonian derivation

## 1. Effective mass

In this section, we derive the effective mass of the vibrational mode. The membrane is a rectangular sus-

pended film of width $L_y = 110$ µm, length $L_z = 140$ µm, and thickness $t = 90$ nm. Including the added mass of the 30 nm-thick aluminum electrode with area 90 µm × 120 µm, the physical mass $M = 5.3$ ng is obtained using the mass densities of silicon nitride $\rho_{SiN} = 3200$ kg/m³ and aluminum $\rho_{Al} = 2700$ kg/m³. The membrane supports a family of out-of-plane vibrational modes, indexed by the number of antinodes $(m, n)$ in the lateral and longitudinal directions, with the displacement profile $u_{m,n}(y, z)$ given by

$$ \lambda\epsilon(y, z) \sin\left(m\pi \frac{y-L_y/2}{L_y}\right) \sin\left(n\pi \frac{z-L_z/2}{L_z}\right). \quad (D1) $$

Here, $\epsilon(y, z)$ is a regularization function which, for thin films under high tensile stress $\sigma_{SiN}$, equals unity everywhere except in an exponentially small region near the clamping edges $y = \pm L_y/2$ and $z = \pm L_z/2$ [36]. While this factor is crucial to capture mechanical damping in the membrane [36], it has negligible effect on the motional mass and will henceforth be set to unity. As the mode profiles (D1) are solutions of the linear elasticity equations, there is a freedom in choosing the normalization constant $\lambda$, leading to an ambiguity in the motional mass $m_{\lambda}$. The kinetic energy of a single (1,2)-mode with amplitude $X$ is given by

$$ T = \frac{1}{2} \dot{X}^2 \iint \rho_s u_{1,2}^2(y, z) dy dz, \quad (D2) $$

where $\rho_s$ is the surface mass density of the membrane. Comparing to a harmonic oscillator, one defines the effective mass

$$ m_{\lambda} = \iint \rho(y, z) u_{1,2}^2(y, z) dy dz = M \frac{\lambda^2}{4}, \quad (D3) $$

which leads to the usual kinetic energy $T = \frac{1}{2}m_{\lambda}\dot{X}^2$. While any normalization factor $\lambda$ is admissible provided the effective mass is redefined according to Eq. D3, a natural choice, given the electrostatic coupling scheme employed in this work, is to choose $\lambda$ such that

$$ \iint_A u_{2,1}(y, z) dy dz = A, \quad (D4) $$

where $A$ denotes the area of overlap in each of the capacitors $C_{m}^{\pm}$. This choice validates the parallel-plate approximation for small displacements $X$, yielding

$$ C_{m}^{\pm}(X) \approx C_{m}(1 \mp X/d). \quad (D5) $$

Indeed,

$$ \begin{aligned} C_{m}^{\pm}(X) &= \iint_A \frac{\epsilon_0}{d \pm X} u_{2,1}(y, z) dy dz \\ &\approx \frac{\epsilon_0 A}{d} \left( 1 \mp \frac{X}{dA} \iint_A u_{2,1}(y, z) dy dz \right). \end{aligned} \quad (D6) $$

With this normalization one finds $\lambda = 1.3$ and hence

$$ m_{\lambda} = M \frac{\lambda^2}{4} = 2.3 \text{ ng.} $$
<figure><img src="image_8.png" /><figcaption>FIG. 7. Equivalent electrical circuit of the flip-chip device. Schematic diagram showing the equivalent circuit and defining the component symbols and node labels used throughout the text.</figcaption></figure>

This effective mass corresponds to the mass of the vibrational mode as seen by the fluxonium electrodes. The second-order term in the expansion (D6) differs from the parallel-plate approximation by 5.5% and will be neglected in further analytical calculations.

## 2. Circuit quantization

Qubit-membrane coupling is achieved via the two capacitors $C_m^+$ and $C_m^-$ [19], formed between the qubit's conducting pad and the metallized membrane (illustrated in Fig. 1(c)). The fluxonium qubit consists of a floating superconducting loop that is symmetrically coupled to ground through two identical capacitors $C_g$. The electrical schematic is shown in Fig. 7, where any additional parasitic capacitances have been neglected. The dynamics are described by the membrane displacement coordinate $X$ and the vector of flux-node coordinates

$$ \vec{\Phi}^T = (\Phi_1, \Phi_2, \Phi_3), \qquad (D7) $$

where $\Phi_1$ is the flux across the qubit, and $\Phi_{2,3}$ denote the fluxes at the coupling nodes (see Fig. 7).

We express the total Lagrangian as the sum of mechanical and electrical contributions:

$$ \mathcal{L} = \mathcal{L}_{\text{mec}}(X, \dot{X}) + \mathcal{L}_{\text{elec}}(\Phi_1, \dot{\vec{\Phi}}, X), \qquad (D8) $$

where overdots indicate time derivatives. The membrane Lagrangian is that of a harmonic oscillator:

$$ \mathcal{L}_{\text{mec}} = \frac{1}{2} m \dot{X}^2 - \frac{1}{2} m \omega_m^2 X^2, \qquad (D9) $$

with $m$ the effective mass of the membrane and $\omega_m$ its resonance frequency. The electrical Lagrangian is split

into kinetic and potential parts, $\mathcal{L}_{\text{elec}} = T_{\text{elec}} - V_{\text{elec}}$:

$$ \begin{aligned} T_{\text{elec}} = & \frac{C_g}{2} (\dot{\Phi}_2 - \frac{1}{2}\dot{\Phi}_1)^2 + \frac{C_g}{2} (\dot{\Phi}_2 + \frac{1}{2}\dot{\Phi}_1)^2 \\ & + \frac{C_m^+}{2} (\dot{\Phi}_3 - \dot{\Phi}_2 - \frac{1}{2}\dot{\Phi}_1)^2 + \frac{C_m^-}{2} (\dot{\Phi}_3 - \dot{\Phi}_2 + \frac{1}{2}\dot{\Phi}_1)^2 \\ & + \frac{C_b}{2} (V_b - \dot{\Phi}_3)^2 + \frac{C_q}{2} \dot{\Phi}_1^2, \end{aligned} \qquad (D10) $$

$$ V_{\text{elec}} = -E_J \cos(2\pi\Phi_q/\Phi_0) + \frac{\Phi_1^2}{2L}. \qquad (D11) $$

Here, $C_g$ is the ground capacitance, $C_b$ the membrane bias capacitance, $C_q$ the qubit shunt capacitance, $V_b$ the bias voltage, $E_J$ the Josephson energy, $\Phi_0 = h/2e$ the flux quantum, and $L$ the inductance of the fluxonium.

For the Legendre transform we rewrite $T_{\text{elec}}$ in matrix form:

$$ T = \frac{1}{2} \sum_{i,j=1}^{3} \dot{\Phi}_i M_{ij}(X) \dot{\Phi}_j - C_b V_b \dot{\Phi}_3, \qquad (D12) $$

with the position-dependent capacitance matrix

$$ M(x) = \begin{pmatrix} \frac{1}{2}(C_g + C_m^\Sigma + 2C_q) & C_m^\Delta & -C_m^\Delta \\ C_m^\Delta & 2(C_g + C_m^\Sigma) & -2C_m^\Sigma \\ -C_m^\Delta & -2C_m^\Sigma & C_b + 2C_m^\Sigma \end{pmatrix}, $$

with $C_m^\Sigma \equiv (C_m^+ + C_m^-)/2$ and $C_m^\Delta \equiv (C_m^+ - C_m^-)/2$ are symmetric and antisymmetric combinations of the membrane-qubit capacitances. The canonical flux-momenta are defined by $Q_k = \partial T/\partial \dot{\Phi}_k$, where $k = 1, 2, 3$. This yields

$$ \dot{\Phi}_k = \sum M_{ki}^{-1} \tilde{Q}_i, \quad \text{where} \quad \tilde{Q}_i = Q_i + C_b V_b \delta_{i3} \qquad (D13) $$

Here $\delta_{ij}$ is the Kronecker delta. The membrane momentum is simply $P_X = m \dot{X}$. Performing the Legendre transform, $H = \sum_{i=1}^3 Q_i \dot{\Phi}_i + P_X \dot{X} - \mathcal{L}$, yields the Hamiltonian

$$ H = \frac{1}{2} \sum_{i,j=1}^{3} \tilde{Q}_i M_{ij}^{-1} \tilde{Q}_j - E_J \cos(2\pi\Phi_1/\Phi_0) + \frac{\Phi_1^2}{2L} \\ + \frac{P_X^2}{2m} + \frac{1}{2} m \omega_m^2 X^2. \qquad (D14) $$

Since $\mathcal{L}$ (Eq. (D8)) does not depend on $\Phi_2$ or $\Phi_3$, these coordinates are cyclic and their conjugate charges are strictly conserved by Noether's theorem. Physically, these conserved quantities correspond to the static charges on the fluxonium and membrane islands, and may therefore be treated as fixed scalar parameters in the subsequent analysis.

We quantize the system by introducing the canonical variables for the qubit phase $\varphi$ and Cooper pair number $n$, along with the dimensionless canonical variables for
the membrane position $x$, momentum $p$ and $a$ ($a^\dagger$) the
mechanical annihilation (creation) operators:

$$
\varphi = 2\pi \frac{\mathbf{\Phi}_1}{\Phi_0}, \quad \mathbf{x} = \frac{\mathbf{X}}{2X_{\text{zpf}}} = \frac{\mathbf{a} + \mathbf{a}^{\dagger}}{2}, \tag{D15}
$$

$$
\mathbf{n} = \frac{Q_1}{2e}, \qquad \mathbf{p} = \frac{\mathbf{P}}{2P_{\text{zpf}}} = i \frac{\mathbf{a}^{\dagger} - \mathbf{a}}{2},
$$

where $X_{zpf} = \sqrt{\hbar/(2m\omega_m)}$ and $P_{zpf} = \sqrt{\hbar m\omega_m/2}$ are
the membrane zero-point fluctuations of position and
momentum, respectively. These variables satisfy the
commutation relations $[\varphi, n] = i$ and $[x, p] = i/2$,
and operators from different sets commute. The ca-
pacitance matrix depends on the membrane position,
$M_{ij}^{-1} = M_{ij}^{-1}(x)$, and must therefore be treated as an
operator. In this basis the Hamiltonian assumes the fa-
miliar form for a fluxonium coupled to a harmonic mode:

$$
\mathbf{H} = 4 \mathbf{E}_C(\mathbf{x}) [\mathbf{n} + \mathbf{n}_g(\mathbf{x})]^2 - E_J \cos \varphi_q + \frac{1}{2} E_L \varphi_q^2 \\
+ \hbar \omega_m (\mathbf{p}^2 + \mathbf{x}^2) + \mathbf{F}_m(\mathbf{x}), \qquad (\text{D16})
$$

where

$$
E_L = \left( \frac{\Phi_0}{2\pi} \right)^2 \frac{1}{L}, \quad E_C(\mathbf{x}) = \frac{e^2}{2} M_{11}^{-1}(\mathbf{x}), \\
\mathbf{n}_g(\mathbf{x}) = \frac{e}{4 E_C(\mathbf{x})} \sum_{i=2}^3 \tilde{Q}_i M_{1i}^{-1}(\mathbf{x}).
$$

This form of the Hamiltonian shows that the coupling
between the fluxonium and the membrane induces an
additional charge **n**<sub>g</sub>(**x**) on the fluxonium. The mem-
brane also experiences an additional electrostatic poten-
tial **H**<sub>ms</sub>(**x**), which, as will be shown later, leads to a
reduction of the membrane frequency. The general form
of **H**<sub>ms</sub> is given by

$$
\mathbf{H}_{\text{ms}} = 4\mathbf{E}_C \mathbf{n}_g^2 + \mathbf{F}_m = \frac{1}{2} \sum_{i,j=2}^{3} \tilde{Q}_i \mathbf{M}_{ij}^{-1} \tilde{Q}_j. \quad (\text{D17})
$$

3. Qubit-membrane coupling

For small membrane displacements we perform a first-
order expansion of the coupling Hamiltonian in the di-
mensionless coordinate $\mathbf{x}$, yielding

$$
\begin{align}
H_{\text{coupling}} &= 8 \mathbf{E}_C(\mathbf{x}) \mathbf{n}_g(\mathbf{x}) \nonumber \\
&\approx 4 E_C \frac{\partial C_m^+}{\partial x} \frac{V_{\text{eff}}(V_b)}{e} \mathbf{n} \mathbf{x}, \tag{D18}
\end{align}
$$

Here, $E_C = e^2/(C_g+C_m+2C_q)$ is the unperturbed capac-
itive energy of the fluxonium (neglecting all terms beyond
linear order in $x$), and

$$
\frac{\partial C_m^{\pm}}{\partial x} = 2 X_{\text{zpf}} \frac{\partial C_m^{\pm}}{\partial X}
$$

relates the derivative with respect to the dimensionless
coordinate **x** back to the physical membrane displace-
ment **X**. The effective bias, $V_{\text{eff}}(V_b)$, is given by:

$$
\begin{equation}
\begin{split}
V_{\text{eff}}(V_b) ={}& \frac{C_g V_b - Q_2/2 + (C_g/C_b) Q_3}{C_g + C_m + 2 (C_g C_m/C_b)} \\
& \rightarrow \frac{C_g V_b - Q_2/2}{C_g + C_m},
\end{split}
\tag{D19}
\end{equation}
$$

where the last expression indicates the limit $C_g/C_b \ll 1$.
This indicates that trapped charges, $Q_2$ and $Q_3$, act as
an effective bias on the qubit. Introducing the dilution
factor $\beta = C_g/(C_g+C_m)$ and the offset voltage $V_{\text{offset}} =
Q_2/2C_g$, we have

$$
V_{\text{eff}}(V_b) = \beta(V_b - V_{\text{offset}}) \quad (\text{D20})
$$

Projecting onto the lowest two levels of the fluxonium
and invoking the rotating-wave approximation yields the
Jaynes–Cummings Hamiltonian,

$$
\mathbf{H}_{JC}/\hbar = \omega_m \mathbf{a}^{\dagger}\mathbf{a} + \frac{\omega_q}{2} \boldsymbol{\sigma}_z + i \frac{\Omega}{2} (\mathbf{a}^{\dagger}\boldsymbol{\sigma}_{-} - \boldsymbol{\sigma}_{+}\mathbf{a}), \quad (\text{D21})
$$

where the effective vacuum Rabi frequency is

$$
\Omega = \omega_q |\langle g | \varphi | e \rangle| \frac{dC_m^+}{dx} \frac{V_{\text{eff}}}{2e}. \quad (\text{D22})
$$

Here, $\omega_q$ is the qubit transition frequency and $\boldsymbol{\sigma}_{\pm}$ are the qubit raising and lowering operators; the particular phase convention for $\varphi_q$ has been chosen to simplify later expressions.

Using the independently determined parameters
$\omega_q/2\pi = 2.35$ MHz (Ramsey experiment), $|\langle g|\phi|e\rangle| = 3.04$ (numerical evaluation using parameters from Sec. 9),
$V_{\text{eff}} = \beta(V_b - V_{\text{offset}}) = 5.6$ V (Ansys simulation for
$\beta$), and $dC_m^+/dx = 2X_{\text{zpf}}C_m/d$, together with $C_m(d = 2.5 \mu\text{m}) = 13.9$ fF (Ansys simulation) and $X_{\text{zpf}} = 0.9$ fm,
we obtain

$$
\Omega^{\text{th}}/2\pi = 1.31(33) \text{ kHz}, \quad (\text{D23})
$$

A value that is in good qualitative agreement with the
value $\Omega/2\pi = 1.50(8)$ kHz extracted directly from the
signal-to-noise ratio of the spectrum analyzer experiment
(see Eq. F44).

4. Electrostatic spring softening

A bias voltage $V_b$ on the membrane produces an elec-
trostatic potential **H**<sub>ms</sub> that effectively softens the me-
chanical spring and lowers the resonant frequency. To
capture this effect, we expand **H**<sub>ms</sub> to second order in
the dimensionless displacement **x**, express it in terms of
the ladder operators and retain only energy-conserving
terms. The membrane Hamiltonian then becomes

$$
\begin{align}
H_{\text{membrane}} &= \hbar \omega_m (\mathbf{a}^{\dagger} \mathbf{a} + \frac{1}{2}) + \mathbf{H}_{\text{ms}} \tag{D24} \\
&= \hbar (\omega_m + \delta\omega_m) (\mathbf{a}^{\dagger} \mathbf{a} + \frac{1}{2}), \nonumber
\end{align}
$$
where the frequency shift is

$$ \delta\omega_m = -\zeta(V_b - V_{\text{offset}})^2 \qquad (\text{D25}) $$

with

$$ \zeta = 2 C_m (C_g + 2C_q) \frac{E_C}{\hbar e^2} \frac{X_{zpf}^2}{d^2} \beta^2. \qquad (\text{D26}) $$

Here we have used the parallel-plate model $C_m^\pm(X) = C_m/(1 \pm X/d)$ to calculate the derivatives $dC_m^\pm/dX$ and $d^2C_m^\pm/dX^2$.

The electrostatic force thus acts as a “negative spring”, lowering the membrane’s resonance frequency. In Fig. 2a we plot this softening and extract a tuning coefficient

$$ \zeta = -0.980 \, \text{Hz/V}^2. $$

By combining Eq. D26 with the capacitance $C_m(d)$ extracted from our ANSYS Q3D simulations and the gap $d$, independently determined (with its uncertainty) in Sec. D6 (Table I), we obtain

$$ \zeta_{\text{th}} = -1.14(42) \, \text{Hz/V}^2. $$

Because $\zeta \propto d^{-3}$, this parameter-free agreement between $\zeta$ and $\zeta_{\text{th}}$ provides a stringent consistency check on our membrane–qubit separation.

## 5. Offset voltage discussion

In our measurements, we consistently observed a residual static charge on the membrane, which offsets the voltage bias required to cancel the electrostatic force by an amount $V_{\text{offset}}$ (see Fig. 2a). This offset arises from trapped charges $Q_2$ and $Q_3$ on the islands between the capacitors (see Eq. D19). Although the qubit remains operational for a DC-bias voltage up to $|V_b| \lesssim 9$ V, the effective coupling strength depends on the difference $V_b - V_{\text{offset}}$. Therefore, the presence of this offset can be used to increase the coupling strength.

We found that $V_{\text{offset}}$ slowly drifts toward the applied bias over a timescale of weeks. To restore the maximal coupling range, we deliberately “reset” the offset by applying a large negative bias of $V_b = -20$ V to the membrane for several days. While this high voltage temporarily disables measurements, it reliably shifts $V_{\text{offset}}$ closer to $-9$ V. After this procedure, we verify the new offset value using Ramsey-type membrane spectroscopy (see Fig. 2a) and ensure that it lies within the desired measurement range.

## 6. Qubit-membrane distance estimation

The value of the vacuum-gap capacitors $C_m^\pm$ depends sensitively on the separation distance $d$ between the fluxonium and the membrane. To estimate $d$ accurately, we measured the same fluxonium device both before

<figure><img src="image_9.png" /><figcaption>FIG. 8. Distance estimation. Empirical approximation of the capacitance found in Ansys Q3D simulations $C_m(d) = \epsilon_0 \times 54 \times 10^{-9}/d^{0.8}$, where both $\epsilon_0$ and $d$ are in SI units.</figcaption></figure>

and after the flip-chip assembly. Plasmon transitions ($\omega_{ef}, \omega_{gh} \propto E_C^{1/2}$) are only weakly dependent on the capacitive energy and are significantly shifted by coupling to nearby array modes. In contrast, the low-frequency qubit transition $\omega_q$ scales exponentially with $E_C$, making it a more precise probe of the capacitor gap. At half a flux quantum bias, $\omega_q$ follows the approximation [47, 48]:

$$ \omega_{ge} = \frac{8 \cdot 2^{\frac{3}{4}}}{\sqrt{\pi}} E_J^{\frac{3}{2}} E_C^{\frac{1}{4}}(d) \exp\left[-\sqrt{\frac{8E_J}{E_C(d)}} + \frac{14\zeta(3)E_L}{\sqrt{8E_J E_C(d)}}\right]. \qquad (\text{D27}) $$

Here, $\zeta(3) \approx 1.2$ is the Riemann zeta function.

We first record the qubit frequency without the membrane as $\omega_q^{\text{BFC}} = 4.93$ MHz. Inverting Eq. D27 for known $E_J$ and $E_L$ yields $E_C(\infty)$. Next, we simulate the qubit–membrane capacitance $C_m(d)$ versus gap $d$ using Ansys Q3D (see Fig. 8(b)) Applying the definition from Eq. D18, we compute

$$ E_C(d) = \frac{e^2}{\frac{e^2}{E_C(\infty)} - C_m(d)}. $$

Substituting $E_C(d)$ back into Eq. D27 gives $\omega_q(d)$ for fixed $E_J$, $E_L$, and $\omega_q^{\text{BFC}}$. While $\omega_q^{\text{BFC}}$ is measured with high precision, we assign a 10 % uncertainty to $E_J$ and $E_L$ and take their values found from the qubit spectroscopy in Sec. E1. Assigning 10 % uncertainties to $E_J$ and $E_L$ and propagating them via a Monte Carlo analysis produces the shaded confidence band in Fig. 8(a). After flip-chip assembly, the qubit frequency drops to $\omega_q^{\text{AFC}} = 2.35$ MHz, from which we extract the gap distance $d = 2.5(3)$ µm.
<figure><img src="image_10.png" /><figcaption>FIG. 9. Two-tone spectroscopy of the heavy-fluxonium qubit: reflected readout amplitude versus probe frequency and external flux $\phi_{ext}$. The primary “diamond” at ~3.46 GHz corresponds to the $\{|g\rangle, |e\rangle\} \rightarrow \{|f\rangle, |h\rangle\}$ transitions. A secondary diamond near ~3.65 GHz, labeled by the transitions $\{|g, 0\rangle, |e, 0\rangle\} \rightarrow \{|g, 1\rangle, |e, 1\rangle\}$ (chain-mode Fock number in the ket), reveals coupling to a single chain mode. Solid lines are the best fit including that mode with parameters $E_J/h = 4.886$ GHz, $E_C/h = 0.408$ GHz, $E_L/h = 0.135$ GHz, $\omega_{chain}/2\pi = 3.650$ GHz, and $g_{qc}/2\pi = 197$ MHz. Dashed lines show an alternative fit fixing $\omega_{ge}/2\pi = 2.32$ MHz, yielding $E_J/h = 4.757$ GHz, $E_C/h = 0.427$ GHz, $E_L/h = 0.121$ GHz, $\omega_{chain}/2\pi = 3.641$ GHz, and $g_{qc}/2\pi = 203$ MHz.</figcaption></figure>

# Appendix E: Characterization of the qubit

## 1. Qubit spectroscopy

The heavy-fluxonium circuit, characterized by $E_J \gg E_C > E_L$, exhibits pronounced anharmonicity: its fundamental transition occurs at $\omega_q/2\pi = 2.2$ MHz, while next excited transitions lie in the GHz range. To probe these higher transitions, we perform two-tone spectroscopy: a variable-frequency probe tone is applied to the fluxonium, followed by a readout pulse that captures the reflected response. In conventional fluxonium devices, transitions between the lower manifold ($|g\rangle, |e\rangle$) and the upper manifold ($|f\rangle, |h\rangle$) produce a characteristic “diamond” pattern composed of four lines in the flux-dependent spectrum, with the $|g\rangle \rightarrow |f\rangle$ and $|e\rangle \rightarrow |h\rangle$ transitions forbidden by symmetry at $\phi_{ext} = \pi$. This pattern appears around 3.46 GHz (see Fig. 9). Our device exhibits, in addition, a second diamond-like structure around 3.7 GHz that the standard fluxonium Hamiltonian does not predict.

We attribute this extra structure to coupling between the fluxonium and a collective mode of the junction chain

employed as a superinductance [49]. This bosonic mode, at $\omega_{chain}$ is strongly coupled to the fluxonium, leading to a rich flux-dependent spectrum near $\phi_{ext} = \pi$. The Hamiltonian of the coupled system, written in terms of the fluxonium flux and charge operators ($\phi, n$) and the chain mode ladder operators ($b, b^\dagger$), is [50]

$$ H = 4E_C n^2 - E_J \cos(\phi - \phi_{ext}) + \frac{1}{2} E_L \phi^2 + \hbar \omega_{chain} b^\dagger b + \hbar g_{qc} (b^\dagger + b) n. \quad (E1) $$

We diagonalize this Hamiltonian numerically to obtain the eigenenergies and fit them to the experimental spectrum.

Beyond reproducing the extra diamond, the chain-mode coupling also reshapes the original $\{|g\rangle, |e\rangle\} \rightarrow \{|f\rangle, |h\rangle\}$ transitions of the fluxonium. A fit to the transitions $\{|g\rangle, |e\rangle\} \rightarrow \{|f\rangle, |h\rangle\}$ with the bare fluxonium model (not represented in the figure) yields a minimum gap of $\omega_q/2\pi = 0.9$ MHz at $\phi_{ext} = \pi$; including the chain-mode (full line fit of Fig. 9) raises this to $\omega_q/2\pi = 1.7$ MHz, in much closer agreement with the observed 2.35 MHz. Although incorporating additional chain modes could further refine the theory-experiment agreement, we have chosen to keep a single chain-mode in the fit to keep the model tractable. We also provide a best fit, where the g-e transition frequency is constrained around 2.35 MHz (dashed lines in Fig. 9). The quoted qubit parameters and uncertainties in the main text, $E_J/h = 4.82(6)$ GHz, $E_C/h = 0.418(9)$ GHz, $E_L/h = 0.128(7)$ GHz are taken as the averages and half-differences of the two fitting results.

## 2. Qubit AC-Stark shift

The lowest $|g\rangle \leftrightarrow |e\rangle$ transition of our fluxonium lies at $\omega_q/2\pi = 2.35$ MHz, while the mechanical mode sits at $\omega_m/2\pi = 4.4$ MHz. By applying a continuous drive at $\omega_d/2\pi = 3.463$ GHz—exactly midway between the $\omega_{ef}$ and $\omega_{gh}$ lines—we induce an AC-Stark shift that pulls the dressed qubit into resonance with the mechanics, at the cost of an increased dephasing rate (see Sec. E 3). To describe this driven system we therefore add a time-dependent charge-drive term to the fluxonium Hamiltonian:

$$ H = H_{\text{fluxonium}} + \hbar \Omega_d \cos(\omega_d t) n, \quad (E2) $$

where $H_{\text{fluxonium}}$ is the undriven fluxonium Hamiltonian and $\Omega_d$ is the coupling of the charge drive to the qubit. We restrict our attention to transitions near $\omega_d$ and retain only the lowest four fluxonium levels, $\{|g\rangle, |h\rangle, |e\rangle, |f\rangle\}$. Projecting the driving term of the Hamiltonian into this four-level subspace, and considering only non-zero matrix elements at $\phi_{ext} = \pi$ yields

$$ H_{\text{Stark}} = \hbar \Omega_d \cos \omega_d t \begin{pmatrix} 0 & n_{\text{gh}} & n_{\text{ge}} & 0 \\ n_{\text{gh}}^* & 0 & 0 & 0 \\ n_{\text{ge}}^* & 0 & 0 & n_{\text{ef}} \\ 0 & 0 & n_{\text{ef}}^* & 0 \end{pmatrix} \quad (E3) $$
<figure><img src="image_11.png" /><figcaption>FIG. 10. **Qubit AC-Stark shift.** (a) Ramsey pulse sequence with a continuous drive applied during the idle interval. (b) Energy-level diagram of the four lowest fluxonium eigenstates. The drive frequency $\omega_d$ is set midway between the $|e\rangle \to |f\rangle$ and $|g\rangle \to |h\rangle$ transitions, inducing level repulsion and opposite AC-Stark shifts on the $\{|g\rangle, |h\rangle\}$ and $\{|e\rangle, |f\rangle\}$ manifolds. (c) The qubit transition frequency as a function of the applied drive amplitude squared (proportional to $\Omega_d^2$). (d) Difference between the $|+\rangle$ and $|-\rangle$ measurement outcome probabilities along the $\sigma_x$ (blue) and $\sigma_y$ (gray) axes at the bare qubit frequency. From these data we reconstruct the complex transverse Bloch vector and fit it with exponentially decaying oscillations; the real and imaginary components are overlaid in yellow.</figcaption></figure>

where $n_{ij} = \langle i|n|j\rangle$ is the charge matrix element. We then transform into the frame rotating at $\omega_d$ via the unitary transformation

$$U = \mathrm{diag}(1 \ e^{i\omega_d t} \ e^{i\omega_{ge} t} \ e^{i(\omega_d + \omega_{ge})t}). \quad (E4)$$

In the rotated frame, dropping all fast rotating terms (beyond or close to $\omega_d$), the Hamiltonian writes

$$H^{\mathrm{int}} = \hbar \begin{pmatrix} 0 & \Omega_d n_{\mathrm{gh}}/2 & 0 & 0 \\ \Omega_d n_{\mathrm{gh}}^*/2 & \Delta_{\mathrm{gh}} & 0 & 0 \\ 0 & 0 & 0 & \Omega_d n_{\mathrm{ef}}/2 \\ 0 & 0 & \Omega_d n_{\mathrm{ef}}^*/2 & \Delta_{\mathrm{ef}} \end{pmatrix}, \quad (E5)$$

with $\Delta_{\mathrm{gh}} = \omega_{\mathrm{gh}} - \omega_d$ and $\Delta_{\mathrm{ef}} = \omega_{\mathrm{ef}} - \omega_d$. This Hamiltonian is block-diagonal, such that it can be diagonalized separately in the two subspaces $\{|g\rangle, |h\rangle\}$ and $\{|e\rangle, |f\rangle\}$. The two lowest energy dressed eigenstates are

$$|\tilde{g}\rangle = \cos \frac{\Theta_{\mathrm{gh}}}{2} |g\rangle - \sin \frac{\Theta_{\mathrm{gh}}}{2} |h\rangle, \quad (E6)$$

$$|\tilde{e}\rangle = \cos \frac{\Theta_{\mathrm{ef}}}{2} |e\rangle - \sin \frac{\Theta_{\mathrm{ef}}}{2} |f\rangle,$$

with mixing angles defined by $\tan \Theta_{\mathrm{gh}} = \Omega_d n_{\mathrm{gh}}/|\Delta_{\mathrm{gh}}|$ and $\tan \Theta_{\mathrm{ef}} = \Omega_d n_{\mathrm{ef}}/|\Delta_{\mathrm{ef}}|$. The eigenergies are given

by

$$\begin{align}
\tilde{E}_g/\hbar &= \frac{1}{2} \mathrm{sign}(\Delta_{\mathrm{gh}}) (|\Delta_{\mathrm{gh}}| - \sqrt{\Delta_{\mathrm{gh}}^2 + \Omega_d^2 n_{\mathrm{gh}}^2}), \tag{E7} \\
\tilde{E}_e/\hbar &= \frac{1}{2} \mathrm{sign}(\Delta_{\mathrm{ef}}) (|\Delta_{\mathrm{ef}}| - \sqrt{\Delta_{\mathrm{ef}}^2 + \Omega_d^2 n_{\mathrm{ef}}^2}).
\end{align}$$

For the circuit at half-flux, $\phi_{\mathrm{ext}} = \pi$, one finds $n_{\mathrm{gh}} \simeq n_{\mathrm{ef}} \equiv n_{\mathrm{up}}$, furthermore, we have chosen a drive frequency precisely midway between $\omega_{\mathrm{gh}}$ and $\omega_{\mathrm{ef}}$ yielding $\Delta_{\mathrm{gh}} = -\Delta_{\mathrm{ef}} \equiv \Delta_{\mathrm{st}}$, and $\Theta_{\mathrm{gh}} \simeq \Theta_{\mathrm{ef}} \equiv \Theta$. In this approximation, the AC-Stark shift is given by

$$\begin{align*}
\delta\omega_q &\equiv \tilde{E}_e - \tilde{E}_g \\
&= \sqrt{\Delta_{\mathrm{st}}^2 + \Omega_d^2 n_{\mathrm{up}}^2} - \Delta_{\mathrm{st}} \\
&\approx \frac{\Omega_d^2 n_{\mathrm{up}}^2}{2\Delta_{\mathrm{st}}} \left( 1 - \frac{1}{4} \frac{\Omega_d^2 n_{\mathrm{up}}^2}{\Delta_{\mathrm{st}}^2} \right).
\end{align*}$$

In the last line, we expanded the expression to fourth order in $\Omega_d/\Delta_{\mathrm{st}}$.

Since the absolute attenuation between source and sample is unknown, we calibrate the AC-Stark shift by Ramsey interferometry Fig. 10a. Beginning at the fluxonium sweet-spot frequency $\omega_q/2\pi \approx 2.2$ MHz, we apply a continuous tone at $\omega_d/2\pi = 3.463$ GHz with varying source power $P$. A $\pi/2$-Ramsey sequence maps the drive-induced phase into $\sigma_z$, from which we extract the dressed qubit frequency. Plotting this frequency versus $P$ reveals a primarily linear Stark pull plus a small quadratic correction (Fig. 10c). Because Eq. E8 fixes the ratio of the linear and quadratic terms, only one free parameter—the scaling factor between $P$ and the on-chip Rabi rate $\Omega_d$—is required to fit all the data. Over the accessible power range the qubit frequency tunes continuously from 2.2 MHz up to ~ 5 MHz (more than one octave), enabling resonance with the 4.4 MHz mechanical mode. The resulting single-parameter fit (solid line, Fig. 10c) quantitatively reproduces both slope and curvature of the Stark shift.

### 3. Decoherence and relaxation times

We measure the qubit relaxation and dephasing rates as functions of the Stark-shifted qubit frequency, $\omega_q^{\mathrm{st}}$. The energy-relaxation rate, $\Gamma_1 = 1/T_1$, is extracted from the exponential decay of the excited-state population. The free-induction-decay rate, $\Gamma_2^* = 1/T_2^*$, is obtained from a Ramsey experiment (Fig. 10a). The resulting rates are plotted in Fig. 11. At the bare qubit frequency we observe $\Gamma_1/\Gamma_2^* \simeq 1.12 < 2$. As the AC-Stark drive increases, the pure dephasing rate $\Gamma_\phi = \Gamma_2^* - \Gamma_1/2$ rises sharply, causing $T_2^*$ to degrade much faster than $T_1$. Both rates exhibit a pronounced peak when $\omega_q^{\mathrm{st}}$ approaches $\omega_m$, due to the qubit transitions induced by the thermally excited mechanical mode.

We attribute the approximately linear rise of $\Gamma_1$ with detuned qubit frequency $\omega_q^{\mathrm{st}}$ to the growing admixture
<figure><img src="image_12.png" /><figcaption>FIG. 11. **Qubit decoherence and relaxation rates.** Energy-relaxation rate $\Gamma_1 = 1/T_1$ (blue markers) and free-induction-decay rate $\Gamma_2^* = 1/T_2^*$ (yellow markers) are plotted versus the AC-Stark-shifted qubit frequency $\omega_q^{st}$. $\Gamma_1$ is extracted from the exponential relaxation of the qubit population toward its thermal equilibrium. $\Gamma_2^*$ is obtained from Ramsey-fringe decay. The dashed blue line shows $\Gamma_1(\omega_q^{st})$ from Eq. (E8), with $\Gamma_1^{ge}$ as the sole fit parameter. The dashed yellow line represents $\Gamma_2^* = \Gamma_1/2 + \Gamma_\phi$, using $\Gamma_\phi(\omega_q^{st})$ from Eq. (E9) and fitting $\delta\Omega_d/\Omega_d$. Vertical black dashed lines indicate the bare qubit frequency $\omega_q$ and the mechanical mode frequency $\omega_m$. Inset: Qubit relaxation at $\omega_m$ following preparation in the state $|e\rangle$. Experimental data (blue) and corresponding exponential fit (orange) are shown.</figcaption></figure>

of the short-lived $\{|f\rangle, |h\rangle\}$ states into the lower dressed states—an effect analogous to dressed-state spontaneous emission in optical tweezers (see Eq. E6). If the bare-manifold decay rates are $\Gamma_1^{ge}$ and $\Gamma_1^{fh}$, then the dressed-state relaxation rate is

$$ \Gamma_1(\delta\omega_q) = \cos^2(\Theta/2) \Gamma_1^{ge} + \sin^2(\Theta/2) \Gamma_1^{fh}. $$

For small mixing angle $\Theta \approx \delta\omega_q/\Delta_{st}$, expanding to first order gives

$$ \Gamma_1(\delta\omega_q) \approx \Gamma_1^{ge} + \frac{\Gamma_1^{fh}}{2\Delta_{st}} \delta\omega_q. \qquad (E8) $$

The dashed blue line in Fig. 11 plots Eq. (E9) using the independently measured value $\Gamma_1^{fh}/2\pi = 139$ kHz, leaving only $\Gamma_1^{ge}$ as a fit parameter.

We attribute the enhanced dephasing rate $\Gamma_2^*(\omega_q^{st})$ to amplitude fluctuations of the Stark drive $\Omega_d$, yielding

$$ \Gamma_{\phi} = \frac{d(\delta\omega_q)}{d\Omega_d} \delta\Omega_d \approx 2\delta\omega_q \frac{\delta\Omega_d}{\Omega_d}, \qquad (E9) $$

where the approximation holds for $\Omega_d \ll \Delta_{st}$. Assuming classical amplitude noise (i.e. constant $\delta\Omega_d/\Omega_d$), fitting this expression to the data (yellow dashed line in Fig. 11) gives $\delta\Omega_d/\Omega_d \approx 0.6\%$, in good qualitative agreement with the ~0.2% specification of our Anapico source.

<figure><img src="image_13.png" /><figcaption>FIG. 12. **Qubit preparation fidelity.** Rabi fringes obtained by driving the $|g\rangle \to |h\rangle$ transition after three different initialization sequences: $|g\rangle$ (blue), $|e\rangle$ (red), and thermal state (gray). Solid lines are simultaneous fits to Eq. E13, yielding the preparation fidelities $\eta_g$ and $\eta_e$.</figcaption></figure>

## 4. Qubit reset

The typical base temperature of our dilution refrigerator is 10 mK, corresponding to a thermal energy of $k_B T/h \approx 0.2$ GHz. Because the qubit transition frequency, $\omega_q/2\pi \approx 2.2$ MHz, lies well below this thermal scale, the ground $|g\rangle$ and excited $|e\rangle$ states are nearly equally populated at equilibrium. To initialize the qubit in a pure state, we exploit a two-photon transition to the second chain mode at 6.42 GHz.

First, to prepare the $|e\rangle$ state, we drive the $|g, 0\rangle \to |e, 1\rangle$ transition with a continuous-wave tone at $(\omega_{chain} + \omega_q)/2 = 3.2118$ GHz. Owing to the short lifetime of the chain-mode, this quickly transfers all qubit population to $|e\rangle$ via optical pumping. We also investigated preparing $|g\rangle$ via the reverse two-photon process $|e, 0\rangle \to |g, 1\rangle$, but found its fidelity to be significantly lower. Instead, we reset the qubit to $|g\rangle$ by first preparing $|e\rangle$ as above and then applying a resonant $\pi$-pulse on the $|e, 0\rangle \to |g, 0\rangle$ transition. To maintain a constant initialization window in our qubit-based spectrum analyzer measurements, we pad the faster $|e\rangle$ preparation with an additional delay equal to the $\pi$-pulse duration, ensuring that both $|g\rangle$ and $|e\rangle$ initializations occupy the same total time.

## 5. Qubit Preparation Fidelity

In this section, we quantify the fidelity of qubit initialization in $|g\rangle$ and $|e\rangle$. We introduce an auxiliary level $|h\rangle$ that is unpopulated in thermal equilibrium at 10 mK.

**Step 1: state preparation.** The system is initialized in one of three target states: ground $|g\rangle$, excited $|e\rangle$, or thermal equilibrium. The true post-initialization reduced density operator in the qubit subspace writes

$$ \rho_0^i = p_i |g\rangle\langle g| + (1-p_i) |e\rangle\langle e|, \qquad (E10) $$

where $i \in \{g, e, \text{thermal}\}$ labels the preparation and $p_i$ is
the population in $|g\rangle$. In particular,

$$p_e = \frac{1 - \eta_e}{2}, \quad p_g = \frac{1 + \eta_g}{2}, \quad p_{\text{thermal}} = \frac{1}{2}, \qquad (\text{E11})$$

with $\eta_e$ and $\eta_g$ the preparation fidelities of $|e\rangle$ and $|g\rangle$, respectively.

**Step 2: transfer to auxiliary state $|h\rangle$.** In order to measure the fidelities $\eta_e$ and $\eta_g$, we apply a controllable pulse of 152-ns, resonant with the $|g\rangle \rightarrow |h\rangle$ transition, varying its amplitude to induce a Rabi rotation by angle $\theta$ between $|g\rangle$ and $|h\rangle$. Restricting attention to the diagonal elements, this pulse leads to

$$|g\rangle\langle g| \rightarrow \frac{1 + \cos \theta}{2} |g\rangle\langle g| + \frac{1 - \cos \theta}{2} |h\rangle\langle h|, \qquad (\text{E12})$$

while $|e\rangle$ remains uncoupled. The system state after the pulse writes

$$\rho_0^i \rightarrow (1 - p_i) |e\rangle\langle e| + p_i \frac{1 + \cos \theta}{2} |g\rangle\langle g| + p_i \frac{1 - \cos \theta}{2} |h\rangle\langle h|. \qquad (\text{E13})$$

**Step 3: dispersive readout.** Finally, for each initial preparation we perform a 2-µs dispersive readout optimized to distinguish the $\{|g\rangle, |e\rangle\}$ manifold from $|h\rangle$. We obtain a set of 3 Rabi-oscillations (Fig. 12) from which we extract contrasts $C^i$ after state preparation $i \in \{\text{th}, g, e\}$. Decay of $|h\rangle$ during readout is neglected, as it uniformly attenuates the contrast. According to (E13), we have

$$C^i = v p_i,$$

where $v$ is a global visibility factor. Since $p_{\text{th}} = \frac{1}{2}$, $p_g = (1 + \eta_g)/2$, and $p_e = (1 - \eta_e)/2$, we have

$$\eta_g = \frac{C^g}{C^{\text{th}}} - 1 = 98.5(2.9)\%, \quad \eta_e = 1 - \frac{C^e}{C^{\text{th}}} = 98.3(1.0)\%.$$

Uncertainties are estimated via bootstrapping, where fits are repeated using subsets of the amplitude points, and the reported errors correspond to the standard deviation of the extracted parameters across the ensemble.

## 6. Qubit readout fidelity

We have developed an optimized readout protocol that enables direct, single-shot discrimination between the $|g\rangle$ and $|e\rangle$ states, without requiring population transfer to higher excited states as in previous implementations [4]. To demonstrate this capability and quantify the readout fidelity, we prepare the qubit in either state and record the resulting I-Q histograms (see Fig. 13[b]). The readout pulse consists of a 7.5 µs waveform with a slowly decaying envelope (see Fig. 13[a]), carefully engineered to maximize contrast in the $|g\rangle-|e\rangle$ manifold. During the readout, the quadratures ($I_m, Q_m$) are sampled every 4 ns. Final values of I and Q are obtained through

<figure><img src="image_14.png" /><figcaption>FIG. 13. Qubit readout fidelity. (a) Amplitude envelope of the 7.5 µs readout pulse. (b) Time-domain integration weights $w_k^I$ and $w_k^Q$ used to compute the I-Q quadratures: these weights are chosen as the optimal linear discriminant (matched filter) between the average $|g\rangle$ and $|e\rangle$ responses. (c) Histogram of single-shot outcomes in the I-Q plane, collected after preparing $|g\rangle$ (blue) or $|e\rangle$ (red). (d) Histograms of the I quadrature for the same data. Colored dashed lines show the fitted Gaussian components (equal variance) that represent the ideal distributions in the absence of preparation and readout errors. The yellow vertical dashed line indicates the discrimination threshold chosen to equalize the contrast of the spectra.</figcaption></figure>

weighted integration using time-dependent weights $w_m^I$ and $w_m^Q$, optimized to enhance the separation between the $|g\rangle$ and $|e\rangle$ outcomes (see Fig. 13[a]).

A fit of the I histogram to a double Gaussian with equal variances is overlaid to highlight deviations arising from finite preparation fidelity and decoherence during readout. The state-discrimination threshold $T_I$, indicated by the dashed gray line, has been chosen to maximize the contrast of Rabi oscillations.

Readout fidelities are defined as $F_g = P(I < T_I | |g\rangle\langle g|)$ and $F_e = P(I > T_I | |e\rangle\langle e|)$. Experimentally, the prepared state are affected by state preparation imperfections, according to (E11). Hence, the measured readout fidelities are

$$\begin{aligned} P(I < T_I | \rho_0^g) &= p_g F_g + (1 - p_g)(1 - F_e) = 83.9(0.9)\%, \\ P(I > T_I | \rho_0^e) &= p_e F_e + (1 - p_e)(1 - F_g) = 68.3(1.0)\%, \end{aligned} \qquad (\text{E14})$$

where $\rho_0^i$ is defined in Eq. E10. From these relations, we
extract intrinsic readout fidelities of

$$F_g = 84.8(1.8)\%, \quad F_e = 69.1(1.2)\%. \qquad (E15)$$

## 7. Flux stabilization

<figure><img src="image_15.png" /><figcaption>FIG. 14. Flux stabilization. (a) I-quadrature of the cavity response (real-part of the reflection coefficient) as a function of the external flux near $\varphi_{\text{ext}} = \pi$. Yellow lines indicate the zoom range shown in panel (b). (b) Dispersive error signal $\mathcal{E}$ generated by dithering the flux around the current value plotted as a function of flux offset $\delta\varphi_{\text{ext}} = \varphi_{\text{ext}} - \pi$.</figcaption></figure>

Fluxonium qubits exhibit enhanced coherence at $\varphi_{\text{ext}} = \pi$, where first-order sensitivity to external flux noise vanishes. In our setup, $\varphi_{\text{ext}}$ is adjusted via both an external coil and an on-chip flux-bias line, but environmental drifts induce slow variations in the required bias current. To maintain $\varphi_{\text{ext}} = \pi$ throughout each experiment presented in the main text, we implement a closed-loop feedback system that actively stabilizes the flux. We exploit the flux-dependent readout cavity response, which exhibits a pronounced dip at $\varphi_{\text{ext}} = \pi$ (see Fig. 14(a)). At each repetition of the experimental sequence, we apply symmetric perturbations of the flux bias $\pm 1.9 \cdot 10^{-4} \times 2\pi$ around the current operating flux point and record the corresponding cavity responses. The difference between these two reflected signals defines the dispersive error signal $\mathcal{E}$, which—near $\varphi_{\text{ext}} = \pi$—is proportional to the flux deviation $\delta\varphi_{\text{ext}} = \varphi_{\text{ext}} - \pi$. Figure 14(b) plots $\mathcal{E}$ versus applied flux, showing an approximately linear response over the range $|\varphi_{\text{ext}} - \pi| < 2 \times 10^{-4} \times 2\pi$. At each correction step, the current through the on-chip flux line is offset by an amount proportional to the measured error signal $\mathcal{E}$, thus realizing an integral controller.

# Appendix F: Spectrum analyzer theory

## 1. Quantum harmonic oscillator spectrum

The mechanical system is modeled as a harmonic oscillator, with annihilation operator $\hat{a}$, coupled to a thermal bath at rate $\kappa_m$ yielding a mean thermal phonon number $n_{th}$. The Hamiltonian and loss operators capturing the system dynamics can be written as

$$\begin{gathered} \mathcal{H}_m/\hbar = \left(\omega_m + \frac{1}{2}\right) \hat{a}^{\dagger}\hat{a}, \\ \mathcal{L}_{m,\downarrow} = \sqrt{\kappa_m(n_{th} + 1)} \hat{a}, \quad \mathcal{L}_{m,\uparrow} = \sqrt{\kappa_m n_{th}} \hat{a}^{\dagger}. \end{gathered}$$

The evolution of the system is modeled by the Lindblad equation $\dot{\hat{\rho}} = \mathcal{L}_m \hat{\rho}$, with

$$\mathcal{L}_m \dot{\hat{\rho}} = -i[\mathcal{H}, \hat{\rho}] + \mathcal{D}[\mathcal{L}_{m,\downarrow}] + \mathcal{D}[\mathcal{L}_{m,\uparrow}]. \qquad (F1)$$

the Lindbladian.

Closely following the reasoning of [25], we'll discuss in this part the quantum noise spectrum. Similarly to classical physics, the noise spectrum of the Heisenberg picture position operator $\hat{x}(t) = \hat{x}e^{\mathcal{L}_m t}$ can be written as

$$S_{xx}(\omega) = \int_{-\infty}^{+\infty} dt \, e^{i\omega t} \langle \mathbf{x}(t)\mathbf{x}(0) \rangle. \qquad (F2)$$

In general, for $t \neq 0$, the commutator $[\hat{x}(t), \hat{x}(0)]$ is a complex non-zero quantity. This reflects the fact that the measurement of $\hat{x}(t)$ is necessarily impacted by the back-action of the measurement of $\hat{x}(0)$. Consequently, the complex-valued auto-correlator $C_{xx}(t) = \langle \hat{x}(t)\hat{x}(0) \rangle$ is a purely mathematical construction that cannot be evaluated by a simple continuous measurement of $\hat{x}(t)$. Indeed, in a practical measurement setting, the added quantum noises ensure that the measured quantities $\hat{x}(t)$ commute with each other at all times $t$, resulting in real-valued correlators. Nevertheless, the complex-nature of $C_{xx}(t)$, which results in a frequency asymmetric spectrum $S_{xx}(\omega)$, has a profound physical meaning. We can gain insights by developing $\hat{x}(t) = (\hat{a}^{\dagger} + \hat{a})/2$, in its elementary components and dropping vanishing terms

$$\langle \mathbf{x}(t)\mathbf{x}(0) \rangle = \frac{1}{4} \left( \langle \hat{a}(t)\hat{a}^{\dagger}(0) \rangle + \langle \hat{a}^{\dagger}(t)\hat{a}(0) \rangle \right). \qquad (F3)$$

The first correlator $\langle \hat{a}(t)\hat{a}^{\dagger}(0) \rangle$ can be understood as the amplitude of probability in a fictitious experiment where an excitation would be recovered at time $t$ after being added to the oscillator at time 0. Therefore, the noise spectrum associated to this correlator is the absorption spectrum of the mechanical mode. Similarly, the noise spectrum associated to $\langle \hat{a}^{\dagger}(t)\hat{a}(0) \rangle$ is the emission spectrum of the mechanical mode, where an excitation is removed from the oscillator at time 0 and re-injected at time $t$. Since the mechanical mode is coupled to a bath, both correlators tend to 0 exponentially as the excitation
can be lost in the interval. Computing this correlators explicitly, we obtain

$$ \langle \mathbf{x}(t)\mathbf{x}(0) \rangle = \frac{1}{4} \left( \langle \mathbf{a}^\dagger(0)\mathbf{a}(0) \rangle e^{-i\omega_m t} + \langle \mathbf{a}(0)\mathbf{a}^\dagger(0) \rangle e^{i\omega_m t} \right) e^{-\frac{\kappa_m}{2}|t|}. \quad (F4) $$

Which in turn yields the spectra

$$ S_{xx}(\omega) = S_{a^\dagger a}(\omega) + S_{aa^\dagger}(\omega) \quad (F5) $$

with

$$ \begin{aligned} S_{a^\dagger a}(\omega) &= \frac{n_{\text{th}}\kappa_m}{(\omega - \omega_m)^2 + (\kappa_m/2)^2}, \\ S_{aa^\dagger}(\omega) &= \frac{(n_{\text{th}} + 1)\kappa_m}{(\omega + \omega_m)^2 + (\kappa_m/2)^2}, \end{aligned} \quad (F6) $$

where the negative and positive frequency parts corresponds to the emission and absorption spectrum, respectively. Crucially, the total spectrum $S_{xx}(\omega)$ is not symmetric-in-frequency and cannot be obtained from the correlator of a real-valued time-dependent observable $\langle \mathbf{A}(t)\mathbf{A}(0) \rangle$. Instead, we separately measure the real part of the autocorrelation functions Re($C_{a^\dagger a}$) and Re($C_{aa^\dagger}$), the Fourier transform of which yields the symmetrized spectra $\overline{S}_{a^\dagger a}$ and $\overline{S}_{aa^\dagger}$, by inspecting the correlation between successive measurements where the qubit is repeatedly prepared in $|g\rangle$ or $|e\rangle$ respectively.

## 2. System definition

In this section we explain how to experimentally measure the emission and absorption spectrum from periodic weak measurements of the membrane with a qubit. The system is composed of a Fluxonium, modeled as a qubit with frequency $\omega_q$ and the membrane mode. The two systems interact with a Jaynes-Cummings Hamiltonian $\mathbf{H}_c/\hbar = i\frac{\Omega}{2}(\mathbf{a} + \mathbf{a}^\dagger)(\mathbf{\sigma} - \mathbf{\sigma}^\dagger)$, yielding the system Hamiltonian

$$ \mathbf{H} = \frac{\hbar\omega_q}{2}\boldsymbol{\sigma}_z + \mathbf{H}_m + \mathbf{H}_c. $$

In the spectrum analyzer experiment, the qubit frequency is regularly Stark-shifted, close to the membrane frequency, to $\omega_m - \Delta$. In the frame rotating at the Stark-shifted qubit frequency for both the membrane and the qubit, after a rotating-wave approximation, the Hamiltonian simplifies to $\mathbf{H}/\hbar = \Delta\mathbf{a}^\dagger\mathbf{a} + i\frac{\Omega}{2}(\mathbf{a}^\dagger\mathbf{\sigma} - \mathbf{a}\mathbf{\sigma}^\dagger)$. When the qubit frequency is not shifted, we can neglect the effect of the coupling. The Lindbladian superoperator of the membrane free evolution is $\mathcal{L}_m\rho = -i[\Delta\mathbf{a}^\dagger\mathbf{a}, \rho] + \mathcal{D}[\mathbf{L}_{m,\downarrow}]\rho + \mathcal{D}[\mathbf{L}_{m,\uparrow}]\rho$. For the sake of simplicity, we will consider for now an ideal qubit and we will treat the qubit imperfections in App. F 5.

## 3. Semi-classical treatment

### a. Autocorrelation function

In the semi-classical hypothesis, we treat classically the mechanical mode degrees of freedom, $\mathbf{x}$ and $\mathbf{p}$. In the frame of the qubit, the mode complex amplitude $\alpha(t)$ is slowly oscillating at rate $\Delta$ with random phase and amplitude fluctuations over a characteristic time-scale $1/\kappa_m$. As the interaction time is short compared to the mechanical mode time-scales, $\alpha(t)$ can be considered constant over one interaction. Moreover, the interaction Hamiltonian can be simplified to $\mathbf{H}_c/\hbar = i\frac{\Omega}{2}(\alpha^*\mathbf{\sigma} - \alpha\mathbf{\sigma}^\dagger)$, which is a classical drive with amplitude $\alpha$ on the qubit. Expanding the terms, we obtain that $\mathbf{H}_c/\hbar = \frac{\Omega}{2}(\text{Im}(\alpha)\mathbf{\sigma}_x + \text{Re}(\alpha)\mathbf{\sigma}_y)$. The unitary evolution after an interaction time $\tau$ can be written as

$$ U_I = \exp\left(-\frac{i}{\hbar}\mathbf{H}_c\tau\right) = \exp\left(-i\frac{\Omega}{2}|\alpha|\vec{n}\cdot\vec{\sigma}\tau\right) \quad (F7) $$

where $\vec{n} = \frac{1}{|\alpha|}(\text{Im}(\alpha), \text{Re}(\alpha), 0)$ and $\vec{\sigma} = (\boldsymbol{\sigma}_x, \boldsymbol{\sigma}_y, \boldsymbol{\sigma}_z)$. Comparing with the standard SU(2) rotation matrix $\exp(-i\frac{\theta}{2}\vec{n}\cdot\vec{\sigma})$ shows that the qubit undergoes a rotation by $\theta = \Omega|\alpha|\tau$ about the axis $\vec{n}$ in the Bloch-sphere equatorial plane, whose complex direction is $i\alpha^*$. Therefore, measuring $\boldsymbol{\sigma}_x$ of the qubit, at time $t_k$, performs a measurement of the position of the mechanical mode as

$$ \langle \boldsymbol{\sigma}_x(t_k) \rangle = \theta_1 x(t_k). \quad (F8) $$

In this semi-classical treatment, where the position $x(t)$ is an external parameter which is not subject to qubit backaction, the measurement outcomes at different times are independent quantum observables. We thus have $\langle \boldsymbol{\sigma}_{x,k}\boldsymbol{\sigma}_{x,k'} \rangle = \langle \boldsymbol{\sigma}_{x,k} \rangle \langle \boldsymbol{\sigma}_{x,k'} \rangle$ for $k \neq k'$. It thus follows from Eq. F8 that

$$ C_k \equiv \overline{\langle \boldsymbol{\sigma}_{x,k} \boldsymbol{\sigma}_{x,0} \rangle} = \theta_1^2 \overline{x(t_k)x(0)}, \quad (F9) $$

where the overline denotes the averaging over the different trajectories of the harmonic oscillator. Moreover, in the case where $k = 0$, the autocorrelation reads $C_0 = 1$ as $\sigma_x^2 = I$. This proves Eq. 4 of the main text.

### b. Wiener-Khintchine theorem

Let us define the discrete Fourier transform of a signal $\{m_k\}_{0 \le k < N}$ sampled at a frequency $1/T$ as $\mathcal{F}[\{m_k\}](\omega_r) = \sum_{n=0}^{N-1} m_n e^{-i\omega_r n}$, where $\omega_r = 2\pi r/NT$ for $0 \le r < N$. These discrete time frequencies maps unambiguously to continuous time frequencies $\omega$ on the first Nyquist zone $0 \le \omega < \pi/T$. We thus drop the subscript $r$ in the following. We define the power spectral density as

$$ S_{mm}(\omega) = T \lim_{N \to \infty} \frac{1}{N} \overline{|\mathcal{F}[\{m_k\}](\omega)|^2}. \quad (F10) $$
We now reintroduce the free evolution of the membrane degrees of freedom. Since this evolution is much slower than all other dynamics in the system, we assume that it commutes with the qubit interaction evolution. To include this contribution, we precede the Kraus map operator with the membrane's free evolution operator, acting over a duration $T = \tau + T_{\text{idle}}$, where $T_{\text{idle}}$ corresponds to the idle time between consecutive weak measurements. The effective total evolution operator when discarding measurement results reads

$$e^{\mathcal{L}_m T} \bar{\mathcal{K}}_i \approx \begin{cases} e^{\mathcal{L}_m T} e^{\frac{\theta_1^2}{4} \mathcal{D}[\mathbf{a}]} & \text{if } i = g, \approx e^{\mathcal{L}'_{m,i} T}, \quad (\text{F19}) \\ e^{\mathcal{L}_m T} e^{\frac{\theta_1^2}{4} \mathcal{D}[\mathbf{a}^\dagger]} & \text{if } i = e, \end{cases}$$

where the effective Lindbladian of the membrane including environmental dissipation and dynamical backaction writes

$$\mathcal{L}'_{m,i} = \mathcal{L}_m + \begin{cases} \kappa \mathcal{D}[\mathbf{a}] & \text{if } i = g, \\ \kappa \mathcal{D}[\mathbf{a}^\dagger] & \text{if } i = e, \end{cases} \quad (\text{F20})$$

where $\kappa = \theta_1^2/4T$ is the dynamical backaction rate. Explicitly, the Lindbladian reads

$$\begin{aligned} \mathcal{L}'_{m,g} &= -i[\mathbf{H}, \cdot] + (\kappa_m(n_{\text{th}}+1)+\kappa) \mathcal{D}[\mathbf{a}] + \kappa_m n_{\text{th}} \mathcal{D}[\mathbf{a}^\dagger]\rho \\ \mathcal{L}'_{m,e} &= -i[\mathbf{H}, \cdot] + \kappa_m(n_{\text{th}}+1)\mathcal{D}[\mathbf{a}] + (\kappa_m n_{\text{th}}+\kappa)\mathcal{D}[\mathbf{a}^\dagger]\rho. \end{aligned} \quad (\text{F21})$$

Inspecting the coefficients in front of $\mathcal{D}[\mathbf{a}]$ and $\mathcal{D}[\mathbf{a}^\dagger]$, one can observe that the membrane's evolution after each measurement cycle is identical to that induced by a coupling at rate $\kappa'_{m,i}$ to an effective environment with thermal occupation $n'_{\text{th},i}$, where

$$\begin{aligned} \kappa'_{m,i} &= \begin{cases} \kappa_m + \kappa & \text{if } i = g, \\ \kappa_m - \kappa & \text{if } i = e; \end{cases} \\ n'_{\text{th}} &= \begin{cases} \frac{\kappa_m n_{\text{th}}}{\kappa'_m} & \text{if } i = g, \\ \frac{\kappa_m n_{\text{th}}}{\kappa'_m} + \frac{\kappa}{\kappa'_m} & \text{if } i = e. \end{cases} \end{aligned} \quad (\text{F22})$$

In other words, the mechanical occupancy is modified by a factor $1/(1 \pm C)$ due to dynamical backaction, where we have introduced the cooperativity $C = \kappa/\kappa_m$ in line with conventions commonly used in the optomechanics community. The last additive term in the occupation $\kappa/\kappa'_m = C/(1+C)$ appearing when the qubit is prepared in state $|e\rangle$ corresponds to the additional fluctuations resulting from quantum backaction. In our experiment where the cooperativity is on the order of 10%, we can safely neglect this contribution in the area of the spectrum $S_{mm}^e(\omega)$.

In the case where the qubit is alternatively prepared in $|g\rangle$ and $|e\rangle$, the mean evolution operator over two cycles reads $e^{\mathcal{L}_m T} \bar{\mathcal{K}}_e e^{\mathcal{L}_m T} \bar{\mathcal{K}}_g \approx e^{2\mathcal{L}''_m T}$, where $\mathcal{L}''_m = \mathcal{L}_m + \frac{\kappa}{2} \mathcal{D}[\mathbf{a}] + \frac{\kappa}{2} \mathcal{D}[\mathbf{a}^\dagger]$. Following the same steps as before, we obtain the effective total mechanical loss rate and steady-state phonon number

$$\kappa_m'' = \kappa_m \quad \text{and} \quad n_{\text{th}}'' = n_{\text{th}} + \frac{\kappa/2}{\kappa_m}. \quad (\text{F23})$$

Once again, the additive term $\frac{\kappa/2}{\kappa_m} = C/2$, corresponding to the quantum backaction is negligible compared to $n_{\text{th}}$. We thus conclude that the alternating preparation protocol effectively cancels the effect of dynamical backaction.

### c. Emission and absorption spectra

The Kraus operators can be written as the product of the mean evolution and the measurement backaction Kraus operator $\mathcal{K}'_{\pm,i}$,

$$\mathcal{K}_{\pm,i} = \frac{\bar{\mathcal{K}}_i}{2} \begin{cases} \mathbb{I} \mp \frac{\theta_1}{2} \mathcal{M}[\mathbf{a}] & \text{if } i = g, \\ \mathbb{I} \pm \frac{\theta_1}{2} \mathcal{M}[\mathbf{a}^\dagger] & \text{if } i = e. \end{cases} = \bar{\mathcal{K}}_i \mathcal{K}'_{\pm,i} \quad (\text{F24})$$

Preparing all the atoms in $|i\rangle \in \{|g\rangle, |e\rangle\}$, the probability to obtain the measurement record $(m_k, \dots, m_0)$, with $m_j \in \{+, -\}$, is

$$\begin{aligned} \mathbb{P}[m_k, \dots, m_0 | \rho_m] &= \\ &= \text{Tr}[e^{\mathcal{L}_m T} \mathcal{K}_{m_k, i} e^{\mathcal{L}_m T} \mathcal{K}_{m_{k-1}, i} \dots e^{\mathcal{L}_m T} \mathcal{K}_{m_0, i} \rho_m]. \end{aligned} \quad (\text{F25})$$

Computing the probability to obtain the 2 results record $(m_k, m_0)$ amounts to summing over all intermediate trajectories

$$\begin{aligned} \mathbb{P}[m_k, m_0 | \rho_m, i] &= \\ &= \sum_{m_{k-1}=\pm} \dots \sum_{m_1=\pm} \mathbb{P}[m_k, \dots, m_0 | \rho_m] \\ &= \text{Tr}[\mathcal{K}'_{m_k, i} (e^{\mathcal{L}_{m,i} T} \bar{\mathcal{K}}_i)^{k-1} e^{\mathcal{L}_m T} \bar{\mathcal{K}}_i \mathcal{K}'_{m_0, i} \rho_m], \end{aligned} \quad (\text{F26})$$

where in the last line we used the linearity of the trace, summing the Kraus operator to the mean ones. The previous equation can be written in terms of the effective Lindbladian operator $\mathcal{L}'_{m,i}$.

$$\mathbb{P}[m_k, m_0 | \rho_m, i] \approx \text{Tr}[\mathcal{K}'_{m_k, i} e^{\mathcal{L}'_{m,i} k T} \mathcal{K}'_{m_0, i} \rho_m]. \quad (\text{F27})$$

In the trivial case where $k=0$, the autocorrelation of the measurement result is $C_0^i = 1$ as $\sigma_x^2 = \mathbb{I}$. In the case where $k \neq 0$, the autocorrelation of the measurement result can be expressed as

$$C_k^i = \mathbb{P}[+k, +0 | \rho_m, i] + \mathbb{P}[-k, -0 | \rho_m, i] \\ - \mathbb{P}[+k, -0 | \rho_m, i] - \mathbb{P}[-k, +0 | \rho_m, i]. \quad (\text{F28})$$

Using the previous result (Eq. F27) and expanding the Kraus operators leads to

$$C_k^g = \frac{\theta_1^2}{4} \text{Tr} \left[ (\mathcal{M}[\mathbf{a}] e^{\mathcal{L}'_m k T} \mathcal{M}[\mathbf{a}]) \rho_m \right], \quad (\text{F29})$$

$$C_k^e = \frac{\theta_1^2}{4} \text{Tr} \left[ (\mathcal{M}[\mathbf{a}^\dagger] e^{\mathcal{L}'_m k T} \mathcal{M}[\mathbf{a}^\dagger]) \rho_m \right].$$

Defining the Heisenberg picture operators $O(t)\rho = Oe^{\mathcal{L}'_m t}\rho = O\rho(t)$, the measurement correlator reads
$$C_k^i = \frac{\theta_1^2}{4} \begin{cases} \text{Tr} \left[ \mathbf{a}(t_k) (\mathbf{a}(0)\boldsymbol{\rho}_m + \boldsymbol{\rho}_m \mathbf{a}^\dagger(0)) + (\mathbf{a}(0)\boldsymbol{\rho}_m + \boldsymbol{\rho}_m \mathbf{a}^\dagger(0)) \mathbf{a}^\dagger(t_k) \right] & \text{if } i = g, \\ \text{Tr} \left[ \mathbf{a}^\dagger(t_k) (\mathbf{a}^\dagger(0)\boldsymbol{\rho}_m + \boldsymbol{\rho}_m \mathbf{a}(0)) + (\mathbf{a}^\dagger(0)\boldsymbol{\rho}_m + \boldsymbol{\rho}_m \mathbf{a}(0)) \mathbf{a}^\dagger(t_k) \right] & \text{if } i = e. \end{cases} \quad (F30)$$

Using the cyclic permutation properties of the trace, the autocorrelation writes

$$C_k^g = \frac{\theta_1^2}{4} \text{Tr} [\mathbf{a}^\dagger(t_k) \mathbf{a}(0) \boldsymbol{\rho}_m + \boldsymbol{\rho}_m \mathbf{a}^\dagger(0) \mathbf{a}(t_k)] \quad (F31)$$

$$C_k^e = \frac{\theta_1^2}{4} \text{Tr} [\mathbf{a}(t_k) \mathbf{a}^\dagger(0) \boldsymbol{\rho}_m + \boldsymbol{\rho}_m \mathbf{a}(0) \mathbf{a}^\dagger(t_k)],$$

where we removed the terms Tr[$\mathbf{a}(t)\mathbf{a}(0)$] = Tr[$\mathbf{a}^\dagger(t)\mathbf{a}^\dagger(0)$] = 0 in the thermal steady-state. We finally find that $C_0^g = C_0^e = 1$ and

$$C_{k|k \neq 0}^i = \frac{\theta_1^2}{4} \begin{cases} \text{Re}(C_{a^\dagger a}(t_k)) & \text{if } i = g, \\ \text{Re}(C_{aa^\dagger}(t_k)) & \text{if } i = e. \end{cases} \quad (F32)$$

Therefore, the spectra measured when the qubit is prepared in $|g\rangle$ (resp. $|e\rangle$) is the emission (resp. absorption) spectrum of the mechanical mode.

$$S_{mm}^g(\omega_k) = T + \frac{\theta_1^2}{4} \bar{S}_{a^\dagger a}(\omega_k + \omega_q), \quad (F33)$$

$$S_{mm}^e(\omega_k) = T + \frac{\theta_1^2}{4} \bar{S}_{aa^\dagger}(\omega_k + \omega_q),$$

where $\omega_k = \frac{2\pi}{T}(k - \frac{N}{2})$ for $k \in [0, N]$ and $\bar{S}_{AB} = \int_{-\infty}^{\infty} dt e^{i\omega t} \text{Re}(\langle A(t)B(0) \rangle)$. Using F6, we find the expression of the measured spectrum computed with N consecutive measurement

$$S_{mm}^g(\omega_k') = T + \frac{\theta_1^2}{4} \frac{n_{th} \kappa'_{m,g}}{(\omega_k' - \omega_m)^2 + (\kappa'_{m,g}/2)^2}, \quad (F34)$$

$$S_{mm}^e(\omega_k') = T + \frac{\theta_1^2}{4} \frac{(n_{th} + 1)\kappa'_{m,e}}{(\omega_k' - \omega_m)^2 + (\kappa'_{m,e}/2)^2},$$

where $\omega_k' = \frac{2\pi}{T}(k - \frac{N}{2}) + \omega_q$ as we measured in the frame of the qubit. As we have set the qubit frequency at a detuning $\Delta$ from the mechanical mode frequency, the frequency part of the denominators simplifies to $\omega_k - \Delta$ where $\omega_k = \frac{2\pi}{T}(k - \frac{N}{2})$. Moreover, the $t = 0$ autocorrelation corresponds to the offset in the final spectrum corresponding to the measurement shot noise. This leads to a signal-to-noise ratio of $4n_{th}\kappa/\kappa_m$.

Performing a similar derivation in the case where the qubit is alternatively prepared in $|g\rangle$ and $|e\rangle$, we obtain

$$S_{mm}^g(\omega_k') = T + \frac{\theta_1^2}{4} \frac{n_{th} \kappa_m}{(\omega_k' - \omega_m)^2 + (\kappa_m/2)^2}, \quad (F35)$$

$$S_{mm}^e(\omega_k') = T + \frac{\theta_1^2}{4} \frac{(n_{th} + 1)\kappa_m}{(\omega_k' - \omega_m)^2 + (\kappa_m/2)^2}.$$

The two spectra only differ by the factor $n_{th}$ and $n_{th} + 1$, corresponding to the fundamental asymmetry between absorption and emission processes.

## 5. Effect of the qubit imperfections

In this section, we account for imperfect state preparation and readout, as well as intrinsic qubit decoherence. We assume the qubit is initialized in $|g\rangle$ and $|e\rangle$ with preparation fidelities $\eta_g$ and $\eta_e$, respectively. The qubit's energy relaxation/excitation (at rate $\kappa_1/2$) and pure dephasing (at rate $\kappa_\phi$) are described by the Lindblad jump operators

$$\mathcal{L}_{q,\downarrow} = \sqrt{\frac{\kappa_1}{2}} \boldsymbol{\sigma}, \quad \mathcal{L}_{q,\uparrow} = \sqrt{\frac{\kappa_1}{2}} \boldsymbol{\sigma}^\dagger, \quad \mathcal{L}_{q,\phi} = \sqrt{\frac{\kappa_\phi}{2}} \boldsymbol{\sigma}_z. \quad (F36)$$

The relevant hierarchy in our experimental regime is $\Delta, \kappa_m n_{th} \ll \Omega \sqrt{n_{th}} \ll \kappa_1, \kappa_\phi$. Therefore, we can perform a perturbation treatment on the losses of the qubit. The Lindblad super-operator, $\mathcal{L}$ and the joint density matrix, $\boldsymbol{\rho}$, can be vectorized using the column-stacking convention on the qubit space. The total Hilbert space is $\mathcal{H}_q \otimes \mathcal{H}_m$, then the vectorized density matrix reads $\vec{\boldsymbol{\rho}} = (\boldsymbol{\rho}_{ee}, \boldsymbol{\rho}_{eg}, \boldsymbol{\rho}_{ge}, \boldsymbol{\rho}_{gg})^T$, where each elements is an operator on $\mathcal{H}_m$. An operator $\mathcal{O}$ acting on the left (resp. right) of the oscillator density matrix is written $\mathcal{O}_-$ (resp. $\mathcal{O}_+$). Using these conventions, the super-operators can be written as, $\mathcal{L} = \mathcal{L}_0 + \mathcal{L}_1 + \mathcal{L}_m$, with

$$\mathcal{L}_0 = \begin{pmatrix} -\kappa_1/2 & 0 & 0 & \kappa_1/2 \\ 0 & -\kappa_2 & 0 & 0 \\ 0 & 0 & -\kappa_2 & 0 \\ \kappa_1/2 & 0 & 0 & -\kappa_1/2 \end{pmatrix},$$

capturing the qubit loss operators and

$$\mathcal{L}_1 = \frac{\Omega}{2} \begin{pmatrix} 0 & \mathbf{a}_{\leftarrow}^{\dagger} & \mathbf{a}_{\rightarrow} & 0 \\ -\mathbf{a}_{\leftarrow} & 0 & 0 & \mathbf{a}_{\rightarrow} \\ -\mathbf{a}_{\rightarrow}^{\dagger} & 0 & 0 & \mathbf{a}_{\leftarrow}^{\dagger} \\ 0 & -\mathbf{a}_{\rightarrow}^{\dagger} & -\mathbf{a}_{\leftarrow} & 0 \end{pmatrix}$$

capturing the qubit membrane coupling.

### a. Perturbative treatment

As $\mathcal{L}_0$ acts only on the qubit and $\mathcal{L}_m$ only on the mechanical mode, $[\mathcal{L}_0, \mathcal{L}_m] = 0$. Under the approximation of a weak coupling, the evolution super-operator can be approximated by the truncated Dyson series
$$ 
\begin{aligned}
e^{\mathcal{L}t} & \approx e^{(\mathcal{L}_0 + \mathcal{L}_m)t} + \int_0^t dt_1 e^{(\mathcal{L}_0 + \mathcal{L}_m)(t-t_1)} \mathcal{L}_1 e^{(\mathcal{L}_0 + \mathcal{L}_m)t_1} + \int_0^t dt_2 \int_0^{t_2} dt_1 e^{(\mathcal{L}_0 + \mathcal{L}_m)(t-t_2)} \mathcal{L}_1 e^{(\mathcal{L}_0 + \mathcal{L}_m)(t_2-t_1)} \mathcal{L}_1 e^{(\mathcal{L}_0 + \mathcal{L}_m)t_1} \\
& \approx e^{\mathcal{L}_m t} \left( e^{\mathcal{L}_0 t} + \int_0^t dt_1 e^{\mathcal{L}_0(t-t_1)} \mathcal{L}_1 e^{\mathcal{L}_0 t_1} + \int_0^t dt_2 \int_0^{t_2} dt_1 e^{\mathcal{L}_0(t-t_2)} \mathcal{L}_1 e^{\mathcal{L}_0(t_2-t_1)} \mathcal{L}_1 e^{\mathcal{L}_0 t_1} \right),
\end{aligned}
\quad (F37) $$

where we used both that $\mathcal{L}_0$ and $\mathcal{L}_m$ commutes and the Baker-Campbell-Hausdorff formula to factorize $e^{\mathcal{L}_m t}$. Upon exponentiation of $\mathcal{L}_0$, we get

$$ e^{\mathcal{L}_0 t} = \begin{pmatrix} \frac{1+e^{-\kappa_1 t}}{2} & 0 & 0 & \frac{1-e^{-\kappa_1 t}}{2} \\ 0 & e^{-\kappa_2 t} & 0 & 0 \\ 0 & 0 & e^{-\kappa_2 t} & 0 \\ \frac{1-e^{-\kappa_1 t}}{2} & 0 & 0 & \frac{1+e^{-\kappa_1 t}}{2} \end{pmatrix}. \quad (F38) $$

The first integral in Eq. F37, can then be computed as

$$ \mathcal{I}_1(t) \equiv \frac{\Omega}{2} \begin{pmatrix} 0 & \mathbf{a}_+^\dagger & \mathbf{a}_- & 0 \\ -\mathbf{a}_+ & 0 & 0 & \mathbf{a}_- \\ -\mathbf{a}_-^\dagger & 0 & 0 & \mathbf{a}_+^\dagger \\ 0 & -\mathbf{a}_-^\dagger & -\mathbf{a}_+ & 0 \end{pmatrix}, \quad (F39) $$

where we defined $\mathbf{a}_\pm = \tau_\Sigma \mathbf{a}_\Sigma \pm \tau_2 \mathbf{a}_\Delta$, with

$$ \begin{aligned} \tau_\Sigma(t) &= \frac{\kappa_1 \tau_1(t) - \kappa_2 \tau_2(t)}{\kappa_1 - \kappa_2}, \\ \tau_1(t) &= \frac{1 - e^{-\kappa_1 t}}{\kappa_1}, \quad \tau_2(t) = \frac{1 - e^{-\kappa_2 t}}{\kappa_2}, \\ \mathbf{a}_\Sigma &= \frac{\mathbf{a}_- + \mathbf{a}_+}{2}, \quad \mathbf{a}_\Delta = \frac{\mathbf{a}_- - \mathbf{a}_+}{2}. \end{aligned} $$

Notably, $\lim_{\kappa_2 \to \kappa_1} \tau_\Sigma(t) = te^{-\kappa_1 t}$ and $\lim_{\kappa_2 \to 0} \tau_1(t) = \tau_2(t) = t$. Due to the matrix structure of this term, it captures the entanglement of the membrane and qubits states, but it doesn't capture the population transfer responsible for heating or cooling of the membrane. In order to capture such effect, we have to consider the second order perturbation term of Eq. F37, that can be written as

$$ \mathcal{I}_2(t) = \frac{\Omega^2}{4} \begin{pmatrix} \mathbf{V}_{1,1} & 0 & 0 & \mathbf{V}_{1,4} \\ 0 & \mathbf{V}_{2,2} & \mathbf{V}_{2,3} & 0 \\ 0 & \mathbf{V}_{3,2} & \mathbf{V}_{3,3} & 0 \\ \mathbf{V}_{4,1} & 0 & 0 & \mathbf{V}_{4,4} \end{pmatrix} \quad (F40) $$

where $\{\mathbf{V}_{i,j}\}_{i,j \in [1,4]}$ terms are not explicitly written for simplicity.

b. *Kraus operators*

The Kraus operators acting on the membrane degree of freedom can be computed as

$$ \begin{aligned} \mathcal{K}_{\pm,i} \rho_m &= \mathrm{Tr}_q[|\pm\rangle \langle\pm| e^{\mathcal{L}\tau} (\rho_m \otimes \rho_{q,i})] \\ \text{with } \rho_{q,i} &= (1 - p_i \quad 0 \quad 0 \quad p_i)^T \end{aligned} \quad (F41) $$

where $p_g = (1 + \eta_g)/2$ and $p_e = (1 - \eta_e)/2$. They can be approximated using Eq. F37 as

$$ \begin{aligned} \mathcal{K}_{\pm,i} = & \frac{e^{\mathcal{L}_m \tau}}{2} \left( \mathbb{I} \pm \frac{\Omega}{4} \left( (2p_i - 1)\tau_\Sigma + \tau_2 \right) \mathcal{M}[\mathbf{a}] \right. \\ & \left. \pm \frac{\Omega}{4} \left( (2p_i - 1)\tau_\Sigma - \tau_2 \right) \mathcal{M}[\mathbf{a}^\dagger] \right. \\ & + \frac{\Omega^2}{4} \left( \frac{\tau - \tau_2}{\kappa_2} + (2p_i - 1) \frac{\tau_2 - \tau_1}{\kappa_1 - \kappa_2} \right) \mathcal{D}[\mathbf{a}] \\ & + \frac{\Omega^2}{4} \left( \frac{\tau - \tau_2}{\kappa_2} - (2p_i - 1) \frac{\tau_2 - \tau_1}{\kappa_1 - \kappa_2} \right) \mathcal{D}[\mathbf{a}^\dagger] \right) \end{aligned} \quad (F42) $$

Due to the limited qubit lifetime, the action of a qubit prepared in $|g\rangle$ (resp. $|e\rangle$) also weakly measures the $\mathbf{a}^\dagger$ (resp. $\mathbf{a}$) operator.

To account for the readout infidelity, we'll write the readout confusion matrix in terms of the single-shot assignment error probabilities $\varepsilon_g = 1 - F_g$ and $\varepsilon_e = 1 - F_e$. As the state is readout in the $|\pm\rangle$ basis with a $\sqrt{X}$ pulse, the readout errors on state $|+\rangle$ (resp. $|-\rangle$) are the one of $|e\rangle$ (resp. $|g\rangle$), defining $\varepsilon_{+/-} = \varepsilon_{e/g}$. The effective Kraus operator acting on the membrane can be written as

$$ \bar{\mathcal{K}}_{\pm,i} = (1 - \varepsilon_{\pm}) \mathcal{K}_{\pm,i} + \varepsilon_{\pm} \mathcal{K}_{\mp,i}. \quad (F43) $$

c. *Emission and absorption spectra*

We can compute all the previous steps to compute the autocorrelation we applied to the lossless qubit, except for the case where $k=0$. Indeed, as the readout of the qubit is asymmetric, the autocorrelation at $k=0$ is modified to $C_k^i = 1 - (\varepsilon_g - \varepsilon_e)^2$. Computing the spectra from this autocorrelations, we obtain that

$$ \begin{aligned} S_{mm}^g(\omega_k') = & (1 - \varepsilon_g - \varepsilon_e)^2 \eta_g \frac{\Omega^2 \tau_\Sigma^2}{4} \left( \frac{\eta_g + \tau_2 / \tau_\Sigma}{2} \bar{S}_{\mathbf{a}^\dagger \mathbf{a}}(\omega_k') + \frac{\eta_g - \tau_2 / \tau_\Sigma}{2} \bar{S}_{\mathbf{a} \mathbf{a}^\dagger}(\omega_k') \right) + T(1 - (\varepsilon_g - \varepsilon_e)^2), \\ S_{mm}^e(\omega_k') = & (1 - \varepsilon_g - \varepsilon_e)^2 \eta_e \frac{\Omega^2 \tau_\Sigma^2}{4} \left( \frac{\eta_e - \tau_2 / \tau_\Sigma}{2} \bar{S}_{\mathbf{a}^\dagger \mathbf{a}}(\omega_k') + \frac{\eta_e + \tau_2 / \tau_\Sigma}{2} \bar{S}_{\mathbf{a} \mathbf{a}^\dagger}(\omega_k') \right) + T(1 - (\varepsilon_g - \varepsilon_e)^2). \end{aligned} \quad (F44) $$

Using the amplitude of the measured thermal peak for the preparation in $|g\rangle$ at 10 mK in Fig. 4 we estimate $\Omega/2\pi = 1.50(8)$ kHz.

The normalized spectra are expressed in units of a calibration peak amplitude resulting from a weak drive on the qubit, the resulting sinus cardinal amplitudes in the spectra are proportionate to $(1 - \varepsilon_g - \varepsilon_e)^2 \eta_e^2$. The area
of the normalized Lorentzian peaks therefore reads

Appendix G: Experimental Setup

$$ 
\begin{aligned}
A_g &= \frac{\Omega^2 \tau_{\Sigma}^2}{4} \left( n_{\text{th}} + \frac{1}{2} - \frac{\tau_2}{2\eta_g \tau_{\Sigma}} \right) \\
A_e &= \frac{\Omega^2 \tau_{\Sigma}^2}{4} \left( n_{\text{th}} + \frac{1}{2} + \frac{\tau_2}{2\eta_e \tau_{\Sigma}} \right)
\end{aligned}
\qquad (\text{F45}) $$

Performing a linear fit of the areas as a function of the cryostat temperature, we are able to express the amplitude in term of phonon number,

$$ 
\begin{aligned}
\bar{n}_g &= n_{\text{th}} + \frac{1}{2} - \frac{\tau_2}{2\eta_g \tau_{\Sigma}} \\
\bar{n}_e &= n_{\text{th}} + \frac{1}{2} + \frac{\tau_2}{2\eta_g \tau_{\Sigma}}
\end{aligned}
\qquad (\text{F46}) $$

The expected asymmetry of the spectrum used in Fig. 4 writes

$$ \bar{n}_e - \bar{n}_g = \frac{\tau_2 \eta_g + \eta_e}{\tau_{\Sigma} 2\eta_g \eta_e}. \qquad (\text{F47}) $$

The experimental setup is shown in Fig. 5. RF signals are controlled and measured using a heterodyne detection scheme, implemented via an FPGA-based controller (OPX, Quantum Machine) combined with a microwave source (Anapico APUASYN20-4). The magnetic flux through the qubit is tuned using a Yokogawa 7651 DC source, which supplies current to a superconducting coil placed directly to the chip enclosure. Fast flux tuning is achieved via an on-chip flux line. Signal amplification is performed using a traveling-wave parametric amplifier (TWPA), a high-electron-mobility transistor (HEMT) amplifier, and a room-temperature amplifier. The DC membrane bias is controlled by a Yokogawa GS200 source and filtered through a commercial 10 kHz low-pass filter (QDevil) to suppress high-frequency noise.

[1] Aspelmeyer, M., Kippenberg, T. J. & Marquardt, F. Cavity optomechanics. Reviews of Modern Physics 86, 1391-1452 (2014).

[2] Manucharyan, V. E., Koch, J., Glazman, L. & Devoret, M. Fluxonium: single Cooper pair circuit free of charge offsets. Science 326, 113-116 (2009).

[3] Zhang, H. et al. Universal fast flux control of a coherent, low-frequency qubit. Physical Review X 11, 011010 (2021).

[4] Najera-Santos, B.-L. et al. High-sensitivity ac-charge detection with a mhz-frequency fluxonium qubit. Physical Review X 14, 011007 (2024).

[5] Diósi, L. Models for universal reduction of macroscopic quantum fluctuations. Phys. Rev. A 40, 1165-1174 (1989). URL https://link.aps.org/doi/10.1103/PhysRevA.40.1165.

[6] Penrose, R. On Gravity's role in Quantum State Reduction. General Relativity and Gravitation 28, 581-600 (1996). URL https://doi.org/10.1007/BF02105068.

[7] Gely, M. F. & Steele, G. A. Superconducting electromechanics to test diósi-penrose effects of general relativity in massive superpositions. AVS Quantum Science 3 (2021).

[8] Rugar, D., Budakian, R., Mamin, H. & Chui, B. Single spin detection by magnetic resonance force microscopy. Nature 430, 329-332 (2004).

[9] Sansa, M. et al. Optomechanical mass spectrometry. Nature Communications 11, 3781 (2020).

[10] Gavartin, E., Verlot, P. & Kippenberg, T. J. A hybrid on-chip optomechanical transducer for ultrasensitive force measurements. Nature nanotechnology 7, 509-514 (2012).

[11] Catalini, L., Tsaturyan, Y. & Schliesser, A. Soft-clamped phononic dimers for mechanical sensing and transduction. Physical Review Applied 14, 014041 (2020).

[12] Jiang, W. et al. Efficient bidirectional piezo-optomechanical transduction between microwave and optical frequency. Nature communications 11, 1166 (2020).

[13] Andrews, R. W. et al. Bidirectional and efficient conversion between microwave and optical light. Nature physics 10, 321-326 (2014).

[14] Bagci, T. et al. Optical detection of radio waves through a nanomechanical transducer. Nature 507, 81-85 (2014).

[15] Hofheinz, M. et al. Synthesizing arbitrary quantum states in a superconducting resonator. Nature 459, 546-549 (2009).

[16] Arragoiz-Arriola, P. et al. Resolving the energy levels of a nanomechanical oscillator. Nature 571, 537-540 (2019).

[17] Bienfait, A. et al. Phonon-mediated quantum state transfer and remote qubit entanglement. Science 364, 368-371 (2019).

[18] Chu, Y. et al. Quantum acoustics with superconducting qubits. Science 358, 199-202 (2017).

[19] Viennot, J. J., Ma, X. & Lehnert, K. W. Phonon-number-sensitive electromechanics. Physical review letters 121, 183601 (2018).

[20] Lee, N. R. et al. Strong dispersive coupling between a mechanical resonator and a fluxonium superconducting qubit. PRX Quantum 4, 040342 (2023).

[21] Peterson, R. et al. Laser cooling of a micromechanical membrane to the quantum backaction limit. Physical review letters 116, 063601 (2016).

[22] Safavi-Naeini, A. H. et al. Observation of quantum motion of a nanomechanical resonator. Physical Review Letters 108, 033602 (2012).

[23] Underwood, M. et al. Measurement of the motional sidebands of a nanogram-scale oscillator in the quantum regime. Physical Review A 92, 061801 (2015).

[24] Xia, Y. et al. Motional sideband asymmetry of a solid-state mechanical resonator at room temperature. arXiv preprint arXiv:2408.06498 (2024).

[25] Clerk, A. A., Devoret, M. H., Girvin, S. M., Marquardt, F. & Schoelkopf, R. J. Introduction to quantum noise,
measurement, and amplification. *Reviews of Modern Physics* **82**, 1155–1208 (2010).

[26] Aguado, R. & Kouwenhoven, L. P. Double quantum dots as detectors of high-frequency quantum noise in mesoscopic conductors. *Phys. Rev. Lett.* **84**, 1986–1989 (2000). URL https://link.aps.org/doi/10.1103/PhysRevLett.84.1986.

[27] Gavish, U., Levinson, Y. & Imry, Y. Detection of quantum noise. *Phys. Rev. B* **62**, R10637–R10640 (2000). URL https://link.aps.org/doi/10.1103/PhysRevB.62.R10637.

[28] Schoelkopf, R., Clerk, A., Girvin, S., Lehnert, K. & Devoret, M. Qubits as spectrometers of quantum noise. *Quantum noise in mesoscopic physics* 175–203 (2003).

[29] This estimate relies on several debated assumptions; for a detailed discussion see App. B, Ref. [7] and references therein.

[30] Yuan, M., Singh, V., Blanter, Y. M. & Steele, G. A. Large cooperativity and microkelvin cooling with a three-dimensional optomechanical cavity. *Nature communications* **6**, 8491 (2015).

[31] Yuan, M., Cohen, M. A. & Steele, G. A. Silicon nitride membrane resonators at millikelvin temperatures with quality factors exceeding 108. *Applied Physics Letters* **107** (2015).

[32] Noguchi, A. et al. Ground state cooling of a quantum electromechanical system with a silicon nitride membrane in a 3d loop-gap cavity. *New Journal of Physics* **18**, 103036 (2016).

[33] Gambetta, J. et al. Qubit-photon interactions in a cavity: Measurement-induced dephasing and number splitting. *Physical Review A* **74**, 042318 (2006).

[34] Koch, J. et al. Charge insensitive qubit design derived from the cooper pair box. *Physical Review A* **76**, 042319 (2007).

[35] Unterreithmeier, Q. P., Weig, E. M. & Kotthaus, J. P. Universal transduction scheme for nanomechanical systems based on dielectric forces. *Nature* **458**, 1001–1004 (2009).

[36] Yu, P.-L., Purdy, T. & Regal, C. Control of material damping in high-q membrane microresonators. *Physical review letters* **108**, 083603 (2012).

[37] Cujia, K., Boss, J. M., Herb, K., Zopes, J. & Degen, C. L. Tracking the precession of single nuclear spins by weak

measurements. *Nature* **571**, 230–233 (2019).

[38] Arcizet, O., Cohadon, P.-F., Briant, T., Pinard, M. & Heidmann, A. Radiation-pressure cooling and optomechanical instability of a micro-mirror. *Nature* **444**, 71–74 (2006).

[39] Schliesser, A., Del'Haye, P., Nooshi, N., Vahala, K. J. & Kippenberg, T. J. Radiation pressure cooling of a micromechanical oscillator using dynamical backaction. *Physical Review Letters* **97**, 243905 (2006).

[40] Tsaturyan, Y., Barg, A., Polzik, E. S. & Schliesser, A. Ultracoherent nanomechanical resonators via soft clamp- ing and dissipation dilution. *Nature nanotechnology* **12**, 776–783 (2017).

[41] Ivanov, E. et al. Edge mode engineering for optimal ul- tracoherent silicon nitride membranes. *Applied Physics Letters* **117** (2020).

[42] Rouchon, P. A tutorial introduction to quantum stochas- tic master equations based on the qubit/photon system. *Annual Reviews in Control* **54**, 252–261 (2022).

[43] Guilmin, P., Rouchon, P. & Tilloy, A. Correlation functions for realistic continuous quantum measurement. *IFAC-PapersOnLine* **56**, 5164–5170 (2023).

[44] Sayrin, C. et al. Real-time quantum feedback prepares and stabilizes photon number states. *Nature* **477**, 73–77 (2011).

[45] Bild, M. et al. Schrödinger cat states of a 16-microgram mechanical oscillator. *Science* **380**, 274–278 (2023).

[46] Patange, H. RF Transduction in a DC-Biased Supercon- ducting Electromechanical System. PhD thesis, Sorbonne Université (2025).

[47] Vakhtel, T. et al. Tunneling of fluxons via a Josephson resonant level. *Physical Review B* **110**, 045404 (2024).

[48] Ardati, W. et al. Using bifluxon tunneling to protect the fluxonium qubit. *Phys. Rev. X* **14**, 041014 (2024).

[49] Masluk, N. A., Pop, I. M., Kamal, A., Minev, Z. K. & Devoret, M. H. Microwave characterization of Joseph- son junction arrays: Implementing<? format?> a low loss superinductance. *Physical review letters* **109**, 137002 (2012).

[50] Viola, G. & Catelani, G. Collective modes in the fluxo- nium qubit. *Physical Review B* **92**, 224511 (2015).