# Supplemental Information:
Measuring High-Order Phonon Correlations in an Optomechanical Resonator

Y. S. S. Patil<EMAIL_ADDRESS>Department of Physics, Yale University,
New Haven, Connecticut 06520, USA  J. Yu  Department of Physics, Yale
University, New Haven, Connecticut 06520, USA  S. Frazier  Department of
Physics, Yale University, New Haven, Connecticut 06520, USA  Y. Wang
Department of Applied Physics, Yale University, New Haven, Connecticut 06520,
USA  K. Johnson  Department of Physics, Yale University, New Haven,
Connecticut 06520, USA  J. Fox  Department of Physics, Yale University, New
Haven, Connecticut 06520, USA  J. Reichel  Laboratoire Kastler Brossel, ENS-
Université PSL, CNRS, Sorbonne Université, Collège de France 24 rue Lhomond,
75005 Paris, France  J. G. E. Harris<EMAIL_ADDRESS>Department of
Physics, Yale University, New Haven, Connecticut 06520, USA  Department of
Applied Physics, Yale University, New Haven, Connecticut 06520, USA  Yale
Quantum Institute, Yale University, New Haven, Connecticut 06520, USA

##

###  .1 Optical Setup

This section describes the optical setup used for the experiments presented in
the main text.

Before describing the setup itself, we first describe the timing sequence used
in these measurements. During an experiment the setup is alternated between
two configurations, each lasting 100 ms. The time sequence is shown in Fig. S1
(the configurations are shown in Fig. S2 and discussed in the following
paragraph). The first of these is the “lock” configuration, during which the
optical path is configured to lock the tunable laser (TL) and filter cavities
to the optomechanical cavity (OMC). The second is the “drive” configuration,
during which the various control voltages used to tune the TL and filter
cavities are held constant. Also, at the start of the drive period, the
optical path is configured so that either a blue- or red-detuned laser drives
the OMC. Photon counting data is collected during the middle 90 ms of the
drive period.

Figure S2 shows the two configurations, with the lock configuration’s optical
path in green and the drive configuration’s optical path in orange. The lock
configuration begins with the TL entering an IQ modulator (IQM) operating in
the single-sideband suppressed-carrier mode. The IQM serves as a frequency
shifter to lock the TL to the OMC. The output light from the IQM is amplified
by an Erbium-doped fiber amplifer (EDFA) and filtered by a broadband
($\Delta\lambda\approx$ 0.3 nm) filter to suppress the EDFA’s amplified
spontaneous emission noise. The beam then passes through an electro-optic
modulator (EOM) which produces a pair of sidebands used for standard Pound-
Drever-Hall (PDH) locking. The beam is then split into three parts, which are
sent to the OMC and the two filter cavities (FC1, FC2). The TL is thus locked
to the OMC via the IQM control voltage, and the filter cavities are locked to
the TL via their piezo and thermal tuning.

Figure S1:  Lock and drive sequence. The time sequence for switching between
the “lock” and the “drive” configurations. Photon counts are acquired only
during the central 90 ms of the “drive” period.

Figure S2:  Optical schematic. This figure shows the two optical
configurations (lock and drive) used in the experiment. The tunable laser (TL)
is locked to the optomechanical cavity during the lock configuration. The blue
(BL) and red (RL) drive lasers are always locked to the TL with a frequency
offset $\sim\pm\omega_{\text{ac}}$. A total of six $1\times 2$ MEMS switches
are used to switch between the lock configuration (green optical path) and the
drive configuration (orange optical path) every 100 ms, such that the optical
components are either locked to their desired frequency, or free running (as
indicated in Fig. S1). PC: polarization controller. VOA: variable optical
attenutator. EDFA: erbium-doped fiber amplifier. EOM: electro-optic modulator.
ILP: in-line polarizer. ISO: isolator. FC: filter cavity.

The blue-/red-detuned drive lasers (BL/RL) are always locked to the TL. This
is accomplished by combining each drive laser with the TL (shown in the upper
left corner of Fig. S2) and monitoring the frequency of the resulting beat
note to ensure that the two drive lasers are detuned from the OMC by $\sim$
$\pm\omega_{\text{ac}}$. Each drive laser passes through a filter cavity
($\kappa/2\pi\approx$ 30 MHz) and is then combined via a 50/50 fiber coupler
to drive the OMC. The part of the optical path that is common to the two drive
beams is shown in orange. A pair of shutters just before the 50/50 coupler are
used to determine whether the BL or the RL drives the OMC. The shutters are
also used to set the desired pulse sequence for pulsed measurements.

Scattered sideband photons that are resonant with the OMC exit the cavity and
pass through filter cavities FC1 and FC2, finally arriving at a pair of
superconducting nanowire single photon detectors (SNSPDs). The variable
optical attenuators (VOAs) are used to adjust the optical power, and the
polarization controllers (PCs) are used to align the polarization of light
entering polarization-sensitive components such as the FCs, the EOM and the
SNSPDs.

###  .2 Characterizating the Spectrum of Photon Counts

Figure 1(d) of the main text shows a typical photon count rate spectrum. Here
we discuss the characterization of such spectra. Fig. S3 shows the count rate
spectra for various values of the drive power $P_{\text{in}}$. Each spectrum
shows a narrow peak (arising from the Stokes- and anti Stokes-scattered
photons) centered at $\omega_{\text{ac}}/2\pi\approx$ 315.3 MHz, a broader
peak centered around 322.3 MHz, and a frequency-independent background.

Figure S3:  Photon count spectrum. Count rates as a function of the detuning
and power of the drive laser. The narrow peak near $315$ MHz corresponds to
sideband photons produced by the creation and annihilation of phonons in the
superfluid-filled cavity. The broad beak near $322$ MHz is due to thermal
fluctuations in the room-temperature optical fibers. The fits are described in
the text.  Figure S4:  Power dependence of the photon count spectrum. (a)
$\Gamma_{\mathrm{bkg}}$ and $\Gamma_{\mathrm{G}}$ as a function of drive
power, and the corresponding linear fits. For $\Gamma_{\text{bkg}}$, the
y-intercepts give background counts arising from SNSPD dark counts and stray
light from the environment. The slopes give the leakage of drive photons
around the filter cavities for the red-detuned and blue-detuned drive each.
$\Gamma_{\mathrm{G}}$ also increases linearly with drive power. (b)
$\Gamma_{\mathrm{res}}$ as a function of drive power. Shown as open squares
are values of $\Gamma_{\mathrm{res}}$ extracted from a 50-frequency
measurement that is fit to Eq. 3. $\Gamma_{\mathrm{res}}$ extracted from a
five-frequency measurement that is fit to Eq. 4 are shown as open circles. The
difference between the two methods is negligible.

The sideband peak is expected to be proportional to the transfer function of
the two cascaded filter cavities FC1 and FC2, given by

$\displaystyle
f_{\text{filter}}(\Delta)=\frac{1}{\bigg{[}1+4\big{(}\frac{|\Delta|-\omega_{\text{ac}}}{\kappa_{\mathrm{FC}1}}\big{)}^{2}\bigg{]}\bigg{[}1+4\big{(}\frac{|\Delta|-\omega_{\text{ac}}}{\kappa_{\mathrm{FC}1}}\big{)}^{2}\bigg{]}}$
(1)

where $\kappa_{\mathrm{FC}1}/2\pi$ = 1.21(5) MHz and
$\kappa_{\mathrm{FC}2}/2\pi$ = 1.71(2) MHz are the filter linewidths. We
identify the broad peak around 322 MHz as the signature of guided-acoustic
wave Brillouin scattering (GAWBS) in the fibers (See Sec. .3 for more
details), which can be characterized by a Lorentzian centered at
$\omega_{\text{G}}$:

$\displaystyle
f_{\text{G}}(\Delta)=\frac{1}{1+4\big{(}\frac{|\Delta|-\omega_{\text{G}}}{\kappa_{\text{G}}}\big{)}^{2}}$
(2)

Thus, for each value of $P_{\text{in}}$ and for each detuning range (i.e.,
$\Delta\sim+\omega_{\text{ac}}$ and $\Delta\sim-\omega_{\text{ac}}$ ), the
count rate spectrum is fit to the form

$\displaystyle\Gamma(\Delta)=\Gamma_{\text{bkg}}+f_{\text{filter}}(\Delta)\Gamma_{\text{res}}+f_{\text{G}}(\Delta)\Gamma_{\text{G}}$
(3)

where $\Gamma_{\text{bkg}}$, $\Gamma_{\text{res}}$, $\Gamma_{\text{G}}$,
$\omega_{\text{ac}}$, $\omega_{\text{G}}$, and $\kappa_{\text{G}}$ are the fit
parameters.

Fig. S4 shows the fit-extracted $\Gamma_{\text{bkg}}$, $\Gamma_{\text{G}}$,
and $\Gamma_{\text{res}}$ as a function of $P_{\text{in}}$. The frequency-
independent background $\Gamma_{\text{bkg}}$ is found to have three sources:
1) the SNSPDs’ darks counts, 2) stray light leaking into the fibers, and 3)
unfiltered drive photons leaking around the two filter cavities in the
detection chain (this last contribution is polarization-dependent, and so
differs between the red-detuned and blue-detuned drives).
$\Gamma_{\text{bkg}}(P_{\text{in}})$ is fit to the form
$\Gamma_{\text{bkg,0}}$ \+ $\Gamma_{\text{bkg,1}}P_{\text{in}}$, which gives
$\Gamma_{\text{bkg,0}}$ = 12.4 $\pm$ 0.8 s-1, and
$\Gamma_{\text{bkg,1}}=9.1\pm 0.2$ s-1µW-1 for the red-detuned drive, and
$\Gamma_{\text{bkg,1}}=12.0\pm 0.4$ s-1µW-1 for the blue-detuned drive. By
acquiring the count rate while blocking the input port directly in front of
the SNSPDs, we find that $\Gamma_{\text{bkg,0}}$ consists of 7 $\pm$ 1 s-1
from dark counts and 5 $\pm$ 1 s-1 from stray light leaking into the
detection-chain fibers. $\Gamma_{\text{bkg,1}}$ quantifies the amount of drive
photons that leak around the filter cavities.

$\Gamma_{\text{G}}(P_{\text{in}})$ is fit to the form
$\Gamma_{\text{G,1}}P_{\text{in}}$. This gives $\Gamma_{\text{G,1}}$ = 28.5
$\pm$ 0.3 s-1µW-1.

A complete characterization of $\Gamma_{\text{res}}$ as a function of drive
power is given in Sec. .4.

The spectra in Fig. S3 (and those used to produce Fig. S4) are measured at
$50$ values of $\Delta/2\pi$. While this is helpful in characterizing the
device, the focus of the main paper is on the Stokes- and anti-Stokes-
scattered photons. To measure the rates of these specific photons, we found
that it was adequate to record the spectrum at just five frequencies close to
$\omega_{\text{ac}}$ (specifically, $\Delta/2\pi=$ 310, 312, 314.9, 315.4,
315.9 MHz), and then fit the results to the form

$\displaystyle\Gamma(\omega)=\Gamma_{\text{bkg}}+f_{\text{filter}}(\omega)\Gamma_{\text{res}}$
(4)

To illustrate the validity of this approach, Fig. S4(b) shows
$\Gamma_{\text{res}}(P_{\text{in}})$ acquired in two ways: by measuring the
count rate at $50$ values of $\Delta/2\pi$ ranging from $\pm$ 310 MHz to $\pm$
330 MHz (as in Fig. S3) and fitting to Eq. 3; and by measuring the count rate
at just the five values of $\Delta/2\pi$ listed above and fitting to Eq. 4.
The difference is negligible, showing the reliability of the five-point
spectrum and the minimal influence of the GAWBS peak (which is not included in
Eq. 4) on this analysis. The data presented in Sec. .4 were acquired using the
“five frequency” approach.

Figure S5:  Wide frequency scan showing several GAWBS peaks. Six individual
Lorentzians are fit (in addition to the tail of $f_{\mathrm{filter}}$), with
the vertical dashed lines marking the measured (best-fit) center frequencies
for the fiber’s transverse acoustic modes.  Figure S6:  The resonant
frequencies of GAWBS. The fitted peak frequencies are shown in solid black
squares and the calculated frequencies of $R_{0m}$ GAWBS resonant modes for a
$125\text{\,}\mathrm{\SIUnitSymbolMicro m}$ single-mode fused silica fiber are
shown in open red circles. The discrepancies are shown on the bottom graph in
solid green circles.

###  .3 Measuring the Wider Photon Count Spectrum and Characterizing GAWBS

This section gives a more detailed characterization of the GAWBS signature
mentioned in Sec. .2. GAWBS is the well-studied interaction between light in a
fiber and the fiber’s thermally populated transverse acoustic modes [1]. In
practice, this interaction produces phase noise in the light with a spectrum
determined by the fiber’s acoustic modes.

Figure S5 shows the photon count rate as the detuning $\Delta/2\pi$ is varied
from 150 to 450 MHz. The overall descending background in Fig. S5 corresponds
to the tail of the filter cavity’s transfer function. Due to the large steps
in the detuning, the sideband peak of the acoustic cavity mode cannot be seen
here, but several other peaks are evident. These peaks are found to be
centered close to the frequencies $f_{m}$ that are expected for the fiber’s
transverse acoustic modes, $R_{0m}$, which vibrate in the radial direction.
These frequencies (labeled as vertical lines in Fig. S5) correspond to the
solutions $y_{m}$ of

$\displaystyle(1-\alpha^{2})J_{0}(y)=\alpha^{2}J_{2}(y)$ (5)

where $J_{x}$ is the $x$-order Bessel function of the first kind, and
$\alpha=0.624$ is the ratio between transverse and longitudinal phonon
velocities for fibers made of fused silica. The resonance frequencies are
given by $f_{m}=\frac{V_{d}}{2\pi a}y_{m}$, where $V_{\text{d}}=5996$ m/s is
the longitudinal phonon velocity and the fiber radius $a$ = 62.5 µm [1, 2].

The data in Fig. S5 is fit to the sum of six Lorentzians (whose center
frequencies, linewidths, and magnitudes are the fit parameters), and the tail
of $f_{\text{filter}}(\omega)$. The best-fit values for the center frequencies
and the theoretical $f_{m}$ calculated from Eq. 5 are shown in Fig. S6.

###  .4 Model and Fits for Drive Power Scans

This section describes measurements of the sideband-photon counting rates as a
function of the mixing chamber temperature $T_{\text{MC}}$ and the laser power
incident on the device $P_{\text{in}}$. These measurements provide information
about the acoustic mode’s temperature and the useful range of $P_{\text{in}}$.
They also provide additional calibration of the device’s optomechanical
coupling rate $g_{0}$.

Figure S7(a) shows $R_{\text{AS}}$ and $R_{\text{S}}$ measured as a function
of $T_{\text{MC}}$ at an incident drive power $P_{\text{in}}=250$ nW. The
values of $R_{\text{AS}}$ and $R_{\text{S}}$ are extracted from fitting five-
point spectra as described in Sec. .2. Both $R_{\text{AS}}$ and $R_{\text{S}}$
increase with $T_{\text{MC}}$ as expected, and the persistent difference
between the two is consistent with the quantum sideband asymmetry (QSA). The
solid lines in Fig. S7(a) show a one-parameter fit to the form
$a(e^{\hbar\omega_{\text{ac}}/k_{B}T_{\text{MC}}}-1)^{-1}$ (for the red data)
and $1+a(e^{\hbar\omega_{\text{ac}}/k_{B}T_{\text{MC}}}-1)^{-1}$ (for the blue
data), where the fit parameter
$a=\gamma_{\text{ac}}C\eta_{\kappa}\eta_{\text{det}}/P_{\mathrm{in}}$ is the
magnitude of the power-normalized QSA. The fit only uses data with
$T_{\text{MC}}>50$ mK, as the calibration of the $\text{RuO}_{2}$ thermometer
used to determine $T_{\text{MC}}$ is uncertain at lower temperatures.

The agreement between data and theory in Fig. S7(a) indicates that the mean
phonon occupancy $n_{\text{ac}}$ is primarily determined by $T_{\text{MC}}$.
However, at higher drive powers, both optical backaction and absorption-
induced heating (in the fibers and mirrors) can alter $n_{\text{ac}}$. To
characterize the role of backaction and heating in these devices, Fig. S7b
shows $R_{\text{AS(S)}}$ as a function of $P_{\text{in}}$. (The data at
$P_{\text{in}}=0$ are measured using pulsed laser excitations and modelling
(and then subtracting) the transient heating.) For sufficiently low
$P_{\text{in}}$ ($\lesssim$ $300\text{\,}\mathrm{nW}$), the values of
$R_{\text{AS(S)}}$ are consistent with $T_{\text{MC}}\approx 20$ mK, and show
a splitting that is dominated by the QSA. As $P_{\text{in}}$ is increased,
both $R_{\text{AS}}$ and $R_{\text{S}}$ increase. However, $R_{\text{S}}$
increases more rapidly with $P_{\text{in}}$ than $R_{\text{AS}}$. This is
consistent with the combined effects of heating from optical absorption (which
increases $n_{\text{ac}}$ independent of $\Delta$) and optomechanical
backaction (which increases (decreases) $n_{\text{ac}}$ when $\Delta$ is
positive (negative)).

Figure S7:  Temperature- and power-dependence of the sidebands.  Normalized
anti-Stokes (red) and Stokes (blue) photon detection rates. (a) Measurements
as a function of $T_{\text{MC}}$ with $P_{\text{in}}=250$ nW. Solid lines: fit
to the form predicted for thermal equilibrium. (b) Measurements as a function
of $P_{\text{in}}$ with $T_{\text{MC}}\approx 20\text{mK}$. Solid lines: a fit
including optomechanical backaction (colored regions) and heating due to
optical absorption (dashed lines). Right-hand axis: the mean phonon occupancy
$n_{\text{ac}}$ ($n_{\text{ac}}+1$) inferred from the detection rates.

To analyze these effects quantitatively, we apply the standard theory of
optomechanical backaction, and model the heating by assuming that the device
is subject to a heatload proportional to $P_{\text{in}}$ and is thermally
linked to the MC with a thermal conductance $\sigma=bT^{k+1}$. Before
describing this approach in detail (see the following paragraphs), we first
summarize the result. The solid lines in Fig. S7(b) show the best fit while
taking $g_{0},T_{\text{MC}},\eta_{\text{det}},b$ and $k$ as fit parameters.
The dashed lines show the change in $n_{\text{ac}}$ attributable to
absorption-induced heating. The best-fit value $k=1.05(3)$ suggests a
predominantly metallic thermal conductance between the device and MC ($k=1$
for metals). The fit also gives a single photon coupling rate
$g_{0}/2\pi=4.58(2)$ kHz and a net detection efficiency of
$\eta_{\text{det}}=0.18(1)$. The best-fit value $T_{\text{MC}}=19(1)$ mK
corresponds to a mean thermal phonon occupancy of $n_{\text{ac}}=0.83(2)$.

In more detail, we start with the power-normalized sideband count rates as
described by standard quantum optomechanics theory:

$\displaystyle\overline{R}_{\text{AS}}$ $\displaystyle=an_{\text{ac}}$ (6)
$\displaystyle\overline{R}_{\text{S}}$ $\displaystyle=a(n_{\text{ac}}+1)$ (7)

where, as defined above,
$a=\gamma_{\text{ac}}C\eta_{\kappa}\eta_{\text{det}}/P_{\text{in}}$, and
$n_{\text{ac}}$ is determined by the joint effects from the thermal bath
($n_{\text{th}}$) and optomechanical backaction [3],

$\displaystyle
n_{\text{ac}}=\frac{\gamma_{\text{opt}}\left(\frac{\kappa_{\text{c}}}{4\omega_{\text{ac}}}\right)^{2}+\gamma_{\text{ac}}n_{\text{th}}}{\gamma_{\text{opt}}+\gamma_{\text{ac}}}$
(8)

where $\gamma_{\text{opt}}\left(=\gamma_{\text{ac}}C=4\cdot
g_{0}^{2}/\kappa_{\text{c}}\cdot n_{\text{c}}\right)$ and $\gamma_{\text{ac}}$
are the optomechanical damping rate and bare mechanical damping rate,
respectively, and $n_{\text{th}}$ accounts for the effect of optical-
absorption-induced heating in the fibers and mirrors.

To model the effect of optical heating, we assume a heat load proportional to
$P_{\text{in}}$, and that the conductivity of the thermal link between the
fiber and the MC is a power law in temperature with exponent $k$, such that
the fiber’s temperature is given by

$\displaystyle
T_{\text{fib}}=(T_{\text{MC}}^{k+1}+\beta^{k+1}P_{\text{in}})^{1/(k+1)}$ (9)

where $\beta$ characterizes the amount of heat generated by a given
$P_{\text{in}}$. In addition, we assume a ballistic transport of heat
(phonons) from the Helium inside the OMC to the MC bath via the Helium channel
connecting them. It follows that

$\displaystyle
n_{\text{th}}=\frac{n_{\text{fib}}\gamma_{\text{ac0}}+n_{\text{MC}}\gamma_{\text{ball}}}{\gamma_{\text{ac0}}+\gamma_{\text{ball}}}$
(10)

where $n_{\text{fib}}$ =
1/($e^{{\hbar\omega_{\text{ac}}}/\left({k_{\text{B}}T_{\text{fib}}}\right)}-1)$
and $n_{\text{MC}}$ =
1/($e^{{\hbar\omega_{\text{ac}}}/\left({k_{\text{B}}T_{\text{MC}}}\right)}-1)$.
The bare mechanical damping rate $\gamma_{\text{ac}}$ = $\gamma_{\text{ac0}}$
\+ $\gamma_{\text{ball}}$ includes: i) the intrinsic mechanical damping rate
$\gamma_{\text{ac0}}$/2$\pi$ = 3.2$\times$103 Hz at zero temperature due to
acoustic loss into the fibers, and ii) the three phonon loss
$\gamma_{\text{ball}}$/2$\pi$ = 2.7$\times$106 $T_{\text{MC}}^{4}$ Hz/K4 at
the base MC temperature.

Figure S8:  Several scans of the drive power. Scattering rates as a function
of drive power, taken over several days. The points are data, and the solid
lines are a fit in which all the fit parameters are global, but the detection
efficiency is allowed to vary between the data sets. For this plot, each data
set is normalized by its fitted detection efficiency.

We have verified the accuracy of this model by using it to analyze drive power
sweep data taken on several different days over the course of ten months. Fig.
S8 shows seven different power sweeps. The solid lines are the best fit of the
data to the model described by Eq. 7 \- 10. The fit uses a single
$T_{\text{MC}}$, $\beta$, $k$, and optomechanical coupling $g_{\text{0}}$ for
all of the sweeps but allows $\eta_{\text{det}}$ to vary between the sweeps,
as the overall detection efficiency is observed to differ from day to day.

The data and the fit are normalized to an ideal detection efficiency (i.e.,
$\eta_{\text{det}}$ = 1) for display in Fig. S8. The fit gives an initial
device temperature of 24.4 $\pm$ 0.4 mK, equivalent to 1.61 $\pm$ 0.03
phonons. The thermal conductivity follows a power law in temperature with an
exponent $k=1.09\pm 0.03$. The fitting parameters are summarized in Table. 1.

$T_{\text{MC}}$ |  24.4 $\pm$ 0.4 mK
---|---
$\beta$ |  0.54 $\pm$ 0.03 K/$\text{W}^{1/(k+1)}$
$k$ |  1.09 $\pm$ 0.03
Table 1:  Fit parameters for thermal transport from fiber to bath The best-fit
in Fig. S8 returns the fit parameters that describe the fiber temperature
according to Eq. 9. $T_{\mathrm{MC}}$ characterizes the initial device
temperature, $k$ the exponent of the power law that describes the thermal
conductivity, and $\beta$ the fraction of heat generated by a given $P_{in}$.

###  .5 Filtering of Counts Detected by the SNSPDs

Figure S9: Afterpulses in the SNSPDs. The second-order correlation function of
a coherent laser shows a sharp feature in the first 50 ns due to the
afterpulsing effect. The flat region before 20 ns is the dead time of the
SNSPDs. For delays $>30$ ns, the correlation decays to 1, as expected for a
coherent state.

Not every count registered by the SNSPDs corresponds to the detection of a
photon. Here we cosider three types of false-counts: (i) afterpulses, (ii)
rapid bursts of counts in a short duration (10 -
$100\text{\,}\mathrm{\SIUnitSymbolMicro s}$), and (iii) SNSPD dark counts.
This section describes the filtering protocols used to prevent false-counts of
types (i) and (ii) from compromising the coherence measurements. False counts
of type (iii) (i.e., dark counts) cannot be identified or filtered, given that
they occur at random times and are uniformly distributed in time; however,
their effect on measured coherences can be accounted for, as is detailed in
Sec. .6.

(i) Afterpulses It is well-known that reflection of the voltage pulses
generated by SNSPDs can cause false counts [4, 5]. Following a detection event
at time $t=0$, an afterpulse can occur at a delay $\tau_{\text{ap}}$ set by
the time it takes for the voltage reflection to reach the SNSPD.

The signature of afterpulses is illustrated in Fig. S9. It shows the second-
order coherence $g^{(2)}(\tau)$ of a power-stabilized laser, measured with a
bin-size of $2$ ns. The peak at delay $\tau_{\text{ap}}\approx 24$ ns is
caused by the afterpulses. For $\tau>50$ ns, $g^{(2)}(\tau)\approx 1$, as
expected for laser light.

To reduce the impact of these afterpulses in the experiments described here,
any count that is registered within $50$ ns of a preceding count (on a given
SNSPD) is tagged as an afterpulse and discarded. This also removes any true
counts that occur within $50$ ns of each other; however the error this
introduces in coherence measurements is to effectively shorten the bin at the
smallest delay by $50$ ns. This can easily be compensated for, but since the
shortest delays used in this work are $\tau\geq$
$1\text{\,}\mathrm{\SIUnitSymbolMicro s}$, we need not perform such a
compensation. For the same reason, we also need not compensate for the
recovery time of the SNSPD, which is specified to be $50$ ns.

(ii) Rapid Bursts There occur rare instances ($<0.008\%$ of all 90 ms DAq
configuration intervals, c.f. Sec. .1) in which the SNSPDs register a periodic
train of counts with an abnormally high rate ($\gtrsim 10^{3}\times$ expected
count rate) for a short duration ($10-100$ µs). While the underlying cause of
these spurious counts is unknown, such events can be identified as outliers
and discarded. This is done by using a statistical model which tags events
having an extremely low probability of occurring, as described below.

Suppose the mean count rate during an experiment is measured to be $I$. The
average number of counts received in a time interval $\Delta t$ is thus
$\lambda=I\Delta t$. The probability of receiving $k$ photons in this interval
can be estimated using a statistical model as $P(k,\lambda)$. Thus, for a
total DAq time $T$, among the $N=T/\Delta t$ number of $\Delta t$\- intervals,
the mean number of intervals in which $k$ counts are recorded is expected to
be $N\times P(k,\lambda)$. To reject outliers, we set a threshold $\epsilon\ll
1$ and search for all $\Delta t$-intervals that receive $\geq k_{\text{thr}}$
counts, where $k_{\text{thr}}$ is the smallest value of $k$ for which $N\times
P(k_{\text{thr}},\lambda)<\epsilon$. For each $\Delta t$-interval that is
identified as an outlier in this way, we discard the entire $90$ ms DAq record
containing this interval. In this work, $N\sim 10^{8}-10^{13}$ and $\epsilon$
is set to $0.1$.

Figure S10: Distribution of counts. Number of intervals as a function of
counts. We binned data into 3 $\mathrm{\mu s}$ intervals. The beige bars are
the histogram of counts in the intervals. The red solid line is the estimated
distribution of counts in the short interval limit $\Delta
t\ll\tau_{\text{c}}$, where $\tau_{\mathrm{c}}$ is the coherence time. We used
a probability threshold $0.1$ (shown as the black dashed line) corresponding
to a count threshold $k_{\mathrm{thr}}$ (shown by the blue dotted line). In
this data, one interval with 5 counts is rejected by this protocol.

Figure S10 illustrates this protocol. It shows the experimentally measured
frequency with which $k-$counts were detected in $\Delta t=$
$3\text{\,}\mathrm{\SIUnitSymbolMicro s}$ intervals for $T\approx 4800$ s
($N\sim 1.6\times 10^{9}$) and $I\approx 2500$ counts/s ($\lambda\sim
1.51\times 10^{-3})$. The red solid line represents the statistical model
$P(k,\lambda)$ (see below). For $k\geq 4$, $N\times P(k,\lambda)<0.1$, so that
$k_{\text{thr}}=4$. As such, the $90$ ms DAq intervals that included any
$\Delta t=$ $3\text{\,}\mathrm{\SIUnitSymbolMicro s}$ intervals that
registered $\geq 4$ counts were discarded from the dataset. In this instance,
$1$ out of $53,421$ DAq intervals was discarded.

We iterated the filtering protocol illustrated above for $\Delta
t=3,10,30,100$ and $300\text{\,}\mathrm{\SIUnitSymbolMicro s}$. The best
apriori statistical guess for $P(k,\lambda)$ depends on the choice of $\Delta
t$, the measured count rate $I$, and the coherence timescale $\tau_{\text{c}}$
of the photon source. For $\Delta t\gg\tau_{\text{c}}$ and
$I^{-1}\ll\tau_{\text{c}}$, the photons detected within the interval $\Delta
t$ can be approximated as uncorrelated, and we model $P(k,\lambda)$ as a
Poissonian, $\text{Pois}(k,\lambda)=\lambda^{k}e^{-\lambda}/k!$ [6]. However,
if $\Delta t\lesssim\tau_{\text{c}}$ or $I^{-1}\gtrsim\tau_{\text{c}}$,
$P(k,\lambda)$ is best modelled by accounting for the coherence of the source.
Given that the acoustic mode generating the sideband photons is expected to be
in thermal equilibrium, $P(k,\lambda)$ is well modelled by
$P_{\text{th}}(k,\lambda)=\lambda^{k}/(1+\lambda)^{k+1}$ for $\Delta
t\ll\tau_{\text{c}}$. For $\Delta t\sim\tau_{\text{c}}$, $P(k,\lambda)$ can be
estimated to be between these two limiting cases of $\Delta
t\gg\tau_{\text{c}}$ and $\Delta t\ll\tau_{\text{c}}$. To avoid making an
apriori guess for $\tau_{\text{c}}$, we chose to model $P(k,\lambda)$ with
$P_{\text{th}}(k,\lambda)$ for all $\Delta t$. This achieves a more
conservative filtering of the data as it allows for more frequent occurrences
of a higher number of counts within any $\Delta t$-interval
($P_{\text{th}}(k,\lambda)>P(k,\lambda)$ for $k>2,\lambda<1$).

Applying this protocol to all of the data acquired used in this work, $133$
out of $\sim 1.86\times 10^{6}$ DAq intervals were tagged and discarded as
outliers.

###  .6 Subtraction of Coincidences due to Background Photons

This section describes the method used to correct the correlation functions
for the presence of background counts in the SNSPD data. As described in Sec.
.2, these background counts result from stray light, leakage through the
filter cavities, and the SNSPD’s dark counts.

We measure (not shown) these background counts to have a fixed mean arrival
rate over the duration of the experiments, and to be uncorrelated among
themselves, i.e. they are Poisson distributed. We also assume them to be
independent of the sideband photons. Therefore, the measured second-order
correlation function $g^{(2)}_{\mathrm{exp}}(\tau)$, including the cross-
correlations between the sideband photons and the background photons, is given
by

$\displaystyle g^{(2)}_{\text{exp}}(\tau)$
$\displaystyle=\frac{\langle[a^{\dagger}+\xi^{\dagger}](0)[a^{\dagger}+\xi^{\dagger}](\tau)[a+\xi](\tau)[a+\xi](0)\rangle}{\langle[a^{\dagger}+\xi^{\dagger}](0)[a+\xi](0)\rangle\langle[a^{\dagger}+\xi^{\dagger}](\tau)[a+\xi](\tau)\rangle}$
$\displaystyle=\frac{g^{(2)}(\tau)+2\epsilon+\epsilon^{2}}{(1+\epsilon)^{2}},$

where $a$ and $\xi$ correspond to the sideband and background photons,
respectively, and we have used $\langle
a^{\dagger}\xi\rangle=0=\langle\xi^{\dagger}a\rangle$, and
$\langle\xi^{\dagger}(0)\xi^{\dagger}(\tau)\xi(\tau)\xi(0)\rangle=\langle\xi^{\dagger}\xi\rangle^{2}$.
$\epsilon=\langle\xi^{\dagger}\xi\rangle/\langle a^{\dagger}a\rangle$ is the
ratio of the mean background and sideband count rates. For the drive powers
used in this work, $0.04\lesssim\epsilon\lesssim 0.2$.

The corrected second-order correlation function $g^{(2)}(\tau)$ is thus given
by

$\displaystyle
g^{(2)}(\tau)=g^{(2)}_{\mathrm{exp}}(\tau)+2(g^{(2)}_{\mathrm{exp}}(\tau)-1)\epsilon+(g^{(2)}_{\mathrm{exp}}(\tau)-1)\epsilon^{2}.$

A similar calculation gives the corrected third-order and fourth-order
correlation functions using

$\displaystyle
g^{(3)}_{\mathrm{exp}}(\tau_{1},\tau_{2})=\dfrac{g^{(3)}(\tau_{1},\tau_{2})+\epsilon\left(g^{(2)}(\tau_{1})+g^{(2)}(\tau_{2})+g^{(2)}(\tau_{1}+\tau_{2})\right)+3\epsilon^{2}+\epsilon^{3}}{(1+\epsilon)^{3}},$
$\displaystyle
g^{(4)}_{\mathrm{exp}}(\tau_{1},\tau_{2},\tau_{3})=\dfrac{g^{(4)}(\tau_{1},\tau_{2},\tau_{3})+\epsilon\sum_{g_{i}^{(3)}\in\mathbb{G}^{(3)}}g_{i}^{(3)}+\epsilon^{2}\sum_{g_{j}^{(2)}\in\mathbb{G}^{(2)}}g_{j}^{(2)}+4\epsilon^{3}+\epsilon^{4}}{(1+\epsilon)^{4}},$

where $\mathbb{G}^{(3)}$ and $\mathbb{G}^{(2)}$ are

$\displaystyle\mathbb{G}^{(3)}$
$\displaystyle=\set{g^{(3)}(\tau_{1},\tau_{2}),g^{(3)}(\tau_{1}+\tau_{2},\tau_{3}),g^{(3)}(\tau_{1},\tau_{2}+\tau_{3}),g^{(3)}(\tau_{2},\tau_{3})},$
$\displaystyle\mathbb{G}^{(2)}$
$\displaystyle=\set{g^{(2)}(\tau_{1}),g^{(2)}(\tau_{2}),g^{(2)}(\tau_{3}),g^{(2)}(\tau_{1}+\tau_{2}),g^{(2)}(\tau_{2}+\tau_{3}),g^{(2)}(\tau_{1}+\tau_{2}+\tau_{3})}.$

###  .7 Correspondence between Photon and Phonon Correlation Functions

In the main text, we introduced the normally ordered phonon coherence
$g_{\text{ac}}^{(n)}$ obtained from the anti-Stokes photon coherence and the
anti-normally ordered phonon coherence $h_{\text{ac}}^{(n)}$ obtained from the
Stokes photon coherence. In this section, we formally demonstrate this one-to-
one correspondence between coherences of the optical fields output by the
cavity and the coherences of the phonons. We also note that optical losses
along the detection chain do not compromise these coherence measurements [7].

In a frame rotating at the optical drive frequency
$\omega_{\text{D}}=\omega_{\text{c}}+\Delta$, the output optical field $o$ of
the optomechanical cavity has a spectrum $S_{o^{\dagger}o}$ that relates to
mechanical spectrum $S_{xx}$ as [8]

$S_{o^{\dagger}o}(\omega)=\kappa_{\text{in}}n_{\text{c}}2\pi\delta(\omega)+\kappa_{\text{in}}g_{0}^{2}n_{\text{c}}|\chi_{\text{c}}(-\omega)|^{2}S_{xx}(\omega),$
(11)

where the first term corresponds to the optical drive field and the second
term corresponds to the optical fields produced by the mechanical motion
through the optomechanical coupling. $\chi_{\text{c}}(\omega)$ is the bare
optical cavity susceptibility
$\chi_{\text{c}}(\omega)=(\kappa_{\text{c}}/2-i(\omega+\Delta))^{-1}$. The
field $p$ detected by the SNSPDs is the cavity output field $o$ passed through
the filter cavities with transmission $f_{\text{filter}}(\omega)$,

$p(t)=\int\limits_{-\infty}^{\infty}\frac{d\omega}{2\pi}e^{-i\omega
t}f_{\text{filter}}(\omega)o(\omega),$

so that the detected normally ordered optical correlations are given by the
Weiner-Khinchin relation

$\displaystyle\langle
p^{\dagger}(t+\tau)p(t)\rangle=\int\limits_{-\infty}^{\infty}\frac{d\omega}{2\pi}\leavevmode\nobreak\
e^{-i\omega\tau}|f_{\text{filter}}(-\omega)|^{2}S_{o^{\dagger}o}(\omega),$

$f_{\text{filter}}(\omega)$ is non-zero only around
$\omega=+\omega_{\text{ac}}$ for the red-detuned drive (and only around
$\omega=-\omega_{\text{ac}}$ for the blue-detuned drive), within the two
filter cavity linewidths
$\kappa_{\text{FC1}},\kappa_{\text{FC2}}\ll\omega_{\text{ac}}$. This ensures
that the contribution of the first term in Eq. 11 to this integral,
corresponding to the contribution of the drive field, is negligible. As such,

$\displaystyle\langle p^{\dagger}(t+\tau)p(t)\rangle=$ (12)
$\displaystyle\kappa_{\text{in}}g_{0}^{2}n_{\text{c}}\int\limits_{-\infty}^{\infty}\frac{d\omega}{2\pi}\leavevmode\nobreak\
e^{-i\omega\tau}|f_{\text{filter}}(-\omega)|^{2}|\chi_{\text{c}}(-\omega)|^{2}S_{xx}(\omega).$

The contribution of $S_{xx}(\omega)$ to these optical correlations can be
decomposed as

$\displaystyle
S_{xx}(\omega)=x_{\text{zpt}}^{2}[S_{b^{\dagger}b}(\omega)+S_{bb^{\dagger}}(\omega)],$

with $x=x_{\text{zpt}}(b+b^{\dagger})$, and $x_{\text{zpt}}$ being the zero-
point motion of the oscillator. $S_{b^{\dagger}b}(\omega)$ is non-zero only
around $\omega=-\omega_{\text{ac}}$, and $S_{bb^{\dagger}}(\omega)$ is non-
zero only around $\omega=+\omega_{\text{ac}}$, corresponding to the anti-
Stokes and Stokes scattering processes respectively, which occur within the
mechanical susceptibility $\chi_{\text{ac}}(\omega)$ of linewidth
$\bar{\gamma}_{\text{ac}}$. This is explicitly seen through the oscillator’s
Langevin equation [9, 10, 11, 12]

$\displaystyle\dot{b}$
$\displaystyle=-\left(\frac{\bar{\gamma}_{\text{ac}}}{2}-i\omega_{\text{ac}}\right)b+\sqrt{\bar{\gamma}_{\text{ac}}}b_{\text{in}},$

where $b_{\text{in}}$ is the thermal noise input with $\langle
b_{\text{in}}\rangle=0$. The correlations for $\tau\geq 0$ are given by
$\langle
b_{\text{in}}^{\dagger}(\tau)b_{\text{in}}(0)\rangle=n_{\text{ac}}\delta(\tau)$
and $\langle
b_{\text{in}}(\tau)b_{\text{in}}^{\dagger}(0)\rangle=(n_{\text{ac}}+1)\delta(\tau)$.
And $n_{\text{ac}}=\langle b(t)^{\dagger}b(t)\rangle$ is assumed to be time-
independent (stationary). The formal integral

$b(\tau)=\sqrt{\bar{\gamma}_{\text{ac}}}\int\limits_{-\infty}^{\infty}dt^{\prime}e^{-(\bar{\gamma}_{\text{ac}}/2-i\omega_{\text{ac}})(\tau-t^{\prime})}b_{\text{in}}(t^{\prime})$

yields the correlations and spectra as

$\langle
b^{\dagger}(\tau)b(0)\rangle=n_{\text{ac}}e^{-(\bar{\gamma}_{\text{ac}}/2+i\omega_{m})\tau},$
(13) $\displaystyle S_{b^{\dagger}b}(\omega)$
$\displaystyle=\int\limits_{-\infty}^{\infty}d\tau\leavevmode\nobreak\
e^{i\omega\tau}\langle b^{\dagger}(\tau)b(0)\rangle$
$\displaystyle=n_{\text{ac}}\bar{\gamma}_{\text{ac}}|\chi_{\text{ac}}(-\omega)|^{2}$
$\displaystyle=\frac{n_{\text{ac}}\bar{\gamma}_{\text{ac}}}{(\bar{\gamma}_{\text{ac}}/2)^{2}+(-\omega-\omega_{\text{ac}})^{2}}.$

And similarly

$\langle
b(\tau)b^{\dagger}(0)\rangle=(n_{\text{ac}}+1)e^{-(\bar{\gamma}_{\text{ac}}/2-i\omega_{m})\tau},$
(14) $\displaystyle S_{bb^{\dagger}}(\omega)$
$\displaystyle=\int\limits_{-\infty}^{\infty}d\tau\leavevmode\nobreak\
e^{i\omega\tau}\langle b(\tau)b^{\dagger}(0)\rangle$
$\displaystyle=(n_{\text{ac}}+1)\bar{\gamma}_{\text{ac}}|\chi_{\text{ac}}(\omega)|^{2}$
$\displaystyle=\frac{(n_{\text{ac}}+1)\bar{\gamma}_{\text{ac}}}{(\bar{\gamma}_{\text{ac}}/2)^{2}+(\omega-\omega_{\text{ac}})^{2}}.$

Thus, for the case of the red-detuned drive $\Delta=-\omega_{\text{ac}}$, only
the $S_{b^{\dagger}b}$ term contributes to the integral around
$\omega=-\omega_{\text{ac}}$ in Eq. 12, corresponding to the anti-Stokes
scattered photons, and yields

$\displaystyle\langle p^{\dagger}(t+\tau)p(t)\rangle=\langle
a_{\text{AS}}^{\dagger}(t+\tau)a_{\text{AS}}(t)\rangle$
$\displaystyle=\kappa_{\text{in}}g_{0}^{2}n_{\text{c}}|\chi_{\text{c}}(\omega_{\text{ac}})|^{2}x_{\text{zpt}}^{2}\int\limits_{-\infty}^{\infty}\frac{d\omega}{2\pi}\leavevmode\nobreak\
e^{-i\omega\tau}S_{b^{\dagger}b}(\omega)$
$\displaystyle=4\kappa_{\text{in}}g_{0}^{2}n_{\text{c}}x_{\text{zpt}}^{2}/\kappa_{\text{c}}^{2}\leavevmode\nobreak\
\langle b^{\dagger}(t+\tau)b(t)\rangle,$

where we have approximated $f_{\text{filter}}(-\omega)\approx 1$ and
$\chi_{\text{c}}(-\omega)\approx\chi_{\text{c}}(\omega_{\text{ac}})$ for
$\omega$ around $-\omega_{\text{ac}}$, as
$\bar{\gamma}_{\text{ac}}\ll\kappa_{\text{FC1}},\kappa_{\text{FC2}}\ll\kappa_{\text{c}}$.

For the case of the blue-detuned drive $\Delta=+\omega_{\text{ac}}$, only the
$S_{bb^{\dagger}}$ term contributes to the integral around
$\omega=+\omega_{\text{ac}}$ in Eq. 12, corresponding to the Stokes scattered
photons, and yields

$\displaystyle\langle p^{\dagger}(t+\tau)p(t)\rangle=\langle
a_{\text{S}}^{\dagger}(t+\tau)a_{\text{S}}(t)\rangle$
$\displaystyle=\kappa_{\text{in}}g_{0}^{2}n_{\text{c}}|\chi_{\text{c}}(-\omega_{\text{ac}})|^{2}x_{\text{zpt}}^{2}\int\limits_{-\infty}^{\infty}\frac{d\omega}{2\pi}\leavevmode\nobreak\
e^{-i\omega\tau}S_{bb^{\dagger}}(\omega)$
$\displaystyle=4\kappa_{\text{in}}g_{0}^{2}n_{\text{c}}x_{\text{zpt}}^{2}/\kappa_{\text{c}}^{2}\leavevmode\nobreak\
\langle b(t+\tau)b^{\dagger}(t)\rangle,$

where we have approximated $f_{\text{filter}}(-\omega)\approx 1$ and
$\chi_{\text{c}}(-\omega)\approx\chi_{\text{c}}(-\omega_{\text{ac}})$ for
$\omega$ around $+\omega_{\text{ac}}$, as
$\bar{\gamma}_{\text{ac}}\ll\kappa_{\text{FC1}},\kappa_{\text{FC2}}\ll\kappa_{\text{c}}$.

This demonstrates the correspondence between the detected optical correlations
and the appropriately ordered phonon correlations. Note that high-order
correlations can be evaluated in terms of such two-time correlations using
Wick’s theorem, for any Gaussian state, thereby implying the correspondence of
the detected optical coherences and the appropriately ordered phonon
coherences. Note that the proportionality constant
$(4\kappa_{\text{in}}g_{0}^{2}n_{\text{c}}x_{\text{zpt}}^{2}/\kappa_{\text{c}}^{2})$
is immaterial for the coherences, which are normalized ratios of such
correlations.

###  .8 Correlation Functions of a Thermal State

This section presents calculations of coherence functions for a thermal state
of an oscillator. These functions are used to fit the data in the main text.

We calculate the coherences for an oscillator’s thermal state using Wick’s
theorem. For any Gaussian state, the theorem implies that high-order moments
can be evaluated in terms of second moments. For example, the numerators in

$g^{(2)}_{\text{ac}}(\tau_{1})=\frac{\langle
b_{0}^{\dagger}b_{1}^{\dagger}b_{1}b_{0}\rangle}{\langle
b_{0}^{\dagger}b_{0}\rangle\langle
b_{1}^{\dagger}b_{1}\rangle},\leavevmode\nobreak\
h^{(2)}_{\text{ac}}(\tau_{1})=\frac{\langle
b_{0}b_{1}b_{1}^{\dagger}b_{0}^{\dagger}\rangle}{\langle
b_{0}b_{0}^{\dagger}\rangle\langle b_{1}b_{1}^{\dagger}\rangle}$

where $b_{0}=b(0)$ and $b_{1}=b(\tau_{1})$, evaluate as

$\displaystyle\langle b_{0}^{\dagger}b_{1}^{\dagger}b_{1}b_{0}\rangle$
$\displaystyle=\langle b_{0}^{\dagger}b_{0}\rangle\langle
b_{1}^{\dagger}b_{1}\rangle+\langle b_{0}^{\dagger}b_{1}\rangle\langle
b_{1}^{\dagger}b_{0}\rangle,$ $\displaystyle\langle
b_{0}b_{1}b_{1}^{\dagger}b_{0}^{\dagger}\rangle$ $\displaystyle=\langle
b_{0}b_{0}^{\dagger}\rangle\langle b_{1}b_{1}^{\dagger}\rangle+\langle
b_{0}b_{1}^{\dagger}\rangle\langle b_{1}b_{0}^{\dagger}\rangle.$

Figure S11:  Fourth-order phonon coherences.
$g^{(4)}_{\text{ac}}(\tau_{1},\tau_{2},\tau_{3})$ and
$h^{(4)}_{\text{ac}}(\tau_{1},\tau_{2},\tau_{3})$ measured for
$P_{\text{in}}\approx$ $5\text{\,}\mathrm{\SIUnitSymbolMicro W}$, with photon
arrival times binned in $10\text{\,}\mathrm{\SIUnitSymbolMicro s}$ bins. Shown
are three sets of representative 2D slices:
$g^{(4)}_{\text{ac}}(\tau_{1},0^{+},\tau_{3})$ and
$h^{(4)}_{\text{ac}}(\tau_{1},0^{+},\tau_{3})$,
$g^{(4)}_{\text{ac}}(\tau_{1},\infty,\tau_{3})$ and
$h^{(4)}_{\text{ac}}(\tau_{1},\infty,\tau_{3})$,
$g^{(4)}_{\text{ac}}(\tau_{1},\tau_{2},\infty)$ and
$h^{(4)}_{\text{ac}}(\tau_{1},\tau_{2},\infty)$, where $\tau=0^{+}$ represents
the bin with $$5\text{\,}\mathrm{\SIUnitSymbolMicro
s}$<\tau<$15\text{\,}\mathrm{\SIUnitSymbolMicro s}$$ and $\tau=\infty$
represents times $\tau>10\bar{\gamma}_{\text{ac}}^{-1}$. Fits are to the
entire 3D ($\tau_{1}$-, $\tau_{2}$-, $\tau_{3}$\- dependent) data. Fit
residuals are shown in black.

For the thermal state, which is a Gaussian state, these second moments are
given for $\tau\geq 0$ by Eq. 13, 14, i.e.

$\displaystyle\langle b^{\dagger}(\tau)b(0)\rangle$
$\displaystyle=n_{\text{ac}}e^{-(\bar{\gamma}_{\text{ac}}/2+i\omega_{\text{ac}})\tau},$
$\displaystyle\langle b(\tau)b^{\dagger}(0)\rangle$
$\displaystyle=(n_{\text{ac}}+1)e^{-(\bar{\gamma}_{\text{ac}}/2-i\omega_{\text{ac}})\tau}.$

The second moments for $\tau\leq 0$ are obtained by conjugation. Combined,
these readily give

$\displaystyle g^{(2)}_{\text{ac}}(\tau)$
$\displaystyle=1+e^{-\bar{\gamma}_{\text{ac}}\tau}=h^{(2)}_{\text{ac}}(\tau).$

To illustrate further, the numerator in

$g^{(3)}_{\text{ac}}(\tau_{1},\tau_{2})=\frac{\langle
b_{0}^{\dagger}b_{1}^{\dagger}b_{2}^{\dagger}b_{2}b_{1}b_{0}\rangle}{\langle
b_{0}^{\dagger}b_{0}\rangle\langle b_{1}^{\dagger}b_{1}\rangle\langle
b_{2}^{\dagger}b_{2}\rangle},$

where $b_{2}=b(\tau_{1}+\tau_{2})$, evaluates as

$\displaystyle\langle
b_{0}^{\dagger}b_{1}^{\dagger}b_{2}^{\dagger}b_{2}b_{1}b_{0}\rangle=\langle
b_{0}^{\dagger}b_{0}\rangle\langle b_{1}^{\dagger}b_{1}\rangle\langle
b_{2}^{\dagger}b_{2}\rangle$ $\displaystyle+\langle
b_{0}^{\dagger}b_{0}\rangle\langle b_{1}^{\dagger}b_{2}\rangle\langle
b_{2}^{\dagger}b_{1}\rangle+\langle b_{0}^{\dagger}b_{1}\rangle\langle
b_{1}^{\dagger}b_{0}\rangle\langle b_{2}^{\dagger}b_{2}\rangle$
$\displaystyle+\langle b_{0}^{\dagger}b_{1}\rangle\langle
b_{1}^{\dagger}b_{2}\rangle\langle b_{2}^{\dagger}b_{0}\rangle+\langle
b_{0}^{\dagger}b_{2}\rangle\langle b_{1}^{\dagger}b_{0}\rangle\langle
b_{2}^{\dagger}b_{1}\rangle$ $\displaystyle+\langle
b_{0}^{\dagger}b_{2}\rangle\langle b_{1}^{\dagger}b_{1}\rangle\langle
b_{2}^{\dagger}b_{0}\rangle.$

This, and evaluations of $h^{(3)}_{\text{ac}},g^{(4)}_{\text{ac}}$, and
$h^{(4)}_{\text{ac}}$ in a similar manner, gives

$\displaystyle
g^{(3)}_{\text{ac}}(\tau_{1},\tau_{2})=h^{(3)}_{\text{ac}}(\tau_{1},\tau_{2})$
$\displaystyle=1+e^{-\bar{\gamma}_{\text{ac}}\tau_{1}}+e^{-\bar{\gamma}_{\text{ac}}\tau_{2}}+3e^{-\bar{\gamma}_{\text{ac}}(\tau_{1}+\tau_{2})},$
$\displaystyle
g_{\text{ac}}^{(4)}(\tau_{1},\tau_{2},\tau_{3})=h_{\text{ac}}^{(4)}(\tau_{1},\tau_{2},\tau_{3})$
$\displaystyle=1+e^{-\bar{\gamma}_{\text{ac}}\tau_{1}}+e^{-\bar{\gamma}_{\text{ac}}\tau_{2}}+e^{-\bar{\gamma}_{\text{ac}}\tau_{3}}$
$\displaystyle+3e^{-\bar{\gamma}_{\text{ac}}(\tau_{1}+\tau_{2})}+3e^{-\bar{\gamma}_{\text{ac}}(\tau_{2}+\tau_{3})}+e^{-\bar{\gamma}_{\text{ac}}(\tau_{1}+\tau_{3})}$
$\displaystyle+9e^{-\bar{\gamma}_{\text{ac}}(\tau_{1}+\tau_{2}+\tau_{3})}+4e^{-\bar{\gamma}_{\text{ac}}(\tau_{1}+2\tau_{2}+\tau_{3})}.$
(15)

Limiting cases amenable to an intuitive interpretation are discussed below.

The fourth-order phonon coherence functions
$g^{(4)}_{\text{ac}}(\tau_{1},\tau_{2},\tau_{3})$ and
$h^{(4)}_{\text{ac}}(\tau_{1},\tau_{2},\tau_{3})$ are three-dimensional data
sets that can be visualized as two-dimensional slices. In Fig. 2(c) of the
main text we show one set of 2D slices
$g^{(4)}_{\text{ac}}(0^{+},\tau_{2},\tau_{3})$ and
$h^{(4)}_{\text{ac}}(0^{+},\tau_{2},\tau_{3})$, where $\tau=0^{+}$ represents
the bin with $$5\text{\,}\mathrm{\SIUnitSymbolMicro
s}$<\tau<$15\text{\,}\mathrm{\SIUnitSymbolMicro s}$$. In Fig. S11, we show
three more sets of 2D slices: $g^{(4)}_{\text{ac}}(\boldsymbol{\tau})$ and
$h^{(4)}_{\text{ac}}(\boldsymbol{\tau})$ measured at
$\boldsymbol{\tau}=(\tau_{1},0^{+},\tau_{3})$ (panel (a)),
$\boldsymbol{\tau}=(\tau_{1},\infty,\tau_{3})$ (panel (b)), and
$\boldsymbol{\tau}=(\tau_{1},\tau_{2},\infty)$ (panel (c)). Here,
$\tau=\infty$ represents times $\tau>10\bar{\gamma}_{\text{ac}}^{-1}$. Fits to
the entire 3D ($\tau_{1}$-, $\tau_{2}$-, $\tau_{3}$\- dependent) data are
shown as 2D surfaces, with fit residuals shown in black.

The 2D slice $g^{(4)}_{\text{ac}}(\tau_{1},\infty,\tau_{3})$ is equivalent to
the product of two one-dimensional $g^{(2)}_{\text{ac}}(\tau)$ functions. If
the third photon arrives with delay $\tau_{2}=\infty$, then the arrivals of
the third and fourth photons are uncorrelated with the arrival of the first
pair of photons. That is, $g^{(4)}_{\text{ac}}(\tau_{1},\infty,\tau_{3})$ is
proportional to the product of the probability of measuring a pair of photons
separated by a delay $\tau_{1}$ and the probability of measuring another pair
of photons separated by a delay $\tau_{3}$. This can be seen from the
expression for $g^{(4)}_{\text{ac}}(\tau_{1},\tau_{2},\tau_{3})$ given in Eq.
15, where setting $\tau_{2}=\infty$ results in

$\displaystyle g_{\mathrm{ac}}^{(4)}(\tau_{1},\infty,\tau_{3})$
$\displaystyle=1+e^{-\bar{\gamma}_{\text{ac}}\tau_{1}}+e^{-\bar{\gamma}_{\text{ac}}\tau_{3}}+e^{-\bar{\gamma}_{\text{ac}}(\tau_{1}+\tau_{3})}$
$\displaystyle=g_{\mathrm{ac}}^{(2)}(\tau_{1})g_{\mathrm{ac}}^{(2)}(\tau_{3}).$

The same equivalence holds for
$h_{\mathrm{ac}}^{(4)}(\tau_{1},\infty,\tau_{3})=h_{\mathrm{ac}}^{(2)}(\tau_{1})h_{\mathrm{ac}}^{(2)}(\tau_{3})$.

Similarly, the slice $g^{(4)}_{\text{ac}}(\tau_{1},\tau_{2},\infty)$ is
equivalent to the two-dimensional $g^{(3)}_{\text{ac}}(\tau_{1},\tau_{2})$. If
the fourth photon arrives with delay $\tau_{3}=\infty$, then its arrival is
uncorrelated with the arrivals of the first three photons. Therefore,
$g^{(4)}_{\text{ac}}(\tau_{1},\tau_{2},\infty)$ is proportional to the
probability of measuring a triplet of photons with a delay $\tau_{1}$ between
the first and second and a delay $\tau_{2}$ between the second and third. We
can also see this by setting $\tau_{3}=\infty$ in Eq. 15, which yields

$\displaystyle g_{\mathrm{ac}}^{(4)}(\tau_{1},\tau_{2},\infty)$
$\displaystyle=1+e^{-\bar{\gamma}_{\text{ac}}\tau_{1}}+e^{-\bar{\gamma}_{\text{ac}}\tau_{2}}+3e^{-\bar{\gamma}_{\text{ac}}(\tau_{1}+\tau_{2})}$
$\displaystyle=g_{\mathrm{ac}}^{(3)}(\tau_{1},\tau_{2}).$

Again, the same is true for
$h_{\mathrm{ac}}^{(4)}(\tau_{1},\tau_{2},\infty)=h_{\mathrm{ac}}^{(3)}(\tau_{1},\tau_{2})$.

###  .9 Coherences of k-quanta- Subtracted/Added Thermal States

In the main text we discussed out-of-equilibrium k-quanta-subtracted/added
thermal states. In this section, we derive the relation between the mean
occupancies and coherences of such states and those of the thermal states they
are generated from.

The subtraction (addition) of a single phonon at time $t=0$ from a state
described by density matrix $\rho$ yields the state with density matrix
$\rho_{-1}$ ($\rho_{+1}$)

$\rho_{-1}=\frac{b(0)\rho b^{\dagger}(0)}{\text{Tr}[b(0)\rho
b^{\dagger}(0)]}=\frac{b(0)\rho b^{\dagger}(0)}{\text{Tr}[\rho
b^{\dagger}(0)b(0)]}=\frac{b(0)\rho b^{\dagger}(0)}{\langle
b^{\dagger}(0)b(0)\rangle_{\rho}}$ (16) $\rho_{+1}=\frac{b^{\dagger}(0)\rho
b(0)}{\text{Tr}[b^{\dagger}(0)\rho b(0)]}=\frac{b^{\dagger}(0)\rho
b(0)}{\text{Tr}[\rho b(0)b^{\dagger}(0)]}=\frac{b^{\dagger}(0)\rho
b(0)}{\langle b(0)b^{\dagger}(0)\rangle_{\rho}}$

where the denominator signifies the normalization, and
$\langle\cdot\rangle_{\rho}$ is the expectation value for state $\rho$. The
coherences and mean phonon number of the single-phonon subtracted (added)
state $\rho_{-1}$ ($\rho_{+1}$) can thus be evaluated in terms of the
coherences of the steady state $\rho$, as follows.

Coherences By way of illustration, consider the second-order coherence of the
single-phonon subtracted state,

$g_{\text{ac}}^{(2)}\big{|}_{-1}(\tau)=\frac{\langle
b^{\dagger}(0)b^{\dagger}(\tau)b(\tau)b(0)\rangle_{\rho_{-1}}}{\langle
b^{\dagger}(0)b(0)\rangle_{\rho_{-1}}\ \langle
b^{\dagger}(\tau)b(\tau)\rangle_{\rho_{-1}}}$

Using Eq. 16, the numerator equates to

$\frac{\langle
b^{\dagger}(0)b^{\dagger}(0)b^{\dagger}(\tau)b(\tau)b(0)b(0)\rangle_{\rho}}{\langle
b^{\dagger}(0)b(0)\rangle_{\rho}}=g_{\text{ac}}^{(3)}(0,\tau)\ \langle
b^{\dagger}(0)b(0)\rangle_{\rho}^{2}$

and the denominator equates to

$\displaystyle\frac{\langle
b^{\dagger}(0)b^{\dagger}(0)b(0)b(0)\rangle_{\rho}}{\langle
b^{\dagger}(0)b(0)\rangle_{\rho}}\ \frac{\langle
b^{\dagger}(0)b^{\dagger}(\tau)b(\tau)b(0)\rangle_{\rho}}{\langle
b^{\dagger}(0)b(0)\rangle_{\rho}}$ $\displaystyle=g_{\text{ac}}^{(2)}(0)\
g_{\text{ac}}^{(2)}(\tau)\ \langle b^{\dagger}(0)b(0)\rangle_{\rho}^{2}$

so that

$g_{\text{ac}}^{(2)}\big{|}_{-1}(\tau)=\frac{g_{\text{ac}}^{(3)}(0,\tau)}{g_{\text{ac}}^{(2)}(0)\
g_{\text{ac}}^{(2)}(\tau)}$ (17)

And similarly,

$h_{\text{ac}}^{(2)}\big{|}_{+1}(\tau)=\frac{h_{\text{ac}}^{(3)}(0,\tau)}{h_{\text{ac}}^{(2)}(0)\
h_{\text{ac}}^{(2)}(\tau)}$ (18)

For the helium acoustic oscillator in this work that is initialized in a
thermal state $\rho_{\text{th}}$, this gives

$g_{\text{ac}}^{(2)}\big{|}_{-1,\text{th}}(\tau)=h_{\text{ac}}^{(2)}\big{|}_{+1,\text{th}}(\tau)=\frac{1+2e^{-\bar{\gamma}_{\text{ac}}t}}{1+e^{-\bar{\gamma}_{\text{ac}}t}}.$

These functions are plotted as solid lines alongside the measurements in Fig.
4(b) of the main text. Note that evaluation of the data points in Fig. 4(b)
makes explicit use of post-selected data. This is because the experimental
realization of Eq. 17, 18,

$\displaystyle\frac{g_{\text{ac}}^{(3)}(0,\tau)}{g_{\text{ac}}^{(2)}(0)\
g_{\text{ac}}^{(2)}(\tau)}$
$\displaystyle=\frac{C^{(3)}_{\text{AS}}(0^{+},\tau)/C^{(3)}_{\text{AS}}(\infty,\infty)}{C^{(2)}_{\text{AS}}(0^{+})/C^{(2)}_{\text{AS}}(\infty)\cdot
C^{(2)}_{\text{AS}}(\tau)/C^{(2)}_{\text{AS}}(\infty)},$
$\displaystyle\frac{h_{\text{ac}}^{(3)}(0,\tau)}{h_{\text{ac}}^{(2)}(0)\
h_{\text{ac}}^{(2)}(\tau)}$
$\displaystyle=\frac{C^{(3)}_{\text{S}}(0^{+},\tau)/C^{(3)}_{\text{S}}(\infty,\infty)}{C^{(2)}_{\text{S}}(0^{+})/C^{(2)}_{\text{S}}(\infty)\cdot
C^{(2)}_{\text{S}}(\tau)/C^{(2)}_{\text{S}}(\infty)}$

involves $C^{(3)}_{\text{AS(S)}}(0^{+},\tau)$ (and
$C^{(2)}_{\text{AS(S)}}(\tau)$) which are photon pairs (and counts)
conditioned on the arrival of an anti-Stokes (Stokes) photon at time $t=0$,
which heralds the subtraction (addition) of one phonon. As described in Sec.
.8 and Fig. S11, each $\infty$ in the normalizations
$C^{(3)}_{\text{AS(S)}}(\infty,\infty)$ and $C^{(2)}_{\text{AS(S)}}(\infty)$
represents times $\tau>10\bar{\gamma}_{\text{ac}}^{-1}$. The time $t=0^{+}$
represents the finite time of $10\text{\,}\mathrm{\SIUnitSymbolMicro s}$ (or
$5\text{\,}\mathrm{\SIUnitSymbolMicro s}$) of the experimental data bin used.
It corresponds to the time of the first click used to measure the
$g^{(2)}_{\text{ac}}$ (or $h^{(2)}_{\text{ac}})$ of the phonon subtracted (or
added) state. This leads to a finite-time correction of the theoretical
predictions, using Eq. 17, 18. This finite-time correction results in a slight
shift of the theoretical curves shown in Fig. 4(b); however since it is
negligible on the scale of Fig. 4(b) we choose not to show it.

More generally, the subtraction (addition) of $k$ phonons at time $t=0$ from a
state described by density matrix $\rho$ yields the state with density matrix
$\rho_{-k}$ ($\rho_{+k}$)

$\rho_{-k}=\frac{(b(0))^{k}\ \rho\
(b^{\dagger}(0))^{k}}{\langle(b^{\dagger}(0))^{k}\ (b(0))^{k}\rangle_{\rho}}$
(19) $\rho_{+k}=\frac{(b^{\dagger}(0))^{k}\ \rho\
(b(0))^{k}}{\langle(b(0))^{k}\ (b^{\dagger}(0))^{k}\rangle_{\rho}}$ (20)

Through an _ab-initio_ evaluation similar to that illustrated above, the
$n^{\text{th}}$-order coherence of a $k$-phonon subtracted (added) state is
evaluated in terms of the coherences of the steady state $\rho$ to be

$g_{\text{ac}}^{(n)}\big{|}_{-k}(\boldsymbol{\tau})=\frac{g_{\text{ac}}^{(k+n)}(\boldsymbol{0}^{\otimes
k},\boldsymbol{\tau})\
\big{(}g_{\text{ac}}^{(k)}(\boldsymbol{0})\big{)}^{n-1}}{g_{\text{ac}}^{(k+1)}(\boldsymbol{0})\big{[}\prod_{p=1}^{n-1}g_{\text{ac}}^{(k+1)}\big{(}\boldsymbol{0}^{\otimes(k-1)},t_{p}\big{)}\big{]}}$
$h_{\text{ac}}^{(n)}\big{|}_{+k}(\boldsymbol{\tau})=\frac{h_{\text{ac}}^{(k+n)}(\boldsymbol{0}^{\otimes
k},\boldsymbol{\tau})\
\big{(}h_{\text{ac}}^{(k)}(\boldsymbol{0})\big{)}^{n-1}}{h_{\text{ac}}^{(k+1)}(\boldsymbol{0})\big{[}\prod_{p=1}^{n-1}h_{\text{ac}}^{(k+1)}\big{(}\boldsymbol{0}^{\otimes(k-1)},t_{p}\big{)}\big{]}}$

where $\boldsymbol{\tau}=(\tau_{1},\tau_{2},...,\tau_{n-1})$,
$\boldsymbol{0}^{\otimes k}=(0,0,...\ k\ \text{times})$, and
$t_{p}=\sum_{j=1}^{p}\tau_{j}$ is the $(p+1)^{\text{th}}$ time. (Recall that
$\tau_{k}$ is the delay between the $k^{\text{th}}$ and $(k+1)^{\text{th}}$
time.)

For the helium acoustic oscillator in this work that is initialized in a
thermal state $\rho_{\text{th}}$, then,

$\displaystyle g^{(2)}_{\text{ac}}|_{-k}(\tau)$
$\displaystyle=\frac{g^{(k+2)}_{\text{ac}}(\textbf{0}^{\otimes
k},\tau)g^{(k)}_{\text{ac}}(\textbf{0})}{g^{(k+1)}_{\text{ac}}(\textbf{0})g^{(k+1)}_{\text{ac}}(\boldsymbol{0}^{\otimes(k-1)},\tau)},$
$\displaystyle g^{(2)}_{\text{ac}}|_{-k,\text{th}}(\tau)$
$\displaystyle=\frac{1+(k+1)e^{-\bar{\gamma}_{\text{ac}}\tau}}{1+ke^{-\bar{\gamma}_{\text{ac}}\tau}},$

and

$\displaystyle h^{(2)}_{\text{ac}}|_{+k}(\tau)$
$\displaystyle=\frac{h^{(k+2)}_{\text{ac}}(\textbf{0}^{\otimes
k},\tau)h^{(k)}_{\text{ac}}(\textbf{0})}{h^{(k+1)}_{\text{ac}}(\textbf{0})h^{(k+1)}_{\text{ac}}(\boldsymbol{0}^{\otimes(k-1)},\tau)},$
$\displaystyle h^{(2)}_{\text{ac}}|_{+k,\text{th}}(\tau)$
$\displaystyle=\frac{1+(k+1)e^{-\bar{\gamma}_{\text{ac}}\tau}}{1+ke^{-\bar{\gamma}_{\text{ac}}\tau}}.$

Mean phonon numbers The mean phonon number of a $k$-phonon subtracted state at
time $\tau$ after its generation is given by (using Eq. 19)

$n_{\text{ac}}^{-k}(\tau)=\langle
b^{\dagger}(\tau)b(\tau)\rangle_{\rho_{-k}}=\frac{\langle(b^{\dagger}(0))^{k}b^{\dagger}(\tau)b(\tau)(b(0))^{k}\rangle_{\rho}}{\langle(b^{\dagger}(0))^{k}\
(b(0))^{k}\rangle_{\rho}}.$

Dividing both sides by $n_{\text{ac}}(0)=\langle
b^{\dagger}(0)b(0)\rangle_{\rho}$ yields

$\frac{n_{\text{ac}}^{-k}(\tau)}{n_{\text{ac}}(0)}=\frac{g_{\text{ac}}^{(k+1)}(\boldsymbol{0}^{\otimes(k-1)},\tau)}{g_{\text{ac}}^{(k)}(\boldsymbol{0})}$
(21)

Eq. 21 has an intuitive interpretation in terms of the scattering rates of
photons by the mechanical oscillator. The numerator
$g_{\text{ac}}^{(k+1)}(\boldsymbol{0}^{\otimes(k-1)},\tau)$ is proportional to
the rate of detecting an anti-Stokes photon at time $t=\tau$ after a
$k-$phonon subtraction has been heralded at time $t=0$ (by the detection of
$k$-photons at time $t=0$). The denominator
$g_{\text{ac}}^{(k)}(\boldsymbol{0})=g_{\text{ac}}^{(k+1)}(\boldsymbol{0}^{\otimes(k-1)},\infty)$
is proportional to the rate of detecting an anti-Stokes photon at time
$t\rightarrow\infty$ after a $k-$phonon subtraction has been heralded at time
$t=0$. The proportionality constant for both is simply the rate of detecting
uncorrelated arrivals of $k+1$ photons. The numerator is also proportional to
$n_{\text{ac}}^{-k}(\tau)$ (as the anti-Stokes scattering rate is proportional
to the mean phonon occupancy – the proportionality constant is simply the
scattering rate $\gamma_{\text{AS}}$, of the main text, times the net
detection efficiency). Whereas the denominator is proportional to
$n_{\text{ac}}^{-k}(\infty)=n_{\text{ac}}(0)$; the latter equality coming from
the fact that for $\tau\rightarrow\infty$, the heralded $k$-phonon subtracted
state equilibrates back to the steady state. Taking a ratio of these rates,
interpreted these two ways, readily yields Eq. 21.

Similarly, the mean phonon number of a $k$-phonon added state at time $\tau$
after its generation is given by (using Eq. 20)

$\displaystyle n_{\text{ac}}^{+k}(\tau)+1$ $\displaystyle=\langle
b^{\dagger}(\tau)b(\tau)\rangle_{\rho_{+k}}+1$ $\displaystyle=\langle
b(\tau)b^{\dagger}(\tau)\rangle_{\rho_{+k}}$
$\displaystyle=\frac{\langle(b(0))^{k}b(\tau)b^{\dagger}(\tau)(b^{\dagger}(0))^{k}\rangle_{\rho}}{\langle(b(0))^{k}\
(b^{\dagger}(0))^{k}\rangle_{\rho}}.$

Dividing both sides by $n_{\text{ac}}(0)+1=\langle
b(0)b^{\dagger}(0)\rangle_{\rho}$ yields

$\frac{n_{\text{ac}}^{+k}(\tau)+1}{n_{\text{ac}}(0)+1}=\frac{h_{\text{ac}}^{(k+1)}(\boldsymbol{0}^{\otimes(k-1)},\tau)}{h_{\text{ac}}^{(k)}(\boldsymbol{0})}.$

This equation has a similar intuitive interpretation as Eq. 21, but with the
rates now being proportional to the Stokes scattering rates, which are
proportional to $n_{\text{ac}}+1$.

Note again that the experimental evaluation of the data points in Fig. 4(a)
involves the explicit use of post-selected data. This is because the
evaluation

$\displaystyle\frac{g_{\text{ac}}^{(k+1)}(\boldsymbol{0}^{\otimes(k-1)},\tau)}{g_{\text{ac}}^{(k+1)}(\boldsymbol{0}^{\otimes(k-1)},\infty)}$
$\displaystyle=\frac{C^{(k+1)}_{\text{AS}}(\boldsymbol{0}^{\otimes(k-1)},\tau)/C^{(k+1)}_{\text{AS}}(\boldsymbol{\infty}^{\otimes
k})}{C^{(k+1)}_{\text{AS}}(\boldsymbol{0}^{\otimes(k-1)},\infty)/C^{(k+1)}_{\text{AS}}(\boldsymbol{\infty}^{\otimes
k})},$
$\displaystyle\frac{h_{\text{ac}}^{(k+1)}(\boldsymbol{0}^{\otimes(k-1)},\tau)}{h_{\text{ac}}^{(k+1)}(\boldsymbol{0}^{\otimes(k-1)},\infty)}$
$\displaystyle=\frac{C^{(k+1)}_{\text{S}}(\boldsymbol{0}^{\otimes(k-1)},\tau)/C^{(k+1)}_{\text{S}}(\boldsymbol{\infty}^{\otimes
k})}{C^{(k+1)}_{\text{S}}(\boldsymbol{0}^{\otimes(k-1)},\infty)/C^{(k+1)}_{\text{S}}(\boldsymbol{\infty}^{\otimes
k})}$

involves $C^{(k+1)}_{\text{AS(S)}}(\boldsymbol{0}^{\otimes(k-1)},\tau)$ which
are photon counts conditioned on the arrival of $k$ anti-Stokes (Stokes)
photons at time $t=0$, which herald the subtraction (addition) of $k$ phonons.
As described earlier, each $\infty$ in the normalizations represents times
$\tau>10\bar{\gamma}_{\text{ac}}^{-1}$.

## References

  * Shelby  _et al._ [1985] R. M. Shelby, M. D. Levenson, and P. W. Bayer, Guided acoustic-wave brillouin scattering, Phys. Rev. B 31, 5244 (1985).
  * Takefushi  _et al._ [2019] N. Takefushi, M. Yoshida, K. Kasai, T. Hirooka, and M. Nakazawa, Gawbs noise characteristics in digital coherent transmission in various optical fibers, in  _24th OptoElectronics and Communications Conference (OECC) and 2019 International Conference on Photonics in Switching and Computing (PSC)_ (2019) pp. 1–3.
  * Marquardt  _et al._ [2007] F. Marquardt, J. P. Chen, A. A. Clerk, and S. M. Girvin, Quantum theory of cavity-assisted sideband cooling of mechanical motion, Phys. Rev. Lett. 99, 093902 (2007).
  * Burenkov  _et al._ [2013] V. Burenkov, H. Xu, B. Qi, R. H. Hadfield, and H.-K. Lo, Investigations of afterpulsing and detection efficiency recovery in superconducting nanowire single-photon detectors, Journal of Applied Physics 113, 213102 (2013).
  * Fujiwara  _et al._ [2011] M. Fujiwara, A. Tanaka, S. Takahashi, K. Yoshino, Y. Nambu, A. Tajima, S. Miki, T. Yamashita, Z. Wang, A. Tomita, _et al._ , Afterpulse-like phenomenon of superconducting single photon detector in high speed quantum key distribution system, Optics express 19, 19562 (2011).
  * Kelley and Kleiner [1964] P. Kelley and W. Kleiner, Theory of electromagnetic field measurement and photoelectron counting, Physical Review 136, A316 (1964).
  * Riedinger  _et al._ [2016] R. Riedinger, S. Hong, R. A. Norte, J. A. Slater, J. Shang, A. G. Krause, V. Anant, M. Aspelmeyer, and S. Gröblacher, Non-classical correlations between single photons and phonons from a mechanical oscillator, Nature (London) 530, 313 (2016).
  * Børkje  _et al._ [2011] K. Børkje, A. Nunnenkamp, and S. M. Girvin, Proposal for Entangling Remote Micromechanical Oscillators via Optical Measurements, Phys. Rev. Lett. 107, 123601 (2011).
  * Lemonde  _et al._ [2014] M.-A. Lemonde, N. Didier, and A. A. Clerk, Antibunching and unconventional photon blockade with gaussian squeezed states, Phys. Rev. A 90, 063824 (2014).
  * Gardiner and Collett [1985] C. W. Gardiner and M. J. Collett, Input and output in damped quantum systems: Quantum stochastic differential equations and the master equation, Phys. Rev. A 31, 3761 (1985).
  * Weiss [2012] U. Weiss,  _Quantum Dissipative Systems_, 4th ed. (World Scientific, 2012).
  * Loudon [2000] R. Loudon,  _The Quantum Theory of Light_ (OUP Oxford, 2000).