# Variable Weights Neural Network For Diabetes Classification

Tanmay Rathi1 and Vipul2 *This work was not supported by any organization and
both authors contributed equally1Tanmay Rathi with the Department of
Electrical Engineering, Indian Institute of Technology Roorkee, Roorkee, INDIA
<EMAIL_ADDRESS>with the Department of Electrical Engineering,
Indian Institute of Technology Roorkee, Roorkee, INDIA<EMAIL_ADDRESS>
###### Abstract

As witnessed in the past year, where the world was brought to the ground by a
pandemic, fighting Life-threatening diseases have found greater focus than
ever. The first step in fighting a disease is to diagnose it at the right
time. Diabetes has been affecting people for a long time, and is growing among
people faster than ever. The number of people who have Diabetes reached 422
million in 2018, as reported by WHO, and the global prevalence of Diabetes
among adults above the age of 18 has risen to 8.5%. Now Diabetes is a disease
that shows no or very few symptoms among the people affected by it for a long
time, and even in some cases, people realize they have it when they have lost
any chance of controlling it. So getting Diabetes diagnosed at an early stage
can make a huge difference in how one can approach curing it. Moving in this
direction in this paper, we have designed a liquid machine learning[2]
approach to detect Diabetes with no cost using deep learning. In this work, we
have used a dataset of 520 instances from [1]. Our approach shows a
significant improvement in the previous state-of-the-art results. Its power to
generalize well on small dataset deals with the critical problem of lesser
data in medical sciences.

##  I Introduction

Diabetes mellitus has been around and affecting the human race since even
before medical science. For a long time, we have tried to conquer Diabetes by
finding new and better treatments and understanding its behavior and patterns.
It is mainly classified as Type 1 and Type 2 Diabetes, where Type 1 occurs due
to the body’s incapability to produce insulin or produce tiny amounts of it.
This insufficiency of insulin occurs mainly in children and adolescents. This
type had no cure for a very long time and would cause fatality within weeks or
months, but all that changed after the invention of artificial insulin. Type 2
is the one most prevalent among adults and occurs due to inefficient use of
the insulin produced inside the body. This type of Diabetes can originate and
grow in a person’s body without showing noticeable symptoms for over seven
years. There is another type of Diabetes, Gestational Diabetes, which mainly
develops in pregnant women due to hormonal changes. We have not yet been able
to discover a perfect single cure for Diabetes, so the best we can do is
reduce the risk of it as early as possible. Due to lack of extreme suggesting
symptoms, it’s medically not possible to confirm the disease. Still, it can be
predicted to happen in the future for a person using the trends seen among
those who suffer from it, even from very nominal symptoms. The text will
introduce a machine learning model to detect Diabetes at no cost from few
simple queries. We used the dataset from [1] containing 520 instances of
diabetic and non-diabetic people to train our model. We provided state-of-the-
art results on this dataset. Despite its small size, we generalize well over
the unseen data.

##  II Architecture

We proposed a liquid machine learning[2] approach for better generalization,
even with a small dataset. We call this layer simply a Variable Weight Layer.
Our layer’s remarkable features are 1) better generalization, 2) easy stacking
with other layers, 3) and full model development. The Variable Weight layer
can provide decent generalizability even on small datasets. It can be stacked
together to build a full architecture like a standard neural network layer or
can be used with a standard neural network and provides an extra variability
in weights for better generalization as in “Fig. 1”. The deriving idea behind
the layer development is that neural networks have stationary weights that
become constant ones trained, making it difficult for them to behave correctly
in changing conditions. Our proposed layer decides its weights based on the
given input and processes input accordingly, hence deciding how to process a
particular input.

The proposed Variable Weight Layer is shown in “Fig. 2”. The inputs from any
previous network layers first pass through the stationary wights to produce
the variable weights. These varying weights are then used to create the final
output with the static bias value as shown in “Fig. 2”. ”f1” and ”f2” are the
two activation functions used. For ”f1,” since predicted weights should be bi-
polar, we experimented with linear and Tanh activation functions. Finally, we
concluded that Tanh works better and also provides weight normalization as an
additional benefit. For ”f2,” we used the ReLU activation function. For more
detailed comparisons between the deriving equation of a neural network layer,
variable weight layer, and variable bias layer, refer to the results section.

Figure 1: Variable Weight Layer can be used with standard neural network layer
providing the additional benefit of Variability and Generalization.

Figure 2: Weight Prediction and Use in Variable Weight Layer.

One worth noting point here is the exponential increase in weights, which we
can deal with by predicting only bias or a particular weight. Indeed,
predicting only a single weight decreases the accuracy compared to the
original layer but reduces the weight drastically and still better than the
standard artificial neural network. As per our final version, we have also
shown that even with a lot fewer parameters, we can generalize well than the
standard artificial neural networks, proving that our changing weights can
generalize well over the unseen data.

##  III Experiments

###  III-A Dataset

The Diabetes dataset[1] contains 520 instances of diabetic and non-diabetic
peoples with 16 attributes. The dataset is a classification dataset where the
result is to classify a person as diabetic or not based on features given. We
decided to go with this dataset as this is a medical dataset and contains very
few training examples to train a big network. Which correctly represents the
current scenario in medical sciences. After preprocessing, our dataset finally
has 500 instances, out of which 314 are non-diabetic people, and 186 are
diabetic. Tables I and II summarise the dataset statistics.The data has been
collected from the patients using a direct questionnaire from Sylhet Diabetes
Hospital of Sylhet, Bangladesh.

TABLE I: Description of Dataset

| Number of Attributes | Number of Instances
---|---|---
| Diabetes Symptom
---
Dataset
16 | 520

TABLE II: Description of Attributes

Attributes | Values
---|---
Age | 1.20-35, 2.36-45, 3.46-55,4.56-65, 6 above 65.
Sex | 1.Male, 2.Female.
Polyuria | 1.Yes, 2.No.
Polydipsia | 1.Yes, 2.No.
sudden weight loss | 1.Yes, 2.No.
Weakness | 1.Yes, 2.No.
Polyphagia | 1.Yes, 2.No.
Genital thrush | 1.Yes, 2.No.
visual blurring | 1.Yes, 2.No.
Itching | 1.Yes, 2.No.
Irritability | 1.Yes, 2.No.
delayed healing | 1.Yes, 2.No.
partial paresis | 1.Yes, 2.No.
muscle stiffness | 1.Yes, 2.No.
Alopecia | 1.Yes, 2.No.
Obesity | 1.Yes, 2.No.
Class | 1.Positive, 2.Negative.

###  III-B Results

The performance of different algorithms on the Diabetes dataset is given in
Tables III and IV. The previous state-of-the-art method on this dataset was
Random Forest Algorithm. Using ten fold cross-validation, 97.4% of instances
were classified correctly, and using the percentage split technique could
classify 99% of the instances correctly. We first tested a standard neural
network and achieved a similar 97.4% cross-validation score and a 99%
percentage split. With our newly designed layer, we achieved a 99.4% cross-
validation score and 100% percentage-split score. Compared with the standard
neural network, our freshly designed layer generalizes well and increases 2%
on cross-validation and 1% percent on the percentage split. To reduce the
weight overhead, we only used our layer at the end of the network. To reduce
the parameters further, we train a network with half of the parameters and
achieve a 1% increase in cross-validation and percentage split metric. To
achieve such a decrease in parameters, we used our newly designed layer
configured to have variable biases rather than variable weights.

The equations behind the three layers we compared above are -
1) Standard Neural Network layer -

In standard Neural Networks, we use the previous layer outputs $l^{i-1}$ with
the current layers stationary weights $W_{s}$ and biases $b_{s}$ to produce
the output of the next layer $l^{i-1}$ as in below-given equations. Here
$l_{j}^{i}$ is the $j^{th}$ output of $i^{th}$ layer. $\widehat{g}$ is the
ReLU activation function.

$l_{1}^{(i)}=\widehat{g}(\,l_{1}^{(i-1)}*w_{s_{11}}+l_{2}^{(i-1)}*w_{s_{21}}+\ldots+l_{n}^{(i-1)}*w_{s_{n1}}+b_{\widehat{s}_{1}}\,)\\\
\\\
l_{2}^{(i)}=\widehat{g}(\,l_{1}^{(i-1)}*w_{s_{12}}+l_{2}^{(i-1)}*w_{s_{22}}+\ldots+l_{n}^{(i-1)}*w_{s_{n2}}+b_{\widehat{s}_{2}}\,)\\\
\\\ \vdots\\\ \\\
l_{n}^{(i)}=\widehat{g}(\,l_{1}^{(i-1)}*w_{s_{1n}}+l_{2}^{(i-1)}*w_{s_{2n}}+\ldots+l_{n}^{(i-1)}*w_{s_{nn}}+b_{\widehat{s}_{n}}\,)\\\
\\\
l^{(i)}=\begin{bmatrix}l_{1}^{(i)}&l_{2}^{(i)}&\ldots&l_{n}^{(i)}\end{bmatrix}$

2) Variable Weight Layer -

In variable weight layer we take the input from any previous layer of the
network $p=l^{k},k\,\epsilon\,[1,i-1]$ where $l^{i}$ is the current layer and
use it to predict the variable weights $W_{v}$ using stationary weights
$W_{s}$ and biases $b_{s}$. The predicted variable weights assist further in
the final output. $f_{1}$ is Tanh activation function and $\widehat{f_{2}}$ is
ReLU. In our results $p=l^{(i-1)}$.

$w_{v_{11}}=f_{1}(\,p_{1}*W_{s_{111}}+p_{2}*W_{s_{112}}+\ldots+p_{n}*W_{s_{11n}}+b_{s_{11}}\,)\\\
\\\
w_{v_{12}}=f_{1}(\,p_{1}*W_{s_{121}}+p_{2}*W_{s_{122}}+\ldots+p_{n}*W_{s_{12n}}+b_{s_{12}}\,)\\\
\\\ \vdots\\\
w_{v_{nn}}=f_{1}(\,p_{1}*W_{s_{nn1}}+p_{2}*W_{s_{nn2}}+\ldots+p_{n}*W_{s_{nnn}}+b_{s_{nn}}\,)\\\
\\\ W_{v}=\begin{bmatrix}w_{v_{11}}&w_{v_{12}}&\ldots&w_{v_{1n}}\\\
w_{v_{21}}&w_{v_{22}}&\ldots&w_{v_{2n}}\\\ \vdots&\vdots&\ddots&\vdots\\\
w_{v_{n1}}&w_{v_{n2}}&\ldots&w_{v_{nn}}\end{bmatrix}\\\ \\\ \\\
l_{1}^{(i)}=\widetilde{f_{2}}(\,p_{1}*w_{v_{11}}+p_{2}*w_{v_{21}}+\ldots+p_{n}*w_{v_{n1}}+b_{\widetilde{s}_{1}}\,)\\\
\\\
l_{2}^{(i)}=\widetilde{f_{2}}(\,p_{1}*w_{v_{12}}+l_{2}*w_{v_{22}}+\ldots+l_{n}*w_{v_{n2}}+b_{\widetilde{s}_{2}}\,)\\\
\\\ \vdots\\\ \\\
l_{n}^{(i)}=\widetilde{f_{2}}(\,p_{1}*w_{v_{1n}}+p_{2}*w_{v_{2n}}+\ldots+p_{n}*w_{v_{nn}}+b_{\widetilde{s}_{n}}\,)\\\
\\\
l^{(i)}=\begin{bmatrix}l_{1}^{(i)}&l_{2}^{(i)}&\ldots&l_{n}^{(i)}\end{bmatrix}\\\
$
3) Variable Bias Layer -

Variable Bias Layer contains a lot lesser parameters as we are only predicting
the bias values. We can also predict $w_{s_{1\,1-n}}$ or any other weights and
make them variable, giving us the freedom to decide the weight overhead. $g$
is linear activation here, and $\widehat{g}$ is ReLU.

$b_{v_{1}}=g(\,p_{1}*W_{s_{11}}+p_{2}*W_{s_{21}}+\ldots+p_{n}*W_{s_{n1}}+b_{s_{1}}\,)\\\
\\\
b_{v_{2}}=g(\,p_{1}*W_{s_{12}}+p_{2}*W_{s_{22}}+\ldots+p_{n}*W_{s_{n2}}+b_{s_{2}}\,)\\\
\\\
b_{v_{n}}=g(\,p_{1}*W_{s_{1n}}+p_{2}*W_{s_{2n}}+\ldots+p_{n}*W_{s_{nn}}+b_{s_{n}}\,)\\\
\\\ b_{v}=\begin{bmatrix}b_{v_{1}}&b_{v_{2}}&\ldots&b_{v_{n}}\end{bmatrix}\\\
l_{1}^{(i)}=\widehat{g}(\,p_{1}*w_{s_{11}}+p_{2}*w_{s_{21}}+\ldots+p_{n}*w_{s_{n1}}+b_{v_{1}}\,)\\\
\\\
l_{2}^{(i)}=\widehat{g}(\,p_{1}*w_{s_{12}}+p_{2}*w_{s_{22}}+\ldots+p_{n}*w_{s_{n2}}+b_{v_{2}}\,)\\\
\\\ \vdots\\\ \\\
l_{n}^{(i)}=\widehat{g}(\,p_{1}*w_{s_{1n}}+p_{2}*w_{s_{2n}}+\ldots+p_{n}*w_{s_{nn}}+b_{v_{n}}\,)\\\
\\\
l^{(i)}=\begin{bmatrix}l_{1}^{(i)}&l_{2}^{(i)}&\ldots&l_{n}^{(i)}\end{bmatrix}\\\
$

Acronyms used in Tables III and IV are NB - Naive Bayes, LR - Logistic
Regression, J48 Decesion Tree, RF - Random forest, NN - Neural Network, VB -
Variable Bias Network, VW - Varible Weight Network

TABLE III: Comparison Of Evaluation Metrices Using 10 Fold Cross-Validation

Evaluation Metrics | Cross-Validation
---|---
NB | LR | J48 | RF | NN | VB | VW
| Total Number
---
of
Instances
500 | 500 | 500 | 500 | 500 | 500 | 500
Correct Classifications | 437 | 462 | 478 | 487 | 487 | 492 | 497
87.4% | 92.4% | 95.6% | 97.4% | 97.4% | 98.4% | 99.4%
Incorrect Classifications | 63 | 38 | 22 | 13 | 13 | 8 | 3
12.6% | 7.6% | 4.4% | 2.6% | 2.6% | 1.6% | 0.6%

TABLE IV: Comparison of Evaluation Metrices Using Percentage Split(80:20)

Evaluation Metrics | Percentage Split
---|---
NB | LR | J48 | RF | NN | VB | VW
| Total Number
---
of
Instances
100 | 100 | 100 | 100 | 100 | 100 | 100
Correct Classifications | 88 | 91 | 95 | 99 | 99 | 100 | 100
88% | 91% | 95% | 99% | 99% | 100% | 100%
Incorrect Classifications | 12 | 9 | 5 | 1 | 1 | 0 | 0
12% | 9% | 5% | 1% | 1% | 0% | 0%

##  IV Conclusion

With the rapid growth in people being affected by Diabetes, it has naturally
become more important and necessary for it to be detected as early as possible
to help cure it, and if figured out at a very early stage, then even some
lifestyle changes can be sufficient, and one may not have to go through a lot
of medications. Early symptoms such as sudden weight loss, obesity, muscle
stiffness can help detect the possibility of having Diabetes through our
modified MLP neural network model, which has achieved fantastic percentage
split accuracy of 100% and has proven to generalize well for new entries with
a cross-validation score of 99.4%. However, this model can be trained and
modified with newly available data to perform even better, and some changes
can be made to the model to have greater efficiency.

##  V Future Directions

In the model architecture number of parameters to be determined during model
training increases exponentially with the number of features or neurons in the
layer preceding the custom or dynamic layer. This will not be feasible for
most classification data as they tend to have a large number of features, and
this could turn out to be highly inefficient in some cases. So several ways
can be tested to reduce this exponential factor to a more linear factor and
still have the advantage of that dynamic layer, such as instead of predicting
all the weights for a layer, it would be more efficient to predict only all
the bias or some weights which could turn out to be more useful in the model
as compared to others. A lot more experiments have to be conducted on many
more diverse kinds of datasets, which would help generalize this approach to a
lot more scenarios. We will soon develop a full report on applications and
provide a library for public use.

## References

  * [1] Islam, M M Faniqul & Ferdousi, Rahatara & Rahman, Sadikur & Bushra, Humayra. (2020). Likelihood Prediction of Diabetes at Early Stage Using Data Mining Techniques. 10.1007/978-981-13-8798-2_12.
  * [2] Ramin Hasani, Mathias Lechner, Alexander Amini, Daniela Rus, & Radu Grosu. (2020). Liquid Time-constant Networks.


* [13] Thomas Jech.  Set theory.  Springer Monographs in Mathematics. Springer-Verlag, Berlin, 2003.  The third millennium edition, revised and expanded.
  * [14] C. A. Johnson.  Distributive ideals and partition relations.  J. Symbolic Logic, 51(3):617–625, 1986.
  * [15] C. A. Johnson.  More on distributive ideals.  Fund. Math., 128(2):113–130, 1987.
  * [16] Akihiro Kanamori.  The higher infinite.  Springer Monographs in Mathematics. Springer-Verlag, Berlin, second edition, 2003.
  * [17] Akihiro Kanamori and Menachem Magidor.  The evolution of large cardinal axioms in set theory.  In Gert H. Müller and Dana S. Scott, editors, Higher Set Theory, pages 99–275, Berlin, Heidelberg, 1978. Springer, Berlin Heidelberg.
  * [18] Kenneth Kunen.  Saturated ideals.  J. Symbolic Logic, 43(1):65––76, 1978.
  * [19] Dan Saattrup Nielsen and Philip D. Welch.  Games and Ramsey-like cardinals.  J. Symbolic Logic, 84(1):408–437, 2019.
  * [20] Marion Scheepers.  Games that involve set theory or topology.  Vietnam J. Math., 23(2):169–220, 1995.
  * [21] Ralf Schindler.  Yet another characterization of remarkable cardinals.  Unpublished note, available at https://ivv5hpp.uni-muenster.de/u/rds/alternative_characterization.pdf, 2017\.
  * [22] Ian Sharpe and Philip D. Welch.  Greatly Erdös cardinals with some generalizations to the Chang and Ramsey properties.  Ann. Pure Appl. Logic, 162(11):863–902, 2011.
  * [23] Saharon Shelah.  Iterated forcing and changing cofinalities.  Israel J. Math., 40(1):1–32, 1981.
  * [24] Boris Šobot.  Games on Boolean algebras.  Ph.D. thesis, University of Novi Sad, 2009.
  * [25] Boban Veličković.  Playful Boolean algebras.  Trans. Am. Math. Soc., 296(2):727–740, 1986.
  * [26] Martin Zeman.  Inner Models and Large Cardinals.  De Gruyter, 2011.


# Sensing force gradients with cavity optomechanics while evading backaction

Elisabet K. Arvidsson  Department of Applied Physics, KTH Royal Institute of
Technology, Hannes Alfvéns väg 12, SE-114 19 Stockholm, Sweden  Ermes Scarano
Department of Applied Physics, KTH Royal Institute of Technology, Hannes
Alfvéns väg 12, SE-114 19 Stockholm, Sweden  August K. Roos  Department of
Applied Physics, KTH Royal Institute of Technology, Hannes Alfvéns väg 12,
SE-114 19 Stockholm, Sweden  Sofia Qvarfort  Nordita, KTH Royal Institute of
Technology and Stockholm University, Hannes Alfvéns väg 12, SE-114 19
Stockholm, Sweden  Department of Physics, Stockholm University, AlbaNova
University Center, SE-106 91 Stockholm, Sweden  David B. Haviland
<EMAIL_ADDRESS>Department of Applied Physics, KTH Royal Institute of
Technology, Hannes Alfvéns väg 12, SE-114 19 Stockholm, Sweden

(May 10, 2024)

###### Abstract

We study force gradient sensing by a coherently driven mechanical resonator
with phase-sensitive detection of motion via the two-tone backaction evading
measurement of cavity optomechanics. The response of the cavity to two
coherent pumps is solved by numerical integration of the classical equations
of motion, showing an extended region of monotonic response. We use Floquet
theory to model the fluctuations, which rise only slightly above that of the
usual backaction evading measurement in the presence of the mechanical drive.
Our analysis indicates that this sensing technique is advantageous for
applications such as Atomic Force Microscopy.

##  I Introduction

Cavity optomechanical systems have proved to be a good platform for sensing
force [1, 2]. Their versatility in size, mechanical frequency and variety of
optomechanical interaction has lead to their use in a wide range of
applications, from gravitational wave detection [3], to micrometer scale
accelerometers [4]. Sensing force gradients with a driven mechanical
oscillation is in most cases superior to sensing force via quasi-static
displacement of the test mass. A high quality factor mechanical resonator that
is coherently driven on resonance experiences a large change in phase for a
small change in the force gradient [5, 6]. Furthermore, the coherently driven
oscillation allows for lockin detection of displacement, which greatly
improves signal to noise ratio.

Cavity optomechanics has achieved remarkable sensitivity, however, no matter
how well the measurements are performed, they will inevitably have uncertainty
from fluctuations, i.e. noise. Two competing sources of noise, imprecision
noise (shot noise) and backaction noise, are minimized at the so-called
standard quantum limit (SQL), where each noise source adds an equivalent of a
quarter of a phonon of noise energy. [7]. There are techniques to avoid
backaction noise, so-called backaction evasion or quantum non-demolition
measurements, where the measured observable is a constant of the motion [8].
These techniques have lead to squeezing of mechanical motion below the zero
point fluctuation [9]. Other techniques have demonstrated improved sensitivity
by exploiting quantum correlations [10] and parametric amplification [11].
Building on backaction evasion, phase-sensitive detection of mechanical motion
has been demonstrated with a coherently driven mechanical mode [12, 13].

In this work we investigate the theoretical sensitivity of force gradient
detection with a backaction-evading measurement of a coherently driven
mechanical mode. Our work is motivated by the goal of improving dynamic-mode
atomic force microscopy (AFM) using a general measurement scheme schematically
described in Figure 1. The paper is structured as follows: In section II we
describe the system and Hamiltonian. Section III describes the pumping scheme
and the classical and quantum dynamics. In Section IV we present the classical
response and noise analysis. Finally we summarize and conclude in Section V.

##  II Model

Figure 1: Schematic of a optical mode with coupling $g_{0}$ to a mechanical
mode perturbed by an external interaction, e.g. the tip-surface force in AFM.
Both modes are driven and experience dissipation. The mechanical mode via
$\hat{b}_{in}$ and $\Gamma$, and the optical mode via $\hat{a}_{in}$ and
$\kappa$. Mechanical displacement is read out via the optical mode through
$\hat{a}_{out}$.

Standard cavity optomechanics couples a high-frequency electromagnetic mode
with resonance frequency $\omega_{c}$, to a mechanical mode with resonance
frequency $\omega_{m}$. We use the term ’optical’ to refer to the high-
frequency mode, which may be microwave or radio frequency. The coupling
results in the modulation of $\omega_{c}$ through mechanical displacement $x$,
described by

$\displaystyle\omega_{c}(x)\approx\omega_{c}+\frac{\partial\omega_{c}}{\partial
x}x+...$ (1)

where the higher order terms can be neglected if the modulation is small. The
total Hamiltonian is

$\displaystyle\hat{H}=\hbar\omega_{c}\hat{a}^{\dagger}\hat{a}+\frac{\hat{p}^{2}}{2m}+\frac{1}{2}m_{\mathrm{eff}}\omega_{m}^{2}\hat{x}^{2}-\frac{\hbar
g_{0}}{x_{\mathrm{zpf}}}\hat{a}^{\dagger}\hat{a}\hat{x},$ (2)

where $m_{\mathrm{eff}}$ is the effective mass,
$x_{\mathrm{zpf}}=\sqrt{\hbar/(2m_{\mathrm{eff}}\omega_{m})}$ the zero point
fluctuation of the mechanical mode, and
$g_{0}=-x_{zpf}\partial\omega_{c}/\partial x$ the single-photon optomechanical
coupling rate. The operators $\hat{a}$ and $\hat{a}^{\dagger}$ are the
annihilation and creation operators of the optical field, obeying the
canonical commutator relation $[\hat{a},\hat{a}^{\dagger}]=1$. The mechanical
position $\hat{x}=x_{\mathrm{zpf}}(\hat{b}+b^{\dagger})$ and momentum
$\hat{p}=-im_{\mathrm{eff}}\omega_{m}x_{\mathrm{zpf}}(\hat{b}-\hat{b}^{\dagger})$
operators obey $[\hat{x},\hat{p}]=i\hbar$.

We consider the example of non-contact atomic force microscopy where the
mechanical resonator is a cantilever beam, clamped at one end and having a
sharp tip at its free end. The tip undergoes small oscillations $x(t)$ about
its equilibrium position. The AFM scanner controls the distance $h$ between
this equilibrium position and the surface. When the tip is close to the
surface it experiences an interaction approximated by a van der Waals
potential between a sphere and a flat surface,

$U(x)=-\frac{HR_{tip}}{6(h+x)}$ (3)

where $H$ is the Hamaker constant and $R_{tip}$ radius of the tip [14]. For
small deviation $x$, we Taylor expand the potential, truncating after the
second order. Promoting $x$ to an operator gives,

$\displaystyle
U(\hat{x})=-\frac{HR_{\text{tip}}}{6h}\left(1-\frac{\hat{x}}{h}+\frac{\hat{x}^{2}}{h^{2}}\right)+\mathcal{O}(\hat{x}^{3}).$
(4)

The first term is simply a static shift which will not impact the equations of
motion and can be disregarded.

The Hamiltonian with the truncated tip-surface force included reads as follows

$\displaystyle\hat{H}$
$\displaystyle=\hbar\omega_{c}\hat{a}^{\dagger}\hat{a}+\frac{\hat{p}^{2}}{2m}+\frac{1}{2}m\omega_{m}^{2}\hat{x}^{2}-\frac{\hbar
g_{0}}{x_{\mathrm{zpf}}}\hat{a}^{\dagger}\hat{a}\hat{x}$ $\displaystyle\qquad-
F_{1}\hat{x}-F_{2}\hat{x}^{2}$ (5)

where $F_{1}=-HR_{\mathrm{tip}}/(6h^{2})$ is a static tip-surface force which
moves the equilibrium position of the mechanical mode. The force gradient
$F_{2}=HR_{\mathrm{tip}}/(6h^{3})$ results in an effective shift of the
mechanical resonance frequency,

$\displaystyle\omega_{\mathrm{eff}}=\sqrt{\omega_{m}^{2}-2F_{2}/m_{\mathrm{eff}}}.$
(6)

Rewriting the Hamiltonian with the effective resonance frequency,

$\displaystyle\hat{H}$
$\displaystyle=\hbar\omega_{c}\hat{a}^{\dagger}\hat{a}+\frac{\hat{p}^{2}}{2m}+\frac{1}{2}m\omega_{\mathrm{eff}}^{2}(F_{2})\hat{x}^{2}$
$\displaystyle\qquad-\frac{\hbar
g_{0}}{x_{\mathrm{zpf}}}\hat{a}^{\dagger}\hat{a}\hat{x}-F_{1}\hat{x}$ (7)

we transform into a frame rotating with respect to an optical frequency which
is detuned from cavity resonance $\omega_{p}=\omega_{c}+\Delta$. In terms of
mechanical annihilation and creation operators the Hamiltonian reads,

$\displaystyle\hat{H}_{rot}$
$\displaystyle=-\hbar\Delta\hat{a}^{\dagger}\hat{a}+\hbar\omega_{\mathrm{eff}}(F_{2})\hat{b}^{\dagger}\hat{b}-\hbar
g_{0}\hat{a}^{\dagger}\hat{a}(\hat{b}+\hat{b}^{\dagger})$
$\displaystyle\qquad-F_{1}x_{\mathrm{zpf}}(\hat{b}+\hat{b}^{\dagger}).$ (8)

Note that the coupling $g_{0}$ depends on the mechanical resonance frequency
through $x_{zpf}$, but here we treat it as a constant, considering only small
changes of $\omega_{m}$.

##  III Force gradient sensing

###  III.1 Optical pumps and mechanical drive

We consider the case where the mechanical mode is coherently driven at a
frequency $\omega_{d}$ which is close to $\omega_{\mathrm{eff}}$. The cavity
is pumped with two tones $\omega_{\pm}=\omega_{p}\pm\omega_{d}$, detuned from
the optical pumps center frequency $\omega_{p}$ by $\omega_{d}$. The resulting
sidebands from the optomechanical interaction interfere at $\omega_{p}$ where
the amplitude of the response depends on the phase of the mechanical motion.
The phase of the mechanical oscillator is thus detected in relation to the
phase of the beating envelope realized by the two optical pumps. We allow for
a slight detuning $\Delta$ of $\omega_{p}$ from the cavity resonance:
$\omega_{p}=\omega_{c}+\Delta$. The optical input field reads

$\displaystyle\hat{a}_{\mathrm{in}}=\bar{a}_{\mathrm{in_{-}}}e^{-i\omega_{-}t}+\bar{a}_{\mathrm{in_{+}}}e^{-i\omega_{+}t}+\hat{d}_{\mathrm{in}}$
(9)

where the first two terms are the classical pump signals, and
$\hat{d}_{\mathrm{in}}$ accounts for any input fluctuations. In a backaction
evading measurement the pumps should be placed with no detunings such that
$\omega_{p}=\omega_{c}$ and $\omega_{d}=\omega_{\mathrm{eff}}$ leading to
$\omega_{\pm}=\omega_{c}\pm\omega_{\mathrm{eff}}$, as any detuning induces
backaction into the measurement.

The mechanical mode is driven by a single tone close to its effective
resonance frequency, $\omega_{d}=\omega_{\mathrm{eff}}+\delta$ where $\delta$
is a small detuning. The mechanical drive reads

$\displaystyle\hat{b}_{\mathrm{in}}$
$\displaystyle=\bar{\beta}_{\mathrm{in}}e^{-i\omega_{d}t}+\hat{c}_{\mathrm{in}}$
(10)

with $\hat{c}_{\mathrm{in}}$ accounting for fluctuations. Here all drive terms
have complex amplitudes denoted with a bar, representing classical signals,
for example
$\bar{\beta}_{\mathrm{in}}=|\bar{\beta}_{\mathrm{in}}|e^{-i\phi_{m}}$.

###  III.2 Quantum and classical dynamics

In the frame rotating with the pump center frequency $\omega_{p}$, we obtain
the following Langevin equations of motion for the optical and mechanical
modes,

$\displaystyle\dot{\hat{a}}$
$\displaystyle=i\Delta\hat{a}+ig_{0}\hat{a}(\hat{b}+\hat{b}^{\dagger})-\frac{\kappa}{2}\hat{a}$
$\displaystyle\qquad-\sqrt{\kappa}(\bar{a}_{in_{-}}e^{i\omega_{d}t}+\bar{a}_{in_{+}}e^{-i\omega_{d}t})-\sqrt{\kappa}\hat{d}_{\mathrm{in}}$
$\displaystyle\dot{\hat{b}}$
$\displaystyle=-i\omega_{eff}(F_{2})\hat{b}+ig_{0}\hat{a}^{\dagger}\hat{a}+\frac{iF_{1}x_{zpf}}{\hbar}-\frac{\Gamma}{2}\hat{b}$
$\displaystyle\qquad-\sqrt{\Gamma}|\bar{\beta}_{\mathrm{in}}|e^{-i\phi_{m}}e^{-i\omega_{d}t}-\sqrt{\Gamma}\hat{c}_{\mathrm{in}}$
(11)

where $\kappa$ and $\Gamma$ are the total optical and mechanical dissipation
rates, respectively. We make the assumption that the internal cavity losses
are small in comparison to the external decay rate, such that
$\kappa_{ext}\approx\kappa$. In other words, the force sensor is designed for
an over-coupled cavity in which damping is dominated by useful signal flowing
to the next stage in the detection chain.

We model the optical and mechanical response as small quantum fluctuations
around a mean,

$\displaystyle\hat{a}$ $\displaystyle=\bar{\alpha}+\hat{d}$
$\displaystyle\hat{b}$ $\displaystyle=\bar{\beta}+\hat{c},$ (12)

where the expectation values $\bar{\alpha}=\braket{\hat{a}}$ and
$\bar{\beta}=\braket{\hat{b}}$ correspond to the classical oscillatory motion
of the optical and mechanical modes respectively. Inserting Eq. (III.2) into
Eq. (11), we find the Langevin equations for the quantum fluctuations.

$\displaystyle\dot{\hat{d}}$
$\displaystyle=i\Delta\hat{d}+ig_{0}(\bar{\alpha}+\hat{d})(\bar{\beta}+\hat{c}+\bar{\beta}^{*}+\hat{c}^{\dagger})-\frac{\kappa}{2}\hat{d}-\sqrt{\kappa}\hat{d}_{\mathrm{in}}$
$\displaystyle\dot{\hat{c}}$
$\displaystyle=-i\omega_{\mathrm{eff}}(F_{2})\hat{c}+ig_{0}(\bar{\alpha}^{*}+\hat{d}^{\dagger})(\bar{\alpha}+\hat{d})-\frac{\Gamma}{2}\hat{c}-\sqrt{\Gamma}\hat{c}_{\mathrm{in}}$
(13)

and the classical motion

$\displaystyle\dot{\bar{\alpha}}$
$\displaystyle=i\Delta\bar{\alpha}+ig_{0}\bar{\alpha}(\bar{\beta}+\bar{\beta}^{*})-\frac{\kappa}{2}\bar{\alpha}$
$\displaystyle\qquad-\sqrt{\kappa}(\bar{a}_{in_{-}}e^{i\omega_{d}t}+\bar{a}_{in_{+}}e^{-i\omega_{d}t})$
$\displaystyle\dot{\bar{\beta}}$
$\displaystyle=-i\omega_{\mathrm{eff}}(F_{2})\bar{\beta}+ig_{0}|\bar{\alpha}|^{2}+\frac{iF_{1}x_{\mathrm{zpf}}}{\hbar}-\frac{\Gamma}{2}\bar{\beta}$
$\displaystyle\qquad-\sqrt{\Gamma}|\bar{\beta}_{\mathrm{in}}|e^{-i\phi_{m}}e^{-i\omega_{d}t}$
(14)

In the usual two-tone backaction evading measurement [15] these equations do
not include the terms with $F_{1}$, $F_{2}$ or $\bar{\beta}_{\mathrm{in}}$,
and can therefore be solved analytically in the frequency domain. In our case,
the introduction of the perturbing force and the mechanical drive leads us to
solve the equations numerically.

##  IV Results

Solving the classical equations of motions (III.2) we show how the force
gradient $F_{2}$ is detected from the cavity field at $\omega_{p}$. Using
(III.2) we model the noise and compare it to the usual backaction evasion
case.

###  IV.1 Modelling the classical response

Figure 2: (a) and (b) Fourier transform of the solutions to the mechanical and
optical classical equations of motion where the frequency components included
in the following noise analysis are labeled. The black dashed lines indicate
the square root of the zero point fluctuations. (c) Change in magnitude of the
optical coherent tone $\bar{\alpha}_{c}$ normalized to a setpoint value
$\bar{\alpha}_{c,stp}$ as a function of mechanical drive phase $\phi_{m}$ and
effective mechanical resonance frequency $\omega_{\mathrm{eff}}$. To the right
a vertical cut is shown detailing the interference pattern and on the bottom a
horizontal cut where the purple region is highlighting the possible region of
operation.

We solve the classical equations of motion (III.2) by numerical integration
using the scipy function ODEint [16]. We simulate a tip-surface interaction
given in Section II with $HR_{\mathrm{tip}}=3.55\cdot 10^{-29}$ $\mathrm{Jm}$,
for example achieved with $H=$ $0.03\text{\,}\mathrm{aJ}$ and $R_{tip}=$
$1.183\text{\,}\mathrm{nm}$, with an optical pump amplitude of
$\bar{a}_{in_{\pm}}=1.62271\cdot 10^{5}$ $\sqrt{\mathrm{photons}/\mathrm{s}}$,
corresponding to the power needed to reach SQL for an optical (microwave) mode
with resonance frequency $\omega_{c}/2\pi=$ $4.5\text{\,}\mathrm{GHz}$ and
total decay rate $\kappa/2\pi=$ $1\text{\,}\mathrm{MHz}$, coupled to a
mechanical mode with resonance frequency $\omega_{m}/2\pi=$
$5.37\text{\,}\mathrm{MHz}$, total decay rate
$\Gamma/2\pi=$$2.3\text{\,}\mathrm{kHz}$, effective mass
$m_{\mathrm{eff}}=$$54\text{\,}\mathrm{pg}$, with an optomechanical coupling
rate $g_{0}/2\pi=$$1\text{\,}\mathrm{kHz}$. The large mechanical linewidth is
chosen to reduce computational time. We apply a mechanical drive to yield an
oscillation amplitude corresponding to 200$x_{\mathrm{zpf}}$. We initially set
the scanner position to $h=$ $0.3\text{\,}\mathrm{nm}$, approaching and
retracting from the surface to obtain a mechanical frequency shift of one
linewidth $\omega_{\mathrm{eff}}\pm\Gamma$. The Fourier transform of the
solutions of the mechanical motion and the optical field are plotted in Figure
2 (a) and (b) for a mechanical drive phase $\phi_{m}=0.86\pi$ and
$\omega_{eff}=\omega_{d}$. The approximated optical and mechanical fields
reads

$\displaystyle\bar{\alpha}(t)$
$\displaystyle\approx\bar{\alpha}_{-}e^{i\omega_{d}t}+\bar{\alpha}_{c}+\bar{\alpha}_{+}e^{-i\omega_{d}t}$
(15) $\displaystyle\bar{\beta}(t)$
$\displaystyle\approx\bar{\beta}_{0}+\bar{\beta}_{1}e^{-i\omega_{d}t}$ (16)

where any contribution less than that of the zero point fluctuations for each
respective mode is neglected. The term $\bar{\beta}_{0}$ is due to the static
force which induces a shift in the optical resonance. The term
$\bar{\beta}_{1}$ is the response to the mechanical drive modulating the
cavity resonance. In the following analysis $2g_{0}|\bar{\beta}_{1}|\ll\kappa$
so that the linearization in Eq. 1 is valid [17].

The quantity of interest is the term $\bar{\alpha}_{c}$ describing the
coherent optical response at frequency $\omega_{p}$, which is a superposition
of the upper mechanical sideband from the red pump, and the lower mechanical
sideband from the blue pump. The magnitude of $\bar{\alpha}_{c}$ depends on
the phase of the mechanical mode in relation to the phase of the beat created
by the two optical pumps. We change the phase of the mechanical mode either
through the mechanical drive phase $\phi_{m}$, or by changing effective
mechanical resonance frequency through the force gradient
$\omega_{\mathrm{eff}}(F_{2})$.

In Fig. 2(b) we plot the change in magnitude
$\Delta\bar{\alpha}_{c}=\bar{\alpha}_{c}-\bar{\alpha}_{c,stp}$ versus
$\phi_{m}$ and $\omega_{\mathrm{eff}}$ where $\bar{\alpha}_{c,stp}$ is a
setpoint for the AFM scanning feedback that controls $h$. A cut along the
vertical dash-dotted line reveals the interference between the two mechanical
sidebands, shown in the rightmost panel. The horizontal dashed line represents
a typical setpoint for feedback, locking to the mechanical phase in the middle
of an interference fringe. We see in the lower panel that the change in
response is monotonic with shift of mechanical resonance frequency
$|\omega_{d}-\omega_{\mathrm{eff}}|\leq 0.2\Gamma$ (purple shaded area). The
extended region of monotonic response to changes in $\omega_{\mathrm{eff}}$
makes $\bar{\alpha}_{c,stp}$ a useful setpoint for scanning feedback.

###  IV.2 Noise analysis

Figure 3: (a) Optical output spectrum detailing the difference between the
pure backaction evasion with no mechanical drive $|\bar{\beta}_{1}|=0$, and
force gradient sensing with a mechanical drive $|\bar{\beta}_{1}|=100$. The
addition of a coherent drive does not significantly add noise. (b) Quadrature
fluctuations of the mechanical mode as a function of the detuning, calculated
for $|\bar{\beta}_{1}|=100$. The range of detuning is chosen to match the
monotonic regime of the classical response. (c) Quadrature fluctuations
normalized to the zero point fluctuation plotted against $|\bar{\beta}_{1}|$
for the case of resonant drive ($\omega_{d}=\omega_{\mathrm{eff}}$) .

In the usual backaction evading measurement the two-tone pumped optical cavity
detects one of the mechanical quadratures without adding fluctuations to that
quadrature. However, in dynamic AFM where the mechanical mode is coherently
driven, the combination of both optical pumping and mechanical driving leads
to couplings which channel backaction fluctuations on to the detected
quadrature. When examining quantum limits to sensing force gradients, it is
important to consider the impact of these additional noise contributions.

Returning to the equations of motions, we insert Eqs. (15) and (16) into Eq.
(III.2) and linearize to obtain

$\displaystyle\dot{\hat{d}}$
$\displaystyle=i\tilde{\Delta}\hat{d}+ig_{0}[(\bar{\alpha}_{-}e^{i\omega_{d}t}+\bar{\alpha}_{c}+\bar{\alpha}_{+}e^{-i\omega_{d}t})\hat{c}+(\bar{\alpha}_{-}e^{i\omega_{d}t}$
$\displaystyle\quad+\bar{\alpha}_{c}+\bar{\alpha}_{+}e^{-i\omega_{d}t})\hat{c}^{\dagger}+\hat{d}(\bar{\beta}_{1}e^{-i\omega_{d}t}+\bar{\beta}^{*}_{1}e^{i\omega_{d}t})]$
$\displaystyle\quad-\frac{\kappa}{2}\hat{d}-\sqrt{\kappa}\hat{d}_{\mathrm{in}}$
(17) $\displaystyle\dot{\hat{c}}$
$\displaystyle=-i\omega_{\mathrm{eff}}(F_{2})\hat{c}+ig_{0}[(\bar{\alpha}_{-}^{*}e^{-i\omega_{d}t}+\bar{\alpha}^{*}_{c}+\bar{\alpha}^{*}_{+}e^{i\omega_{d}t})\hat{d}$
$\displaystyle\quad+(\bar{\alpha}_{-}e^{i\omega_{d}t}+\bar{\alpha}_{c}+\bar{\alpha}_{+}e^{-i\omega_{d}t})\hat{d}^{\dagger}]-\frac{\Gamma}{2}\hat{c}-\sqrt{\Gamma}\hat{c}_{\mathrm{in}}$
(18)

where $\tilde{\Delta}=\Delta+g_{0}(\bar{\beta}_{0}+\bar{\beta}_{0}^{*})$. At
the setpoint of operation the placement of the pump center frequency can be
adjusted such that $\tilde{\Delta}=0$. As the AFM scanner moves and the static
displacement $\bar{\beta}_{0}$ changes, a small non-zero detuning
$\tilde{\Delta}$ needs to be taken into account. We see that both the optical
(Eq. IV.2) and mechanical (Eq. IV.2) fluctuations are periodic with
$\omega_{d}$, complicating the analysis. To solve for these fluctuations we
turn to Floquet theory using the same approach as in [18, 19] where we expand
according to

$\displaystyle\hat{d}(t)=\sum_{n=-\infty}^{\infty}e^{in\omega_{d}t}\hat{d}^{(n)}(t),\quad\hat{d}^{{\dagger}}(t)=\sum_{n=-\infty}^{\infty}e^{in\omega_{d}t}\hat{d}^{(n){\dagger}}(t).$
(19)

where $\hat{d}^{(n)}$ denotes a Fourier operator (in either time or frequency
domain) rotating at frequency $n\omega_{d}$. Rewriting Eqs. (IV.2) and (IV.2)
in the form of stationary Langevin equations of motion and Fourier
transforming the equations gives

$\displaystyle-i\omega\hat{d}^{(n)}$
$\displaystyle=(i\tilde{\Delta}-in\omega_{d})\hat{d}^{(n)}+ig_{0}(\bar{\alpha}_{-}\hat{x}^{(n-1)}+\bar{\alpha}_{+}\hat{x}^{(n+1)})$
$\displaystyle\quad+ig_{0}\bar{\alpha}_{c}\hat{x}^{(n)}+ig_{0}(\bar{\beta}_{1}\hat{d}^{(n+1)}+\bar{\beta}_{1}^{*}\hat{d}^{(n-1)})$
$\displaystyle\quad-\frac{\kappa}{2}\hat{d}^{(n)}-\sqrt{\kappa}\hat{d}^{(n)}_{\mathrm{in}}$
(20) $\displaystyle-i\omega\hat{c}^{(n)}$
$\displaystyle=-i[\omega_{\mathrm{eff}}(F_{2})+n\omega_{d}]\hat{c}^{(n)}+ig_{0}(\bar{\alpha}_{-}\hat{d}^{(n-1){\dagger}}$
$\displaystyle\quad+\bar{\alpha}_{+}\hat{d}^{(n+1){\dagger}})+ig_{0}(\bar{\alpha}_{-}^{*}\hat{d}^{(n+1)}+\bar{\alpha}_{+}^{*}\hat{d}^{(n-1)})$
$\displaystyle\quad+ig_{0}(\bar{\alpha}_{c}^{*}\hat{d}^{(n)}+\bar{\alpha}_{c}\hat{d}^{(n){\dagger}})-\frac{\Gamma}{2}\hat{c}^{(n)}-\sqrt{\Gamma}\hat{c}_{\mathrm{in}}^{(n)}$
(21)

where $\hat{x}^{(n)}=\hat{c}^{(n)}+\hat{c}^{(n){\dagger}}$. For the full
calculation see Appendix A.

To understand how the mechanical drive and force gradient impacts
fluctuations, we first calculate the fluctuations of the mechanical mode and
later the spectrum of the output optical field.

We define the mechanical quadrature rotating at frequency $\omega_{d}$,
coupled to the optical field in the back action evading measurement [1] as

$\displaystyle\hat{X}=\frac{1}{\sqrt{2}}\left(\hat{c}e^{i\omega_{d}t}+\hat{c}^{\dagger}e^{-i\omega_{d}t}\right).$
(22)

The variance of this quadrature is (see Appendix C)

$\displaystyle\braket{\hat{X}^{2}}=\int\frac{d\omega}{2\pi}S^{(0)}_{XX}(\omega)=\frac{1}{2}\int\frac{d\omega}{2\pi}[S^{(-2)}_{cc}(\omega+\omega_{d})$
$\displaystyle+S^{(0)}_{cc^{\dagger}}(\omega+\omega_{d})+S^{(2)}_{c^{\dagger}c^{\dagger}}(\omega-\omega_{d})+S^{(0)}_{c^{\dagger}c}(\omega-\omega_{d})]$
(23)

where the power spectrum is defined as

$\displaystyle
S_{OO^{\prime}}^{(m)}(\omega)=\sum_{n}\int\frac{d\omega^{\prime}}{2\pi}\braket{\hat{O}^{(n)}(\omega+n\omega_{d})\hat{O}^{\prime(m-n)}(\omega^{\prime})}.$
(24)

To calculate the optical spectrum we first write the equations for the terms
at $n\pm 1$

$\displaystyle\hat{d}^{(n+1)}$
$\displaystyle=ig_{0}\chi_{c}\boldsymbol{(}\omega-\omega_{d}(n+1)\boldsymbol{)}(\bar{\alpha}_{-}\hat{x}^{(n)}+\bar{\alpha}_{+}\hat{x}^{(n+2)}$
$\displaystyle\quad+\bar{\alpha}_{c}\hat{x}^{(n+1)}+\bar{\beta}_{1}\hat{d}^{(n+2)}+\bar{\beta}_{1}^{*}\hat{d}^{(n)})$
$\displaystyle\hat{d}^{(n-1)}$
$\displaystyle=ig_{0}\chi_{c}\boldsymbol{(}\omega-\omega_{d}(n-1)\boldsymbol{)}(\bar{\alpha}_{-}\hat{x}^{(n-2)}+\bar{\alpha}_{+}\hat{x}^{(n)}$
$\displaystyle\quad+\bar{\alpha}_{c}\hat{x}^{(n-1)}+\bar{\beta}_{1}\hat{d}^{(n)}+\bar{\beta}_{1}^{*}\hat{d}^{(n-2)})$
(25)

where we have introduced the optical susceptibility
$\chi_{c}(\omega)=[\kappa/2-i(\omega+\tilde{\Delta})]^{-1}$. Discarding the
terms at $n\pm 2$, valid for fast periodic modulation
$\omega_{d}\gg\kappa,\Gamma$ [20], we use Eq. (25) to arrive at an expression
for the Fourier operators of the optical mode

$\displaystyle\hat{d}^{(n)}=\chi_{c}^{\prime}(\omega)ig_{0}(\bar{\alpha}_{-}^{\prime}x^{(n-1)}+\bar{\alpha}_{+}^{\prime}x^{(n+1)}+\bar{\alpha}_{c}^{\prime}x^{(n)})$
(26)

where the modified optical susceptibility is

$\displaystyle\chi_{c}^{\prime}(\omega)$
$\displaystyle=[\chi_{c}(\omega-\omega_{d}n)^{-1}+g_{0}^{2}|\bar{\beta}_{1}|^{2}\chi_{c}\boldsymbol{(}\omega-\omega_{d}(n+1)\boldsymbol{)}+$
$\displaystyle\quad+g_{0}^{2}|\bar{\beta}_{1}|^{2}\chi_{c}\boldsymbol{(}\omega-\omega_{d}(n-1)\boldsymbol{)}]^{-1}$
(27)

and where we have defined

$\displaystyle\bar{\alpha}_{-}^{\prime}=\bar{\alpha}_{-}+ig_{0}\bar{\beta}_{1}^{*}\chi_{c}\boldsymbol{(}\omega-\omega_{d}(n-1)\boldsymbol{)}\bar{\alpha}_{c}$
(28)
$\displaystyle\bar{\alpha}_{+}^{\prime}=\bar{\alpha}_{+}+ig_{0}\bar{\beta}_{1}\chi_{c}\boldsymbol{(}\omega-\omega_{d}(n+1)\boldsymbol{)}\bar{\alpha}_{c}$
(29)
$\displaystyle\bar{\alpha}_{c}^{\prime}=\bar{\alpha}_{c}+ig_{0}\bar{\beta}_{1}\chi_{c}\boldsymbol{(}\omega-\omega_{d}(n+1)\boldsymbol{)}\bar{\alpha}_{-}$
$\displaystyle\qquad+ig_{0}\bar{\beta}_{1}^{*}\chi_{c}\boldsymbol{(}\omega-\omega_{d}(n-1)\boldsymbol{)}\bar{\alpha}_{+}.$
(30)

The measured power spectral density is given by the $m=0$ component [18] of
Eq. 24 (see Appendix B)

$\displaystyle S_{d^{\dagger}d}^{(0)}(\omega)$
$\displaystyle=|\chi_{c}^{\prime}(\omega)|^{2}g_{0}^{2}[|\bar{\alpha}_{-}^{\prime}|^{2}S_{xx}^{(0)}(\omega+\omega_{d})$
$\displaystyle+|\bar{\alpha}_{+}^{\prime}|^{2}S_{xx}^{(0)}(\omega-\omega_{d})+\bar{\alpha}_{-}^{*\prime}\bar{\alpha}_{+}^{\prime}S_{xx}^{(-2)}(\omega+\omega_{d})$
$\displaystyle+\bar{\alpha}_{+}^{*\prime}\bar{\alpha}_{-}^{\prime}S_{xx}^{(2)}(\omega-\omega_{d})+|\bar{\alpha}_{c}^{\prime}|^{2}S_{xx}^{(0)}(\omega)$
$\displaystyle+\bar{\alpha}_{-}^{*\prime}\bar{\alpha}_{c}^{\prime}S_{xx}^{(-1)}(\omega+\omega_{d})+\bar{\alpha}_{+}^{*\prime}\bar{\alpha}_{c}^{\prime}S_{xx}^{(1)}(\omega-\omega_{d})$
$\displaystyle+\bar{\alpha}_{c}^{*\prime}\bar{\alpha}_{-}^{\prime}S_{xx}^{(1)}(\omega)+\bar{\alpha}_{c}^{*\prime}\bar{\alpha}_{+}^{\prime}S_{xx}^{(-1)}(\omega)]$
(31)

If $\bar{\alpha}_{-}=\bar{\alpha}_{+}$ the first four terms in Eqs. (IV.2)
constitute the standard backaction evasion spectrum, where we have not used
the rotating wave approximation.

Figure 3(a) shows a comparison of the optical output power spectral density
$S_{d_{out}^{\dagger}d_{out}}^{(0)}(\omega)=\kappa
S_{d^{\dagger}d}^{(0)}(\omega)$ for two cases: a coherent mechanical state
driven at $\omega_{d}=\omega_{\mathrm{eff}}$ with $|\bar{\beta}_{1}|=100$
corresponding to a oscillation magnitude of $200x_{zpf}$, and the usual
backaction evading measurement with $|\bar{\beta}_{1}|=0$. In both cases the
thermal mechanical phonon occupancy $\bar{n}_{m}^{th}=0$ and thermal cavity
photon occupancy $\bar{n}_{c}^{th}=0$. We see that the spectrum in the driven
case shows additional noise peaks around the two optical pump frequencies,
stemming from the mixing of the coherent optical response $\bar{\alpha}_{c}$
with the mechanical thermal noise. The peaks at $\pm 2\omega_{\mathrm{eff}}$
are due to backaction usually neglected in the rotating wave approximation.

To further see the effect of the mechanical resonance frequency shift from the
influence of the tip-surface force gradient, we calculate the variance of the
measured mechanical quadrature. The mechanical resonance frequency shift
results in an effective detuning of the two optical pumps from the optimal
pumping configuration of the backaction evading measurement
($\omega_{\pm}=\omega_{c}\pm\omega_{eff}$). The effect of this detuning on the
quadrature fluctuations has previously been investigated in [21] but without a
mechanical drive. Figure 3(b) shows the increase in the variance of the
mechanical quadrature due to backaction noise as a function of the effective
mechanical resonance frequency. For no detuning the quadrature fluctuation is
very close to that of the unmeasured mechanical state. As the tip-surface
force gradient introduces a change in the effective resonance frequency,
additional backaction noise is introduced into the quadrature. This leads to
an increase of the noise around the center frequency component
$\bar{\alpha}_{c}$ in the output optical field.

Figure 3(c) shows the fluctuations at zero detuning for several coherent
oscillation amplitudes $|\bar{\beta}_{1}|$ of the mechanical mode. We see that
the additional backaction noise increases with the amplitude of the mechanical
motion, however, the increase in quadrature fluctuations from the detuning is
much larger than that from the mechanical drive. The upper x-axis gives the
corresponding values of $2g_{0}|\bar{\beta}_{1}|/\kappa\ll 1$ where the
linearized optomechanical interaction is valid.

##  V Summary and Conclusions

We examined a scheme for force gradient sensing with a driven mechanical mode
where motion is detected with the two-tone backaction evading measurement of
cavity optomechanics. We envision this technique being useful in for example
low temperature atomic force microscopy but also for other types of force
gradient detection. The detection principle gives a monotonic response to
changes in force gradient which cause frequency shifts on the order of 20% of
the mechanical linewidth. A smaller mechanical linewidth improves the force
sensitivity, but limits the region of monotonic response for feedback control
of Atomic Force Microscopy.

Using Floquet theory to analyze the equations of motion of fluctuating terms,
we showed that quadrature fluctuations of the mechanical mode, as detected in
the optical output spectrum, are only slightly larger than that achieved with
the usual backaction evading measurement (see Fig. 3(b) and (c)). The majority
of the increase in fluctuations comes from the change of the effective
mechanical resonance due to the force gradient, resulting in a detuning of the
optical pumps from the optimal placement for backaction evasion. The
contribution of the mechanical drive to the increase in fluctuations is much
less than that of the change of the mechanical resonance frequency.

To further decrease the noise of the mechanical quadrature and the noise
around the center pump frequency, one can adjust the amplitudes of the two
optical pumps such that $|\bar{a}_{in-}|>|\bar{a}_{in+}|$ to produce a
squeezed mechanical state. For a mechanical resonator close to its ground
state, this squeezing can result in a quadrature fluctuation less than the
zero point fluctuations. However, un-equal pump amplitudes changes the region
of monotonic response in the classical interference pattern, as complete
destructive and constructive interference is no longer achieved.

## Acknowledgements

The authors would like to thank Clara Wanjura as well as the QAFM team for
fruitful discussions: T. Glatzel, M. Zutter, E. Tholén, D. Forchheimer, I.
Ignat, M. Kwon, and D. Platz. E.K.A., E.S., A.K.R., are supported by an EIC
Pathfinder (FET Open) grant 828966 – QAFM and the Swedish SSF grant
ITM17-0343. S.Q. is funded in part by the Wallenberg Initiative on Networks
and Quantum Information (WINQ) and in part by the Marie Skłodowska–Curie
Action IF programme Nonlinear optomechanics for verification, utility, and
sensing (NOVUS) – Grant-Number 101027183. Nordita is supported in part by
NordForsk.

## Data Availability

The data that support the findings of this study are openly available in
Zenodo at https://doi.org/10.5281/zenodo.11175476, reference number 11175476.

## References

  * Aspelmeyer et al. [2014] Markus Aspelmeyer, Tobias J. Kippenberg, and Florian Marquardt.  Cavity optomechanics.  _Rev. Mod. Phys._ , 86:1391–1452, Dec 2014.  doi: 10.1103/RevModPhys.86.1391.  URL https://link.aps.org/doi/10.1103/RevModPhys.86.1391.
  * Li et al. [2021] Bei-Bei Li, Lingfeng Ou, Yuechen Lei, and Yong-Chun Liu.  Cavity optomechanical sensing.  _Nanophotonics_ , 10(11):2799–2832, 2021.  doi: doi:10.1515/nanoph-2021-0256.  URL https://doi.org/10.1515/nanoph-2021-0256.
  * Abbott et al. [2016] B. P. Abbott et al.  Observation of gravitational waves from a binary black hole merger.  _Phys. Rev. Lett._ , 116:061102, Feb 2016.  doi: 10.1103/PhysRevLett.116.061102.  URL https://link.aps.org/doi/10.1103/PhysRevLett.116.061102.
  * Krause et al. [2012] Alexander G. Krause, Martin Winger, Tim D. Blasius, Qiang Lin, and Oskar Painter.  A high-resolution microchip optomechanical accelerometer.  _Nature Photonics_ , 6(11):768–772, Nov 2012\.  ISSN 1749-4893.  doi: 10.1038/nphoton.2012.245.  URL https://doi.org/10.1038/nphoton.2012.245.
  * Smith [1995] D. P. E. Smith.  Limits of force microscopy.  _Review of Scientific Instruments_ , 66(5):3191–3195, 05 1995.  ISSN 0034-6748.  doi: 10.1063/1.1145550.  URL https://doi.org/10.1063/1.1145550.
  * Albrecht et al. [1991] T. R. Albrecht, P. Grütter, D. Horne, and D. Rugar.  Frequency modulation detection using high‐Q cantilevers for enhanced force microscope sensitivity.  _Journal of Applied Physics_ , 69(2):668–673, 01 1991.  ISSN 0021-8979.  doi: 10.1063/1.347347.  URL https://doi.org/10.1063/1.347347.
  * Clerk et al. [2010] A. A. Clerk, M. H. Devoret, S. M. Girvin, Florian Marquardt, and R. J. Schoelkopf.  Introduction to quantum noise, measurement, and amplification.  _Rev. Mod. Phys._ , 82:1155–1208, Apr 2010.  doi: 10.1103/RevModPhys.82.1155.  URL https://link.aps.org/doi/10.1103/RevModPhys.82.1155.
  * Braginsky et al. [1980] Vladimir B. Braginsky, Yuri I. Vorontsov, and Kip S. Thorne.  Quantum nondemolition measurements.  _Science_ , 209(4456):547–557, Aug 1980.  doi: 10.1126/science.209.4456.547.  URL https://doi.org/10.1126/science.209.4456.547.
  * Wollman et al. [2015] E. E. Wollman, C. U. Lei, A. J. Weinstein, J. Suh, A. Kronwald, F. Marquardt, A. A. Clerk, and K. C. Schwab.  Quantum squeezing of motion in a mechanical resonator.  _Science_ , 349(6251):952–955, 2015.  doi: 10.1126/science.aac5138.  URL https://www.science.org/doi/abs/10.1126/science.aac5138.
  * Mason et al. [2019] David Mason, Junxin Chen, Massimiliano Rossi, Yeghishe Tsaturyan, and Albert Schliesser.  Continuous force and displacement measurement below the standard quantum limit.  _Nature Physics_ , 15(8):745–749, Aug 2019.  ISSN 1745-2481.  doi: 10.1038/s41567-019-0533-5.  URL https://doi.org/10.1038/s41567-019-0533-5.
  * Motazedifard et al. [2019] Ali Motazedifard, A. Dalafi, F. Bemani, and M. H. Naderi.  Force sensing in hybrid bose-einstein-condensate optomechanics based on parametric amplification.  _Phys. Rev. A_ , 100:023815, Aug 2019.  doi: 10.1103/PhysRevA.100.023815.  URL https://link.aps.org/doi/10.1103/PhysRevA.100.023815.
  * Hertzberg et al. [2009] J. B. Hertzberg, T. Rocheleau, T. Ndukum, M. Savva, A. A. Clerk, and K. C. Schwab.  Back-action-evading measurements of nanomechanical motion.  _Nature Physics_ , 6(3):213–217, December 2009\.  doi: 10.1038/nphys1479.  URL https://doi.org/10.1038/nphys1479.
  * Roos et al. [2023] August K. Roos, Ermes Scarano, Elisabet K. Arvidsson, Erik Holmgren, and David B. Haviland.  Kinetic inductive electromechanical transduction for nanoscale force sensing.  _Phys. Rev. Appl._ , 20:024022, Aug 2023.  doi: 10.1103/PhysRevApplied.20.024022.  URL https://link.aps.org/doi/10.1103/PhysRevApplied.20.024022.
  * Giessibl [2003] Franz J. Giessibl.  Advances in atomic force microscopy.  _Rev. Mod. Phys._ , 75:949–983, Jul 2003.  doi: 10.1103/RevModPhys.75.949.  URL https://link.aps.org/doi/10.1103/RevModPhys.75.949.
  * Bowen and Milburn [2015] W.P. Bowen and G.J. Milburn.  _Quantum Optomechanics_.  Taylor & Francis, Boca Raton, first edition, 2015.  ISBN 9781482259155.  URL https://doi.org/10.1201/b19379.
  * Virtanen et al. [2020] Pauli Virtanen, Ralf Gommers, Travis E. Oliphant, Matt Haberland, Tyler Reddy, David Cournapeau, Evgeni Burovski, Pearu Peterson, Warren Weckesser, Jonathan Bright, Stéfan J. van der Walt, Matthew Brett, Joshua Wilson, K. Jarrod Millman, Nikolay Mayorov, Andrew R. J. Nelson, Eric Jones, Robert Kern, Eric Larson, C J Carey, İlhan Polat, Yu Feng, Eric W. Moore, Jake VanderPlas, Denis Laxalde, Josef Perktold, Robert Cimrman, Ian Henriksen, E. A. Quintero, Charles R. Harris, Anne M. Archibald, Antônio H. Ribeiro, Fabian Pedregosa, Paul van Mulbregt, and SciPy 1.0 Contributors.  SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python.  _Nature Methods_ , 17:261–272, 2020.  doi: 10.1038/s41592-019-0686-2.
  * Law [1995] C. K. Law.  Interaction between a moving mirror and radiation pressure: A hamiltonian formulation.  _Phys. Rev. A_ , 51:2537–2541, Mar 1995.  doi: 10.1103/PhysRevA.51.2537.  URL https://link.aps.org/doi/10.1103/PhysRevA.51.2537.
  * Malz and Nunnenkamp [2016a] Daniel Malz and Andreas Nunnenkamp.  Floquet approach to bichromatically driven cavity-optomechanical systems.  _Physical Review A_ , 94:023803, Aug 2016a.  doi: 10.1103/PhysRevA.94.023803.  URL https://link.aps.org/doi/10.1103/PhysRevA.94.023803.
  * Malz and Nunnenkamp [2016b] Daniel Malz and Andreas Nunnenkamp.  Optomechanical dual-beam backaction-evading measurement beyond the rotating-wave approximation.  _Phys. Rev. A_ , 94:053820, Nov 2016b.  doi: 10.1103/PhysRevA.94.053820.  URL https://link.aps.org/doi/10.1103/PhysRevA.94.053820.
  * Malz [2019] Daniel Hendrik Malz.  _Periodic driving and nonreciprocity in cavity optomechanics_.  PhD thesis, University of Cambridge, 2019.  URL https://www.repository.cam.ac.uk/handle/1810/283253.
  * Wollman [2015] Emma Edwina Wollman.  _Quantum Squeezing of Motion in a Mechanical Resonator_.  PhD thesis, California Institute of Technology, 2015.

##  Appendix A Solving the noise equation of motions

We start with the equations of motion, where the optics is in a rotating frame
and the mechanics in the lab frame.

$\displaystyle\dot{\hat{d}}$
$\displaystyle=i\Delta\hat{d}+-ig_{0}(\bar{\alpha}_{-}e^{i\omega_{d}t}+\bar{\alpha}_{c}+\bar{\alpha}_{+}e^{-i\omega_{d}t}+\hat{d})(\bar{\beta}_{0}+\bar{\beta}_{1}e^{-i\omega_{d}t}+\hat{c}+\bar{\beta}_{0}+\bar{\beta}^{*}_{1}e^{i\omega_{d}t}+\hat{c}^{\dagger})-\frac{\kappa}{2}\hat{d}-\sqrt{\kappa}\hat{d}_{\mathrm{in}}$
(32) $\displaystyle\dot{\hat{c}}$
$\displaystyle=-i\omega_{\mathrm{eff}}(F_{2})\hat{c}+ig_{0}(\bar{\alpha}_{-}^{*}e^{-i\omega_{d}t}+\bar{\alpha}^{*}_{c}+\bar{\alpha}^{*}_{+}e^{i\omega_{d}t}+\hat{d}^{\dagger})(\bar{\alpha}_{-}e^{i\omega_{d}t}+\bar{\alpha}_{c}+\bar{\alpha}_{+}e^{-i\omega_{d}t}+\hat{d})-\frac{\Gamma}{2}\hat{c}-\sqrt{\Gamma}\hat{c}_{\mathrm{in}}$
(33)

Linearizing these equations gives

$\displaystyle\dot{\hat{d}}$
$\displaystyle=i\tilde{\Delta}\hat{d}+ig_{0}(\bar{\alpha}_{-}e^{i\omega_{d}t}\hat{c}+\bar{\alpha}_{c}\hat{c}+\bar{\alpha}_{+}e^{-i\omega_{d}t}\hat{c}+\bar{\alpha}_{-}e^{i\omega_{d}t}\hat{c}^{\dagger}+\bar{\alpha}_{c}\hat{c}^{\dagger}+\bar{\alpha}_{+}e^{-i\omega_{d}t}\hat{c}^{\dagger}+\hat{d}\bar{\beta}_{1}e^{-i\omega_{d}t}$
$\displaystyle+\hat{d}\bar{\beta}^{*}_{1}e^{i\omega_{d}t})-\frac{\kappa}{2}\hat{d}-\sqrt{\kappa}\hat{d}_{\mathrm{in}}$
(34) $\displaystyle\dot{\hat{c}}$
$\displaystyle=-i\omega_{\mathrm{eff}}(F_{2})\hat{c}+ig_{0}(\bar{\alpha}_{-}^{*}e^{-i\omega_{d}t}\hat{d}+\bar{\alpha}^{*}_{c}\hat{d}+\bar{\alpha}^{*}_{+}e^{i\omega_{d}t}\hat{d}+\bar{\alpha}_{-}e^{i\omega_{d}t}\hat{d}^{\dagger}+\bar{\alpha}_{c}\hat{d}^{\dagger}+\bar{\alpha}_{+}e^{-i\omega_{d}t}\hat{d}^{\dagger})-\frac{\Gamma}{2}\hat{c}-\sqrt{\Gamma}\hat{c}_{\mathrm{in}}$
(35)

and the conjugate

$\displaystyle\dot{\hat{d}}^{\dagger}$
$\displaystyle=-i\tilde{\Delta}\hat{d}^{\dagger}-ig_{0}(\bar{\alpha}_{-}^{*}e^{-i\omega_{d}t}\hat{c}^{\dagger}+\bar{\alpha}_{c}^{*}\hat{c}^{\dagger}+\bar{\alpha}_{+}^{*}e^{i\omega_{d}t}\hat{c}^{\dagger}+\bar{\alpha}_{-}^{*}e^{-i\omega_{d}t}\hat{c}+\bar{\alpha}_{c}^{*}\hat{c}+\bar{\alpha}_{+}^{*}e^{i\omega_{d}t}\hat{c}+\hat{d}^{\dagger}\bar{\beta}_{1}^{*}e^{i\omega_{d}t}$
$\displaystyle+\hat{d}^{\dagger}\bar{\beta}_{1}e^{-i\omega_{d}t})-\frac{\kappa}{2}\hat{d}^{\dagger}-\sqrt{\kappa}\hat{d}^{\dagger}_{\mathrm{in}}$
(36) $\displaystyle\dot{\hat{c}}^{\dagger}$
$\displaystyle=i\omega_{\mathrm{eff}}(F_{2})\hat{c}^{\dagger}-ig_{0}(\bar{\alpha}_{-}e^{i\omega_{d}t}\hat{d}^{\dagger}+\bar{\alpha}_{c}\hat{d}^{\dagger}+\bar{\alpha}_{+}e^{-i\omega_{d}t}\hat{d}^{\dagger}+\bar{\alpha}^{*}_{-}e^{-i\omega_{d}t}\hat{d}+\bar{\alpha}_{c}^{*}\hat{d}+\bar{\alpha}_{+}^{*}e^{i\omega_{d}t}\hat{d})-\frac{\Gamma}{2}\hat{c}^{\dagger}-\sqrt{\Gamma}\hat{c}^{\dagger}_{\mathrm{in}}$
(37)

where $\tilde{\Delta}=\Delta+g_{0}(\bar{\beta}_{0}+\bar{\beta}_{0}^{*})$. We
note that both in the optical equation of motion as well as in the mechanical
one there are rotating terms. By going into a rotating frame of the mechanics
some of the terms could become stationary. However, without neglecting some of
the contributions further analysis is needed. Following [18] the Floquet
theory stipulates that we expand the operators in terms of a Fourier
components:

$\displaystyle\hat{d}(t)$
$\displaystyle=\sum_{n=-\infty}^{\infty}e^{in\omega_{d}t}\hat{d}^{(n)}(t)$
(38) $\displaystyle\hat{d}^{\dagger}(t)$
$\displaystyle=\sum_{n=-\infty}^{\infty}e^{in\omega_{d}t}\hat{d}^{(n){\dagger}}(t)$
(39) $\displaystyle\hat{c}(t)$
$\displaystyle=\sum_{n=-\infty}^{\infty}e^{in\omega_{d}t}\hat{c}^{(n)}(t)$
(40) $\displaystyle\hat{c}^{\dagger}(t)$
$\displaystyle=\sum_{n=-\infty}^{\infty}e^{in\omega_{d}t}\hat{c}^{(n){\dagger}}(t)$
(41)

and then also define

$\displaystyle\hat{d}^{(n)}(\omega)$
$\displaystyle=\int_{-\infty}^{\infty}\mathrm{d}te^{i\omega
t}\hat{d}^{(n)}(t)$ (42) $\displaystyle\hat{d}^{(n){\dagger}}(\omega)$
$\displaystyle=\int_{-\infty}^{\infty}\mathrm{d}te^{i\omega
t}\hat{d}^{(n){\dagger}}(t)$ (43)

We note that this convention uses
$[\hat{d}^{(n)}(\omega)]^{\dagger}=\hat{d}^{(-n){\dagger}}(-\omega)$. Moving
along, when we take the derivative of one of the Fourier components, we must
invoke the chain law:

$\displaystyle\dot{\hat{d}}=\frac{d}{dt}\sum_{n=-\infty}^{\infty}e^{in\omega_{d}t}\hat{d}^{(n)}(t)=(in\omega_{d})\sum_{n=-\infty}^{\infty}e^{in\omega_{d}t}\hat{d}^{(n)}(t)+\sum_{n=-\infty}^{\infty}e^{in\omega_{d}t}\dot{\hat{d}}^{(n)}(t)$
(44)

We now insert these expressions into the equations of motion to find

$\displaystyle\sum_{n=-\infty}^{\infty}e^{in\omega_{d}t}\dot{\hat{d}}^{(n)}(t)$
$\displaystyle=\sum_{n=-\infty}^{\infty}e^{in\omega_{d}t}(i\tilde{\Delta}-in\omega_{d})\hat{d}^{(n)}(t)+ig_{0}\left(\bar{\alpha}_{-}e^{i\omega_{d}t}\sum_{n=-\infty}^{\infty}e^{in\omega_{d}t}\hat{c}^{(n)}(t)+\bar{\alpha}_{+}e^{-i\omega_{d}t}\sum_{n=-\infty}^{\infty}e^{in\omega_{d}t}\hat{c}^{(n){\dagger}}(t)\right)$
$\displaystyle+ig_{0}\left(\bar{\alpha}_{-}e^{i\omega_{d}t}\sum_{n=-\infty}^{\infty}e^{in\omega_{d}t}\hat{c}^{(n){\dagger}}(t)+\bar{\alpha}_{+}e^{-i\omega_{d}t}\sum_{n=-\infty}^{\infty}e^{in\omega_{d}t}\hat{c}^{(n)}(t)\right)$
$\displaystyle+ig_{0}\bar{\alpha}_{c}\left(\sum_{n=-\infty}^{\infty}e^{in\omega_{d}t}\hat{c}^{(n)}(t)+\sum_{n=-\infty}^{\infty}e^{in\omega_{d}t}\hat{c}^{(n){\dagger}}(t)\right)$
$\displaystyle+ig_{0}\left(\bar{\beta}_{1}e^{-i\omega_{d}t}\sum_{n=-\infty}^{\infty}e^{in\omega_{d}t}\hat{d}^{(n)}(t)+\bar{\beta}_{1}^{*}e^{i\omega_{d}t}\sum_{n=-\infty}^{\infty}e^{in\omega_{d}t}\hat{d}^{(n)}(t)\right)$
$\displaystyle-\frac{\kappa}{2}\sum_{n=-\infty}^{\infty}e^{in\omega_{d}t}\hat{d}^{(n)}(t)-\sqrt{\kappa}\sum_{n=-\infty}^{\infty}e^{in\omega_{d}t}\hat{d}^{(n)}_{\mathrm{in}}(t)$
$\displaystyle\sum_{n=-\infty}^{\infty}e^{in\omega_{d}t}\dot{\hat{c}}^{(n)}(t)$
$\displaystyle=\sum_{n=-\infty}^{\infty}e^{in\omega_{d}t}(-i\omega_{\mathrm{eff}}(F_{2})-in\omega_{d})\hat{c}^{(n)}(t)-$
$\displaystyle+ig_{0}\left(\bar{\alpha}_{-}e^{i\omega_{d}t}\sum_{n=-\infty}^{\infty}e^{in\omega_{d}t}\hat{d}^{(n)}(t)+\bar{\alpha}_{+}e^{-i\omega_{d}t}\sum_{n=-\infty}^{\infty}e^{in\omega_{d}t}\hat{d}^{(n){\dagger}}(t)\right)$
$\displaystyle+ig_{0}\left(\bar{\alpha}_{-}^{*}e^{-i\omega_{d}t}\sum_{n=-\infty}^{\infty}e^{in\omega_{d}t}\hat{d}^{(n)}(t)+\bar{\alpha}_{+}^{*}e^{i\omega_{d}t}\sum_{n=-\infty}^{\infty}e^{in\omega_{d}t}\hat{d}^{(n){\dagger}}(t)\right)$
(45)
$\displaystyle+ig_{0}\left(\bar{\alpha}_{c}^{*}\sum_{n=-\infty}^{\infty}e^{in\omega_{d}t}\hat{d}^{(n)}(t)+\bar{\alpha}_{c}\sum_{n=-\infty}^{\infty}e^{in\omega_{d}t}\hat{d}^{(n){\dagger}}(t)\right)-\frac{\Gamma}{2}\sum_{n=-\infty}^{\infty}e^{in\omega_{d}t}\hat{c}^{(n)}(t)$
$\displaystyle-\sqrt{\Gamma}\sum_{n=-\infty}^{\infty}e^{in\omega_{d}t}\hat{c}_{\mathrm{in}}^{(n)}(t).$
(46)

If we multiply in the rotating terms into the sums and then shift the indices
we can remove the sums such that

$\displaystyle\dot{\hat{d}}^{(n)}$
$\displaystyle=(i\tilde{\Delta}-in\omega_{d})\hat{d}^{(n)}+ig_{0}\left(\bar{\alpha}_{-}\hat{c}^{(n-1)}+\bar{\alpha}_{+}\hat{c}^{(n+1){\dagger}}\right)+ig_{0}\left(\bar{\alpha}_{-}\hat{c}^{(n-1){\dagger}}+\bar{\alpha}_{+}\hat{c}^{(n+1)}\right)+ig_{0}\bar{\alpha}_{c}\left(\hat{c}^{(n)}+\hat{c}^{(n){\dagger}}\right)$
$\displaystyle\qquad+ig_{0}\left(\bar{\beta}_{1}\hat{d}^{(n+1)}+\bar{\beta}_{1}^{*}\hat{d}^{(n-1)}\right)-\frac{\kappa}{2}\hat{d}^{(n)}-\sqrt{\kappa}\hat{d}^{(n)}_{\mathrm{in}}$
$\displaystyle\dot{\hat{c}}^{(n)}$
$\displaystyle=(-i\omega_{\mathrm{eff}}(F_{2})-in\omega_{d})\hat{c}^{(n)}+ig_{0}\left(\bar{\alpha}_{-}\hat{d}^{(n-1){\dagger}}+\bar{\alpha}_{+}\hat{d}^{(n+1){\dagger}}\right)+ig_{0}\left(\bar{\alpha}_{-}^{*}\hat{d}^{(n+1)}+\bar{\alpha}_{+}^{*}\hat{d}^{(n-1)}\right)+ig_{0}\left(\bar{\alpha}_{c}^{*}\hat{d}^{(n)}+\bar{\alpha}_{c}\hat{d}^{(n){\dagger}}\right)$
$\displaystyle\qquad-\frac{\Gamma}{2}\hat{c}^{(n)}-\sqrt{\Gamma}\hat{c}_{\mathrm{in}}^{(n)}.$
(47)

The equivalent conjugate equations are

$\displaystyle\dot{\hat{d}}^{(n){\dagger}}$
$\displaystyle=\left(-i\tilde{\Delta}-in\omega_{d}\right)\hat{d}^{(n){\dagger}}-ig_{0}\left(\bar{\alpha}_{-}^{*}\hat{c}^{(n+1){\dagger}}+\bar{\alpha}_{+}^{*}\hat{c}^{(n-1)}\right)-ig_{0}\left(\bar{\alpha}_{-}^{*}\hat{c}^{(n+1)}+\bar{\alpha}_{+}^{*}\hat{c}^{(n-1){\dagger}}\right)$
$\displaystyle\qquad-
ig_{0}\bar{\alpha}_{c}^{*}\left(\hat{c}^{(n){\dagger}}+\hat{c}^{(n)}\right)-ig_{0}\left(\bar{\beta}_{1}^{*}\hat{d}^{(n-1){\dagger}}+\bar{\beta}_{1}\hat{d}^{(n+1){\dagger}}\right)-\frac{\kappa}{2}\hat{d}^{(n){\dagger}}-\sqrt{\kappa}\hat{d}_{\mathrm{in}}^{(n){\dagger}}$
$\displaystyle\dot{\hat{c}}^{(n){\dagger}}$
$\displaystyle=\left(i\omega_{\mathrm{eff}}(F_{2})-in\omega_{d}\right)\hat{c}^{(n){\dagger}}-ig_{0}\left(\bar{\alpha}_{-}^{*}\hat{d}^{(n+1)}+\bar{\alpha}_{+}^{*}\hat{d}^{(n-1)}\right)-ig_{0}\left(\bar{\alpha}_{-}\hat{d}^{(n-1){\dagger}}+\bar{\alpha}_{+}\hat{d}^{(n+1){\dagger}}\right)-ig_{0}\left(\bar{\alpha}_{c}\hat{d}^{(n){\dagger}}+\bar{\alpha}_{c}^{*}\hat{d}^{(n)}\right)$
$\displaystyle\qquad-\frac{\Gamma}{2}\hat{c}^{(n){\dagger}}-\sqrt{\Gamma}\hat{c}_{\mathrm{in}}^{(n){\dagger}}.$
(48)

Casting this into a matrix equation we get

$\displaystyle\begin{pmatrix}i(\omega+\omega_{d})+A^{(0)}&A^{(-1)}&0\\\
A^{(1)}&i\omega+A^{(0)}&A^{(-1)}\\\
0&A^{(1)}&i(\omega-\omega_{d})+A^{(0)}\end{pmatrix}\crossproduct\begin{pmatrix}\hat{\mathbf{x}}^{(-1)}\\\
\hat{\mathbf{x}}^{(0)}\\\ \hat{\mathbf{x}}^{(1)}\\\
\end{pmatrix}=\begin{pmatrix}0\\\ \hat{\mathbf{x}}_{\mathrm{in}}\\\ 0\\\
\end{pmatrix}$ (49)

where

$\displaystyle
A^{(0)}=\begin{pmatrix}-i\tilde{\Delta}+\frac{\kappa}{2}&-ig_{0}\bar{\alpha}_{c}&0&-ig_{0}\bar{\alpha}_{c}\\\
-ig_{0}\bar{\alpha}_{c}^{*}&i\omega_{\mathrm{eff}}+\frac{\Gamma}{2}&-ig_{0}\bar{\alpha}_{c}&0\\\
0&ig_{0}\bar{\alpha}_{c}^{*}&i\tilde{\Delta}+\frac{\kappa}{2}&ig_{0}\bar{\alpha}_{c}^{*}\\\
ig_{0}\bar{\alpha}_{c}^{*}&0&ig_{0}\bar{\alpha}_{c}&-i\omega_{\mathrm{eff}}+\frac{\Gamma}{2}\end{pmatrix}$
(50)

and

$\displaystyle A^{(-1)}$
$\displaystyle=-ig_{0}\begin{pmatrix}\bar{\beta}_{1}&\bar{\alpha}_{+}&0&\bar{\alpha}_{+}\\\
\bar{\alpha}_{-}^{*}&0&\bar{\alpha}_{+}&0\\\
0&-\bar{\alpha}_{-}^{*}&-\bar{\beta}_{1}&-\bar{\alpha}_{-}^{*}\\\
-\bar{\alpha}_{-}^{*}&0&-\bar{\alpha}_{+}&0\end{pmatrix}$ (51) $\displaystyle
A^{(1)}$
$\displaystyle=-ig_{0}\begin{pmatrix}\bar{\beta}_{1}^{*}&\bar{\alpha}_{-}&0&\bar{\alpha}_{-}\\\
\bar{\alpha}_{+}^{*}&0&\bar{\alpha}_{-}&0\\\
0&-\bar{\alpha}_{+}^{*}&-\bar{\beta}_{1}^{*}&-\bar{\alpha}_{+}^{*}\\\
-\bar{\alpha}_{+}^{*}&0&-\bar{\alpha}_{-}&0\end{pmatrix}$ (52)

and

$\hat{\mathbf{x}}_{\mathrm{in}}=\begin{pmatrix}-\sqrt{\kappa}\hat{d}_{\mathrm{in}}&-\sqrt{\Gamma}\hat{d}_{\mathrm{in}}&-\sqrt{\kappa}\hat{d}_{\mathrm{in}}^{\dagger}&-\sqrt{\Gamma}\hat{c}_{\mathrm{in}}^{\dagger}\end{pmatrix}^{T}.$
(53)

where we’ve truncated at $n\pm 2$.

##  Appendix B Optical spectrum

The optical equation of motion is

$\displaystyle\chi_{c}^{\prime-1}(\omega)\hat{d}^{(n)}=ig_{0}(\bar{\alpha}_{-}^{\prime}x^{(n-1)}+\bar{\alpha}_{+}^{\prime}x^{(n+1)}+\bar{\alpha}_{c}^{\prime}x^{(n)})$
(54)

and for the conjugate

$\displaystyle(\chi_{c}^{\prime-1})^{*}(\omega)\hat{d}^{(n){\dagger}}=-ig_{0}(\bar{\alpha}_{-}^{\prime*}x^{(n+1)}+\bar{\alpha}_{+}^{\prime*}x^{(n-1)}+\bar{\alpha}_{c}^{\prime*}x^{(n)})$
(55)

To calculate the mechanical contributions to the optical spectrum we return to
the time domain and calculate the spectrum of two operators rotating at
different frequencies,
$\hat{x}_{a}=\hat{x}(t)e^{ia\omega_{d}t}=\hat{x}^{(n-a)}(t)$ and
$\hat{x}_{b}=\hat{x}(t)e^{ib\omega_{d}t}=\hat{x}^{(n-b)}(t)$. In the equation
of motion for the optical mode $a,b=-1,0,1$. We have

$\displaystyle S_{x_{a}x_{b}}(\omega,t)$
$\displaystyle=\int_{-\infty}^{\infty}d\tau
e^{i\omega\tau}\braket{\hat{x}(t+\tau)e^{ia\omega_{d}(t+\tau)}\hat{x}(t)e^{ib\omega_{d}t}}$
$\displaystyle=\sum_{n,m}\int_{-\infty}^{\infty}d\tau
e^{i\omega\tau}\braket{\hat{x}^{(n)}(t+\tau)e^{in\omega_{d}(t+\tau)}e^{ia\omega_{d}(t+\tau)}\hat{x}^{(m)}(t)e^{im\omega_{d}t}e^{ib\omega_{d}t}}$
$\displaystyle=\sum_{n,m}e^{i(n+m)\omega_{d}t}e^{i(a+b)\omega_{d}t}\int_{-\infty}^{\infty}d\tau
e^{i\omega\tau}e^{ia\omega_{d}\tau}e^{in\omega_{d}\tau}\braket{\hat{x}^{(n)}(t+\tau)\hat{x}^{(m)}(t)}$
$\displaystyle=\sum_{m}S_{x_{a}x_{b}}^{(m)}(\omega)e^{im\omega_{d}t}$ (56)

If we let $m=m^{\prime}-n$ we get

$\displaystyle S_{x_{a}x_{b}}(\omega,t)$
$\displaystyle=\sum_{n,m^{\prime}}e^{im^{\prime}\omega_{d}t}e^{i(a+b)\omega_{d}t}\int_{-\infty}^{\infty}d\tau
e^{i\omega\tau}e^{ia\omega_{d}\tau}e^{in\omega_{d}\tau}\braket{\hat{x}^{(n)}(t+\tau)\hat{x}^{(m^{\prime}-n)}(t)}$
(57)

The last integral we can instead write as

$\int_{-\infty}^{\infty}d\tau
e^{i\omega\tau}e^{ia\omega_{d}\tau}e^{in\omega_{d}\tau}\braket{\hat{x}^{(n)}(t+\tau)\hat{x}^{(m)}(t)}=\int_{-\infty}^{\infty}\frac{d\omega^{\prime}}{2\pi}\braket{\hat{x}^{(n)}(\omega+a\omega_{d}+n\omega_{d})\hat{x}^{(m^{\prime}-n)}(\omega^{\prime})}=S_{xx}^{(m^{\prime})}(\omega+a\omega_{d})$
(58)

We therefore have

$\displaystyle S_{x_{a}x_{b}}(\omega,t)$
$\displaystyle=\sum_{m}S^{(m)}_{x_{a}x_{b}}(\omega)e^{im\omega_{d}t}=\sum_{m}e^{im\omega_{d}t}e^{(a+b)i\omega_{d}t}S_{xx}^{(m)}(\omega+a\omega_{d})$
(59)

For the steady state spectrum we require that $m=0$, or $m^{\prime}=-a-b$.
This means that the contribution to the optical spectrum becomes

$\displaystyle
S_{x_{a}x_{b}}^{(0)}(\omega)=S_{xx}^{(-a-b)}(\omega+a\omega_{d})$ (60)

This means that the spectrum of the optical mode becomes

$\displaystyle S^{(0)}_{d^{\dagger}d}(\omega)$
$\displaystyle=g_{0}^{2}|\chi_{c}^{\prime}(\omega)|^{2}(|\bar{\alpha}_{-}^{\prime}|^{2}S_{xx}^{(0)}(\omega+\omega_{d})+|\bar{\alpha}_{+}^{\prime}|^{2}S_{xx}^{(0)}(\omega-\omega_{d})+\bar{\alpha}_{-}^{\prime*}\bar{\alpha}_{+}^{\prime}S_{xx}^{(-2)}(\omega+\omega_{d})+\bar{\alpha}_{+}^{\prime*}\bar{\alpha}_{-}^{\prime}S_{xx}^{(2)}(\omega-\omega_{d})+$
$\displaystyle+|\bar{\alpha}_{c}^{\prime}|^{2}S_{xx}^{(0)}(\omega)+\bar{\alpha}_{-}^{\prime*}\bar{\alpha}_{c}^{\prime}S_{xx}^{(-1)}(\omega+\omega_{d})+\bar{\alpha}_{+}^{\prime*}\bar{\alpha}_{c}^{\prime}S_{xx}^{(1)}(\omega-\omega_{d})+\bar{\alpha}_{c}^{\prime*}\bar{\alpha}_{-}^{\prime}S_{xx}^{(1)}(\omega)+\bar{\alpha}_{c}^{\prime*}\bar{\alpha}_{+}^{\prime}S_{xx}^{(-1)}(\omega)$
(61)

where

$\displaystyle S_{xx}^{(m)}(\omega)$
$\displaystyle=\sum_{n}\int\frac{d\omega^{\prime}}{2\pi}\braket{\hat{x}^{(n)}(\omega+n\omega_{d})\hat{x}^{(m-n)}(\omega^{\prime})}=$
$\displaystyle=\sum_{n}\int\frac{d\omega^{\prime}}{2\pi}\braket{(\hat{c}^{(n)}(\omega+n\omega_{d})+\hat{c}^{(n){\dagger}}(\omega+n\omega_{d}))(\hat{c}^{(m-n)}(\omega^{\prime})+\hat{c}^{(m-n){\dagger}})(\omega^{\prime})}$
(62)

##  Appendix C Mechanical fluctuations

We define the mechanical quadrature as
$\hat{X}=\frac{1}{\sqrt{2}}\left(\hat{c}e^{i\omega_{d}t+i\theta}+\hat{c}^{\dagger}e^{-i\omega_{d}t-i\theta}\right)$.
Following the same steps as in Appendix B we get for the mechanical quadrature

$\displaystyle S_{XX}(\omega,t)$
$\displaystyle=\frac{1}{2}\int_{-\infty}^{\infty}d\tau
e^{i\omega\tau}\braket{(\hat{c}(t+\tau)e^{i\omega_{d}(t+\tau)+i\theta}+\hat{c}^{\dagger}(t+\tau)e^{-i\omega_{d}(t+\tau)-i\theta})(\hat{c}(t)e^{i\omega_{d}t+i\theta}+\hat{c}^{\dagger}(t)e^{-i\omega_{d}t-i\theta})}$
$\displaystyle=\frac{1}{2}\sum_{n,m}\int_{-\infty}^{\infty}d\tau
e^{i\omega\tau}\braket{\hat{c}^{(n)}(t+\tau)e^{in\omega_{d}(t+\tau)}e^{i\omega_{d}(t+\tau)+i\theta}\hat{c}^{(m)}(t)e^{im\omega_{d}t}e^{i\omega_{d}t+i\theta}}$
$\displaystyle+\frac{1}{2}\sum_{n,m}\int_{-\infty}^{\infty}d\tau
e^{i\omega\tau}\braket{\hat{c}^{(n)}(t+\tau)e^{in\omega_{d}(t+\tau)}e^{i\omega_{d}(t+\tau)+i\theta}\hat{c}^{(m){\dagger}}(t)e^{im\omega_{d}t}e^{-i\omega_{d}t-i\theta}}$
$\displaystyle+\frac{1}{2}\sum_{n,m}\int_{-\infty}^{\infty}d\tau
e^{i\omega\tau}\braket{\hat{c}^{(n){\dagger}}(t+\tau)e^{in\omega_{d}(t+\tau)}e^{-i\omega_{d}(t+\tau)-i\theta}\hat{c}^{(m)}(t)e^{im\omega_{d}t}e^{i\omega_{d}t+i\theta}}$
$\displaystyle+\frac{1}{2}\sum_{n,m}\int_{-\infty}^{\infty}d\tau
e^{i\omega\tau}\braket{\hat{c}^{(n){\dagger}}(t+\tau)e^{in\omega_{d}(t+\tau)}e^{-i\omega_{d}(t+\tau)-i\theta}\hat{c}^{(m){\dagger}}(t)e^{im\omega_{d}t}e^{-i\omega_{d}t-i\theta}}$

For each term individually we have

$\displaystyle\sum_{n,m}\int_{-\infty}^{\infty}d\tau
e^{i\omega\tau}\braket{\hat{c}^{(n)}(t+\tau)e^{in\omega_{d}(t+\tau)}e^{i\omega_{d}(t+\tau)+i\theta}\hat{c}^{(m)}(t)e^{im\omega_{d}t}e^{i\omega_{d}t+i\theta}}$
$\displaystyle=\sum_{n,m^{\prime}}e^{i(n+m)\omega_{d}t}e^{2i\omega_{d}t}\int_{-\infty}^{\infty}d\tau
e^{i\omega\tau}e^{i\omega_{d}\tau}e^{in\omega_{d}\tau}\braket{\hat{c}^{(n)}(t+\tau)\hat{c}^{(m)}(t)}e^{2i\theta}$
(64)

and

$\displaystyle\sum_{n,m}\int_{-\infty}^{\infty}d\tau
e^{i\omega\tau}\braket{\hat{c}^{(n)}(t+\tau)e^{in\omega_{d}(t+\tau)}e^{i\omega_{d}(t+\tau)+i\theta}\hat{c}^{(m){\dagger}}(t)e^{im\omega_{d}t}e^{-i\omega_{d}t-i\theta}}$
$\displaystyle=\sum_{n,m}e^{i(n+m)\omega_{d}t}\int_{-\infty}^{\infty}d\tau
e^{i\omega\tau}e^{i\omega_{d}\tau}e^{in\omega_{d}\tau}\braket{\hat{c}^{(n)}(t+\tau)\hat{c}^{(m){\dagger}}(t)}$
(65)

and

$\displaystyle\sum_{n,m}\int_{-\infty}^{\infty}d\tau
e^{i\omega\tau}\braket{\hat{c}^{(n){\dagger}}(t+\tau)e^{in\omega_{d}(t+\tau)}e^{-i\omega_{d}(t+\tau)-i\theta}\hat{c}^{(m)}(t)e^{im\omega_{d}t}e^{i\omega_{d}t+i\theta}}$
$\displaystyle=\sum_{n,m}e^{i(n+m)\omega_{d}t}\int_{-\infty}^{\infty}d\tau
e^{i\omega\tau}e^{-i\omega_{d}\tau}e^{in\omega_{d}\tau}\braket{\hat{c}^{(n){\dagger}}(t+\tau)\hat{c}^{(m)}(t)}$
(66)

and

$\displaystyle\sum_{n,m}\int_{-\infty}^{\infty}d\tau
e^{i\omega\tau}\braket{\hat{c}^{(n){\dagger}}(t+\tau)e^{in\omega_{d}(t+\tau)}e^{-i\omega_{d}(t+\tau)-i\theta}\hat{c}^{(m){\dagger}}(t)e^{im\omega_{d}t}e^{-i\omega_{d}t-i\theta}}$
$\displaystyle=\sum_{n,m}e^{i(n+m)\omega_{d}t}e^{-2i\omega_{d}t}\int_{-\infty}^{\infty}d\tau
e^{i\omega\tau}e^{-i\omega_{d}\tau}e^{in\omega_{d}\tau}\braket{\hat{c}^{(n){\dagger}}(t+\tau)\hat{c}^{(m){\dagger}}(t)}e^{-2i\theta}.$
(67)

The power spectral density of the quadrature fluctuations will be given by the
time-independant contribution for each of the terms above:

$\displaystyle
S_{XX}^{(0)}(\omega)=\frac{1}{2}\boldsymbol{(}S_{cc}^{(-2)}(\omega+\omega_{d})e^{2i\theta}+S_{cc^{\dagger}}^{(0)}(\omega+\omega_{d})+S_{c^{\dagger}c}^{(0)}(\omega-\omega_{d})+S_{c^{\dagger}c^{\dagger}}^{(2)}(\omega-\omega_{d})e^{-2i\theta}\boldsymbol{)}$
(68)