# Thermodynamics of a single mesoscopic phononic mode

Ilya Golokolenov, Arpit Ranadive, Luca Planat, Martina Esposito, Nicolas

Roch, Xin Zhou, Andrew Fefferman, Eddy Collin

► To cite this version:

Ilya Golokolenov, Arpit Ranadive, Luca Planat, Martina Esposito, Nicolas Roch, et al.. Thermodynamics of a single mesoscopic phononic mode. Physical Review Research, 2023, 5 (1), pp.013046. 10.1103/PhysRevResearch.5.013046. [hal-03765514]

HAL Id: hal-03765514

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

Submitted on 28 Feb 2023

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 CC BY 4.0 - Attribution - International License
Thermodynamics of a single mesoscopic phononic mode

Ilya Golokolenov,¹ Arpit Ranadive,¹ Luca Planat,¹,² Martina Esposito,¹,¹ Nicolas Roch,¹ Xin Zhou,¹,³
Andrew Fefferman,¹ and Eddy Collin¹,*,

¹Institut Néel UGA, Centre National de la Recherche Scientifique, 25 Rue des Martyrs, 38042 Grenoble, France

²Silent Waves, 25 Rue Ponsard, 38100 Grenoble, France

³IEMN, Université de Lille, Centre National de la Recherche Scientifique, Avenue Henri Poincaré,
59650 Villeneuve d'Ascq, France

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

(Received 19 August 2022; accepted 13 January 2023; published 25 January 2023)

In recent decades, the laws of thermodynamics have been pushed down to smaller and smaller scales, within the theoretical field of stochastic thermodynamics and state-of-the-art experiments performed on microfabricated mesoscopic systems. These measurements concern thermal properties of electrons, photons, and mesoscopic mechanical objects. Here we report on the measurements of thermal fluctuations of a *single mechanical mode* in equilibrium with a heat reservoir. The device under study is a nanomechanical beam with a first flexural mode resonating at 3.8 MHz, cooled down to temperatures in the range from 100 to 400 mK. The technique is constructed around a microwave optomechanical setup using a cryogenic high electron mobility transistor, and is based on *two parametric amplifications* implemented in series: an in-built optomechanical “blue detuned” pumping plus a traveling wave parametric amplifier stage. We demonstrate our ability to resolve energy fluctuations of the mechanical mode in real time up to the fastest relevant speed given by the mechanical relaxation rate. The energy probability distribution is then exponential, matching the expected Boltzmann distribution. The variance of fluctuations is found to be $(k_B T)^2$ with no free parameters. Our microwave detection floor is about three times the standard quantum limit at 6 GHz; the resolution of our fastest acquisition tracks reached about 100 phonons, and is directly related to the rather poor optomechanical coupling of the device ($g_0/2\pi \approx 0.5$ Hz). This result is deeply in the classical regime, but shall be extended to the quantum case in the future with systems presenting a much larger $g_0$ (up to $2\pi \times 250$ Hz), potentially reaching the resolution of a *single* mechanical quantum. We believe that it will open an experimental field, *phonon-based quantum stochastic thermodynamics*, with fundamental implications for quantum heat transport and macroscopic mechanical quantum coherence.

DOI: 10.1103/PhysRevResearch.5.013046

I. INTRODUCTION

Statistical physics, and by induction thermodynamics,
are the basis of our understanding of macroscopic properties
from the microscopic entities and laws that structure matter.
One of the key results is the second law of thermodynamics,
which explains the arrow of time from purely reversible mi-
croscopic processes [1]. Fluctuations $\delta X$ of a quantity $X$ are
then Gaussian and vanishingly small, leading to a well-defined
mean value $\langle X \rangle$.

But many of our intuitive understandings break down at
small scales: fluctuations can become as large as mean values,
and a specific class of theories known as fluctuation theo-
rems has been developed to describe them [2,3]. With today's
technologies, these concepts (and their related paradoxes) can
even be probed experimentally using mesoscale and nanoscale
devices. For instance, a Maxwell demon has been realized

by monitoring the charge in a single electron box, and feed-
ing back this information through a gate voltage controlled
by a computer; work is thus extracted [4]. Such electronic
systems are extremely promising, since one can cool them
down low enough (tens of millidegrees Kelvin) so that they
behave according to the laws of quantum mechanics. It should
then (at least in principle) be possible to probe the impact of
quantum coherence on thermodynamic concepts, which is the
new exciting field of quantum thermodynamics [5,6].

Beyond electric circuits, thermodynamics is conveying
concepts which are at the intersection of physics, chemistry,
and biology: after all, motion is a key ingredient there. Indeed,
the Landauer erasure principle has for instance been tested
using soft cantilevers and trapped colloids [7,8], demonstrat-
ing that erasing one bit of information produces a minimum
$k_B T$ ln(2) amount of heat. Similar stochastic thermodynamics
implementations have been realized on DNA molecules, e.g.,
monitoring their folding/unfolding and extracting work from
it [9]. Motion is thus at the core of the definition of heat: after
all, phonons are elementary (quasi-)particles constructed from
the (quantized) collective motion of atoms [10]. The quantum
limit of heat fluctuations [11] and phonon thermal conduc-
tance [12,13] are still subjects of debate today [14,15], with
very few experiments available in the literature [16,17]. Be-
sides, center-of-mass motion of mesoscopic objects is thought

<footnote>*eddy.collin@neel.cnrs.fr</footnote>
to be sensitive to *quantum aspects of gravity* (or any other fluctuating fields that might be at the source of wave-function collapse) [18]. Such mechanisms predict an imprint on mechanical fluctuations that might be measurable [19,20]. But of course, having a large mechanical object cold enough to host very few thermal excitations (population $n_{th} < 1$), while being in equilibrium with a heat reservoir, is a technological challenge. This has been recently demonstrated with a 15-µm aluminum drumhead device cooled down to 500 µK [21].

Studying quantum fluctuations at equilibrium of a macroscopic mechanical object down to the single quantum might thus be within reach in the near future [15]. We present in this paper a specific scheme enabling such kinds of measurements, based on microwave optomechanics. We focus here on a strict classical description of the experiment, which is a mandatory preliminary step towards the quantum realization, which we discuss in the Conclusion of the present paper. In Sec. II we describe the apparatus around which the experiment is constructed. In Sec. III we present the electric circuit modeling leading to the detected signal definition, while in Sec. IV the measurement protocol is explained and mathematically analyzed. The results are finally discussed in Sec. V. We separate what is directly the expression of expected properties of a single phononic mode in contact with a heat bath, from features (certainly material dependent) which are *not* expected. The former is an energy power spectrum typical of an *Ornstein-Uhlenbeck* process [1,22], with exponentially distributed fluctuations. The latter are visible as 1/f-type contributions to the spectra and out-of-equilibrium signatures, which will be discussed in the framework of the *two level systems* (TLSs) theory [23,24] (Sec. VI).

## II. EXPERIMENTAL SETUP

The optomechanical device we use has been presented in Ref. [25]. It consists of a 50-µm-long beam of 300-nm width and about 100-nm thickness, embedded in a microwave cavity (gap about 100 nm). The beam is a bilayer, made of high-stress silicon-nitride (SiN) covered with a thin layer of aluminum. The cavity is patterned with a 100-nm layer of niobium. The first in-plane flexural mode of the beam resonates at approximately $\Omega_m/(2\pi) = 3.8$ MHz, while the cavity resonates at $\omega_c/(2\pi) = 5.988$ GHz. The motion $x$ of the beam modulates the cavity's mode effective capacitance $C(x)$, leading to a frequency change characterized by the (first order) coupling strength $G = d\omega_c/dx$, which is measured to be about $G/(2\pi) \approx 1.8 \times 10^{13}$ Hz/m. The cavity is coupled evanescently (with an effective capacitance $C_c$) to a transmission line which enables one to connect the device to the drive/measurement circuitry.

A schematic of the setup is presented in Fig. 1. A first microwave pump tone is used to drive the optomechanics. This signal is split in order to create a “cancellation line” that opposes whatever remains from this pump at the input of the detection amplifying stage. This cancellation is performed by a computer that checks periodically the pump amplitude on the output, and adjusts a voltage-controlled attenuator and phase shifter on the cancellation line. We can suppress this signal by at least 60 dB. Three distinct microwave amplifiers are in use: first a traveling wave parametric amplifier (TWPA)

<figure><img src="image_3.png" /><figcaption>FIG. 1. Experimental optomechanical setup. A microwave pump tone drives the mechanical mode, while the sideband signal is amplified by a TWPA, a cryogenic HEMT, plus a room temperature HEMT. The TWPA is powered by a separate pump tone, and is protected from saturation due to the optomechanical drive by means of a cancellation tone. The measurement is performed through down-mixing followed by a lock-in amplifier (see text for details). *Bottom inset*: Measurement scheme; a pump tone (of power $P_{in}$) is applied at $\omega_p = \omega_c + \Omega_m$ (“blue pumping” arrow), with $\omega_c$ the cavity resonance frequency (whose susceptibility is shown in violet, arb. units). *Top inset*: Mechanical sideband spectrum measured at $\omega_c$, imprinted by the amplified Brownian NEMS motion at $\Omega_m$ (Lorentz curve of half-height width $\Gamma_{eff}$, arb. height; see text).</figcaption></figure>

that is powered by a separate pump tone [26]. The characteristics of this device are explicitly given in Appendix C, and lead to an equivalent noise at its input of 0.8 K (±0.1 K) for the whole chain (noise figures being quoted at the readout frequency). Cancellation and TWPA pump tones are added to the signal line by means of a power combiner and a directional coupler (light blue rectangles). The ambient noise coming from the pumps’ injection lines is attenuated by about 50 dB (orange zigzag blocks in Fig. 1). Besides, both the Nano Electro Mechanical System (NEMS) cell and the TWPA are protected by (two-stage) circulators (yellow disks; the orange squares are $Z_0 = 50 \, \Omega$ loads). The signal is then further amplified by a cryogenic high electron mobility transistor (HEMT) amplifier from Low Noise Factory (LNF) with noise temperature 2.5 K, and then a room temperature HEMT. We finally mix down the signal with a local oscillator (LO, a microwave tone shifted from $\omega_c$ by a fixed $\Delta\omega$ frequency) and a Marki mixer (brown circle). This megahertz signal at $\Delta\omega$ is finally fed into a Zurich Instruments (ZI) high frequency lock-in amplifier that is used to demodulate and digitize the data.

The setup is mounted on a commercial BlueFors cryostat, and experiments are performed while regulating the mixing chamber base plate from 100 to 400 mK. At higher temperatures, the TWPA amplifier stops working properly, while at lower temperatures an internal optomechanical instability known in the community as “spikes” corrupts the results [25]. These points shall be commented upon in more details in Appendices C and D, respectively.
<figure><img src="image_4.png" /><figcaption>FIG. 2. Electric circuit model describing the microwave cavity with a movable capacitor (the boxed “NEMS + CAVITY” in Fig. 1; see text).</figcaption></figure>

### III. DETECTED SIGNAL

Our detection scheme is based on “blue detuned” optomechanical pumping (bottom inset Fig. 1) [25,27]. A microwave pump tone (of power $P_{in}$) is applied at frequency $\omega_p = \omega_c + \Omega_m$. The (dynamical) back-action of the light field onto the mechanics leads to antidamping, transforming the mechanical relaxation rate $\Gamma_m$ into $\Gamma_{eff}$ (see Appendix B for details):

$$ \Gamma_{\text{eff}} = \Gamma_{\text{m}} + \Gamma_{\text{opt}}, \quad (1) $$

$$ \Gamma_{\text{opt}} = -2G^2 \left( \frac{\kappa_{\text{ext}}/2}{\kappa_{\text{tot}}} \right) \frac{P_{\text{in}}}{\Omega_m \omega_c k_0}, \quad (2) $$

with $\kappa_{ext}$ the coupling rate of the cavity to the transmission line, and $\kappa_{tot} = \kappa_{ext} + \kappa_{in}$ the total cavity relaxation rate ($\kappa_{in}$ stands for internal losses). Since $\Gamma_{eff} < \Gamma_m$, the motion is (parametrically) amplified. We fit about 100 kHz for $\kappa_{ext}$ and 200–150 kHz for $\kappa_{tot}$ (decreasing with increasing $P_{in}$) on the susceptibility curves (Fig. 1 bottom inset, violet line) [25]. This decrease of losses with microwave power is presumably due to the presence of TLSs (see discussion below). Note the factor 1/2 that arises in Eq. (2) because of evanescent coupling. We write $k_0$ and $m_0$ for the effective spring constant and mass of the mechanical mode, respectively, having $\Omega_m^2 = k_0/m_0$.

In the following, we describe the measurement in classical terms, following the electric circuit modeling of Refs. [28,29]; only when appropriate shall we make the link to quantum theory. For the modeling, we shall assume $\omega \sim \omega_c$, and $Z_0 C_c \omega_c \ll 1$. The voltage $V_{out}$ to be detected is written

$$ V_{\text{out}}(t) = -\omega_c^2 C_c \frac{Z_0}{2} \phi(t) + V_{\text{amp}}(t), \quad (3) $$

with $V_{amp}$ the voltage noise of the detection port, and $\phi$ the generalized flux (primitive integral of the voltage drop) across the microwave cavity. This cavity is modeled as an RLC circuit, with $C_t = C(0)$ the effective capacitance of the mode in the absence of motion and $L$ its effective inductance, such that $\omega_c^2 = 1/(LC_t)$ (see Fig. 2). The microwave damping is modeled through two effective resistors in parallel, an internal $R_{in}$ for all material-dependent losses, and a $R_{ext}$ corresponding to the leakage towards the outside (beyond capacitance $C_c$, not shown). The total losses simply write $1/R_t = 1/R_{in} + 1/R_{ext}$, with $\kappa_{in} = 1/(R_{in}C_t)$ and $\kappa_{ext} = 1/(R_{ext}C_t)$. One defines

$$ \kappa_{\text{ext}} = \omega_c^2 \frac{Z_0}{2} \frac{C_c^2}{C_t}, \quad (4) $$

from circuit theory, considering that the end load of the output line is equal to $Z_0$ [28]. The linear coupling between mechanics and microwaves generates a comb in the output signal, that presents components at $\omega_n = \omega_c + n\Omega_m$ with $n \in \mathbb{Z}$. This leads to the expressions [28]

$$ V_{\text{out}}(t) = \sum_{n=-\infty}^{\infty} \frac{V_{M,n}(t) e^{-i\omega_n t} + V_{M,n}(t)^* e^{+i\omega_n t}}{2}, \quad (5) $$

$$ \phi(t) = \sum_{n=-\infty}^{\infty} \frac{\mu_n(t) e^{-i\omega_n t} + \mu_n(t)^* e^{+i\omega_n t}}{2}, \quad (6) $$

$$ V_{\text{amp}}(t) = \sum_{n=-\infty}^{\infty} \frac{V_{N,n}(t) e^{-i\omega_n t} + V_{N,n}(t)^* e^{+i\omega_n t}}{2}, \quad (7) $$

with $V_{M,n}$, $\mu_n$, and $V_{N,n}$ the respective (complex) amplitudes in the frames rotating at $\omega_n$. The motion $x$ is itself expressed in a frame rotating at $\Omega_m$:

$$ x(t) = \frac{x_0(t) e^{-i\Omega_m t} + x_0(t)^* e^{+i\Omega_m t}}{2}, \quad (8) $$

with $x_0$ the complex motion amplitude. The drive voltage $V_d$ created by the microwave generator can be defined as

$$ V_d(t) = \frac{1}{2} V_p e^{-i\omega_p t} + \frac{1}{2} V_p^* e^{+i\omega_p t} + \sum_{n=-\infty}^{\infty} \frac{V_{P,n}(t) e^{-i\omega_n t} + V_{P,n}(t)^* e^{+i\omega_n t}}{2}, \quad (9) $$

with $V_{P,n}$ the complex noise amplitude of the drive field around frequency $\omega_n$. The injected power is thus $P_{in} = |V_p|^2/(2Z_0)$.

In the sideband-resolved limit $\Omega_m \gg \kappa_{tot}/2$, only three components of the comb are relevant: the pump tone at $\omega_p$, and the two sidebands at $\omega_p \pm \Omega_m$ (i.e., $n = \pm 1$). For a blue detuned scheme, $\omega_p = \omega_c + \Omega_m$. In this case, only the sideband at $\omega_p - \Omega_m = \omega_c$ is measurable, the other one being strongly suppressed. The corresponding voltage amplitude $V_{M,n=-1}$ is found to be [28]

$$ V_{M,-1}(t) \approx -G x_0^*(t) \frac{\kappa_{\text{ext}}/2}{\kappa_{\text{tot}}} \frac{V_p}{\Omega_m} + \frac{\kappa_{\text{ext}}}{\kappa_{\text{tot}}} V_{P,-1}(t) + V_{N,-1}(t), \quad (10) $$

keeping only noise terms inside the cavity. The reverse scheme is “red detuned” pumping (pump tone applied at $\omega_c - \Omega_m$), where the sign of $\Gamma_{opt}$ in Eq. (2) is opposite (– should read +, leading to attenuation instead of amplification). The measurable sideband is again the one at $\omega_c$, but it corresponds now to $n = +1$. The voltage amplitude $V_{M,n=+1}$ is then similar to Eq. (10), with a change of sign in front of $G$ and a replacement $x_0^* \to x_0$. The two blue and red detuned schemes are explicitly compared in the following in order to validate the data analysis.

The mixing process can be formally written as

$$ \alpha[V_{\text{out}}(t) \times \cos(\omega_d t)]_{\text{filter}} = V_{\text{meas}}(t), \quad (11) $$

with $\omega_d$ the frequency of the LO in Fig. 1; we define $\omega_n - \omega_d = \Delta\omega$ the demodulation frequency ($n = \pm 1$ depending on the scheme, blue or red detuned pumping). In Eq. (11), the term “filter” means that only the component at $\omega_n - \omega_d$
is processed, while the one at $\omega_n + \omega_d$ is filtered out. The coefficient $\alpha$ conveniently contains all calibration from the detection chain, which is discussed in more detail in Appendix C. One subtlety arises concerning the detection noise background: the component at $\omega'_n = \omega_d - \Delta\omega = \omega_n - 2\Delta\omega$ is actually mixed down equally well as $\omega_n$ in this process, and adds up to the initial noise background $V_{N,n}$ appearing in Eq. (10). The voltage digitized by the lock-in amplifier therefore reads

$$V_{\text{meas}}(t) = \frac{\alpha}{2} \left[ -G \bar{x}(t) \frac{\kappa_{\text{ext}}/2}{\kappa_{\text{tot}}} \frac{V_p}{\Omega_m} + \frac{\kappa_{\text{ext}}}{\kappa_{\text{tot}}} V_p(t) + V_N(t) \right], \quad (12)$$

having defined for blue detuned pumping

$$V_P(t) = \frac{V_{P,-1}(t)e^{-i\Delta\omega t} + V_{P,-1}^*(t)e^{+i\Delta\omega t}}{2}, \quad (13)$$

$$V_N(t) = \frac{1}{2}[V_{N,-1}(t) + V_{N,n'}^*(t)]e^{-i\Delta\omega t} + \frac{1}{2}[V_{N,-1}(t) + V_{N,n'}^*(t)]^*e^{+i\Delta\omega t}, \quad (14)$$

$$\bar{x}(t) = \frac{x_0^*(t)e^{-i\Delta\omega t} + x_0(t)e^{+i\Delta\omega t}}{2}, \quad (15)$$

the noise component due to the pump tone, the amplification chain noise background, and the “effective motion” signal, respectively. Note that the latter corresponds exactly to Eq. (8) under the replacement $\Omega_m \to -\Delta\omega$. In Eq. (14), $n'$ refers to the component at $\omega_c - 2\Delta\omega$; no such term exists for the input noise, in the limit $\kappa_{\text{tot}} \ll \Delta\omega$ which is taken experimentally [we chose arbitrarily $\Delta\omega/(2\pi) = \Omega_m/(2\pi) + 2$ MHz, within the lock-in bandwidth]. A similar writing holds for red detuned pumping. For more details on the classical circuit theory, we refer the interested reader to Ref. [28].

Up to this point, the time-dependent variables introduced above ($V_P$, $V_N$ and $\bar{x}$) correspond mathematically to one realization of the stochastic processes they correspond to. Let us define $V_{\text{meas}}(\omega) = \mathcal{F}\mathcal{T}[V_{\text{meas}}(t)](\omega)$ the voltage in frequency space, where $\mathcal{F}\mathcal{T}[f(t)](\omega) = \int_{-\infty}^{+\infty} f(t)e^{-i\omega t} dt$. Using Kubo’s notations [30], we define the instantaneous (i.e., before ensemble averaging) voltage power spectrum as $2\pi S_{V_{\text{meas}}}(\omega) \delta_0(\omega' - \omega) = V_{\text{meas}}(\omega) V_{\text{meas}}(\omega')^*$ [and similar expressions for the constitutive random variables Eqs. (13)–(15)], with $\omega \in [-\infty; +\infty]$. From Eq. (12), we obtain

$$S_{V_{\text{meas}}}(\omega) = \frac{|\alpha|^2}{4} \left[ G^2 \left( \frac{\kappa_{\text{ext}}/2}{\kappa_{\text{tot}}} \right)^2 \frac{|V_p|^2}{\Omega_m^2} S_{\bar{x}}(\omega) + \left( \frac{\kappa_{\text{ext}}}{\kappa_{\text{tot}}} \right)^2 S_{V_p}(\omega) + S_{V_N}(\omega) + \text{cross terms} \right], \quad (16)$$

with “cross terms” referring to all cross-correlation spectra. Knowing that input and output noises are uncorrelated, and that correlations between $\bar{x}$ and $V_P$ and $V_N$ (which are responsible for sideband asymmetry [29]) are negligible here, these terms shall vanish when computing statistical properties in Sec. IV.

The voltage power spectral densities $S_{V_P}$ and $S_{V_N}$ are reasonably flat over the width of the microwave cavity resonance $\kappa_{\text{tot}}$: we can therefore treat them as being white. Besides, since

the voltage noise amplitudes $V_{N,-1}$ and $V_{N,n'}$ are uncorrelated, and essentially of equal intensity, the power spectral density $S_{V_N}$ is twice the level measured before the mixer. This is the price to pay in the down-mixing process (see Appendix C). Dividing Eq. (16) by $Z_0$, one obtains the (double-sided) power spectral density of detected power [in W/(rad/s), therefore joule]. Further dividing by $\hbar \omega_c$ one converts it into a photon flux power spectral density [in (photons/s)/(rad/s), therefore photons]:

$$S_{\bar{\phi}}(\omega) = 2G^2 \left( \frac{\kappa_{\text{ext}}/2}{\kappa_{\text{tot}}} \right)^2 \frac{P_{\text{in}}}{\Omega_m \omega_c k_0} S_n(\omega) + \left[ \left( \frac{\kappa_{\text{ext}}}{\kappa_{\text{tot}}} \right)^2 S_{\text{in}\phi}(\omega) + S_{\text{out}\phi}(\omega) \right], \quad (17)$$

where we dropped the calibration factor $|\alpha|^2/4$ for simplicity. Note that the formula reads the same for red detuned pumping. Expressing input and output noises $S_{\text{in}\phi}$ and $S_{\text{out}\phi}$ in terms of photons enables one to evaluate the technique for future quantum measurements [31]: we reach about three photons, which is the state of the art [26] (see discussion in Appendix C). Explicitly,

$$S_n(\omega) = \frac{k_0 S_{\bar{x}}(\omega)}{\hbar \Omega_m}, \quad (18)$$

which corresponds to the instantaneous mechanical energy power spectral density (expressed in phonons), peaked around $\Delta\omega$ (instead of $\Omega_m$). All of these classical spectra are even, therefore experimentally what is presented for a quantity $X$ is the single-sided $2S_X(f > 0)$, with $f = \omega/(2\pi)$ in Hz.

## IV. MEASUREMENT PROTOCOL

In practice, each measured power spectrum is acquired over a finite time $\delta t$ around time $t$: $\langle S_{\bar{\phi}}(\omega) \rangle_{\delta t}(t)$. Neglecting for now the photon background noise in Eq. (17), this is simply proportional to the mechanical energy spectrum $\langle S_n(\omega) \rangle_{\delta t}(t)$, within a coefficient $(\kappa_{\text{ext}}/2)|\Gamma_{\text{opt}}| \propto P_{\text{in}}$. If the time span $\delta t$ is infinitely long, this quantity should become $t$ independent and reproduce the well-known Lorentzian mechanical spectrum, with a half-height width $\Gamma_{\text{eff}}$ and an area $k_B T/(\hbar\Omega_m)$ (schematic in Fig. 1, top inset) [27]. This is not perfectly true experimentally because of $1/f$ drifts (see the discussion below). On the other hand if $\delta t$ becomes infinitely short, one is supposed to obtain (mathematically) a Dirac peak (essentially, the motion is a well-defined oscillation at $\Omega_m$ for timescales $\ll 1/\Gamma_{\text{eff}}$), which fluctuates over longer times $t$. Again this suffers from experimental limitations: the frequency resolution is inversely proportional to the acquisition speed, which means that the peak is “blurred” over a frequency span $1/\delta t$. This aspect is explicitly discussed in Appendix F. The actual experimental dependence of measured spectra for not-too-long and not-too-short speeds is shown in Fig. 4, top insets.

From the photon flux spectra $\langle S_{\bar{\phi}}(\omega) \rangle_{\delta t}$, we define the sideband peak power (in photons/s),

$$\langle \dot{E}_{\bar{\phi}} \rangle_{\delta t}(t) = \frac{1}{2\pi} \int_{-\infty}^{+\infty} \langle S_{\bar{\phi}}(\omega) \rangle_{\delta t}(t) d\omega, \quad (19)$$
which is then proportional to the mechanical energy (in phonons):

$$ \langle E_n \rangle_{\delta t}(t) = \frac{1}{2\pi} \int_{-\infty}^{+\infty} \langle S_n(\omega) \rangle_{\delta t}(t) d\omega. \quad (20) $$

Technically, the area of the sideband is not obtained through integration (which would be very noisy), but rather with a Lorentz fit (from PYTHON routines SciPy optimize .curve fit [32] and lmfit .minimize [33], black line top inset on left, Fig. 4). Indeed, the background noise that we neglected up to now impacts strongly the quality of the numerical analysis, especially at the fastest acquisition speeds (smallest $\delta t$; see Appendix F). Besides, at not-too-fast speeds, it enables us to fit also both the peak width (extracting thus $\Gamma_m$ from $\Gamma_{eff}$) and the peak position $\Omega_m$ (defined from the demodulation reference). It turns out that these parameters are not constant, and fluctuate over time; this is explicitly discussed in Sec. VI.

Data acquisition is performed over a time $\Delta T$, that we arbitrarily chose as being 1000 $\delta t$ for convenience (each of our sets is made of $N = 1000$ points exactly). We construct

$$ \langle \dot{E}_{\varphi} \rangle = \langle \langle \dot{E}_{\varphi} \rangle_{\delta t}(t) \rangle_{\Delta T}, \quad (21) $$

$$ C_{\dot{E}_{\varphi}}(\tau) = \langle \langle \dot{E}_{\varphi} \rangle_{\delta t}(t) \langle \dot{E}_{\varphi} \rangle_{\delta t}(t-\tau) \rangle_{\Delta T}, \quad (22) $$

the photon flux mean value and the corresponding autocorrelation function, respectively. In Eq. (22), $\tau$ takes discrete values from -1000 $\delta t$ to +1000 $\delta t$. Similar expressions hold for the mechanical energy with $\langle E_n \rangle$ and $C_{E_n}(\tau)$. For each set, the whole procedure is repeated from 2 to 100 times (depending on acquisition speed) in order to improve the quality of the data and assess the impact of 1/f drifts on quoted values [indeed, ideally Eqs. (21) and (22) should be t independent]. As always, taking the experimental averaging over $\Delta T$ for an ensemble average is based on a fundamental hypothesis: ergodicity. This assumption is not that straightforward here, precisely because of the 1/f detected features; this point shall be specifically discussed in the following.

## V. PHONON MEAN POPULATION AND POWER SPECTRAL DENSITY

In order to interpret the experiment, we first remind the reader of basic classical statistical physics results [1,22]. We model the mechanical mode as being in contact with a thermodynamic bath at temperature $T$ (the cryostat), and an optical bath at an effective temperature $T_{opt}$ [29]. We define for the mechanical energy $E_m = \hbar \Omega_m E_n$

$$ \Delta E = E_m - \langle E_m \rangle, \quad (23) $$

the amplitude of fluctuations around the mean $\langle E_m \rangle$. This quantity follows the dynamics equation

$$ \frac{d\Delta E(t)}{dt} = -(\Gamma_m + \Gamma_{\text{opt}})\Delta E(t) + R_m(t) + R_{\text{opt}}(t), \quad (24) $$

with $R_m$ and $R_{opt}$ the two random energy flows associated with each bath. They verify

$$ C_m(\tau) = \langle R_m(t)R_m(t-\tau) \rangle = 2\Gamma_m(k_B T)^2 \delta_0(\tau), \quad (25) $$

$$ C_{\text{opt}}(\tau) = \langle R_{\text{opt}}(t)R_{\text{opt}}(t-\tau) \rangle = 2|\Gamma_{\text{opt}}|(k_B T_{\text{opt}})^2 \delta_0(\tau), \quad (26) $$

$$ 0 = \langle R_m(t)R_{\text{opt}}(t-\tau) \rangle, \quad (27) $$

<figure><img src="image_5.png" /><figcaption>FIG. 3. Phonon population $\langle E_n \rangle / \mathcal{G}$ as a function of applied pump power $P_{in}$ measured at 200 mK. Labels stand for different measuring schemes. The dashed line is a simple guide for the eye (see text; note the log-log scale spanning two orders of magnitude in power). Inset: Intrinsic phonon population (value extrapolated at zero pump power) as a function of temperature; error bars are from the reproducibility scatter of the main graph; the line is calculated from theory ($k_B T / [\hbar \Omega_m]$), see text).</figcaption></figure>

meaning that they have no intrinsic finite correlation time, with the last line stipulating that the two baths are uncorrelated. One obtains

$$ \langle E_m \rangle = \frac{\Gamma_m(k_B T) + |\Gamma_{\text{opt}}|(k_B T_{\text{opt}})}{\Gamma_m + \Gamma_{\text{opt}}}, \quad (28) $$

$$ S_{\Delta E}(\omega) = \frac{2[\Gamma_m(k_B T)^2 + |\Gamma_{\text{opt}}|(k_B T_{\text{opt}})^2]}{(\Gamma_m + \Gamma_{\text{opt}})^2 + \omega^2}, \quad (29) $$

for the mean energy and power spectral density. Consider $\Gamma_{opt} = 0$; then one recovers the simple case of a canonical ensemble with a bath at temperature $T$, a situation which has been studied experimentally with macroscopic objects [34]. Energy fluctuations are Gaussian [35], a simple consequence of the central limit theorem because of the large number of degrees of freedom involved. However, the single mode statistics is different: it is described by the Boltzmann distribution $p(E) = e^{-E/(k_B T)}/Z$, where in the classical limit the partition function $Z = k_B T$ [ensuring $\int_0^{+\infty} p(E) dE = 1$]. One can recover these results from pure (classical) mechanical arguments, starting from the fluctuation-dissipation theorem and its associated Langevin force (see Appendix A).

Consider now $\Gamma_{opt} \neq 0$ but $T_{opt} \approx 0$. Equation (28) reads

$$ \langle E_m \rangle = \mathcal{G} k_B T, \quad (30) $$

$$ \mathcal{G} = \frac{\Gamma_m}{\Gamma_m + \Gamma_{\text{opt}}}, \quad (31) $$

with $\mathcal{G}$ the optomechanics gain. We plot in Fig. 3 the mean population of the mechanical mode, recalculated from the mean measured photon flux Eq. (21). The measurements have been performed with both red detuned ($\mathcal{G} < 1$) and blue
<figure><img src="image_6.png" /><figcaption>FIG. 4. The central part shows the power spectral density (PSD) of the photon flux $S_{\phi}$ measured at $T = 200$ mK, with a pump power $P_{in} = 10$ nW. Different colors correspond to different acquisition speeds ($\delta t$) and durations ($\Delta T = 1000$ $\delta t$); see labels $\delta t = \{t_1, t_2, t_3, t_4\}$ in the top of the figure. The black solid line is a fit, which corresponds to the expected spectrum of fluctuations [demonstrating the high-frequency cutoff at $\Gamma_{eff}/(2\pi)$], but also reveals an unexpected $1/f$ addendum. Note the overlap between different spectra ("stitching"; see text). Top inset: Typical raw signals for different $\delta t$ (see legend $t_1 - t_4$ and color code), which are fitted with a Lorentzian function (black line demonstrated on the left panel) to define sideband peak area $(\tilde{E}_{\phi})_{\delta t}$, position [shift from reference value $\Delta\omega/(2\pi)$], and linewidth $\Gamma_{eff}/(2\pi)$; for slow acquisitions, the peak is well defined (left graph; see black line fit), while we lose resolution as the acquisition speed is increased; note the different lines plotted in the two right panels, taken from the same respective statistical batches; at the fastest, the linewidth is essentially given by the sampling (see text). Bottom inset: Corresponding probability distribution functions (PDF) vs amplitude normalized to mean (equivalent to $E_n/(\langle E_n \rangle)$), demonstrating the change of shape as the acquisition speed increases (from Gaussian to exponential; see black lines and discussion in text).</figcaption></figure>

detuned ($G > 1$) schemes, with the TWPA amplifier on and off (see labels). All the data are in very good agreement, and the remaining scatter is due to the reproducibility of the measurement ($1/f$ drifts). Note that $\langle E_n \rangle$ (which in the quantum language is nothing but the mode's thermal population $n_{th}$) verifies $\langle E_n \rangle \gg 1$: we are deeply in the classical limit at all studied temperatures.

In Fig. 3, we see that the measured mean phonon population (normalized to gain $G$, ranging from 0.5 to 3) is actually increasing as we increase the injected microwave power $P_{in}$, a phenomenon known as *technical heating* in the community

(see, e.g., Refs. [21,25]). A very natural guess is to assume that the effective temperature of the optical bath $T_{opt} \neq 0$ increases and becomes relevant [see Eq. (28)]. This is actually inconsistent. In the classical formalism (with a large enough cavity photon population), $T_{opt} = T_{cav} \Omega_m / \omega_c$ from back-action; we safely neglect sideband asymmetry, which would induce a reversed-in-sign correction for red and blue pumping which is not observed here [29]. $T_{cav}$ is the effective temperature of the cavity, which could be due to both a real physical heating of the chip (microwave absorption) or to out-of-equilibrium photons arising from the generator noise (see Appendix C).
At the highest powers (around 15–20 nW), $T_{cav}$ would reach about 300 K (about 1000 photons), which is unphysical (and obviously not observed when measuring directly the output spectrum of the cavity). We thus have to conclude that this effect has another (unknown) origin, with no clear temperature and power dependencies (the line in Fig. 3 is above all a guide to the eye, here a simple linear law). In the inset of Fig. 3, we plot the mean population extrapolated at zero power. The line is the theoretical prediction from Eq. (30) with $\Gamma_{opt} = 0$: $k_B T/(\hbar\Omega_m)$. The agreement is fairly good, with a slight (unexplained) deviation at low temperatures. Note that the scatter is essentially due to the reproducibility of the measurement, impacted by the $1/f$ noise. We shall comment on the features which are not understood in the following.

From the photon flux correlation function Eq. (22), we compute the fast Fourier transform (FFT, NumPy.correlate autocorrelation and NumPy FFT algorithms [36] in PYTHON), leading to the experimental spectrum definition

$$2S_{\dot{E}_{\psi}}(f) = \frac{\delta t}{N^2} \text{FFT}[C_{\dot{E}_{\psi}}(\tau)](f), \quad (32)$$

for different acquisition speeds $\delta t$, with the factor 2 on the left-hand side due to the experimental convention $f > 0$. These spectra are plotted in Fig. 4, main graph (see color code for $\delta t$, top insets). The normalization factor in Eq. (32) takes into account both the number of points $N$ of the discretized acquisition and the bandwidth $1/\delta t$. For not-too-fast acquisition speeds, the data overlap very well, demonstrating “stitching.” We conclude that ergodicity is well verified even at the slowest speeds, where $1/f$ drifts are non-negligible. However, the fastest tracks (only red data) should be rescaled because of the acquisition finite bandwidth (see Appendix G).

The full spectrum displayed in Fig. 4 is fit by the expression (black full line)

$$2S_{\dot{E}_{\psi}}(f) = \frac{A_f}{f} + \frac{S_0}{1 + \left(\frac{f}{\Gamma_{\text{eff}}/(2\pi)}\right)^2}. \quad (33)$$

The $\Gamma_{eff}$ is actually not fitted, but obtained from the known power dependence of the measured peak width (see Appendix B), demonstrating very good agreement with data: the mode cannot exchange energy with its environment at speeds exceeding its relaxation rate. The impact of detection noise and of the fitting routine is analyzed in Appendix F. $S_0$ then gives us the level of energy fluctuations while $A_f/f$ corresponds to an unexpected contribution that shall be discussed in Sec. VI.

Applying now the assumption $T_{opt} \approx 0$ to Eq. (29), one can write

$$S_{\Delta E}(|\omega| \ll \Gamma_{\text{eff}}) \approx G^2 \frac{2(k_B T)^2}{\Gamma_m}, \quad (34)$$

with $2(k_B T)^2/\Gamma_m$ the $\omega \to 0$ value that characterizes an unpumped mode ($\Gamma_{opt} = 0$). Reversing this expression and making use of the transduction coefficient between the optical and the mechanical fields, we can therefore recalculate from the best fit value of $S_0$ the actual thermodynamical mechanical fluctuation level. Comparison between red detuned ($G < 1$) and blue detuned ($G > 1$) schemes is discussed in Appendix E. Making this experiment at various temperatures, we

<figure><img src="image_7.png" /><figcaption>FIG. 5. Normalized spectrum level (no units) $\sqrt{\frac{\Gamma_m S_0}{G^2}}$ at different temperatures and pump powers (size of dots from small/low to big/large $P_{in}$). The line is the theoretical calculation with no free parameters (see text). Inset: Time trace of population fluctuations at the fastest acquisition speed (200 mK, power 10 nW); the line is the mean value; resolution is about 100 phonons; see text for discussion.</figcaption></figure>

present this quantity (normalized to the phonon energy $\hbar\Omega_m$) as a function of $T$ in Fig. 5. The black line is the theoretical prediction, with no free parameters, demonstrating very good agreement with data. The scatter seems to be due to our reproducibility and fitting capability (see Appendix F), with no specific link to the microwave power $P_{in}$. It confirms the magnitude of the variance $\langle\Delta E^2\rangle = 1/(2\pi) \int_{-\infty}^{+\infty} S_{\Delta E}(\omega)d\omega = (k_B T)^2$ in this canonical ensemble, the subtlety being that the specific heat associated with the single mode is precisely $k_B$ [34].

Finally, from the acquired time tracks of $\langle\dot{E}_{\psi}\rangle_{\delta t}(t)$ we can build histograms; this is done for each acquisition speed $\delta t$ (see bottom insets in Fig. 4). We plot them with an area normalized to 1 [directly compatible with a probability distribution function (PDF)], and an energy amplitude axis normalized such that the mean is also 1. At the fastest speed the shape is clearly exponential (right plot; the line is a theoretical function with no fit parameter). The standard deviation is then also 1. However, as we slow down the acquisition, the distribution becomes gradually Gaussian (see line middle-left plot of Fig. 4 bottom insets), with a smaller standard deviation, as it should. Note that the first PDF (on the left) does not display a fit, the histogram being quite impacted by the slow drifts of the (nonstationary) $1/f$ noise.

An attempt had been made in Ref. [21] to develop the measurement method producing Fig. 4. However, the resolution was very far from enough for this purpose (no TWPA was used), and the heavy averaging was essentially filtering out the thermodynamic contribution, leaving only what should have been the equivalent of our $1/f$ term (see discussion in Appendix H). On the other hand, we reach here at the fastest tracks a resolution of about 100 phonons (see real-time data in inset Fig. 5), being limited only by our relatively poor
<figure><img src="image_8.png" /><figcaption>FIG. 6. $A_f$ constant of $1/f$ noise contribution fitted at different powers and temperatures, in units of photons on chip. Dashed lines are exponential fits: $A_f = A_0(T)e^{P/P_0(T)}$ (see text). Inset: $P_0$ fit parameter from the main graph as a function of temperature $T$.</figcaption></figure>

optomechanical coupling $G$. In quantum mechanics terms with the zero point fluctuation $x_{ZPF} \approx 27$ fm, we have $g_0 = Gx_{ZPF} \approx 2\pi \ 0.5$ rad/s. Using drumhead aluminum devices in the future, one can reach couplings as high as $g_0 \approx 2\pi \ 250$ rad/s, winning therefore a factor $\approx 10^5$ on the detected signal (all other parameters being kept equivalent) [29].

## VI. 1/f NOISES

A striking unexpected feature observed in our measurements is the $1/f$ contribution to the photon flux fluctuations, main plot in Fig. 4. We show the fit parameter $A_f$ in Fig. 6 as a function of both injected power $P_{in}$ and temperature $T$. We observe that this coefficient can be fit by an exponential input power dependence (note the y axis of the figure), with a smooth temperature dependence (see inset).

The origin of this effect remains unknown. It is not even clear if it originates in the phonon or in the photon field (see discussion in Appendix E, comparing red and blue detuned pumping), which is why we characterize it in terms of photons. Besides, since $1/f$ drifts are responsible for very slow (close to $\omega \to 0$) dynamics, one could wonder whether this signature has to be linked to the *technical heating* of Fig. 3 (a supposedly true dc effect). Again, this remains an open question.

Besides, one should also keep in mind that the mechanical parameters $\Omega_m$ and $\Gamma_m$ are also fluctuating; this had been already reported in Ref. [21] in the framework of low-temperature optomechanics, but also in more conventional experiments [37,38]. Since we fit the mechanical response (on not-too-fast tracks), we can extract these parameters and compute their statistical properties. This is summarized in Fig. 7. The probability distributions look reasonably Gaussian, and the power spectral densities present a clear $1/f$ trend (see example in the inset). We find out that the damping noise is essentially constant in temperature $T$, which is consistent with findings from Ref. [38] taken at slightly higher temperatures.

<figure><img src="image_9.png" /><figcaption>FIG. 7. Allan variance of mechanical frequency $\Omega_m/(2\pi)$ (blue) and of mechanical damping $\Gamma_m/(2\pi)$ (green) as a function of temperature. Tracks are acquired with our slowest acquisition speed, over about 14 h. Sizes of dots are shown for different microwave powers (from small/low to big/large as in Fig. 5). Dashed lines are guides for the eyes. Inset: Typical power spectral density, here computed for frequency (at 200 mK, with drive power 10 nW); the line is a $1/f$ fit (see text); the doubling of the line comes from the folding of the negative frequency axis onto the positive one.</figcaption></figure>

On the other hand, the frequency noise *grows* as we cool down, a feature also seen in Ref. [21] down to much lower temperatures. Both damping and frequency fluctuations are of the same order as reported values for SiN beam devices cooled at cryogenic temperatures [38]. The scatter in Fig. 7 is rather large (as is usually the case when measuring $1/f$), but no specific drive power dependence can be seen. Again, the mechanism behind these features might be linked to the previous properties impacting energy fluctuations, but no microscopic theory has been formulated yet.

It is nonetheless tempting to imagine that TLSs are responsible for these facts. Indeed, low temperature mechanical properties of microfabricated structures are analyzed in the framework of this model [39,40], while experiments on microwave cavities also reported a growing frequency variance with lowering temperature, which was interpreted as a signature of interacting TLSs [41]. Furthermore, an internal instability of microwave optomechanics has been reported below typically 150 mK [25]. This feature, nicknamed “spikes,” is discussed in Appendix D. The deviation at low $T$ of the mechanical mode energy with respect to theory (Fig. 3 inset) might be related to this; the fact that the most impacting $1/f$ signatures grow as we cool down is a rather intriguing fact, that might suggest that a profound link exists between them.

## VII. CONCLUSION

We report on a technology that enables us to measure in real time the energy fluctuations of a mesoscopic mechanical mode. The setup is constructed around a state-of-the-art microwave optomechanical cryogenic platform, presenting a background noise of three times the standard quantum
limit (three photons at 6 GHz). The resolution is then about 100 phonons at our fastest acquisition rates (about 20 ms), with a 4-MHz mode. The limiting parameters are the intrinsic losses of the TWPA [26], and mostly the weak optomechanical coupling of the device we used here [25]. We believe that both aspects can be greatly improved, leading to an effective resolution on the detected sideband spectrum equivalent to a single mechanical quantum.

This means not only that these experimental capabilities surpass the best microfabricated calorimetric setups to date (zepto-joule calorimetry [42]), but also that the phonons themselves become the quantum bath being probed, performing quantum calorimetry with phonons, opening thus an experimental field [43,44].

The next experimental step is therefore to mount this measurement setup onto a cryostat enabling one to “brute-force” cool down to the quantum regime such MHz mechanical modes, as demonstrated in Ref. [21]. Quantum stochastic thermodynamics experiments would be within reach [15], but this requires one to further analyze the setup in quantum mechanics terms (the theoretical treatment presented here being purely classical).

The measurement is constructed around the observable $\hat{x} \propto \hat{a} + \hat{a}^{\dagger}$ (motion amplitude), not energy $\hbar\Omega_m \hat{n} \propto \hat{a}^{\dagger}\hat{a}$, which means that we shall not detect single-phonon tunnelings per se, but their (dispersive) imprint onto the optical field. This should nonetheless enable one to study single phonon events, transposing to mechanics what has been beautifully achieved for electrons. As an example, one should be able to demonstrate, at extremely low temperatures where the mechanical thermal population $n_{th} < 1$, how the system can be absolutely free of excitations over macroscopic timescales (similarly to electrons in a superconductor) [45]: a rather counterintuitive possibility, which essentially means that the system could be said to be $T = 0$ K exactly for a short period of time, while obviously on average $T > 0$ K is always recovered.

## ACKNOWLEDGMENTS

We wish to acknowledge the use of the Néel Institute Nanofab facility. The authors are grateful to J. L. Garden, O. Bourgeois, O. Maillet, and A. D. Armour for very useful discussions. We acknowledge support from the European Research Council under CoG ULT-NEMS Grant No. 647917 (E.C.) and StG UNIGLASS Grant No. 714692 (A.F.), and by the European Union’s Horizon 2020 Research and Innovation program under Grant No. 899561 (N.R.). The research leading to these results has received funding from the European Union’s Horizon 2020 Research and Innovation program, under Grant No. 824109, the European Microkelvin Platform.

## APPENDIX A: FROM LANGEVIN FORCE TO BOLTZMANN DISTRIBUTION

One can easily construct the energetics description of a mode from its motion, at least in the so-called high-Q limit. Consider the dynamics equation of a harmonic oscillator ($m_0$ being its mass, and $k_0$ its spring constant, $\Omega_m^2 = k_0/m_0$):

$$ \ddot{x} + \Gamma_m \dot{x} + \Omega_m^2 x = \delta F/m_0, \quad (A1) $$

with $\delta F(t)$ the Langevin force linked to the damping $\Gamma_m$ through the fluctuation-dissipation theorem. Both originate from a thermal bath at temperature $T$, and the stochastic force is by definition described by a centered $(\langle \delta F \rangle = 0)$ Gaussian probability distribution with correlator:

$$ C_{\delta F}(\tau) = \langle \delta F(t) \delta F(t - \tau) \rangle = 2m_0 \Gamma_m k_B T \delta_0(\tau). \quad (A2) $$

The presence in the above equation of the Dirac distribution simply means that there is no finite correlation time characterizing the bath (the associated spectrum is white). This obviously poses a mathematical problem for our definitions: the variance of this noise which defines the width of the Gaussian probability distribution is infinite (since it is related to the integral of the fluctuation spectrum). It shall not impact the final result of the modeling, which is cut off at high frequencies by the mechanical relaxation rate. One should therefore clarify that $\delta F$ fluctuations are Gaussian for any bandwidth $\Delta\omega$ cut in the white noise spectrum, around any frequency $\omega_0$.

Let us now transpose the dynamics into the rotating frame associated to the mode (at frequency $\Omega_m$),

$$ \delta F(t) = F_X(t) \cos(\Omega_m t) + F_Y(t) \sin(\Omega_m t), \quad (A3) $$

$$ x(t) = X(t) \cos(\Omega_m t) + Y(t) \sin(\Omega_m t), \quad (A4) $$

having introduced the two quadratures of force and motion. Equation (A1) can be rewritten, in matrix form,

$$ \begin{pmatrix} 1 & +\frac{1}{2Q} \\ -\frac{1}{2Q} & 1 \end{pmatrix} \begin{pmatrix} -\dot{X} \\ +\dot{Y} \end{pmatrix} = -\frac{\Gamma_m}{2} \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix} \begin{pmatrix} -X \\ +Y \end{pmatrix} + \frac{\Omega_m}{2k_0} \begin{pmatrix} F_Y \\ F_X \end{pmatrix}, \quad (A5) $$

with $Q = \Omega_m/\Gamma_m$ the quality factor, having neglected the slow components $\ddot{X}, \ddot{Y}$ (rotating wave approximation, valid for $Q \gg 1$). Similarly, we write

$$ \dot{x}(t) = \Omega_m [-X(t) \sin(\Omega_m t) + Y(t) \cos(\Omega_m t)], \quad (A6) $$

for the velocity, neglecting $\dot{X}$ and $\dot{Y}$. From the definitions of kinetic energy $E_c = m_0 \dot{x}^2/2$ and potential energy $E_p = k_0 x^2/2$, we simply obtain for the total energy $E_m = E_c + E_p$

$$ E_m(t) = k_0 \frac{X(t)^2 + Y(t)^2}{2}. \quad (A7) $$

Let us take the limit $Q \to +\infty$ in Eq. (A5); the $X$ and $Y$ equations then separate. Multiplying the first one by $k_0 X$, and the second one by $k_0 Y$, we write

$$ \begin{aligned} k_0 X \dot{X} &= -\frac{\Gamma_m}{2} k_0 X^2 - \frac{\Omega_m}{2} X F_Y, \\ k_0 Y \dot{Y} &= -\frac{\Gamma_m}{2} k_0 Y^2 + \frac{\Omega_m}{2} Y F_X, \end{aligned} \quad (A8) $$

which after adding up leads to the result

$$ \frac{dE_m}{dt} = -\Gamma_m E_m + \frac{\Omega_m}{2}[Y F_X - X F_Y]. \quad (A9) $$

This equation can be recast into the form of Eq. (24) by introducing the energy difference $\Delta E(t) = E_m(t) - \langle E_m \rangle$ and the bath stochastic energy flow $R_m(t)$:

$$ R_m(t) = \frac{\Omega_m}{2} [Y(t)F_X(t) - X(t)F_Y(t)] - \Gamma_m \langle E_m \rangle, \quad (A10) $$
obtained here in the high-Q limit; a more generic discussion can be found in Ref. [46]. The mean energy can be inferred from the equipartition result (see below): $\langle E_m \rangle = k_B T$. Equation (A9) is finally easily solved in frequency space as

$$S_{\Delta E}(\omega) = \frac{S_{R_m}(\omega)}{\Gamma_m^2 + \omega^2}, \quad (A11)$$

with $S_{R_m}$ the spectrum associated to $R_m$. We should now construct the statistical properties of this variable, from the initial properties of $\delta F$.

To do so, Eq. (A10) is rewritten as

$$R_m(t) = R_X(t) + R_Y(t) - \Gamma_m \langle E_m \rangle, \quad (A12)$$

$$R_X(t) = \frac{\Omega_m}{2}(\chi * F_X)(t) F_X(t), \quad (A13)$$

$$R_Y(t) = \frac{\Omega_m}{2}(\chi * F_Y)(t) F_Y(t), \quad (A14)$$

where we introduced the mechanical susceptibility $\chi(t)$ Fourier transform (in the rotating frame), and $*$ designates the convolution product $[(f * g)(t) = \int_{-\infty}^{+\infty} f(t-x)g(x) dx]$. We have

$$\chi(t) = \frac{\Omega_m}{2k_0} e^{-\frac{\Gamma_m}{2}t} \Theta(t), \quad (A15)$$

$$\langle F_X(t)F_X(t - \tau) \rangle = 4m_0\Gamma_m k_B T \delta_0(\tau), \quad (A16)$$

$$\langle F_Y(t)F_Y(t - \tau) \rangle = \langle F_X(t)F_X(t - \tau) \rangle, \quad (A17)$$

$$\langle F_X(t)F_Y(t - \tau) \rangle = 0, \quad (A18)$$

with $\Theta(t)$ the Heaviside function $[0 \text{ for } t < 0$, and $1$ for $t > 0$; for the time being, we only require $\Theta(0)$ to be finite]. Equations (A17) and (A18) simply state that the phase of the random force is irrelevant, which would not be the case in the presence of a squeezed noise. Equation (A15) solves Eq. (A5) in the limit $1/Q \approx 0$, while Eq. (A16) is deduced from Eqs. (A2) and (A3) (note the extra factor of 2 in the rotating frame noise amplitude). $F_X$ and $F_Y$ have by construction the same probability distribution as $\delta F$ (namely, Gaussian). The mean values verify

$$\langle R_X \rangle = \langle R_Y \rangle = \Gamma_m k_B T \int_{-\infty}^{+\infty} e^{-\frac{\Gamma_m}{2}(t-x)} \Theta(t-x) \delta_0(x-t) dx \\ = \Gamma_m k_B T \Theta(0), \quad (A19)$$

which introduces the value $\Theta(0)$ which has not been defined yet. In order to impose $\langle R_m \rangle = 0$, we have to take $\Theta(0) = 1/2$. Consider now the correlation functions of the type

$$\langle R_A(t)R_B(t') \rangle = \left(\frac{\Omega_m^2}{4k_0}\right)^2 \times \int_{-\infty}^{+\infty} \int_{-\infty}^{+\infty} e^{-\frac{\Gamma_m}{2}(t-x)} \\ \Theta(t-x) e^{-\frac{\Gamma_m}{2}(t'-x')} \Theta(t'-x') \\ \times \langle F_A(x)F_A(t)F_B(x')F_B(t') \rangle dx dx', \quad (A20)$$

with $A, B = X, Y$ in all possible combinations. The second order force correlator can be decomposed using Wick's

theorem:

$$\langle F_A(x)F_A(t)F_B(x')F_B(t') \rangle = +\langle F_A(x)F_A(t) \rangle \langle F_B(x')F_B(t') \rangle \\ + \langle F_A(x)F_B(x') \rangle \langle F_A(t)F_B(t') \rangle \\ + \langle F_A(x)F_B(t') \rangle \langle F_B(x')F_A(t) \rangle. \quad (A21)$$

Reinjecting this result into Eq. (A20), one obtains

$$\langle R_X(t)R_X(t') \rangle = \langle R_Y(t)R_Y(t') \rangle \\ = (\Gamma_m k_B T)^2 \Theta(0)^2 + \Gamma_m (k_B T)^2 [\delta_0(t - t') \\ + \Gamma_m \Theta(t - t') \Theta(t' - t)], \quad (A22)$$

$$\langle R_X(t)R_Y(t') \rangle = \langle R_Y(t)R_X(t') \rangle = (\Gamma_m k_B T)^2 \Theta(0)^2. \quad (A23)$$

In Eq. (A22), the product $\Theta(t - t')\Theta(t' - t)$ can obviously be dropped when compared to the $\delta_0(t - t')$ term. Finally, making use of all of these findings, we deduce from Eq. (A12)

$$\langle R_m(t)R_m(t') \rangle = 2\Gamma_m(k_B T)^2 \delta_0(t - t'), \quad (A24)$$

which leads to the corresponding (white) spectrum $S_{R_m}(\omega) = 2\Gamma_m(k_B T)^2$; Eq. (A11) therefore reproduces Eq. (29), as it should.

To conclude this Appendix, let us focus on the distributions of the random variables. Since the forces $F_X$ and $F_Y$ are Gaussian distributed, their corresponding motions $X$ and $Y$ are also Gaussian (linear response). Their PDFs are written

$$p(X) = \frac{1}{\sqrt{2\pi \sigma_X^2}} e^{-\frac{x^2}{2\sigma_X^2}}, \quad (A25)$$

$$p(Y) = \frac{1}{\sqrt{2\pi \sigma_Y^2}} e^{-\frac{y^2}{2\sigma_Y^2}}, \quad (A26)$$

with $\sigma_X^2$ and $\sigma_Y^2$ the corresponding variances, which are defined as

$$\sigma_X^2 = \sigma_Y^2 = \frac{1}{2\pi} \int_{-\infty}^{+\infty} |\chi(\omega)|^2 S_{F_X}(\omega) d\omega = \frac{k_B T}{k_0}, \quad (A27)$$

in which we introduced

$$\chi(\omega) = \frac{\Omega_m}{2k_0} \frac{1}{\Gamma_m/2 + i\omega}, \quad (A28)$$

$$S_{F_X}(\omega) = 4m_0\Gamma_m k_B T, \quad (A29)$$

the mechanical susceptibility $\chi(\omega)$ in frequency space and the force noise spectrum $S_{F_X}(\omega)$. The variances are finite, as they should be.

Equation (A7) tells us that energy is the sum of two uncorrelated squared Gaussian variables: this is known as a $(\chi_2)^2$ law. It results in an exponential distribution function:

$$p(E_m) = \frac{1}{\sigma_E} e^{-\frac{E_m}{\sigma_E}} \Theta(E_m), \quad (A30)$$

which verifies $\langle E_m \rangle = \sigma_E$ and $\langle E_m^2 \rangle = 2\sigma_E^2$, leading to an energy variance of $\sigma_E^2$. Since $\langle E_m \rangle = k_0(\sigma_X^2 + \sigma_Y^2)/2$, one infers
<figure><img src="image_10.png" /><figcaption>FIG. 8. Power dependence of the effective damping (peak half-height width) $\Gamma_{\text{eff}}$, at different temperatures. The lines are fits leading to $\Gamma_m$ ($P_{\text{in}} = 0$ values) and coupling $G$ (slope); see text. *Inset*: Temperature dependence of mechanical damping $\Gamma_m$; the line is a linear fit [25].</figcaption></figure>

immediately that $\sigma_E = k_B T$ which matches the well-known equipartition result.

Equation (A30) is nothing but the classical version of the Boltzmann distribution. The final message is then that an exponential energy distribution is equivalent to a Gaussian motion distribution; with the magnitudes of the associated white spectra related to $T$, the bath characteristic temperature.

## APPENDIX B: EFFECTIVE DAMPING $\Gamma_{\text{eff}}$

The microwave setup calibration is discussed in Appendix C below; on the other hand, the optomechanics coupling $G$ requires a specific measurement that we present here. It is based on the exploitation of the mean sideband peak characteristics (averaging together all the data measured during the period $\Delta T$, using the blue detuned pumping). This is done for all acquisition speeds $\delta t$, except the fastest one (red curves in Fig. 4) which distorts the measured line, an effect commented upon in Appendix G. The Lorentzian fit enables one to extract area (leading to mean energy, Fig. 3), peak position, and half-height width $\Gamma_{\text{eff}}$. The latter can be fit to Eqs. (1) and (2) as a function of injected power $P_{\text{in}}$, for each temperature $T$. This is shown in Fig. 8.

The slopes of these lines define the optomechanical coupling $G/(2\pi) \approx 1.8 \cdot 10^{13}$ Hz/m, and the $P_{\text{in}} \to 0$ extrapolation gives us the mechanical damping rate $\Gamma_m$. It is found to be linear in temperature in this range (see inset in Fig. 8), in accordance with Ref. [25].

The same calibration can be done in red detuned pumping, with a change of sign in the slope. Note that these slopes are independent of $T$, as they should be since the coupling is a pure geometrical effect. The scatter in Fig. 8 is genuine, and comes from the fluctuations of mechanical parameters (see Sec. VI). Finally, the measurement of $\Gamma_{\text{eff}}$ enables one to compute the optomechanical gain $G$ for any $T$ and $P_{\text{in}}$ (and any of the two red or blue detuned schemes).

## APPENDIX C: TWPA CHARACTERIZATION

Prior to the run, the microwave setup is characterized. The gains and losses are measured carefully using a vector network analyzer. For the injection lines, room temperature noise is suppressed by 50-dB attenuation affixed at different stages of the cryostat (zigzag blocks in Fig. 1).

On the detection side, the cryogenic HEMT provides a gain of 40 dB while the room temperature one has 30 dB. The microwave noise background around 6 GHz, referenced to the input of the cryo-HEMT, is then about 2.5 K. This is what is obtained from measurements performed before the mixer, with a spectrum analyzer; when using the mixing technique, we obtain twice this background, as explained in the main part of the paper. Besides the calibration of our (passive) microwave elements, the HEMT noise figure has been verified, and validated, in a run using the cold/hot load technique: comparing the measured noise generated by a 50-Ω termination located on the 3-K plate to the one of a similar load bolted onto the mixing chamber stage (at 10 mK). A microwave switch mounted on the same mixing chamber plate enabled us to connect one or the other loads while keeping the cryostat cold.

At the lowest temperatures, the TWPA on/off gain around 6 GHz is 10 dB, and the insertion loss of “TWPA + directional coupler + circulator” is about 4 dB. This is not optimal since we work on the higher side of the bandpass of the amplifier, which was designed for a 15-dB gain at slightly lower frequencies. Optimizing the center frequency and the insertion losses, we infer that we could potentially win about 5–10 dB [26]. The gain degrades quickly as we increase the temperature of the parametric amplifier above 300 mK; besides, the intrinsic noise of the TWPA itself increases dramatically above this limit (because of free quasiparticles thermally excited in the aluminum layer; see Fig. 9). Typically, it can still be used at 400 mK with a marginal gain, and essentially ruins the measurement chain above.

The measured noise of the full setup with the TWPA amplifier “on” is shown in Fig. 9 (the three amplifier gains are taken into account, the noise being referenced to the input of the TWPA). As we increase the mixing chamber plate temperature $T$, we see a thermal increase in this noise which is due to all the components present on the mixing chamber plate: our 50-Ω load, the microwave components (circulators and couplers), and the TWPA itself (see Fig. 1). As such, modeling properly the $T$-dependent losses becomes quite involved, and is far outside of the scope of our paper. We therefore simply linearly fit the data in Fig. 9 in order to extract a reasonable estimate of our $T \to 0$ background noise. The extrapolation slightly underestimates the real limiting noise figure, since for $k_B T/(\hbar\omega_c) \ll 1$ it should flatten out below typically 100 mK. We therefore retain a conservative value of about 0.8 K ($\pm$0.1 K) at $T = 0$ K, corresponding thus to about three photons, which matches expectations [26].

In the inset of Fig. 9, we show the input power dependence of the measured noise: the background level increases with $P_{\text{in}}$, and similarly for the two pumping schemes (red and blue detuned). By measuring the cavity spectrum, we observe that this is actually due to an out-of-equilibrium photon population that increases with the pump power [25]. Classically, we
<figure><img src="image_11.png" /></figure>

FIG. 9. Microwave noise amplitude around 6 GHz (referenced to the input of the TWPA amplifier), for low injection powers $P_{in}$. The line is a linear guide showing a temperature-dependent contribution when the mixing chamber plate temperature is regulated, and the $T \to 0$ limit of $0.8 \pm 0.1$ K (three photons). Note the fast increase in noise above 0.35 K (see text). *Inset*: Noise dependence to power $P_{in}$ measured at 200 mK, for the two different pumping schemes; the line is a linear guide to the eye, demonstrating the impact of input noise (see text).

model it with the voltage noise amplitudes $V_{P,\pm 1}$ appearing in Eq. (10), the sign $\pm$ referring to red or blue detuned. By using different microwave sources, or adding a notch filter, we can demonstrate that this noise is directly related to the quality of the pump signal. The actual power dependence seems to depend on the setup; it looks reasonably linear in Fig. 9 which is different from the findings of Ref. [25]. Besides, within our resolution it seems that the so-called technical heating of Fig. 3, and the $1/f$ contribution to the spectrum of Fig. 4, *do not* depend on this pump noise.

## APPENDIX D: “SPIKES” INSTABILITY

In Ref. [25], an instability in the dynamics of the beam was reported for temperatures lower than about 150 mK. It is visible as large amplitude peaks appearing in the sideband spectrum, which were nicknamed “spikes.” The origin of this phenomenon is still unknown, and we presume that it should also impact our measurements, the extent of which is the point of the present Appendix.

In Fig. 10 we plot two slow time tracks taken at slightly different temperatures: 100 and 200 mK. Strikingly, we see that the amplitude of fluctuations is *much larger* on the colder data set, which is a clear signature of spikes.

Let us comment in more detail upon the data. When computing the mean energy $\langle E_n \rangle$, we see that the value extrapolated at zero injected power is actually *smaller* than expected (see inset of Fig. 3). On the other hand, the relative importance of technical heating grows as we cool down; if this is not taken into account properly (which is particularly difficult at low temperatures without TWPA), the inferred mode energy

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

FIG. 10. Slow time tracks normalized to mean for two temperatures, at microwave power around 7 nW (see legend). The colder set exhibits much *larger* excursions than the hotter one (see discussion in text).

at a finite power $P_{in}$ would then be much *larger* than the thermodynamic value.

Interestingly, while the $1/f$ contribution to energy fluctuations increases as we cool down (Fig. 6), the flat part still seems to reproduce very well the thermodynamic value, even at 100 mK (see Fig. 5). It is thus very tempting to suggest that spikes are inherently linked to the $1/f$ fluctuations, whatever might be the microscopic mechanism behind this. How to understand a smaller mean energy compared to the thermodynamic temperature $T$ remains also an open question.

## APPENDIX E: “RED” AND “BLUE” PUMPING SCHEMES COMPARISON

In order to validate our analysis, we compare red and blue detuned schemes. For the mean mechanical energy, this is done in Fig. 3. We see that, indeed, correcting for the optomechanical gain $\mathcal{G}$ produces the same result. In this Appendix, we shall concentrate on fluctuations.

In Fig. 11 we plot the normalized phonon spectrum level $\sqrt{\frac{\Gamma_m S_0}{\mathcal{G}^2}} / (\frac{\kappa_{\text{ext}}}{\kappa_{\text{tot}}} |\Gamma_{\text{opt}}|)$, recalculated from the photon flux spectrum fit (right axis, 200-mK data). We present both blue detuned data (which can be found also in Fig. 5) and red detuned data. As for the mean mechanical energy, the agreement between the two pumping methods is very good, the difference being obviously that it is not possible to follow the fastest tracks with a red detuned scheme, since we need the optomechanical gain $\mathcal{G} > 1$ to do so.

In Fig. 11, we also show the corresponding $1/f$ components (left axis). However, we kept them in units of photons because transforming the graph into phonons following the same procedure as for the flat spectrum $S_0$ does not produce a much better plot [the fit of $A_f$ in Eq. (33) is not that good]. Determining whether the origin of this effect is in the mechanics or the optics remains thus an open question. However, we clearly see that red and blue detuned pumping data follow
<figure><img src="image_13.png" /></figure>

FIG. 11. Spectral contributions as a function of pump power $P_{\text{in}}$ measured at 200 mK with the two pumping schemes: red detuned (red circled dots) or blue detuned (blue circled). Left axis (log. scale): $1/f$ contribution $A_f$, in photons$^2$/s$^2$ (the dashed line is the exponential fit of Fig. 6). Right axis (lin. scale): Recalculated normalized phonon spectrum level (same parameter as in Fig. 5, no units). The horizontal line is obtained from the fit value of Fig. 4.

the same trend as a function of $P_{\text{in}}$: the $1/f$ term grows very quickly with increasing power (see exponential fit).

## APPENDIX F: FASTEST TRACKS FIT

The most difficult measurements are obviously the ones realized at the fastest acquisition speed. For the slower tracks, the peaks are sufficiently well defined that the Lorentz fit error is small compared to the reproducibility; this is not true anymore for the fastest tracks, where both an increased error (discussed in this Appendix) and a bias (presented in Appendix G) exist. When opening the bandwidth, the background noise increases as well, and we can resolve the sideband peak only at the highest powers $P_{\text{in}}$, with the largest gains $G$. Besides, with an acquisition bandwidth $1/\delta t$ larger than the peak width $\Gamma_{\text{eff}}$, we lose information on the shape of the sideband: the imprint of the motion is visible as only one to three points higher than the background (see right top inset in Fig. 4).

The fitting procedure is then as follows: we first average together all the data sets taken over the period $\Delta T = 1000 \delta t$. This produces a sideband peak which looks reasonably Lorentzian, with a width essentially given by the acquisition bandwidth. Then, in the fitting routine that infers the sideband parameters of each data file, we fix the Lorentz peak width to its mean value, and constrain the peak position to be around 0 Hz within only a few frequency steps $1/\delta t$. As such, the fit peak frequency position distribution looks like a (centered) truncated Gaussian, and our main fitting parameter is the height of the peak, or equivalently its area. Obviously, acquiring data faster than $\Gamma_{\text{eff}}$ should impact quantitatively the extracted area value: this point is explicitly discussed in Appendix G. The optomechanics gain $G$ is then computed from the known power (and temperature) dependence of $\Gamma_{\text{eff}}$;

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

FIG. 12. The 200 mK fastest track power spectral densities for peak (red data) and background (gray data). The full line corresponds to Eq. (33) without the $1/f$ term, but with a constant noise contribution attributed to fitting noise (which has been subtracted in Fig. 4). The dashed line marks the level of background noise (see text for details). Inset: Corresponding probability distribution functions; the line is the exponential curve (see text).

note that at these acquisition speeds, $1/f$ noise in damping $\Gamma_m$ is irrelevant.

The great capability of this technique is that we can easily separate what is genuinely characteristic of the sideband, from what is simply due to the noise background. This is illustrated in Fig. 12, where we compare the power spectral density calculated from the previous fitting, and the one obtained when constraining the fit position of the Lorentz peak far from the central value of 0 Hz. In the latter case, the obtained spectrum is white, and more than one order of magnitude smaller than what is obtained when fitting on the sideband (main graph). The probability distribution function is centered on zero, and clearly distinct from the exponential tail obtained with the sideband peak data (inset). Note the slight negativity which comes from a nonconstrained fit that also captures the background noise when no signal is to be seen; this has been truncated in Fig. 4 for clarity.

Finally, only the sideband processed data show the cutoff at $\Gamma_{\text{eff}}$ in the spectrum, but the computed power spectral density does not fall to zero above this value (see fit in Fig. 12). This is presumably due to the fit error, which is distinct from the background noise (which is subtracted in Fig. 4 for clarity); see Appendix G for details of fast tracks fit corrections. As a matter of fact, the final scatter in Fig. 5 corresponds to our ability to fit the flat region of the FFTs of the type of Fig. 4 (main graph), obtained for different temperatures. This finite error bar (which can be understood as our capability of defining $k_B$) corresponds here to 100 spectra averaged together, producing the typical scatter seen in Figs. 4 and 12.

## APPENDIX G: FAST ACQUISITION RATE CORRECTIONS

At the fastest acquisition rate, we saw in the preceding Appendix that the shape of the sideband peak is altered: it
<figure><img src="image_15.png" /></figure>

FIG. 13. Mean area of the sideband peak as a function of acquisition bandwidth $1/\delta t$ (only the two fastest settings), in normalized units ($A_0$ being the value corresponding to the slowest track; data are taken at different powers and different temperatures). The line is a simple fit function: $f(x) = 1 + x/4$. Inset: Corresponding normalized flat spectrum level as a function of bandwidth ($S_0$ being the reference value for slower acquisition speeds); the fit is $1/f(x)^2$.

becomes broader, its width being defined by the sampling rate. This is not the only impact of the fast tracking. Comparing the mean area obtained at different speeds $\delta t$, we also find out that it is overestimated. In contrast, when stitching the fastest spectrum to the others, we realize that we underestimate fluctuations. This is summarized in Fig. 13, in a universal plot with normalized axes.

The x axis corresponds to the sampling bandwidth normalized to $\Gamma_{\text{eff}}$. The y axis is the mean area normalized to its value obtained at slow acquisitions (main), or the fit plateau in the fluctuation spectrum normalized to the value extracted with slow tracking (inset). These can be fit by very simple empirical laws (see caption of Fig. 13).

In practice, with our settings only the fastest tracks ($\delta t = 22$ ms, Fig. 4) need a rescaling. Note that it does not impact the fit of the plateau $S_0$ in power spectral density, which is very clearly defined by slower acquisition rates. It is only needed for display purposes, when plotting the full-range data from $1/f$ to cutoff $\Gamma_{\text{eff}}$. The fitting routine itself also impacts Fig. 13, and its inherent bias is contained therein within our empirical dependencies. For a more profound analysis of fit biases, please see Ref. [47].

## APPENDIX H: SLIDING AVERAGING

A first attempt to produce energy fluctuations spectra had been made in Ref. [21], with measurements performed down

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

FIG. 14. Power spectral density obtained from data acquired at 200 mK (slowest rate, same set as Fig. 4). We present the original spectrum, and the one obtained when applying a “sliding average” on the measurement. Lines are fits (see text).

to the quantum regime. While the idea was clearly defined, the resolution of the experiment was very far from the requirements needed to produce the results we describe here. To obtain fittable data, the authors had to process a “sliding averaging” on the acquired measurements (averaging together $\#n$ neighboring files, while shifting this window through the whole set of data), and the extracted spectrum characteristics did not present the expected thermodynamic behavior. The nature of these slow fluctuations remained thus an open question in this paper.

In the present Appendix we study the effect of sliding averaging on our own data. In Fig. 14, we plot the power spectral density obtained with our raw slow data measured at 200 mK, together with the one obtained when processing a sliding averaging (with a $\#n = 10$ file averaging window). We clearly see that the averaging acts as a low pass filter, which transforms the initial $1/f$ component of our data into a $1/f^2$ (see fit lines); it completely suppresses the thermodynamic plateau $S_0$.

We therefore conclude that sliding averaging essentially preserves only the $1/f$ component of energy fluctuations. There is thus no particular information in the shape of the spectra obtained this way, since they are characteristic only of the filtering method. Especially, what looked like a very low frequency cutoff with a plateau is nothing but an artifact of filtering + FFT method. However, the $\sigma_E \propto \sqrt{T}$ law observed in Ref. [21] contains genuine information, which is characteristic of the (unknown) mechanism causing these slow fluctuations.

<footnote>[1] L. Landau and E. Lifshitz, Statistical Physics (Pergamon Press, London, 1959), Vol. 5.</footnote>

<footnote>[2] C. Jarzynski, Nonequilibrium Equality for Free Energy Differences, Phys. Rev. Lett. 78, 2690 (1997).</footnote>
[3] C. Jarzynski, Equalities and inequalities: Irreversibility and the second law of thermodynamics at the nanoscale, Annu. Rev. Condens. Matter Phys. **2**, 329 (2011).

[4] J. V. Koski, V. F. Maisi, J. P. Pekola, and D. V. Averin, Experimental realization of a Szilard engine with a single electron, Proc. Natl. Acad. Sci. USA **111**, 13786 (2014).

[5] J. P. Pekola, Towards quantum thermodynamics in electronic circuits, Nat. Phys. **11**, 118 (2015).

[6] M. Campisi, P. Hänggi, and P. Talkner, Colloquium: Quantum fluctuation relations: Foundations and applications, Rev. Mod. Phys. **83**, 771 (2011).

[7] A. Bérut, A. Arakelyan, A. Petrosyan, S. Ciliberto, R. Dillenschneider, and E. Lutz, Experimental verification of Landauer's principle linking, Nature (London) **483**, 187 (2012).

[8] S. Dago, J. Pereda, N. Barros, S. Ciliberto, and L. Bellon, Information and Thermodynamics: Fast and Precise Approach to Landauer's Bound in an Underdamped Micromechanical Oscillator, Phys. Rev. Lett. **126**, 170601 (2021).

[9] M. Ribezzi-Crivellari and F. Ritort, Large work extraction and the Landauer limit in a continuous Maxwell demon, Nat. Phys. **15**, 660 (2019).

[10] A. N. Cleland, *Foundations of Nanomechanics: From Solid-State Theory to Device Applications* (Springer, New York, 2013).

[11] D. V. Averin and J. P. Pekola, Violation of the Fluctuation-Dissipation Theorem in Time-Dependent Mesoscopic Heat Transport, Phys. Rev. Lett. **104**, 220601 (2010).

[12] L. G. C. Rego and G. Kirczenow, Quantized Thermal Conductance of Dielectric Quantum Wires, Phys. Rev. Lett. **81**, 232 (1998).

[13] M. P. Blencowe and V. Vitelli, Universal quantum limits on single-channel information, entropy, and heat flow, Phys. Rev. A **62**, 052104 (2000).

[14] J. P. Pekola and B. Karimi, Colloquium: Quantum heat transport in condensed matter systems, Rev. Mod. Phys. **93**, 041001 (2021).

[15] E. Collin, Mesoscopic quantum thermo-mechanics: A new frontier of experimental physics, AVS Quantum Science **4**, 020501 (2022).

[16] K. Schwab, E. Henriksen, J. Worlock, and M. Roukes, Measurement of the quantum of thermal conductance, Nature (London) **404**, 974 (2000).

[17] A. Tavakoli, K. Lulla, T. Crozes, N. Mingo, E. Collin, and O. Bourgeois, Heat conduction measurements in ballistic 1d phonon waveguides indicate breakdown of the thermal conductance quantization, Nat. Commun. **9**, 4287 (2018).

[18] A. Bassi, K. Lochan, S. Satin, T. P. Singh, and H. Ulbricht, Models of wave-function collapse, underlying theories, and experimental tests, Rev. Mod. Phys. **85**, 471 (2013).

[19] L. Diósi, Testing Spontaneous Wave-Function Collapse Models on Classical Mechanical Oscillators, Phys. Rev. Lett. **114**, 050403 (2015).

[20] A. Vinante and H. Ulbricht, Gravity-related collapse of the wave function and spontaneous heating: Revisiting the experimental bounds, AVS Quantum Science **3**, 045602 (2021).

[21] D. Cattiaux, I. Golokolenov, S. Kumar, M. Sillanpää, L. Mercier de Lépinay, R. R. Gazizulin, X. Zhou, A. D. Armour, O. Bourgeois, A. Fefferman, and E. Collin, A macroscopic object passively cooled into its quantum ground state of motion beyond single-mode cooling, Nat. Commun. **12**, 6182 (2021).

[22] G. E. Uhlenbeck and L. S. Ornstein, On the theory of the Brownian motion, Phys. Rev. **36**, 823 (1930).

[23] W. A. Phillips, Tunneling states in amorphous solids, J. Low Temp. Phys. **7**, 351 (1972).

[24] P. W. Anderson, B. I. Halperin, and C. M. Varma, Anomalous low-temperature thermal properties of glasses and spin glasses, Philos. Mag. **25**, 1 (1972).

[25] X. Zhou, D. Cattiaux, R. R. Gazizulin, A. Luck, O. Maillet, T. Crozes, J.-F. Motte, O. Bourgeois, A. Fefferman, and E. Collin, On-Chip Thermometry for Microwave Optomechanics Implemented in a Nuclear Demagnetization Cryostat, Phys. Rev. Appl. **12**, 044066 (2019).

[26] L. Planat, A. Ranadive, R. Dassonneville, J. Puertas Martínez, S. Léger, C. Naud, O. Buisson, W. Hasch-Guichard, D. M. Basko, and N. Roch, Photonic-Crystal Josephson Traveling-Wave Parametric Amplifier, Phys. Rev. X **10**, 021021 (2020).

[27] M. Aspelmeyer, T. J. Kippenberg, and F. Marquardt, Cavity optomechanics, Rev. Mod. Phys. **86**, 1391 (2014).

[28] X. Zhou, D. Cattiaux, D. Theron, and E. Collin, Electric circuit model of microwave optomechanics, J. Appl. Phys. **129**, 114502 (2021).

[29] I. Golokolenov, D. Cattiaux, S. Kumar, M. Sillanpää, L. M. de Lépinay, A. Fefferman, and E. Collin, Microwave single-tone optomechanics in the classical regime, New J. Phys. **23**, 053008 (2021).

[30] R. Kubo, The fluctuation-dissipation theorem, Rep. Prog. Phys. **29**, 255 (1966).

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

[32] P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, S. J. van der Walt, M. Brett, J. Wilson, K. J. Millman, N. Mayorov, A. R. J. Nelson, E. Jones, R. Kern, E. Larson, C. J. Carey et al., SciPy 1.0: Fundamental algorithms for scientific computing in Python, Nature Methods **17**, 261 (2020).

[33] M. Newville, T. Stensitzki, D. B. Allen, and A. Ingargiola, LMFIT: Non-linear least-square minimization and curve-fitting for Python (2014), doi:10.5281/zenodo.11813.

[34] T. C. P. Chui, D. R. Swanson, M. J. Adriaans, J. A. Nissen, and J. A. Lipa, Temperature fluctuations in the canonical ensemble, Phys. Rev. Lett. **69**, 3005 (1992).

[35] P. K. Day, I. Hahn, T. C. P. Chui, A. W. Harter, D. Rowe, and J. A. Lipa, The fluctuation-imposed limit for temperature measurement, J. Low Temp. Phys. **107**, 359 (1997).

[36] C. R. Harris, K. J. Millman, S. J. van der Walt, R. Gommers, P. Virtanen, D. Cournapeau, E. Wieser, J. Taylor, S. Berg, N. J. Smith, R. Kern, M. Picus, S. Hoyer, M. H. van Kerkwijk, M. Brett, A. Haldane, J. F. del Río, M. Wiebe, P. Peterson, P. Gérard-Marchant et al., Array programming with NumPy, Nature (London) **585**, 357 (2020).

[37] M. Sansa, E. Sage, E. C. Bullard, M. Gély, T. Alava, E. Colinet, A. K. Naik, L. G. Villanueva, L. Duraffourg, M. L. Roukes, G. Jourdan, and S. Hentz, Frequency fluctuations in silicon nanoresonators, Nat. Nanotechnol. **11**, 552 (2016).

[38] O. Maillet, X. Zhou, R. R. Gazizulin, R. Ilic, J. M. Parpia, O. Bourgeois, A. D. Fefferman, and E. Collin, Measuring frequency fluctuations in nonlinear nanomechanical resonators, ACS Nano **12**, 5753 (2018).
[39] B. D. Hauer, P. H. Kim, C. Doolin, F. Souris, and J. P. Davis, Two-level system damping in a quasi-one-dimensional optomechanical resonator, Phys. Rev. B 98, 214303 (2018).

[40] O. Maillet, D. Cattiaux, X. Zhou, R. R. Gazizulin, O. Bourgeois, A. D. Fefferman, and E. Collin, Nanomechanical damping via electron-assisted relaxation of two-level systems, arXiv:2009.03804v1.

[41] J. Burnett, L. Faoro, I. Wisby, V. L. Gurtovoi, A. V. Chernykh, G. M. Mikhailov, V. A. Tulin, R. Shaikhaidarov, V. Antonov, P. J. Meeson, A. Y. Tzalenchuk, and T. Lindström, Evidence for interacting two-level systems from the 1/f noise of a superconducting resonator, Nat. Commun. 5, 4119 (2014).

[42] A. Tavakoli, K. J. Lulla, T. Puurtinen, I. Maasilta, E. Collin, L. Saminadayar, and O. Bourgeois, Specific heat of thin phonon cavities at low temperature: Very high values revealed by zep-tojoule calorimetry, Phys. Rev. B 105, 224313 (2022).

[43] C. K. Stahle, D. McCammon, and K. D. Irwin, Quantum calorimetry, Phys. Today 52(8), 32 (1999).

[44] B. Karimi, F. Brange, P. Samuelsson, and J. P. Pekola, Reaching the ultimate energy resolution of a quantum detector, Nat. Commun. 11, 367 (2020).

[45] E. T. Mannila, P. Samuelsson, S. Simbierowicz, J. T. Peltonen, V. Vesterinen, L. Grönberg, J. Hassel, V. F. Maisi, and J. P. Pekola, A superconductor free of quasiparticles for seconds, Nat. Phys. 18, 145 (2022).

[46] Y. Mishin and J. Hickman, Energy spectrum of a Langevin oscillator, Phys. Rev. E 94, 062151 (2016).

[47] S. F. Nørrelykke and H. Flyvbjerg, Power spectrum analysis with least-squares fitting: Amplitude bias and its elimination, with application to optical tweezers and atomic force microscope cantilevers, Rev. Sci. Instrum. 81, 075103 (2010).