Quantum estimation of coupling strengths in driven-dissipative optomechanics

We exploit local quantum estimation theory to investigate the measurement of linear and quadratic coupling strengths in a driven-dissipative optomechanical system. For experimentally realistic values of the model parameters, we find that the linear coupling strength is considerably easier to estimate than the quadratic one. Our analysis also reveals that the majority of information about these parameters is encoded in the reduced state of the mechanical element, and that the best estimation strategy for both coupling parameters is well approximated by a direct measurement of the mechanical position quadrature. Interestingly, we also show that temperature does not always have a detrimental effect on the estimation precision, and that the effects of temperature are more pronounced in the case of the quadratic coupling parameter.


I. INTRODUCTION
Quantum optomechanics focuses on the interaction between the electromagnetic radiation and motional degrees of freedom of mechanical oscillators [1][2][3]. The simplest optomechanical system consists of a single cavity mode interacting with a single mechanical mode and is realised, for example, in an optical cavity with a movable mirror. In this case the mechanism responsible for the interaction is radiation pressure, which entails momentum exchange between light and matter. The presence of a cavity boosts the otherwise weak radiation pressure force, enhancing the light-matter interaction. The quantum effects of radiation pressure forces and the associated limits they set on the precision of mirrordisplacement measurements are of great importance for many applications including gravitational wave detectors, scanning probe microscopy and force sensing [1][2][3][4].
Although the radiation pressure interaction is intrinsically non-linear [5], approximate models are usually used [1,2] which assume a linear dependence of the cavity frequency ω(X b ) on the dimensionless position of the movable mirror, X b . So far these "linear" models have proved extremely successful, aided by the fact that the bare (or 'singlephoton') optomechanical coupling strength is usually very small [1,2]. However, researchers are continuously exploring ways of enhancing the optomechanical coupling, as well as the potential of optomechanics for ultra-high-accuracy applications such as Planck physics [6]. Hence, extensions to the linear model are becoming a necessity. The next step beyond the linear approach is to expand the cavity frequency up to and including second order in X b , leading to what we will call the "quadratic model" [7]. Accurate knowledge of all the relevant optomechanical coupling parameters will indeed be crucial for virtually any application of these systems. With such motivation in mind, this paper exploits local quantum estimation theory [8] (QET) to investigate how precisely the linear and quadratic coupling strengths may be measured in a model quantum optomechanical system. In a nutshell, QET looks for the best strategy for estimating unknown parameters encoded in the density matrix of a quantum system (i.e. a quantum statistical model ) [9,10]. The ultimate limits to the precision with which the desired parameters can be estimated may be quantified via the quantum Cramér-Rao bounds [8,11,12].
Within such framework, we consider an optomechanical set-up featuring a driven (and lossy) cavity, whose dynamics is described by a Lindblad master equation. As is typical in recent optomechanics experiments, we assume sufficiently strong driving to approximate the dynamics via a master equation that is bilinear in the canonical operators. This leads to a Gaussian steady state [13], whose first and second moments we characterise via a combination of matrix algebra and numerical methods. This state, with its explicit dependence on all the model parameters (in particular the unknown coupling strengths), will embody our quantum statistical model. This in turn may be attacked via general closed-form expressions that are available for QET in Gaussian models [14]. The visual representation of the parameter estimation methodology for driven-dissipative optomechanics is shown in Fig. 1.
Using this approach, and for typical values of the model parameters (inspired by recent experiments), we found that it is significantly easier to estimate FIG. 1: Schematic of the parameter estimation methodology for driven-dissipative optomechanics. We consider a driven-dissipative optomechanical system featuring a driven (by an external laser) and lossy (photons escaping the cavity) cavity and a damped mechanical oscillator. The mechanical support has low but finite temperature (leading to a non-zero thermal occupation number). The optomechanical coupling arises due to the radiation pressure on the movable mirror. Once the system has reached a steady state we measure an observable. We repeat the measurement many times to get the statistics. Finally, we process the data to find the best guess for the coupling parameters of interest. the linear coupling strength than the quadratic one. This result is indeed to be expected in typical experimental conditions, given the significant difference in the relative strength of the two constants. Our analysis also reveals that the majority of information about the parameters is encoded in the reduced state of the mechanical element.
We then investigate how well some specific measurements perform when compared to the fundamental limits imposed by QET. We in particular focus on measurements of the mechanical position X b , field amplitude Q, mechanical momentum P b and field phase P, which can all be analyzed within the Gaussian formalism [15]. Among those we find that the best strategy to estimating the coupling parameters is through a direct measurement of the mechanical position X b .
We additionally explored the influence of temperature on the estimation precision of the coupling strengths. In the case of the linear coupling parameter we found that the effect of temperature is mostly significant at lower intracavity photon num-bers (i.e. lower driving strengths), where it improves the estimation precision. At higher intracavity photon numbers, the zero temperature scenario predicts a better estimation precision instead. Interestingly, in the case of the quadratic coupling parameter we found that a hotter mechanical bath gives a higher estimation precision at all intracavity photon numbers in our range.
We note that the application of QET in closely related optomechanical set-ups was previously considered in works by Bernad, Sanavio and Xuereb [10,16], who adopted the linear model of optomechanics. In Ref. [10], purely Hamiltonian (non-dissipative) dynamics was assumed, and it was found that larger intracavity photon numbers would facilitate the estimation of the linear coupling strength. Analogously, our results to follow show that a similar conclusion holds in the strong driving regime. However, we also find that for weaker drive strengths the picture is more complicated when considering finitetemperature effects. We note that very recently the same authors went on to consider a driven-damped system [16] though using a somewhat different approach to ours and without considering quadratic couplings. In particular, Ref. [16] neglects the contribution of the steady state's first moments (i.e. the averages of the canonical operators) to the quantum Fisher information (QFI). We checked that this is a well-justified assumption for the model parameters adopted therein. However, our results show clearly that there are experimentally accessible parameter regimes where the picture changes dramatically and the first moments can come to dominate the QFI. Finally, Schneiter et al. [17] have applied local QET to a time-dependent, purely Hamiltonian and quadratic optomechanical system. Such a model, however, is sufficiently different from ours to prevent a simple and direct comparison between the two. This paper is organised as follows. In Sec. II we introduce our model of driven-dissipative optomechanics, including both linear and quadratic coupling terms, and outline how the dynamics may be approximated via a bilinear master equation in the limit of strong driving. In Sec. III we calculate the steady state of the system, which is Gaussian within the considered approximations, and hence is fully characterized by its first and second moments. We then develop the necessary QET tools to investigate the optimal estimation of linear and quadratic coupling parameters. In Sec. IV we present and discuss the findings of our research. Finally, in Sec. V we summarise our results.

II. MODEL
We consider a simple optomechanical system consisting of two quantum harmonic oscillators, describing a single-mode cavity field and a single mechanical mode, respectively. The two modes are coupled non-linearly via radiation pressure [1][2][3]. We thus assume that our system is described by the Hamiltonian withQ andP the amplitude and phase quadratures for the cavity mode,X b andP b the dimensionless position and momentum operators of the movable mirror, whilst m and ω m are the effective mass and frequency of the mechanics [5,7,10,18]. The only nontrivial commutators read [Q,P] = [X b ,P b ] = i. The optomechanical coupling arises from a parametric dependence of the cavity frequency on the mechanical position ω(X b ).
In the widely explored "linear" regime of optomechanics, the mechanical motion is assumed to be very small and ω(X b ) is approximated by an expansion to linear order inX b . For the "quadratic model" of optomechanics, instead, terms up to and includinĝ X 2 b are retained: where ω 0 is the bare cavity frequency [1,2,7]. The strength of the optomechanical interaction can be quantified with the linear and quadratic coupling strengths, which for a generic set-up are defined as respectively [2]. Note that we can always ensure that g 1 is positive by a redefinition of the positive direction of X b , and that the linear model is recovered by setting g 2 to zero. A purely Hamiltonian description of the system is however not sufficient for our purposes, since we aim to describe a (more realistic) driven-dissipative optomechanical system featuring a driven and lossy cavity, and a damped mechanical oscillator. In order to conveniently introduce coherent driving in the model, we shall move to a frame rotating at the frequency of the driving laser, ω L . In this frame the Hamiltonian of the driven system may be written as where ∆ 0 = ω 0 − ω L is the detuning between the cavity and driving laser and E is the drive amplitude. In our model we will include cavity decay at a rate κ and mechanical damping at a rate Γ m , assuming that the thermal occupation of the cavity mode is negligible. We assume that the corresponding master equation describing the dynamics of the system is of the general Lindblad form [19,20]: where we defined the vector of quadrature operatorŝ while γ is the damping matrix: n m = 1/ e ωm/k B T − 1 is the mean occupancy of the mechanical oscillator, k B is the Boltzmann constant and T is the temperature of the mechanical reservoir [1,2]. We note that, in choosing a Lindblad form, we automatically excluded the use of the standard Brownian motion master equation (SB-MME) [20] to describe mechanical damping. Indeed, a Lindblad form greatly simplifies our analysis, since it avoids non-positivity issues that are known to occur in the SBMME [21]. The main effect of the drive is to displace the steady states of both the cavity field and the mechanical position [2]. We assume that the cavity is driven sufficiently strongly (and that the optomechanical couplings are weak enough) so that the system dynamics can be approximated via a bilinear master equation description, where only small quantum fluctuations around the semi-classical steady state are considered [2,16]. In detail, we start by displacing our canonical operators as perR →R + R 0 , where is the vector of steady-state quadrature averages.
Here, x 0 and p 0 are the average position and momentum of the mechanics in the steady state, while Q 0 and P 0 are the steady state displacements of the amplitude and phase quadratures, respectively [2]. Of course, the steady-state expectation values of the transformed operators will now vanish. This results in the following equations for the steady-state values of the system's first moments: where The non-linearity of Eqs. (10-13) suggests that multiple steady state solutions are possible [1,2]. This is known as dynamical multistability [1,2]. In detail, depending on the driving strength up to five (quadratic model) or three (linear model) different steady states solutions can exist. In this work, we only focus on parameters regimes where the system is stable, i.e., where a unique real solution to Eqs. (10)(11)(12)(13) exists. This in turn places an upper bound to the drive strength |E| 2 [10].
After the displacement has been implemented, we neglect terms that are beyond quadratic order in the transformed canonical operators [2]. The corresponding master equation readṡ where we note that the Lindblad operators remain unchanged, while the Hamiltonian now takes the bilinear form with ω eff = ω m + 2g 2 |α| 2 the effective mechanical frequency and g eff = − √ 2g 1 + 2g 2 x 0 the effective coupling strength. Note that the assumption of strong cavity driving translates into the condition that the intracavity photon number is large, i.e., |α| 2 ≡ (Q 2 0 + P 2 0 )/2 1.

A. Covariance Matrix Formalism
Due to its bilinear form, the master equation (14) admits a Gaussian steady state that can, in general, be fully characterised by its first and second moments of the quadrature operatorsR [13]. After having determined our steady state, we will be able to exploit general closed-form expressions that are available for QET in Gaussian models [14].
As anticipated in the previous section, the first moments of our Gaussian steady state are given by R 0 = (Q 0 , P 0 , x 0 , p 0 ), and are found by solving Eqs. (10-13) -recall also that we will only consider parameter regimes in which such solution is unique. The second moments are instead encoded in the steady state covariance matrixσ, which in our displaced frame of reference is given by [13] where {Â,B} ≡ÂB +BÂ is the anticommutator. As detailed in Appendix A, master equation (14) implies the following Lyapunov equation for the steady state covariance matrix: where We note that Eq. (17) can be solved analytically in terms of the model parameters and the vector of averages R 0 . The latter, however, may in general not admit an analytical expression in terms of the model parameters, as we recall it is the solution to the nonlinear system of equations (10)- (13). In the next section we shall show how to develop a comprehensive QET analysis of the coupling parameters solely from the knowledge of the first and second moments of our Gaussian steady state.

B. Quantum Estimation Theory for Gaussian States
The aim of quantum estimation theory (QET) is to identify the best strategy for estimating one or more parameters encoded in the density matrix of a quantum system [8][9][10]. Here we focus on local QET, which seeks a strategy that maximises the Fisher information over all possible measurements, and implicitly assumes that a rough estimate of the parameter value is known in advance [8].
In our model of driven-dissipative optomechanics, the parameters to be estimated shall be the coupling strengths g 1 and g 2 . As anticipated, all of the information about these parameters will be contained in the steady-state averages, R 0 , as well as in the steady state covariance matrix,σ. Specifically, for our coupling parameters (g 1 , g 2 ) the elements of the quantum Fisher information matrix (QFIM) are given by where Lσ(A) =σAσ, L W (A) = WAW. Note also that the term (4Lσ + L W ) −1 refers to the pseudoinverse if the term inside the bracket is singular [14,22]. The first term is the contribution due to the averages, while the second term is the contribution due to the variances and covariances towards the total QFI [22]. This terminology will be convenient later on as we seek to unravel how the different terms contribute across different parameter regimes. We note, however, that our terminology only describes the origin of the dependence of the gradients with respect to the coupling parameters. Hence whilst the first term in eqn. (24) only contains gradients of the averages with respect to the coupling constants, and we will therefore call it the contribution of the averages, it does also depend on the covariance matrix.
The ultimate limit to parameter estimation in this context is set by the QCRB [8,11,12]. In multiparameter estimation theory both coupling parameters are assumed to be unknown (or only known with low precision). In this case, the QCRB relates the covariance matrix of any pair of unbiased estimators for the parameters (g 1 , g 2 ) to the QFIM. For M experimental runs, the corresponding QCRB reads [12,14] Cov (g 1 , The limiting case of single-parameter estimation theory can be reached if we assume that only one parameter is unknown, say g i . In this case the QCRB relates the variance Var (g i ) of an unbiased estimator of the parameter g i to the corresponding diagonal element of the QFIM. For M experimental runs, the corresponding bound reads [14,15,[22][23][24][25] Var In other words, the diagonal elements of the QFI matrix quantify the "best-case-scenario" performance for the estimation of each individual parameter.
Hence, in what follows we shall pick Eq. (26) as our benchmark in evaluating the performance of various measurements (see below). Note that in single-parameter estimation theory the saturation of the QCRB is guaranteed, at least in the limit M → ∞, and assuming that every mathematically allowed quantum measurement can be implemented [12,14,22]. This, however, is not true for the estimation of multiple parameters: in this case the optimal measurements for different parameters may not be compatible [12]. While the QFI quantifies the ultimate quantum limit to parameter estimation [8,12,26], the estimation performance of specific measurement strategies may be quantified via the classical Fisher information (FI) matrix [8,27]. In our context, the FI measures the amount of information that a classical random variable s (the outcome of a quantum measurement) contains about the parameters g 1 , g 2 [28]. The FI matrix elements take the form where p g1,g2 (s) is the probability distribution of the measurement outcome s, assumed to be a smooth function of (g 1 , g 2 ) [14,28]. Depending on the chosen observables, analytical solutions to the integral in Eq. (27) may exist. This is particularly true for quadrature measurements (i.e. a measurement of Q,P,X b ,P b or a linear combination thereof), provided that the measured state is Gaussian [15].
In the case of optomechanics, it is well known that one can use a homodyne detection scheme to measure the light quadraturesQ,P [15,29]. However, we shall also consider a direct measurement of the mechanical quadratures,X b andP b , for completeness. In practice this could potentially be achieved using e.g. another optical mode of the cavity. In this scenario, the probability distribution associated with a measurement ofŝ ∈ {Q,P,X b ,P b } has the following expression [15]: where s 0 (g 1 , g 2 ) is the steady state average of the chosen quadrature, appropriately chosen from the set {Q 0 , P 0 , x 0 , p 0 }, whileσ kk (g 1 , g 2 ) is the corresponding diagonal element of the steady state covariance matrix (σ 11 forŝ =Q,σ 22 forŝ =P and so on). In this setting an analytical solution to the integral in Eq. (27) exists and is given by Note that the choice of a strategy to estimate the parameters (g 1 , g 2 ) is optimal if the FI and QFI matrices are equal, i.e. J = I. As anticipated, we shall focus solely on the diagonal elements of the QFIM (Eq. (24)): I 1,1 and I 2,2 , respectively. As noted above, the diagonal elements are indeed the "most optimistic" quantifiers of estimation precision of the coupling strengths. In general, however, the combined precision of the two parameter estimations will be worse than what the diagonal elements suggest.
Both, the definitions of the QFI (Eq. (24)) and the FI (Eq. (29)) rely on the derivatives of the steady state covariance matrix and the averages with respect to the coupling strengths. Since we have seen that bothσ and R 0 are determined by the nonlinear system of equations (10)- (13), which in general can only be solved numerically, we use implicit differentiation to calculate the derivatives in question. This allows us to express all our quantities of interest in terms of the numerical solution to the above nonlinear equations, and allows us to avoid numerical differentiation altogether.

IV. RESULTS
For simplicity we consider a specific geometry, that of a Fabry-Perot cavity [1], in which one mirror is fixed and the other is mounted on the mechanical oscillator. In this case, assuming an ideal onedimensional cavity field, the cavity frequency takes the specific form with L the bare cavity length and x zp = /2mω m the ground state position uncertainty of the mechanical oscillator [2]. Making the standard assumption that the mechanical motion is very small on the scale of the cavity lengthX b /L 1, for the quadratic model of optomechanics ω(X b ) is approximated as where g 1 = ω 0 x zp /L and g 2 = 2ω 0 x 2 zp /L 2 are the linear and quadratic coupling strengths in accordance with Eqs. (3) and (4), respectively [1,2,7].
Here, we examine three scenarios: zero temperature (T = 0 K), low temperature (T = 1 mK) and "high" temperature (T = 80 mK) scenario. In each case we are looking for the best strategy to estimate the linear, g 1 , and quadratic, g 2 , coupling strengths. First, we establish the fundamental quantum limits on the estimation precision, which, in accordance with the quantum Cramér-Rao bound (QCRB), are quantified with the "global" QFIs (i.e. the QFIs calculated from the bipartite state of light plus mechanics). Additionally, by tracing out the mechanical (light) mode, we can calculate "local" QFIs that are relevant when only the light (mechanical) mode is directly measurable. Comparing these local QFIs with the global ones will also reveal how much information about the coupling parameters is contained in the reduced states of light and mechanics. Finally, we compare the QFI limits to the performance of a small selection of "realistic" measurements (quantified with the respective FI), including those ofQ,P, X b andP b . This can help us discern which of the experimentally common measurements constitute the best strategy to parameter estimation in each scenario. In many cases we have chosen to measure the FIG. 2: (a) Log-log plot of the relative error bound on g 1 (as implied by QFI) against the intracavity photon number, |α| 2 as predicted by the linear and quadratic models in the zero temperature (green circles and orange dashed line, respectively), low temperature (red triangles and black line, respectively) and high temperature (blue squares and purple dot-dashed line, respectively) scenarios. (b) Log-log plot of the relative error bound on g 2 (as implied by QFI) against the intracavity photon number, |α| 2 as predicted by the quadratic model in the zero temperature (purple dot-dashed line), low temperature (black line) and high temperature (orange dashed line) scenarios. estimation precision with the relative error Our choice of parameter values is motivated by recent experiments where the ground state of a mechanical oscillator was approached via back-action cooling arising from a red-detuned laser drive, specifically [31]. Correspondingly, we adopt the following parameter values ω m = 1.1×10 7 Hz, m = 4.8×10 −14 kg, Γ m = 32 Hz, ∆ 0 = ω m , κ = 10 5 Hz, g 1 = 2 × 10 2 Hz and g 2 = 2g 2 1 /ω 0 ≈ 1.1 × 10 −5 Hz [31]. In order to ensure that the driving is strong enough for the Gaussian approximation to hold and we do not encounter any stability issues we consider a region 10 8 ≤ E ≤ 3.8 × 10 9 Hz in all three scenarios. In terms of the intracavity photon number, |α| 2 , this corresponds to a region: 80 |α| 2 1.2 × 10 5 (or 1.9 log 10 (| α | 2 ) 5.1).
In Fig. 2(a) we investigate the effects of the higher order g 2 term, temperature and driving on the estimation precision of the linear coupling strength, g 1 . Clearly, in all scenarios the effect of the higher order corrections due to the g 2 term is minimal: the linear and quadratic models show a very good agreement at all |α| 2 . This is to be expected as g 2 is several orders of magnitude below g 1 in our example. As an interesting aside, in a membrane-in-the-middle optomechanical system it is possible to engineer a purely quadratic coupling via the position of the membrane, in which case g 2 would clearly play the key role. [2,30]. Figure 2(a) reveals a surprisingly complex dependence on temperature. There is a crossover around log 10 (|α| 2 ) ∼ 4.6 (or |α| 2 ∼ 4 × 10 4 ): below this value the high temperature scenario offers the best precision for estimating g 1 , but above it the best precision is found at lower temperatures. As discussed in Appendix B (see in particular Fig. 5), much of this behavior can be understood by looking at the relative contributions of the variances and averages to the QFI and how these change with temperature. The contribution of the averages to the QFI always increases monotonically with the intracavity photon number and hence it always eventually dominates. This, taken together with the fact that the contribution of the averages to the QFI is reduced by in- creasing the temperature of the mechanical reservoir, means that the best precision is eventually expected at zero-temperature. However, for non-zero temperatures the contribution to the QFI from the variances is important and it is dominant at sufficiently low intracavity photon numbers. This is not surprising given the strong cooling effect that the driven cavity can have on the mechanics [1,2], leading to a strong dependence of the corresponding variances on the coupling strength (and the intracavity photon number). The impact of the cooling effect gets stronger at higher temperatures. In contrast, at zero temperature there is a very weak dependence of the variances on g 1 , so that in that case the averages always dominate. Away from the zero-temperature limit, the contribution of the variances to the QFI develops a peak at a particular intra-cavity photon number, reflected as the maxima in the relative error for g 1 seen at finite temperatures in Fig. 2(a).
In Fig. 2(b) we explore the effect of cavity driving and temperature on the estimation precision of the quadratic coupling strength, g 2 . Overall the relative estimation errors are much higher for g 2 than g 1 , which is unsurprising as the former is several orders of magnitude smaller.
Interestingly, at all |α| within our allowed range the high temperature scenario predicts the lowest relative errors bound on g 2 : a hotter mechanical bath gives a better estimation precision for all driving strengths below the instability threshold. This can be traced back to the fact that, also in the estimation of g 2 , the information content of the variances is again higher than that of the averages at lower driving strengths, and in the high temperature case it remains so for all driving strengths up until the instability threshold. The effect gets weaker as the temperature is reduced, and in the T = 0 case a crossover is seen with the contribution of the averages eventually becoming dominant (see Fig. 6 in Appendix B). The overall result is that the relative error bounds on g 2 decrease monotonically with increasing drive across the parameter range studied.
In Fig. 3 we compare global and local QFIs for g 1 (Fig. 3(a)) and g 2 (Fig. 3(b)). In a nutshell, we find that the majority of information about the coupling parameters is contained in the reduced state of the mechanics. Note that, in standard optomechanical experiments, measurements are typically performed on the light mode. Nevertheless, our results suggest that significantly more information about the couplings might be available by probing the mechanical motion more directly.
For g 1 the uncertainties found from either just the mechanical subsystem, or just the optical subsystem only drop monotonically with drive strength at T = 0, matching what happens with the full system. In this case, even at 80mK the uncertainties obtained from the reduced state of the mechanics almost match those from the full system. In contrast, for g 2 the uncertainty obtained from the state of just the mechanics only reaches that achieved with the full state in the zero temperature limit.
Finally, in Fig. 4 we show how some realistic measurements perform in comparison to the ultimate limits given by the QFI. The figure shows that, out of the measurements considered, the mechanical po-sition almost always does best at estimating the coupling parameters. The ultimate limits to estimation precision of the coupling strengths can only be approached at low and intermediate drive strengths via measurement ofX b at zero temperature. For higher temperatures this limit is approached for g 1 at very high intracavity photon numbers, whilst for g 2 it is never achieved.

V. CONCLUSIONS
We employed local QET to the problem of estimating linear and quadratic coupling parameters in driven-dissipative optomechanics. For experimentally realistic values of the model parameters, inspired by Ref. [31], we have found that it is considerably easier to estimate the linear coupling strength than the quadratic one. Our analysis has also showed that the best strategy for estimating the coupling parameters can be well approximated by a direct measurement of the mechanical positionX b .
Exploring the effect of temperature on the estimation precision of the coupling strengths, we found that higher temperatures are not always detrimental to the estimation performance. The effect of temperature is particularly striking when analyzing the estimation of the quadratic coupling parameter: in this case we found that a hotter mechanical bath (T = 80mK) resulted in a higher estimation precision for all drive strengths below the instability threshold. In contrast, in the case of the linear coupling strength the effect of temperature is most significant at lower driving. Past a certain drive strength, better estimation precision for the linear coupling parameter is instead achieved at lower temperatures. is the dissipator [20]. Using the same convention, we find that the vector of first moments evolves according to [13,15]: where γ is a matrix with elements γ ij that has been conveniently decomposed into its symmetric and antisymmetric parts: Instead, the covariance matrix obeys the following equation of motioṅ The steady state covariance matrix,σ, can thus be found by solving the Lyapunov equation where In order for the steady state covariance matrix to describe a physical state ρ(t) it must be a real, symmetric and a positive semi-definite matrix [13,15]. The semi-positivity requirement is satisfied provided that B is stable, i.e. the real parts of the eigenvalues of B are all negative [13]. Additionally, the Robertson-Schrödinger uncertainty relation must be satisfied, i.e.
The inequality (A15) is in fact a necessary and sufficient condition forσ to represent the steady state covariance matrix of a Gaussian state [13,15]. In our case, assuming that B is indeed stable, Eq. (A15) can be further simplified to Wγ T W ≥ 0, which is equivalent to the simpler condition γ ≥ 0. It is easy to check that the latter is always satisfied for the choice of γ adopted in the main text. the different contributions to the QFI (those due to the averages and the variances) behave in the zero and high temperature cases.
In our set-up, the coupling strengths have dimensions of frequency, whilst the averages and the steady state covariance matrix are dimensionless. Therefore, as the QFI is a dimensionful quantity that involves derivatives with respect to these parameters in its definition, it will have dimensions of the inverse frequency squared. We can remove the complication with units by switching to the dimensionless linear and quadratic coupling strengths defined asg 1 = g 1 /ω m andg 2 = g 2 /ω m , respectively. It is straightforward to show that the QFIs for the dimensionless and the original coupling strengths are related viaĨ ii = ω 2 m I ii for i = {1, 2}. Recall that non-monotonic behavior of the QFI was only observed in the high (and low) temperature scenarios in the case ofg 1 . In order to understand this behavior we compare the contributions from the variances and averages (as defined below eqn. (24)) to the QFI in the cases ofg 1 andg 2 , separately. Figure 5 shows the evolution of the QFI contributions with drive strength forg 1 . In the zero temperature limit (Fig. 5(a)) the contribution of the averages dominates at all |α|. Thus, here the majority of information about g 1 is always encoded in the averages. In contrast, in the high temper-ature scenario we observe a crossover of the contributions from the variances and the averages at log 10 (|α| 2 ) ∼ 4.7. The variances encode most of the QFI at lower drive strengths, but the contribution of the averages grows monotonically until it eventually becomes dominant. The contribution of the variances increases with drive before reaching a maximum and then declining (albeit rather slowly).
The important role of the variances at lower drive strengths for the high temperature case can be associated with the cavity back-action cooling effect [1,2] whereby the mechanical variances are progressively reduced below their values in thermal equilibrium as the drive strength is increased. Thus we naturally expect the variances to encode more information about the coupling as the temperature is increased and the impact of the back-action on the mechanical variances becomes more important. Furthermore, the back-action cooling effect saturates at very strong drive strengths [1,31] and hence it is perhaps no surprise that the contribution of the variances eventually starts to decline.
In Fig. 6 we investigate the QFI contributions for g 2 . Here, a crossover is only observed in the zero temperature scenario at log 10 (|α| 2 ) ∼ 4. However, unlike in theg 1 case, this crossover of the variances and averages contributions does not give rise to any non-monotonic behavior. In the high temperature scenario the majority of information about the pa-rameter is contained in the variances over the whole range of |α| studied.