Simultaneous retrodiction of multi-mode optomechanical systems using matched filters

TL;DR

This paper develops a method for retrodicting the initial states of multiple optomechanical oscillators using matched filters, achieving quantum-limited sensitivity and demonstrating the approach experimentally for two modes.

Contribution

It introduces a linear matched filter approach for simultaneous multi-mode state retrodiction in optomechanical systems, reaching the Standard Quantum Limit.

Findings

Optimal estimators are derived from a detailed model of the homodyne signal.

High-cooperativity measurements enable retrodiction at the SQL.

Experimental two-mode retrodiction demonstrates practical feasibility.

Abstract

Generation and manipulation of many-body entangled states is of considerable interest, for applications in quantum simulation or sensing, for example. Measurement and verification of the resulting many-body state presents a formidable challenge, however, which can be simplified by multiplexed readout using shared measurement resources. In this work, we analyze and demonstrate state retrodiction for a system of optomechanical oscillators coupled to a single-mode optical cavity. Coupling to the shared cavity field facilitates simultaneous optical measurement of the oscillators' transient dynamics at distinct frequencies. Optimal estimators for the oscillators' initial state can be defined as a set of linear matched filters, derived from a detailed model for the detected homodyne signal. We find that the optimal state estimate for optomechanical retrodiction is obtained from high-cooperativity measurements, reaching estimate sensitivity at the Standard Quantum Limit (SQL). Simultaneous estimation of the state of multiple oscillators places additional limits on the estimate precision, due to the diffusive noise each oscillator adds to the optomechanical signal. However, we show that the sensitivity of simultaneous multi-mode state retrodiction reaches the SQL for sufficiently well-resolved oscillators. Finally, an experimental demonstration of two-mode retrodiction is presented, which requires further accounting for technical fluctuations of the oscillator frequency.

Equations

$$\displaystyle\langle\hat{X}_{1}(t)\rangle=e^{-\Gamma_{1}t/2}\Bigl{(}\langle\hat{X}_{1}(0)\rangle\cos\omega_{1}t+\langle\hat{P}_{1}(0)\rangle\sin\omega_{1}t\Bigr{)}\text{,}$$

$$⟨\Delta {\hat{X}}_i^2⟩=⟨\Delta {\hat{P}}_i^2⟩={\Delta n}_i\ge \frac{1}{2}\text{,}$$

$${\hat{c}}^{\text{AM}}$$

$${\hat{X}}_i$$

$$\mathcal{H}=-\hslash \Delta {\hat{c}}^{†}\hat{c}+\sum _i\hslash {\omega }_i{\hat{a}}_i^{†}{\hat{a}}_i+2\sum _i\hslash \sqrt{\bar{n}}{g}_i{\hat{c}}^{\text{AM}}{\hat{X}}_i$$

$$\dot{\hat{c}}=-i\sum _i\sqrt{2\bar{n}}{g}_i{\hat{X}}_i-\kappa \hat{c}+\sqrt{2\kappa }\hat{\xi }$$

$${\hat{c}}^{\text{AM}}(t)$$

$${\hat{c}}^{\text{PM}}(t)$$

$$\hat{S}(t)=\sqrt{2}\sum _i{g}_i{\hat{X}}_i(t)+\sqrt{P_{\text{SN}}}{\hat{\xi }}_{\text{SN}}(t)\text{.}$$

$${\dot{\hat{a}}}_i=(-{\Gamma }_i/2-i{\omega }_i){\hat{a}}_i-i\sqrt{C_i{\Gamma }_i}{\hat{\xi }}^{\text{AM}}+\sqrt{{\Gamma }_i}{\hat{\eta }}_i\text{,}$$

$$C_i=4\bar{n}{g}_i^2/\kappa {\Gamma }_i\text{,}$$

$${\hat{X}}_i(t)={\hat{\mathbf{Q}}}_i^{\mathsf{T}}{\mathbf{r}}_i(t)+{\hat{D}}_i(t)\text{,}$$

$${\mathbf{r}}_i(t)=e^{-{\Gamma }_{i}t/2}\left(\begin{array}{c}\cos {\omega }_{i}t\\ \sin {\omega }_{i}t\end{array}\right)\Theta (t)\text{,}$$

$${\hat{D}}_i(t)={\int }_0^{\infty }d\tau {\hat{\mathbf{d}}}_i^{\mathsf{T}}(\tau ){\mathbf{r}}_i(t-\tau )\text{,}$$

$${\hat{\mathbf{d}}}_i(\tau )=\sqrt{{\Gamma }_i}\left(\begin{array}{c}{{\hat{\eta }}^{\text{AM}}}_i(\tau )\\ {{\hat{\eta }}^{\text{PM}}}_i(\tau )-\sqrt{2C_i}{\hat{\xi }}^{\text{AM}}(\tau )\end{array}\right)\text{.}$$

$$\hat{S}(t)=\sum _i\sqrt{2}{g}_i\left[{\hat{\mathbf{Q}}}_i^{\mathsf{T}}{\mathbf{r}}_i(t)+{\hat{D}}_i(t)\right]+\sqrt{P_{\text{SN}}}{\hat{\xi }}_{\text{SN}}(t)\text{,}$$

$$\breve{\mathbf{q}}={\int }_0^{t_f}dt\mathbf{m}(t)\hat{S}(t)\text{,}$$

$$⟨\breve{\mathbf{q}}⟩$$

$$\hat{\mathbf{Q}}$$

$${\mathrm{J}}_{ij}=\sqrt{2}{g}_j{\int }_0^{t_f}dt{\mathbf{m}}_i(t){\mathbf{r}}_j^{\mathsf{T}}(t)$$

$${\mathbf{m}}_i^{\text{OLS}}(t)\equiv {\mathbf{r}}_i(t)\text{.}$$

$$\breve{\mathbf{Q}}=\left(\begin{array}{c}{\breve{\mathbf{Q}}}_1\\ {\breve{\mathbf{Q}}}_2\\ ⋮\\ {\breve{\mathbf{Q}}}_N\end{array}\right)\equiv {\mathrm{J}}^{-1}{\int }_0^{t_f}dt\mathbf{m}(t)\hat{S}(t)\text{,}$$

$$\displaystyle\Phi^{\text{OLS}}[{\bm{m}}(t)]\equiv\int_{0}^{{t_{f}}}dt\,\Bigl{[}\hat{S}(t)-\sqrt{2}\sum_{j}g_{j}\breve{{\bm{Q}}}^{\mathsf{T}}_{j}{\bm{r}}_{j}(t)\Bigr{]}^{2}\text{,}$$

$${\mathrm{\Omega }}_{nm}=2\sum _{kl}{g}_k{g}_l⟨{\hat{D}}_k(t_n){\hat{D}}_l(t_m)⟩+P_{\text{SN}}{f}_s{\delta }_{nm}\text{,}$$

$${\mathbf{m}}^{\text{GLS}}(t_n)\equiv \sum _{m=0}^{N_t-1}[{\mathrm{\Omega }}^{-1}{]}_{nm}\mathbf{r}(t_m)\text{.}$$

$$\mathrm{cov}[\hat{\mathbf{Q}}]\equiv \mathrm{Re}[⟨\hat{\mathbf{Q}}{\hat{\mathbf{Q}}}^{\mathsf{T}}⟩]-⟨\hat{\mathbf{Q}}⟩⟨{\hat{\mathbf{Q}}}^{\mathsf{T}}⟩\text{.}$$

$$\mathrm{\Sigma }\equiv \mathrm{cov}[\breve{\mathbf{Q}}]=\mathrm{cov}[\hat{\mathbf{Q}}]+\mathrm{T}+\mathrm{B}+\mathrm{M}\text{,}$$

$$\mathrm{M}$$

$${\mathrm{N}}_k(\tau )=\sqrt{2{\Gamma }_k}{g}_k{\mathrm{J}}^{-1}{\int }_0^{t_f}dt\mathbf{m}(t){\mathbf{r}}_k^{\mathsf{T}}(t-\tau )\text{,}$$

$$\mathrm{T}$$

$$\mathrm{B}$$

Full text

Building many-body quantum systems by assembling ensembles of well-controlled quantum modes with tunable interactions is a promising path toward quantum simulation and quantum information processing. The increased dimensionality of many-body systems, however, makes measurement of entangled states challenging, because of the large number of observables required to fully characterize the quantum state. The physical resources necessary to perform these measurements can be reduced by using a shared measurement ‘bus,’ coupled to multiple quantum degrees of freedom, facilitating multiplexed measurement of their quantum states, such as demonstrated with arrays of super-conducting qubits coupled to a common strip-line resonator Jerger et al. (2012); Chen et al. (2012). Each mode can be independently measured by being sequentially coupled to the measurement bus, reading out each of their states with independent temporal modes of the output field Palomaki et al. (2013a). Alternatively, if the dynamics of each mode are spectrally resolved at distinct frequencies, the many-body state of a system can be simultaneously measured through a continuous weak measurement.

High-finesse optical resonators provide a particularly powerful tool for measuring and controlling the dynamics of diverse systems, demonstrated in cavity optomechanics Aspelmeyer et al. (2014), collective atomic spin optodynamics Leroux et al. (2010); Vasilakis et al. (2015); Kohler et al. (2017), and in the emerging field of cavity optomagnonics Zhang et al. (2016); Haigh et al. (2015); Osada et al. (2016); Liu et al. (2016); Viola Kusminskiy et al. (2016). Multiple modes of diverse systems can be simultaneously coupled to a common cavity field, facilitating simultaneous optical measurement Botter et al. (2013), in addition to long-range optically mediated interactions Shkarin et al. (2014); Spethmann et al. (2016); Kohler et al. (2018).

The optical field leaking out of the cavity provides a continuous measurement of the system dynamics. Retrodiction of past quantum states from continuous measurements has been theoretically described through back-propagation of an effect matrix Gammelmark et al. (2013), closely related to the quantum theory of smoothing Tsang (2009). For Gaussian dynamics, these estimates allow a particularly compact description by evolving the phase space mean and covariance Zhang and Mølmer (2017); Huang and Sarovar (2018); Lammers (2018); Lammers and Hammerer , recently demonstrated on an optomechanical system Rossi et al. (2019). Related analysis employing linear filters applied to the output of an optical interferometer has been proposed for conditional quantum state preparation Müller-Ebhardt et al. (2009) and verification Miao et al. (2010), applicable for tests of macroscopic quantum mechanics in gravitational-wave detectors such as LIGO. For a multi-mode optomechanical system, estimation of a collective quadrature has also been proposed through temporal modulation of the measurement strength Moore et al. (2016).

In this work, we consider estimation of the initial state, at time $t=0$, of a multi-mode optomechanical system of $N$ harmonic oscillators, illustrated in Fig. 1a, retrodicted using matched filters applied to continuous measurement of its subsequent free evolution. The oscillators are dispersively coupled to a common single-mode optical cavity, which is driven on resonance with a coherent probe. The reflected optical field performs a continuous weak measurement of the sum of oscillator displacements, recorded using a balanced homodyne detector. Knowledge of the coherent system evolution allows the initial state to be inferred from the observed transient dynamics. However, measurement backaction, arising from quantum fluctuations of the cavity field, perturbs the trajectory of each oscillator’s evolution, adding incoherent noise to the subsequent measurement record, which must be considered in obtaining the optimal state estimate. We focus in particular on the experimentally relevant system of linear cavity optomechanics, allowing derivation of analytic results for the optimal state estimators which can be directly applied to experimental measurements. However, the formalism developed in this work can be directly generalized for measurement and retrodiction of any linearizable system undergoing transient or non-steady state dynamics driven by Markovian noise.

Retrodiction can be illustrated by considering continuous measurement of the position $\hat{X}_{1}(t)$ of a single harmonic oscillator. Estimates for the average position quadrature $\langle\hat{X}_{1}(0)\rangle$ and momentum quadrature $\langle\hat{P}_{1}(0)\rangle$ of the oscillator’s initial state can be recovered from measurement of its subsequent coherent evolution

[TABLE]

where $\omega_{1}$ and $\Gamma_{1}$ are the oscillators frequency and energy damping rate, respectively. Incoherent noise from measurement backaction and the oscillator’s intrinsic thermal bath also perturbs the oscillator’s trajectory during measurement, as simulated in Fig. 1c. The accumulated diffusion from these noise baths reduces the relative signal-to-noise in the measurement record at later times.

Quadrature estimators for the oscillator’s initial state can be defined as linear filters of the recorded homodyne signal $\hat{S}(t)$, which appropriately weight the measured signal at each subsequent time $t$. The optimally ‘matched’ filters must appropriately balance the coherent evolution and incoherent diffusion, in order to minimize the total estimate error. The use of matched filters to recover signals of a known form from stationary additive Gaussian noise is well described in standard textbooks Wainstein and Zubakov (1970) and commonly employed, for instance, for gravitational wave detection Abbott et al. (2004). However, for retrodiction of optomechanical systems, diffusion of the oscillators’ state, driven by quantum backaction and the thermal bath, introduces non-stationary noise which accumulates throughout the observed transient signals, requiring careful accounting of the full two-time correlation of the signal noise.

For high-quality oscillators, where $\omega_{i}\gg\Gamma_{i}$, the initial state undergoes multiple coherent oscillations during the subsequent ringdown. Retrodiction from the observed trajectory, therefore, obtains approximately equal information about each initial quadrature amplitude, recovering an estimate of the oscillator’s initial state that is independent of the oscillator’s phase. The Heisenberg uncertainty principle $\langle\Delta\hat{X}_{i}^{2}\rangle\langle\Delta\hat{P}_{i}^{2}\rangle\geq 1/4$ establishes a fundamental bound, known as the Standard Quantum Limit (SQL) Caves (1980); Braginsky et al. (1992), for the minimum noise added by such phase-independent measurements

[TABLE]

quantified here as an added effective thermal phonon occupation ${\Delta n}_{i}$.

In this work, we consider retrodiction of Gaussian states and demonstrate inference of two-mode squeezed states in a multi-mode optomechanical system. Full tomography of a general many-body state involves estimation of the entire density matrix, which contains further information about all higher-order moments of the oscillator quadratures. The following analysis could be extended to estimate these higher-order moments of the multi-mode state. Though beyond the scope of this work, it would be worth considering which features of multi-mode quantum states can be retrodicted from such phase-independent optomechanical measurements Nha et al. (2010) or how measurements beyond the SQL Miao et al. (2010); Lei et al. (2016); Ockeloen-Korppi et al. (2016); Møller et al. (2017) can be performed simultaneously on multi-mode systems.

.1 Summary of main results

We summarize here the primary conclusions of this work and give an overview of the following sections. In Sec. I, we describe a general model for homodyne measurement of a multi-mode optomechanical system. Measurement backaction arises from radiation pressure shot noise, which appears as a shared noise bath driving correlated diffusion of the oscillators. A general set of quadrature estimators is defined in Sec. II, in terms of linear filters applied to the recorded homodyne signal. Optimal quadrature filters, which minimize the estimator variance, are derived using linear regression, accounting for temporal correlations introduced by the oscillators’ diffusive motion.

In Sec. III, we derive the estimator covariance matrix, describing the imprecision added to the quadrature estimates by each source of noise in the recorded signal. The estimator covariance can be measured using an ensemble of repeated measurements, which allows estimation of the covariance of a squeezed state, after correcting for the added noise covariance.

We derive an analytic approximation for the optimal single-oscillator filter in Sec. IV, which minimizes the total added covariance. The measurement strength, parameterized by the cooperativity $C_{i}$, quantifies the relative rate that information is gained from the system. The optimal estimate, with sensitivity reaching the SQL, is obtained from a high-cooperativity measurement $\nu_{i}+1\ll C_{i}\ll\omega_{i}/\Gamma_{i}$, which is bounded below by the thermal bath occupation $\nu_{i}$ and above by the oscillator quality factor. This condition ensures the measurement rate far exceeds the loss of information to the thermal environment. The corresponding measurement backaction drives rapid diffusion of the oscillators, which is suppressed in the quadrature estimates by the appropriately optimized filters.

We demonstrate these matched-filter estimators on a simulated two-mode system in Sec. V. Simultaneous estimation of the state of multiple oscillators places additional constraints on the optimal estimate sensitivity, explored in Sec. VI. We show that the state of multiple oscillators can be retrodicted from the measurement record with precision at the SQL, if their frequencies are resolved by many linewidths.

Experimental results of matched filter estimates are presented in Sec. VII, obtained from a recent demonstration of the negative-mass instability between collective atomic spin and motion Kohler et al. (2018). Additional experimental complications from shot-to-shot fluctuations of system parameters had to be included in the model to recover accurate retrodicted estimates. Finally, conclusions and outlook are summarized in Sec. VIII.

For clarity of notation throughout, vectors will be notated in bold (e.g. ${\bm{v}}$), and matrices in Roman typeface (e.g. $\mathrm{M}$). Hermitian amplitude and phase quadratures of bosonic operators, such as $\hat{c}$, are defined according to

[TABLE]

respectively. In particular, the quadratures of the optomechanical oscillators will be notated as a generalized position and momentum,

[TABLE]

respectively.

I Simultaneous optomechanical measurement

Consider an ensemble of $N$ harmonic oscillators, illustrated in Fig. 1a, described by bosonic operators $\hat{a}_{i}$ evolving at frequencies $\omega_{i}$, which are dispersively coupled to a driven single-mode optical cavity with independent linear optomechanical coupling strengths $g_{i}$ Aspelmeyer et al. (2014). The cavity resonance frequency is shifted by the sum of the oscillators’ displacements, which can be continuously measured by driving the cavity on resonance, such that the oscillators’ motion modulates the phase quadrature of reflected light and is recorded using optical homodyne detection.

For small displacements, the dispersive shift of the cavity frequency is small relative to the cavity linewidth ${\sum_{i}g_{i}\langle\hat{a}^{\dagger}_{i}+\hat{a}_{i}\rangle\ll\kappa}$, and the dynamics of the cavity field can be linearized in terms of fluctuations $\hat{c}$ around an average cavity photon number $\bar{n}$, in a frame rotating at the cavity drive frequency ${\omega_{\text{p}}}$, yielding a multi-mode generalization of the linearized optomechanical Hamiltonian Aspelmeyer et al. (2014)

[TABLE]

where $\Delta={\omega_{\text{p}}}-{\omega_{\text{c}}}$ is the detuning between the drive and cavity resonance frequency ${\omega_{\text{c}}}$, having dropped constant energy terms.

In order to simultaneously measure the intrinsic dynamics of multiple independent oscillators – without introducing optically mediated coupling Shkarin et al. (2014); Spethmann et al. (2016), spring shifts Sheard et al. (2004); Corbitt et al. (2006), or damping Arcizet et al. (2006); Gigan et al. (2006); Schliesser et al. (2006) – we always consider a resonantly driven cavity, with $\Delta=0$. The Heisenberg-Langevin equation of motion for the state of the cavity field $\hat{c}(t)$

[TABLE]

is obtained from Eq. (5), with the addition of input and output terms Gardiner and Collett (1985) introducing vacuum fluctuations $\hat{\xi}$ from optical coupling to the environment, parametrized by the cavity half-linewidth $\kappa$.

For simultaneous measurement of multiple oscillators, it is advantageous to work in the fast-cavity (unresolved-sideband) regime defined by $\kappa\gg\omega_{i}$, such that the cavity field is sensitive across a wide bandwidth to dynamics of many oscillators at well-resolved frequencies. Eq. (6) can then be solved under the adiabatic approximation $\dot{\hat{c}}\approx 0$, assuming the cavity field equilibrates to the oscillators’ motion nearly instantaneously, yielding solutions for the amplitude and phase quadratures

[TABLE]

respectively. The cavity input fluctuations are assumed to be in the vacuum state, described by the two-time correlation $\langle\hat{\xi}(t){\hat{\xi}^{\dagger}}(t^{\prime})\rangle=\delta(t-t^{\prime})$.

The optical field leaking out of the cavity, determined by the boundary condition $\sqrt{2\kappa}\hat{c}={\hat{c}_{\text{out}}}-\hat{\xi}$, carries information about the oscillator dynamics in its phase quadrature. The optical phase is recorded using a balanced homodyne detector, resulting in a signal proportional to the instantaneous shift of the cavity frequency due to the displacement of each oscillator

[TABLE]

The second term describes the added measurement shot noise, due to vacuum fluctuations of the optical probe, with two-time correlation $\langle\hat{\xi}_{\text{SN}}(t)\hat{\xi}_{\text{SN}}(t^{\prime})\rangle=\delta(t-t^{\prime})$ and normalized shot noise PSD ${P_{\text{SN}}}=\kappa/(8\epsilon\bar{n})$ in terms of the total cavity photon detection efficiency $\epsilon$.

I.1 Derivation of oscillator trajectories

Accurate retrodiction of the initial state of the oscillators from their subsequent evolution requires knowledge of the system’s coherent dynamics, in addition to a complete stochastic model for all incoherent noise sources. The oscillator equations of motion derived from Eq. (5) are given by 111 assuming $g_{i}>0$, without loss of generality.

[TABLE]

assuming each oscillator is intrinsically coupled with strength $\Gamma_{i}$ to an independent Markovian bath $\hat{\eta}_{i}$, with two-time correlation $\langle\hat{\eta}_{i}(t){\hat{\eta}^{\dagger}}_{i}(t^{\prime})\rangle=(\nu_{i}+1)\delta(t-t^{\prime})$ parametrized by the equilibrium thermal occupation $\nu_{i}$.

Radiation-pressure forces introduce measurement backaction, by coupling each oscillator to the cavity amplitude fluctuations ${\hat{c}^{\text{AM}}}$ described by Eq. (7a). The optomechanical cooperativity

[TABLE]

which parametrizes the measurement strength, quantifies here the added equilibrium occupation due to diffusion from measurement backaction. This backaction noise represents a common-mode bath, driving correlated diffusion of each oscillator during measurement Spethmann et al. (2016).

Each oscillator’s trajectory is found by solving Eq. (9), simulated numerically in Fig. 1c, and is readily separated into two parts:

[TABLE]

coherent evolution of the initial phase space quadratures, summarized by the two-element vector $\hat{{\bm{Q}}}_{i}=\begin{pmatrix}\hat{X}_{i}(0),\hat{P}_{i}(0)\end{pmatrix}^{\mathsf{T}}$, and accumulated incoherent diffusion $\hat{D}_{i}(t)$. The coherent state evolution is described by a vector of quadrature impulse response functions

[TABLE]

in terms of the Heaviside step function $\Theta(t)$. The accumulated oscillator diffusion is given by the convolution

[TABLE]

of the oscillator’s response with the stochastic input noise, summarized by the input quadrature vector

[TABLE]

Substituting the oscillator trajectory given by Eq. (11) into Eq. (8), the quantum mechanical model for the measured homodyne signal can be written as

[TABLE]

which is a sum of the coherent ringdown of the initial states $\hat{{\bm{Q}}}_{i}$, the accumulated diffusion of each oscillator during measurement, and measurement shot noise. The incoherent part of Eq. (15) determines the total noise PSD of the homodyne signal, displayed for simulated signals in Fig. 1b, and contains thermal- and backaction-driven optomechanical responses in addition to the broadband shot noise floor.

Although these results are derived in the unresolved-sideband limit, Eq. (15) can be generalized for any resonantly driven cavity, provided only that $\kappa\gg\Gamma_{i}$, by accounting for a reduced effective coupling strength $g_{i}\rightarrow g_{i}\kappa/\sqrt{\kappa^{2}+\omega_{i}^{2}}$, due to suppression of shot noise fluctuations and the optomechanical response by the cavity susceptibility, and an effective phase delay of the optically measured oscillator amplitude $\hat{a}_{i}\rightarrow e^{i\phi_{i}}\hat{a}_{i}$ with $\tan\phi_{i}=\omega_{i}/\kappa$.

II Linear filter estimation

The goal of retrodiction, considered here, is to estimate the state of all oscillators at time $t=0$, represented by the quadrature vectors $\hat{{\bm{Q}}}_{i}$, from continuous measurement of their subsequent dynamics Gammelmark et al. (2013). The homodyne photocurrent is amplified electronically and then digitally sampled, resulting in a classical recorded signal that contains noise arising from measurement shotnoise and quantum backaction. For the linear systems considered in this work, such state estimation from the recorded signals can be approached as an essentially classical signal-filtering problem Müller-Ebhardt et al. (2009).

A general set of linear filters applied to the observed homodyne signal is defined as

[TABLE]

in terms of a vector of real-valued temporal weight functions ${\bm{m}}(t)$. Here, the vector of filter outputs $\breve{{\bm{q}}}$ represent projections from the infinite-dimensional space of the continuous signal $\hat{S}(t)$ onto temporal modes defined by the filter functions ${\bm{m}}(t)$.

Assuming all noise sources in Eq. (15) have zero mean, the average filter outputs can be directly evaluated and expressed as a matrix equation

[TABLE]

in terms of the $2N$-element vector of initial quadrature amplitudes $\hat{{\bm{Q}}}$ and $2N\times 2N$ normalization matrix $\mathrm{J}$

[TABLE]

defined by concatenation of the individual 2-dimensional quadrature vectors and $2\times 2$ block matrices

[TABLE]

for 2-element vectors of filter functions ${\bm{m}}_{i}(t)$ defined for each oscillator $i$.

As an initial example, the filters can be defined according to each oscillator’s quadrature response functions

[TABLE]

For this choice of filters functions, the estimator results represent the projection of the signal onto each quadrature’s coherent response function. The diagonal elements of the matrix $\mathrm{J}$ describe the normalization of each filter and the off-diagonal elements reflect the non-orthogonality between quadrature responses, due to spectral overlap from the finite oscillator linewidths.

This overlap between quadrature filters introduces spurious correlations in the raw filter outputs $\breve{q}_{i}$. Provided there are $2N$ linearly independent filter functions ${\bm{m}}(t)$ that span the space of the quadrature response functions ${\bm{r}}_{i}(t)$, then Eq. (17) can be solved to recover the average initial quadrature amplitudes by inverting the normalization matrix $\mathrm{J}$.

By extension, a complete set of unbiased quadrature estimators can, therefore, be defined as

[TABLE]

satisfying $\langle\breve{{\bm{Q}}}_{i}\rangle=\langle\hat{{\bm{Q}}}_{i}\rangle$.

Results obtained by applying these filter estimators to simulated measurements of a single oscillator’s trajectory are shown in Fig. 2a. The notation $\breve{{\bm{Q}}}$ is used here to indicate an estimator for the vector of quadrature operators $\hat{{\bm{Q}}}$, corresponding to a temporal mode of the detected optical field. The distribution of measured samples, obtained by application of Eq. (25) to the recorded homodyne traces, can be described by the statistics of the thermal and quantum noise contained in the estimator model defined by Eqs. (15) and (25).

II.1 Generalized least-squares optimization

The optimal set of filter functions ${\bm{m}}_{i}(t)$, which provide a minimum-variance unbiased estimate for any linear combination of quadratures $\hat{{\bm{Q}}}$, can be derived using the method of Least Squares, such as commonly used for linear curve-fitting. The filters defined by Eq. (24) are obtained from an Ordinary Least Squares (OLS) linear regression, by minimizing the sum of square residuals between the measured signal and the coherent model

[TABLE]

parameterized by estimators $\breve{{\bm{Q}}}_{j}$.

However, the Gauss-Markov theorem McElroy (1967) proves that these estimators are optimal only when the signal noise is temporally uncorrelated (white noise). Diffusive motion of the oscillators during the measurement, driven by thermal and backaction fluctuations, generates temporal correlations in the signal noise, indicated by structure in the PSD shown in Fig. 1b. In the presence of temporally correlated noise, a minimum variance unbiased estimator can be constructed through linear regression using the Generalized Least Squares (GLS) method Aitken (1936). This method can be understood conceptually as de-correlating the temporal signal by inverting the known two-time noise correlation function prior to performing linear regression.

The inverted noise correlation function is more easily defined for a signal sampled at discrete times $t_{n}=n/f_{s}$, with sample frequency $f_{s}$ and count ${N_{t}}=f_{s}{t_{f}}$, which is typical for most experimental applications. The discrete two-time correlation function of the added diffusive and measurement noise defines the ${N_{t}}\times{N_{t}}$ square matrix

[TABLE]

which can be inverted numerically to derive the discrete GLS filter functions

[TABLE]

These filters implicitly transform the signal to de-correlate the noise, recovering conditions to satisfy the Gauss-Markov theorem and, therefore, providing the minimum-variance unbiased estimate for any general optomechanical measurement.

The effect of diffusive noise on the optimal filters is further illustrated in Sec. IV, where an analytic formula is derived that approximates the optimal filters for a single oscillator. The numerical method defined in Eqs. (27) and (28), however, facilitates calculation of optimal filters in any conditions and is fully general for state estimation with multiple oscillators, considered in Sec. V.

III Linear estimator statistics

In addition to providing an estimate of the mean quadrature amplitudes $\langle\hat{{\bm{Q}}}_{i}\rangle$, an ensemble of estimates from repeated measurements of identically prepared states can be used to estimate the multi-mode quadrature covariance of the initial state

[TABLE]

Noise in the measured signal introduces errors into the quadrature estimates, increasing the uncertainty of the retrodicted state. For linear filters, the additive noise in the signal described by Eq. (15) introduces a systematic bias to the covariance of the observed quadrature estimates

[TABLE]

which is a sum of the actual state covariance $\operatorname{cov}[\hat{{\bm{Q}}}]$ and the additive covariance from each independent noise source – thermal diffusion $\mathrm{T}$, quantum backaction $\mathrm{B}$, and measurement shot noise $\mathrm{M}$ – indicated by shaded regions in Fig. 2.

The estimate covariance added by measurement shot noise is given by

[TABLE]

in terms of the overlap matrix between each pair of normalized quadrature filters. The filter functions of interest are generally not orthogonal, therefore measurement shot noise will induce correlated errors in the quadrature estimates, due to the spectral overlap between pairs of filters, described by non-zero off-diagonal components of $\mathrm{M}$.

Diffusion of the oscillators’ states during the measurement, driven by their intrinsic thermal baths as well as quantum backaction, also adds to both the variances and covariances of the quadrature estimates. The response of oscillator $k$ to a generic bath fluctuation at time $\tau$ has an integrated effect on the quadrature estimates, described by the $2N\times 2$ matrix-valued function

[TABLE]

which arises from applying the estimator defined by Eq. (25) to the diffusion term of the full signal model in Eq. (15). The thermal baths of each oscillator are assumed to be independent, so that the total thermal noise covariance is simply given by a sum over the variance of induced estimate perturbations, weighted by each oscillator’s bath occupation

[TABLE]

The quantum backaction, however, induces correlated diffusion of the oscillators during the measurement, since the oscillators all respond to the same amplitude fluctuations of the cavity field and have finite spectral overlap of their susceptibilities (assuming non-zero oscillator linewidths). Diffusive motion from this common optical bath induces correlated errors in the quadrature estimates, with covariance given by a sum over all oscillator pairs

[TABLE]

assuming in the last line that $\omega_{k}+\omega_{l}\gg\Gamma_{k}+\Gamma_{l}$.

As defined above, the added noise covariance matrices $\mathrm{T}$, $\mathrm{B}$, and $\mathrm{M}$ are expressed in units of an equivalent thermal phonon occupation and represent the measurement uncertainty for any single estimate obtained from these filters. For a given measurement configuration, the optimal filters minimize this added noise covariance and provide estimates of the quadrature amplitudes with the least uncertainty.

The quadrature covariance of the initial multi-mode state can be inferred from an ensemble of repeated measurements, by inverting Eq. (30)

[TABLE]

assuming identical preparation of the initial state for each measurement. If the system and bath parameters are independently calibrated, then the bias matrices $\mathrm{T}$, $\mathrm{B}$, and $\mathrm{M}$ can be precisely calculated and subtracted to recover the inferred multi-mode state covariance. The statistical uncertainty of the inferred covariance is then limited by the uncertainty of the estimator covariance $\mathrm{\Sigma}$, which can be reduced by minimizing the total added noise covariance $\mathrm{T}+\mathrm{B}+\mathrm{M}$ for each sample, in addition to increasing the sample size ${n_{s}}$ (see Appendix A).

Fig. 2b demonstrates retrodiction of an oscillator initially prepared in a squeezed state $\ket{\zeta}=\hat{\mathcal{S}}(\zeta)\ket{0}$, defined by the single-mode squeezing operator Walls and Milburn (2008)

[TABLE]

Squeezing of the quadrature variance below the ground-state zero-point-motion cannot be directly observed with the phase-independent quadrature estimates considered in this work, since they are constrained by the SQL. However, quadrature squeezing can be inferred from the covariance of an ensemble of estimates, demonstrated in Fig. 2b, only after subtracting the added noise covariance matrices, which are numerically evaluated from Eqs. (31), (33), and (34) as described in Appendix B.

IV Standard quantum limit for retrodiction

Realizing optimal state retrodiction of a given ensemble of oscillators involves two separate choices. First, the intra-cavity intensity $\bar{n}$ must be chosen at the time of measurement, determining the measurement cooperativity $C_{i}$ for each oscillator. Second, for a given measurement strength, the optimal filter must be derived to obtain the best state estimates from each recorded trace, provided by Eq. (28). In this section, we consider retrodiction of the state of a single oscillator and determine the measurement conditions for achieving sensitivity at the SQL.

The OLS filters defined in Eq. (24) weight the filter estimate according to the coherent decay of the initial state. However, diffusion of the oscillator state after $t=0$ adds noise to the observed trajectories, illustrated in Fig. 1c, accumulating at a rate proportional to the measurement cooperativity $C_{1}$ and the thermal bath occupation $\nu_{1}$. An optimal filter should appropriately weight the relative signal to noise at each subsequent time $t$, implying that the optimal filter envelope must decay faster than the coherent response functions ${\bm{r}}_{1}(t)$.

The GLS filters derived from Eq. (28) can be used to calculate the optimal filter for any particular system parameters, however an analytic model yields further intuition about the optimal measurement cooperativity. For estimation of a single oscillator, the GLS filters are well-approximated by a set of exponentially damped sinusoidal filters

[TABLE]

parametrized by an arbitrary exponential decay rate $\gamma_{1}$, for sufficiently long observation intervals $\Gamma_{1},\gamma_{1}\gg{t_{f}}^{-1}$. In general, the total noise added to the quadrature estimates for any oscillator $i$ can be quantified by the average variance

[TABLE]

where $\mathrm{T}_{ii}$ refers to the $2\times 2$ diagonal block matrix corresponding to oscillator $i$. For estimation of a single oscillator with the exponential filters, this added noise occupation can be numerically evaluated for a given measurement cooperativity $C_{1}$ and filter decay rate $\gamma_{1}$, with results summarized in Fig. 3a.

The individual noise covariance matrices given by Eqs. (31), (33), and (34) can also be evaluated analytically and shown to be approximately proportional to the identity matrix, for filters that decay slower than the oscillation frequency $\gamma_{1}\ll\omega_{1}$ and assuming a measurement duration that is sufficiently long to capture the full transient response $\Gamma_{1},\gamma_{1}\gg{t_{f}}^{-1}$. The uncertainty contributed from each noise source, therefore, can be fully described by the added quadrature variance, expressed as equivalent added thermal occupation

[TABLE]

These three expressions sum to give the total added noise, as displayed in Fig. 3b-c.

For measurements performed at a particular cooperativity $C_{1}$, the uncertainty of estimates obtained from the recorded traces can be evaluated individually for each noise source as a function of the filter decay rate $\gamma_{1}$, as shown in Fig. 3b. As predicted by the Gauss-Markov theorem Rao (1976), the OLS estimators defined by Eq. (24), equivalently $\gamma_{1}=\Gamma_{1}$, only minimize the added variance due to the temporally uncorrelated measurement noise ${\Delta n}_{1}(\mathrm{M})$. The OLS optimization does not minimize the estimator variance for temporally correlated noise Rao (1976), driven by the thermal bath and quantum backaction, and the optimal filter decay rate is increased $\gamma_{1}>\Gamma_{1}$, because a shorter temporal filter envelope captures less of the accumulated diffusion at later times.

Minimizing the total added variance, the optimal exponential filter decay rate for a given measurement cooperativity $C_{1}$ is

[TABLE]

The added noise for this optimized exponential filter is plotted in Fig. 3c as a function of measurement cooperativity $C_{1}$. The optimal measurement condition for retrodiction is achieved in the limit of high cooperativity, corresponding to backaction-dominated diffusion of the oscillator during the measurement.

In this high-cooperativity regime, where diffusion from backaction exceeds that from the oscillator’s thermal motion (satisfying $C_{1}\gg\nu_{1}+1/2$), and where the effect of backaction in the measured signal exceeds that of measurement shot noise (satisfying $C_{1}\gg\epsilon^{-1/2}$), the optimal filter decay rate is approximately

[TABLE]

with corresponding minimized total added noise occupation

[TABLE]

The SQL for state retrodiction, therefore, is reached in this backaction-dominated regime, assuming ideal detection efficiency $\epsilon=1$. The result in Eq. (42) also defines the optimal filter envelope for an oscillator with negligible intrinsic damping $\Gamma_{1}\rightarrow 0$. In this limit, the optimal measurement sensitivity is realized at any finite measurement strength, provided an observation period much longer than the backaction diffusion timescale ${t_{f}}\gg 2\pi\kappa/(4\bar{n}g_{1}^{2})$.

In the strong-measurement limit $C_{1}\rightarrow\infty$, this best exponential filter evolves toward a delta function, which would describe a fully projective, instantaneous measurement of the oscillator’s position. However, as observed in Fig. 3a,c, the optimal measurement cooperativity is bounded from above by the oscillator’s quality factor $2\sqrt{\epsilon}C_{1}\ll\omega_{1}/\Gamma_{1}$. Beyond this bound, the filter envelope decays within an oscillation period, violating the approximation $\gamma_{1}\ll\omega_{1}$ above, and retrodiction no longer provides a phase-independent estimate of both oscillator quadratures.

It is noteworthy to contrast these imprecision limits for state retrodiction with those for the well-demonstrated limits for continuous displacement Teufel et al. (2009); Anetsberger et al. (2010); Rossi et al. (2018) and force Schreppler et al. (2014); Mason et al. (2019) detection. In each case, the optimal measurement is obtained with an equal imprecision added by measurement noise and quantum backaction. When an optomechanical oscillator is employed as a sensor for external forces or displacements, for instance from gravitational waves Abbott et al. (2016), then diffusion driven by quantum backaction increases the measurement imprecision at later times, and the optimal sensitivity on mechanical resonance is achieved with cooperativity $C_{1}=(2\sqrt{\epsilon})^{-1}$ Schreppler et al. (2014), typically of unity order.

However, for the case of retrodiction, the results described here indicate that the optimal sensitivity is reached in the high-cooperativity regime $C_{1}\gg\nu_{1}+1/2$, where information about the oscillator’s initial state is rapidly extracted. Measurement of the oscillator’s state by the cavity mode inherently results in backaction noise added to the oscillator. Nevertheless, it is preferable to increase the measurement cooperativity such that the measurement rate and the backaction-induced diffusion far exceed the loss of state information to the unmeasured modes of the oscillator’s thermal bath. The additional diffusive noise added to the oscillator’s trajectory during measurement is suppressed in the estimate, by using an appropriately short filter profile.

V Two-mode state estimation

The formalism developed in Sec. I and II is already fully general for estimation of multi-mode states of $N$ oscillators. Similar multi-mode estimation has been applied in experiments to obtain estimates for two-mode states Palomaki et al. (2013b); Spethmann et al. (2016); Kohler et al. (2018) and to demonstrate entanglement Palomaki et al. (2013a) from correlations observed in subsequent measurements of two modes.

As an example, consider simultaneous retrodiction of two oscillators prepared in a two-mode squeezed state (TMSS) $\ket{z}=\hat{\mathcal{S}}_{2}(z)\ket{0,0}$, generated from ground-state oscillators through the action of the two-mode squeezing operator Walls and Milburn (2008)

[TABLE]

After the squeezing interaction is turned off, estimates of each quadrature of the two-mode system are obtained by applying optimized filters, calculated from Eqs. (27) and (28), to the subsequently observed free transient decay.

The resulting $2N$-dimensional Gaussian phase space distribution defines an ellipsoid in phase space, fully characterized by its mean and covariance. The distribution can be visualized in terms of orthogonal 2-d projections, as shown in Fig. 4 for estimates obtained from simulated homodyne measurements of a TMSS. Projections onto the 2-d phase space for each individual oscillator, shown in Fig. 4a, are equivalent to tracing over the other oscillator’s state and reflect an effective thermal occupation $\langle\hat{a}^{\dagger}_{i}\hat{a}_{i}\rangle=\sinh^{2}|z|$.

The presence of two-mode squeezing is revealed by correlations between quadratures of the different oscillators, which are displayed in Fig. 4b. ‘In-phase’ correlations, between the positions of the two oscillators $\hat{X}_{1}(0)$ and $\hat{X}_{2}(0)$, are produced by the real part of $z$, while the imaginary part generates ‘out-of-phase’ correlations, between the position of one oscillator and the momentum of the other, as shown in Fig. 4b. The other two orthogonal projections of the 4-d phase space are qualitatively similar, but not shown. Once again, recovering the actual state covariance requires subtracting the covariances added by each noise source, which are numerically evaluated from independent knowledge of the system parameters.

V.1 Verification of estimated correlation

The two-mode covariance inferred from the matched-filter estimates can be experimentally validated through comparison to model-independent statistics of the recorded signals. For a multi-mode system, the mean-squared homodyne signal also reveals information about the initial correlations Kohler et al. (2018)

[TABLE]

again assuming $\omega_{i}+\omega_{j}\gg\Gamma_{i}+\Gamma_{j}$, where $f_{BW}$ is the bandwidth of the recorded signal. This expression reflects coherent time evolution of the mean and covariance of the initial quadrature amplitudes, in addition to relaxation given by

[TABLE]

to an equilibrium signal variance determined by the thermal baths and measurement backaction.

Correlations between oscillators appear in this signal as transient beatnotes at the sum and difference frequencies, unperturbed by the thermal, backaction, and measurement noise that biases the filter covariances. This signal, therefore, serves as an independent statistic for comparison to the matched-filter estimates, through reconstructing the predicted beatnote for the inferred two-mode state, shown in Fig. 4c. This comparison indicates how faithfully the filter model matches the system dynamics and serves as a check for calibration of the added noise covariance that must be subtracted from the filter estimates.

The ‘out-of-phase’ correlations for the TMSS simulated in Fig. 4 are characteristic of those generated by the negative-mass instability observed in Ref. Kohler et al. (2018), produced by resonant coupling between positive- and negative-mass oscillators. For an effective negative-mass oscillator, the coherent state evolution corresponds to an opposite rotation in phase space, corresponding to $\omega_{i}\rightarrow-\omega_{i}$ in Eq. (9) and Eq. (12). By consequence, for the TMSS represented in Fig. 4c, the 2nd-order coherence $\langle\hat{a}_{1}\hat{a}_{2}\rangle$ evolves at the frequency difference, with amplitude and phase directly reflecting the magnitude and phase of correlations between the two oscillators.

VI Limits of multi-mode estimate precision

For a given multi-mode measurement record, with independently calibrated system and noise parameters, the GLS method facilitates numerical calculation of the optimal filter to recover the multi-mode state estimates. For a single oscillator, the optimal estimate imprecision approaches the SQL in the limit of high cooperativity, as demonstrated in Sec. IV. However, diffusive motion of each oscillator in the measurement record introduces additional imprecision to the estimate results. Here we explore the optimal two-mode measurement strength and additional limits to the estimate precision due to the presence of a second oscillator

Consider state retrodiction for simultaneous observation of two oscillators that differ only in their resonance frequency, with frequency separation $\delta=\omega_{2}-\omega_{1}$. The distinguishability of the two oscillator responses is parametrized by their spectral resolution $\delta/\Gamma_{1}$ and the measurement strength by the cooperativity $C_{1}=C_{2}$.

The total estimate uncertainty for the first oscillator ${\Delta n}_{1}$, defined by Eq. (39), was numerically computed as a function of measurement cooperativity $C_{1}$ for a few different oscillator detunings, assuming perfect detection efficiency $\epsilon=1$, with results shown in Fig. 5a. For sufficiently low-cooperativity measurements, the single-oscillator exponential filters defined by Eqs. (38) and (41) achieve optimal results, with the measurement imprecision decreasing with increasing measurement strength $C_{1}$. When the optimal single-oscillator filter decay rate ${\gamma_{\text{opt}}}$ approaches the oscillators’ frequency difference $\delta$, the simple exponential filters fail to optimally distinguish the response of each oscillator, resulting in additional estimate imprecision.

The GLS filters defined by Eq. (28), achieve reduced imprecision at high cooperativity. These filters optimally distinguish the motion of each oscillator by acquiring a notch in the filter spectrum, shown in Fig. 5b, that suppresses signal components in the frequency band of the other oscillator. However, for larger cooperativity, the back-action broadened filters for each oscillator become increasingly indistinguishable. The filter normalization matrix $\mathrm{J}$, defined by Eq. (23), becomes nearly singular, resulting in growth of the normalized estimate imprecision.

The optimal measurement cooperativity for two-mode retrodiction at a given detuning $\delta$ was found by numerically minimizing the total estimate imprecision ${\Delta n}_{i}$, summarized in Fig. 5c. When $\Gamma_{1}\ll\delta\ll\omega_{1}$, the optimal cooperativity is approximated by ${C_{\text{opt}}}=\delta/2\Gamma_{1}$, the threshold where the corresponding single-oscillator exponential filter linewidth defined by Eq. (42) matches the oscillator detuning $\delta$.

At the optimal measurement cooperativity, the minimum estimate imprecision for each oscillator, shown in Fig. 5d, approaches the SQL in the limit of well-resolved oscillators $\delta/\Gamma_{i}\gg 1$. When the oscillators are not well resolved $\delta\lesssim\Gamma_{1}$, the filters cannot distinguish the motion of each oscillator. The homodyne signal is always sensitive to the center of mass motion of the oscillators, but contains negligible information about their relative motion, resulting in strong correlated errors in the individual oscillator estimates, as displayed in Fig. 5d, described by the off-diagonal blocks of the noise covariance matrices

[TABLE]

When $\delta\gtrsim\omega_{i}$, the optimal measurement cooperativity ${C_{\text{opt}}}$ is once more limited by the single-oscillator quality factor, as seen in Fig. 2c.

VII Experimental demonstration

We previously employed this matched-filter analysis for retrodicting correlated states of a multi-mode optodynamical system, reported in Ref. Kohler et al. (2018). In that work, correlations were generated through a negative-mass instability from collective coupling between the motion and spin of an atomic ensemble, resulting in resonant pair-creation analogous to a parametric amplifier.

Real-world measurements invariably bring additional complications, such as variations in system parameters and additional incoherent dynamics. These must all be included in the preceding statistical analysis to accurately estimate the state covariance. In particular, for Ref. Kohler et al. (2018), inhomogeneous broadening of the observed average homodyne PSD, due to shot-to-shot fluctuations of the oscillator frequencies $\omega_{i}$, caused dephasing of the ensemble-averaged signals and complicates the matched filter analysis.

To include this inhomogeneous broadening in the filter statistics, the oscillator frequency $\omega_{i}$ was treated as a classical random variable with variance $\sigma_{i}$. The quadrature response functions ${\bm{r}}_{i}(t)$ are then also stochastic, with mean trajectory

[TABLE]

which includes an additional decay envelope with dephasing rate $\sigma_{i}$, assuming a Gaussian frequency distribution. The intrinsic damping rate $\Gamma_{i}$ and inhomogeneous dephasing rate $\sigma_{i}$ were experimentally calibrated by fitting the observed homodyne PSD with a Voigt profile.

This inhomogeneous broadening of the oscillator response motivates definition of modified quadrature filters

[TABLE]

which are spectrally broadened to match the average PSD, in order to better capture the oscillator response across the distribution of frequencies.

The filter normalization matrix $\mathrm{J}$, defined in Eq. (23), must also be considered as a matrix-valued random variable. Assuming here that fluctuations of the system parameters are uncorrelated with the initial quadrature operators, its expectation can be independently computed in Eq. (17) as

[TABLE]

Individual quadrature estimates from data reported in Ref. Kohler et al. (2018) are displayed in Fig. 6, obtained using filters defined by Eq. (49) and normalized by $\langle\mathrm{J}\rangle^{-1}$.

Estimation of the quadrature covariance is further complicated. The 2nd-order moment of any pair of the unnormalized filter estimates, defined by Eq. (16), can be written

[TABLE]

Here, the indices $i$, $j$, $k$, and $l$ run over all $2N$ quadrature elements, unlike the block-matrix notation used above. The covariances added from each noise source in Eq. (51) are given by the matrices

[TABLE]

These expressions can be evaluated numerically as described in Appendix B, based on independent calibrations of the system parameters. The system of $(2N)^{2}$ equations defined by Eq. (51) can then be inverted to recover the 2nd-moments of the quadrature estimates and the added noise covariances, which define the covariance ellipses shown in Fig. 6a-b.

Finally, incoherent coupling in the experimental system between the spin and motion of individual atoms also resulted in additional thermal diffusion of the collective spin, which was included as an additional noise drive in Eq. (14). Diffusion from this interaction increased with $\bar{n}$, proportional to the measurement strength, and therefore limited the experimental system to low-cooperativity measurements, which prevented achieving retrodiction at the SQL.

VIII Conclusion

In summary, we have derived a detailed model for linear state estimation from optical measurements of multi-mode optomechanical systems. We have demonstrated that retrodiction of the past state of a single oscillator from high-cooperativity measurements approaches the SQL, when estimated with a matched filter that optimally suppresses the accumulated backaction noise. Furthermore, we have described a general numerical method to derive optimal filters, which directly facilitates optimal estimation from simultaneous measurements of multiple oscillators. We explored the optimal measurement strength for retrodiction of a multi-mode state and demonstrate additional constraints on the oscillator detuning in order to reach estimate imprecision at the SQL. This work lays an experimentally motivated framework for simultaneous measurements of multi-mode systems, which provides an essential tool toward efficient measurements of many-body systems.

Appendix A Convergence of sample covariance matrix

Assuming the noise covariance matrices can be calculated to arbitrary precision, from independent calibration of the system and bath parameters, then the statistical uncertainty of the state covariance matrix inferred from Eq. (30) is determined solely by the uncertainty of the estimator covariance $\mathrm{\Sigma}$.

This covariance matrix can be estimated from an ensemble of ${n_{s}}$ identical measurements $\{\breve{{\bm{Q}}}[i]\}$ according to the sample covariance

[TABLE]

where $\breve{{\bm{\mu}}}=(1/{n_{s}})\sum_{i}\breve{{\bm{Q}}}[i]$ is the sample mean. The covariance estimator $\breve{\mathrm{\Sigma}}$ is itself a random variable, which, assuming Gaussian initial states and input noise, is distributed according to the $2N$-dimensional central Wishart distribution Chatfield and Collins (1980)

[TABLE]

with $({n_{s}}-1)$ degrees of freedom. This distribution is the multi-dimensional generalization of the $\chi^{2}$ distribution, and allows calculation of estimate uncertainties from the variance of the matrix elements

[TABLE]

The uncertainty for estimating the multi-mode state covariance $\breve{\Sigma}$ from an ensemble of measurements, therefore, is reduced by minimizing the total systematic noise bias $\mathrm{T}+\mathrm{B}+\mathrm{M}$ added to the retrodicted estimate covariance $\mathrm{\Sigma}$ in Eq. (30).

Appendix B Evaluation of estimate covariance from diffusive noise

The covariance matrices of estimator noise added from thermal- and backaction-driven diffusion, defined in Eq. (33) and Eq. (34), respectively, contain a triple integral over rapidly oscillatory integrands. To aid in numerical evaluation, the integral which appears in each of these equations, describing the estimate covariance for correlated diffusion of oscillator $k$ and $l$ driven by a common bath, can be rewritten as

[TABLE]

in terms of the generic oscillator two-time correlation for response to a common noise bath defined in Eq. (46).

Generalizing to a complex response function

[TABLE]

in terms of which ${\bm{r}}_{i}(t)=\begin{pmatrix}\operatorname{Re}[\rho_{i}(t)],&-\operatorname{Im}[\rho_{i}(t)]\end{pmatrix}^{\mathsf{T}}$ then Eq. (46) can be rewritten as

[TABLE]

This integral can be evaluated analytically, giving

[TABLE]

in terms of the complex oscillator responses. This result can be used for efficient numerical evaluation of the double integral in Eq. (59) for discretely sampled signals and filter functions.

When considering fluctuations of the oscillator frequencies in Sec. VII, it is necessary to evaluate the expectation of this two-time response product. Assuming the variances of the oscillator frequencies are small compared to their difference, $\sigma_{k}+\sigma_{l}\ll|\omega_{k}-\omega_{l}|$ for $k\neq l$, this expectation can be approximated by

[TABLE]

Acknowledgements.

The authors acknowledge helpful discussions with Mohan Sarovar, Zhishen Huang, and Klemmens Hammerer about related work on quantum state estimation, retrodiction, and smoothing. Thanks also to Johannes Zeiher for feedback on the manuscript. This work was supported by the Air Force Office of Scientific Research, through grants FA9550-14-1-0257 and FA9550-19-1-0328. J.K. was supported by the U.S. Department of Defense through the National Defense Science and Engineering Graduate Fellowship program, and J.G. and E.D. by the National Science Foundation Graduate Fellowship.