# Macroscopic quantum entanglement between an optomechanical cavity and a
continuous field in presence of non-Markovian noise

S. Direkci<EMAIL_ADDRESS>Theoretical Astrophysics 350-17, California
Institute of Technology, Pasadena, California 91125, USA  K. Winkler  Vienna
Center for Quantum Science and Technology (VCQ), Faculty of Physics & Vienna
Doctoral School in Physics, University of Vienna, A-1090 Vienna, Austria  C.
Gut  Vienna Center for Quantum Science and Technology (VCQ), Faculty of
Physics & Vienna Doctoral School in Physics, University of Vienna, A-1090
Vienna, Austria  K. Hammerer  Institute for Theoretical Physics and Institute
for Gravitational Physics (Albert-Einstein-Institute), Leibniz University
Hannover, Appelstrasse 2, 30167 Hannover, Germany  M. Aspelmeyer  Vienna
Center for Quantum Science and Technology (VCQ), Faculty of Physics,
University of Vienna, A-1090 Vienna, Austria  Institute for Quantum Optics and
Quantum Information (IQOQI) Vienna, Austrian Academy of Sciences,
Boltzmanngasse 3, 1090 Vienna, Austria  Y. Chen  Theoretical Astrophysics
350-17, California Institute of Technology, Pasadena, California 91125, USA

###### Abstract

Probing quantum entanglement with macroscopic objects allows to test quantum
mechanics in new regimes. One way to realize such behavior is to couple a
macroscopic mechanical oscillator to a continuous light field via radiation
pressure. In view of this, the system that is discussed comprises an
optomechanical cavity driven by a coherent optical field in the unresolved
sideband regime where we assume Gaussian states and dynamics. We develop a
framework to quantify the amount of entanglement in the system numerically.
Different from previous work, we treat non-Markovian noise and take into
account both the continuous optical field and the cavity mode. We apply our
framework to the case of the Advanced Laser Interferometer Gravitational-Wave
Observatory (Advanced LIGO) and discuss the parameter regimes where
entanglement exists, even in the presence of quantum and classical noises.

††preprint: APS/123-QED

##  I Introduction

Entanglement is one of the hallmarks of the “quantumness” of physical systems.
Ideally, it is possible for macroscopic objects, massive and/or containing a
high number of degrees of freedom, to be entangled with each other. Yet in
practice, such macroscopic entanglement can be very delicate in presence of
decoherence. It is an intriguing challenge to create and verify macroscopic
entanglement, which is often viewed as expanding the limits of the quantum
regime.

Optomechanical systems are promising candidates for experimental demonstration
of macroscopic entanglement, partly due to their theoretical robustness
against temperature effects [1]. They can also be used to engineer the quantum
state of the mechanical system [2], where the entanglement is generated by the
momentum exchange between the light reflecting from the mechanical oscillator
– a phenomenon known as radiation pressure. It is theoretically well
understood and broadly discussed in the literature [3, 4, 5, 6, 7, 8], see for
example [1, 9] for a review.

Entanglement in optomechanical devices has been widely studied and there have
been several successful experimental realizations: stationary entanglement
between simultaneous light tones mediated by an optomechanical device [10,
11], generation of entanglement between spaced mechanical oscillators both in
the micro and macro regime via radiation pressure [12, 13, 14, 15], and
optomechanical entanglement between the light field and the mechanical
oscillator in a pulsed scheme [16] are examples of such demonstrations. There
also exist many proposals in the literature to further study macroscopic
quantum phenomena in optomechanical systems [17, 18, 19, 20, 21] and
entanglement between coupled oscillators in presence of non-Markovian baths
[22]. In this work, we consider stationary optomechanical entanglement, where
the system parameters (e.g. the driving), and statistical behavior thereof,
are not changing over time. Schemes to verify stationary optomechanical
entanglement were proposed in [6, 4, 23, 24, 25], whereas an experimental
demonstration has, to our knowledge, not been performed yet.

Our system consists of a single mechanical mode interacting with an optical
cavity mode and the quadratures of the light field exiting the cavity. At any
time $t$, we study the bipartite entanglement that is present in the joint
quantum state between mechanical mode, optical mode, and the light that has
exited the system during $t^{\prime}\leq t$. See Fig. 1 of Ref. [6], which
includes a space-time diagram that illustrates the configuration. In the
regime where the dynamics are linear, the state is Gaussian, and the noise
processes are Markovian (white), the open-system optomechanical dynamics is
solvable analytically and the state of the system can be known exactly [3].
The white-noise model describes well devices with high-frequency oscillators,
where only thermal excitations are expected, and in the limit of large bath
temperature where $k_{B}T\gg\hbar\omega_{m}$, $T$ is the temperature of the
bath and $\omega_{m}$ is the resonance frequency of the oscillator [26].

In this work, we extend the description to non-Markovian Gaussian noise
processes where analytical results are, to our knowledge, not available, thus
requiring numerical methods. This approach is applicable whenever non-white
noises processes, such as structural damping [27, 28], are relevant. We extend
on the methods developed in [6], incorporating a cavity, and, more
importantly, non-Markovian noise processes. The technique consists in
computing the minimal symplectic eigenvalue of the partially transposed
covariance matrix of the system, constructed with numerical methods, which
provides a measure of appropriate bipartite entanglement.

We first investigate entanglement in a generalized setting, considering a
heavy suspended oscillator with a low mechanical resonance frequency. This
corresponds to the free-mass limit, where the mechanical resonance frequency
is much smaller than the other characteristic frequencies of the system. Then,
we focus our attention on Advanced LIGO (aLIGO) [29] and use it as a case
study. It has been recently shown that by injecting squeezed vacuum, the
detector’s quantum noise can in principle surpass the free-mass standard
quantum limit (SQL) by 3 dB [30]. It is natural to ask whether this can
already imply that aLIGO has built quantum entanglement between the mirrors
and the light field. The answer to this question is non-trivial. First, from
[30], we see that the level of classical noise is not yet below the SQL 111The
quantum noise’s beating of the SQL is only inferred by subtracting a classical
noise floor that was obtained through calibration. Second, the strict
definition of entanglement we use here requires integrating over all
frequencies: it remains uncertain whether having noise below the SQL within a
certain finite frequency band automatically leads to entanglement. Therefore,
we parametrize aLIGO’s noise curves to investigate regimes where entanglement,
according to its strict definition, exists.

This paper is organized as follows: In Sec. II, we introduce the dynamics of
the system and its equations of motion. In Sec. III, we state our entanglement
criterion and the covariance matrix of the system for two partitions of
interest. To show the usefulness of our technique for systems with low
mechanical resonance frequencies, we investigate entanglement in a generalized
setting in Sec. IV. In Sec. V, we give details about aLIGO’s noise budget, and
talk about how we model it in our calculations. Finally, in Sec. VI, we
investigate whether there is entanglement between the mechanical oscillator
and the light field at aLIGO for the partitions of interest, given different
parametrizations of the classical noise curves.

##  II System Dynamics

Let us consider an optical cavity with a movable mirror, driven by a laser
with frequency $\omega_{0}$ close to one of the resonant frequencies of the
cavity, $\omega_{0}+\Delta$ [32]. The quantity $\Delta$ is often referred to
as the detuning frequency of the cavity. For such a system, the linearized
Hamiltonian in the interaction picture with the rotating-wave approximation
(RWA) is given by [33]

$\begin{gathered}H=\hbar\omega_{m}\hat{B}^{\dagger}\hat{B}+\hbar\Delta\hat{A}^{\dagger}\hat{A}-\hbar
G\hat{x}(\hat{A}^{\dagger}+\hat{A})\\\
+i\hbar\sqrt{2\gamma}\int_{-\infty}^{\infty}\frac{d\Omega}{2\pi}\left[\hat{A}^{\dagger}\hat{c}(\omega_{0}+\Omega)-\hat{A}\hat{c}^{\dagger}(\omega_{0}+\Omega)\right]\\\
+\int_{-\infty}^{\infty}\frac{d\Omega}{2\pi}\left[\hbar\Omega\hat{c}^{\dagger}(\omega_{0}+\Omega)\hat{c}(\omega_{0}+\Omega)\right]\,,\end{gathered}$
(1)

where $\hat{B}$ and $\hat{B}^{\dagger}$ are the annihilation and creation
operators of the mechanical mode (center of mass motion of the mirror),
$\hat{x}$ is the position of the center of mass of the mirror, $\omega_{m}$ is
the mechanical resonance frequency, $\hat{A}$ and $\hat{A}^{\dagger}$ are the
annihilation and creation operators of the cavity mode,
$\hat{c}(\omega_{0}+\Omega)$ and $\hat{c}^{\dagger}(\omega_{0}+\Omega)$ are
the annihilation and creation operators of the external vacuum light field at
frequency $\omega_{0}+\Omega$, $G$ is the linear optomechanical coupling
constant, and $\gamma$ is the decay rate of the cavity mode. The position and
momentum operators of the mirror are related to the creation and annihilation
operators of the mechanical mode by

$\displaystyle\hat{x}$
$\displaystyle=\sqrt{\frac{\hbar}{M\omega_{m}}}\,\frac{\left(\hat{B}+\hat{B}^{\dagger}\right)}{\sqrt{2}}\,,$
(2a) $\displaystyle\hat{p}$ $\displaystyle=\sqrt{\hbar
M\omega_{m}}\,\frac{\left(\hat{B}-\hat{B}^{\dagger}\right)}{\sqrt{2}i}\,.$
(2b)

Note that for the sake of convenience we chose a displaced frame where all
operators have zero mean. The mode operators satisfy the canonical commutation
relations,

$\begin{gathered}[\hat{A},\hat{A}^{\dagger}]=[\hat{B},\hat{B}^{\dagger}]=1.\end{gathered}$
(3)

aLIGO detectors are power- and signal-recycled Fabry-Perot Michelson
interferometers, which contain a high number of degrees of freedom. However,
the core optomechanics can still be studied by the Hamiltonian given above;
this reduction manifests itself in the “scaling-law” relations governing
aLIGO’s sensitivity as parameters of the signal-recycling cavity are modified
[32]. From the scaling-law, the coupling constant $G$ is related to the
parameters of the interferometer by,

$\begin{gathered}G=\sqrt{\frac{2\omega_{0}P_{c}}{\hbar Lc}}\,,\end{gathered}$
(4)

where $L$ is the arm length of the interferometer (i.e. the cavity length),
$P_{c}$ is the power circulating inside the cavity, and $c$ is the speed of
light in vacuum.

We can transform the Hamiltonian such that the cavity mode $(A,A^{\dagger})$
couples with the traveling wave at $z=0$ (where the point-wise cavity
interface is located). We derive this transformation in Appendix A. We use
$\hat{u}$ and $\hat{v}$ to label the field right before entering and right
after exiting the cavity, respectively.

In this paper, we restrict ourselves to $\Delta=0$. In this resonant case, the
system is unconditionally stable and it reaches a steady state, in which the
Heisenberg equations can be solved using Fourier transformation 222We use the
convention $\mathcal{F}\\{f(t)\\}=\int_{-\infty}^{\infty}f(t)e^{i\omega t}dt$.
To write down and solve the Heisenberg equations, instead of annihilation and
creation operators we use the Caves-Schumaker quadrature operators [35, 36]:

$\displaystyle\hat{u}_{1}(\Omega)=\frac{\hat{u}(\omega_{0}+\Omega)+\hat{u}^{\dagger}(\omega_{0}-\Omega)}{\sqrt{2}}\,,$
(5a)
$\displaystyle\hat{u}_{2}(\Omega)=\frac{\hat{u}(\omega_{0}+\Omega)-\hat{u}^{\dagger}(\omega_{0}-\Omega)}{\sqrt{2}i}\,,$
(5b)

where $\hat{u}_{j}^{\dagger}(\Omega)=\hat{u}_{j}(-\Omega)$. Quadratures
$\hat{v}_{1}(\Omega)$ and $\hat{v}_{2}(\Omega)$ are defined from
$\hat{v}(\omega_{0}+\Omega)$ and $\hat{v}(\omega_{0}-\Omega)$ in a similar
fashion. After setting $c=1$, their commutation relations are [37]

$\displaystyle[\hat{u}_{1}(\Omega),\hat{u}_{2}(\Omega^{\prime})]$
$\displaystyle=[\hat{v}_{1}(\Omega),\hat{v}_{2}(\Omega^{\prime})]=2\pi\delta(\Omega+\Omega^{\prime})\,,$
(6a) $\displaystyle[\hat{u}_{j}(\Omega),\hat{u}_{j}(\Omega^{\prime})]$
$\displaystyle=[\hat{v}_{j}(\Omega),\hat{v}_{j}(\Omega^{\prime})]=0\,,$ (6b)

for $j=1,2$. Then, in the time domain, we have

$\displaystyle\hat{u}_{j}(t)$
$\displaystyle=\int_{-\infty}^{\infty}\frac{d\Omega}{2\pi}\hat{u}_{j}(\Omega)e^{-i\Omega
t}\,,$ (7a) $\displaystyle[\hat{u}_{1}(t),\hat{u}_{2}(t^{\prime})]$
$\displaystyle=i\delta(t-t^{\prime})\,,$ (7b)
$\displaystyle[\hat{u}_{1}(t),\hat{u}_{1}(t^{\prime})]$
$\displaystyle=[\hat{u}_{2}(t),\hat{u}_{2}(t^{\prime})]=0\,,$ (7c)

for $j=1,2$. We similarly define quadrature operators $\hat{A}_{1,2}$ and
$\hat{B}_{1,2}$, in the time domain, with

$\hat{A}_{1}(t)=\frac{\hat{A}(t)+\hat{A}^{\dagger}(t)}{\sqrt{2}}\,,\quad\hat{A}_{2}(t)=\frac{\hat{A}(t)-\hat{A}^{\dagger}(t)}{\sqrt{2}i}\,,$
(8)

and similarly for $\hat{B}_{1,2}$. We also have

$\displaystyle[\hat{A}_{1}(t),\hat{A}_{1}(t)]$
$\displaystyle=[\hat{A}_{2}(t),\hat{A}_{2}(t)]=0\,,$ (9)
$\displaystyle[\hat{A}_{1}(t),\hat{A}_{2}(t)]$ $\displaystyle=i\,,$ (10)

and the same for $\hat{B}_{1,2}$. Note here that the commutators are for same-
time operators.

We include two classes of “classical” noises [37] in our system: a force noise
$\hat{n}_{F}$ and a sensing noise $\hat{n}_{X}$, arising from a quantum
treatment of the interaction of the system with its environment. Force noise
affects the center-of-mass motion of the mechanical oscillator by introducing
fluctuations in its momentum. We also introduce a velocity damping of the
oscillator, with a damping rate $\gamma_{m}$. $\gamma_{m}$ and $n_{F}$ are
associated with the heat bath(s) the mass is coupled to, with the value of
$\gamma_{m}$ and the spectrum of $n_{F}$ related by the fluctuation-
dissipation theorem [38]. Sensing noise affects how the position is measured
by the light field. In our model below it arises from fluctuations of the
reflecting surface that introduces noise in the cavity field.

In the Heisenberg picture, the dynamics are given by the Langevin equations of
motion. In the Fourier domain they are written as

$\displaystyle-i\Omega\hat{A}_{1}$
$\displaystyle=-\gamma\hat{A}_{1}+\sqrt{2\gamma}\hat{u}_{1}\,,$ (11a)
$\displaystyle-i\Omega\hat{A}_{2}$
$\displaystyle=-\gamma\hat{A}_{2}+\sqrt{2\gamma}\hat{u}_{2}+\sqrt{2}G(\hat{x}+\hat{n}_{X})\,,$
(11b) $\displaystyle-i\Omega\hat{x}$ $\displaystyle=\hat{p}/M\,,$ (11c)
$\displaystyle-i\Omega\hat{p}$
$\displaystyle=-\gamma_{m}\hat{p}-M\omega_{m}^{2}\hat{x}+\sqrt{2}\hbar
G\hat{A}_{1}+\hat{n}_{F}\,,$ (11d) $\displaystyle\hat{v}_{1}$
$\displaystyle=\hat{u}_{1}-\sqrt{2\gamma}\hat{A}_{1}\,,$ (11e)
$\displaystyle\hat{v}_{2}$
$\displaystyle=\hat{u}_{2}-\sqrt{2\gamma}\hat{A}_{2}\,.$ (11f)

We refer to Eqs. (11) as the Heisenberg equations for the rest of the article.
It is straightforward to solve them to obtain
$(\hat{x},\hat{p},\hat{A}_{1,2},\hat{v}_{1,2})$ in terms of the input fields,
$(\hat{u}_{1,2},\hat{n}_{X},\hat{n}_{F})$, referred as the input-output
relations of the system. More specifically, quantum fluctuations in the in-
going quadratures $\hat{u}_{1,2}(\Omega)$ drive the system’s quantum noise
[39]. From Eq. (11e) and (11f), reading the out-going field quadratures are
subject to noises in $\hat{u}_{1}$ and $\hat{u}_{2}$, giving rise to the shot
noise (SN) for that readout strategy [39]. On the other hand, from Eq. (11a),
we see that $\hat{u}_{1}$ drives $\hat{A}_{1}$, which in Eq. (11d) drives the
momentum of the test mass, which then shows up in the position of the test
mass via Eq. (11c), giving rise to quantum radiation pressure noise (QRPN),
also known as backaction noise in the literature. In general, the power
spectrum of the SN is inversely proportional to circulating power in the
cavity, while that of the QRPN is proportional to circulating power.

##  III Entanglement Criteria and Partitions

The canonical commutation relations imply that
$\mathbf{V}+\frac{1}{2}\mathbf{K}$ is positive semidefinite, where
$\mathbf{V}$ is the covariance matrix with
$\mathbf{V}_{ij}=\langle\\{\hat{X}_{i}-\langle\hat{X}_{i}\rangle,\hat{X}_{j}-\langle\hat{X}_{j}\rangle\\}\rangle/2$
and $\mathbf{K}_{ij}=[\hat{X}_{i},\hat{X}_{j}]$ is the commutator matrix of
the quadratures in the system. This relation can be stated as

$\mathbf{V}+\frac{1}{2}\mathbf{K}\geq 0\,.$ (12)

Here for an $N$-partite system containing $N$ harmonic oscillators, the
matrices $\mathbf{V}$ and $\mathbf{K}$ are $2N$-dimensional.

Figure 1:  Figurative representation of the two different partitions that is
used while testing for entanglement, which are partitioning by tracing over
(top row) and not tracing over (bottom row) the cavity. Note that the system
configuration is not changed, i.e. cavity is still present for both
partitions.

To test for bipartite entanglement in a multimode system, we use the
positivity of the partial transpose (PPT) criterion, which is necessary and
sufficient to test for the separability of one of the modes from the rest for
Gaussian systems [40, 41, 42]. To use the PPT criterion in this context, one
obtains the partial-transposed covariance matrix $\mathbf{V}_{\rm pt}$ by
reverting the momentum of that one mode (which puts a minus sign on the column
and the row that contains the momentum in question) [43]. The PPT criterion
for separability is expressed as

$\mathbf{V}_{\rm pt}+\frac{1}{2}\mathbf{K}\geq 0\Leftrightarrow{\rm
Separability}\,.$ (13)

The amount of entanglement is quantified by the logarithmic negativity,
$E_{\mathcal{N}}$ [44]. For a Gaussian state of N modes, it is defined as

$E_{\mathcal{N}}=\sum_{j=1}^{N}\text{max}\\{0,-\text{log}_{2}(\tilde{\nu}_{j})\\}\,,$
(14)

where $\tilde{\nu}_{j}$, $j=1,\dots N$ are the symplectic eigenvalues of the
partially transposed covariance matrix, $\mathbf{V}_{\rm pt}$, which are given
by the absolute values of the eigenvalues of $\mathbf{K}^{-1}\mathbf{V}_{\rm
pt}$. For 1 vs. $N-1$ mode partitions, only one of the symplectic eigenvalues
of $\mathbf{V}_{\rm pt}$ can have a magnitude smaller than 1 [45], therefore
there can be at most one negative eigenvalue of $\mathbf{V}_{\rm
pt}+\frac{1}{2}\mathbf{K}$. We label the corresponding symplectic eigenvalue
as $\tilde{\nu}_{min}$.

Using the PPT criterion, we test for entanglement between the mechanical
oscillator and the optical field in two ways: first, we construct the
covariance matrix $\mathbf{V}$ with the mechanical mode and the modes of the
light field, essentially tracing out the cavity mode. Here, we perform the
partial transpose operation with respect to the mechanical oscillator. Second,
we include the cavity mode in the covariance matrix while still taking the
partial transpose with respect to the mechanical oscillator, which corresponds
to measuring the entanglement between the oscillator and the joint system of
the cavity plus external light field. The two ways of partitioning are
depicted in Fig. 1. The elements of the covariance matrix $\mathbf{V}$ for
both types of partitions, as well as the discretization of $\mathbf{V}$ can be
found in Appendix D.

##  IV Entanglement in Presence of Non-Markovian Noises

Due to the numerical nature of the algorithm, we can tackle any noise spectral
density associated with $\hat{u}_{1}(\Omega)$, $\hat{u}_{2}(\Omega)$,
$\hat{n}_{F}(\Omega)$ and $\hat{n}_{X}(\Omega)$ using the PPT criterion
defined in Eq. (13) to determine whether entanglement is present for a given
partition. Conversely, this problem is analytically solvable only for some
simplified noise models to our knowledge, such as assuming all the noise
sources to have a white spectrum [6, 23].

To show the usefulness of the method, we investigate entanglement in heavy
suspended oscillators with relatively low mechanical resonance frequencies.
For such systems, the mechanical resonance frequency, $\omega_{m}$, is much
smaller than the other frequencies of the system, which is referred to as the
free-mass limit. In this setting, $\omega_{m}$ essentially does not affect the
dynamics. Furthermore, in this limit where $\Omega\gg\omega_{m},\gamma_{m}$,
the tradeoff between shot noise and QRPN give rise to the SQL [46], given by

$S_{\rm SQL}(\Omega)=\frac{2\hbar}{M\Omega^{2}}.$ (15)

In the context of suspended oscillators, $\hat{n}_{F}(\Omega)$ is the force
that gives rise to the suspension thermal noise, whereas $\hat{n}_{X}(\Omega)$
is an effective displacement that gives rise to coating thermal noise. When
thermal noise is due to the internal friction of the suspension or the
oscillator, the noise spectrum of $\hat{n}_{F}(\Omega)$ and
$\hat{n}_{X}(\Omega)$ decrease as $1/\Omega$ above internal resonances, which
is referred to as structural damping [47, 48, 49]. Evidently, structural
damping gives rise to non-Markovian noises, and the position-referred noise
spectral densities of $\hat{n}_{F}(\Omega)$ and $\hat{n}_{X}(\Omega)$ are
given, in the free-mass limit, by

$\displaystyle S_{F}(\Omega)$
$\displaystyle=\frac{2\hbar}{M}\frac{\Omega_{F}^{3}}{|\Omega|^{5}}\,,$ (16a)
$\displaystyle S_{X}(\Omega)$
$\displaystyle=\frac{2\hbar}{M}\frac{1}{\Omega_{X}|\Omega|}\,,$ (16b)

where $\Omega_{F}$ and $\Omega_{X}$ are the frequencies where the respective
noise curves cross the SQL, given in Eq. (15). Accordingly, they encode the
strength of the noise processes $n_{F}$ and $n_{X}$, relative to the SQL
level. For $\hat{n}_{F}(\Omega)$, the position-referred spectrum is related to
the noise spectral density, labeled as $S_{n_{F}}(\Omega)$, with
$S_{F}(\Omega)=S_{n_{F}}(\Omega)/M^{2}\Omega^{4}$. Whereas for
$\hat{n}_{X}(\Omega)$, the position-referred spectrum $S_{X}(\Omega)$ is also
the noise spectral density $S_{n_{X}}(\Omega)$ 333To find spectra of the form
of Eqs. (16), one assumes the free-mass and high-$Q$ limits where the
mechanical susceptibility behaves as $1/\Omega^{2}$; additionally, one must be
in frequency ranges where $S_{n_{F}}$ and $S_{n_{X}}$ are constant polynomial
roll-off – typical at frequencies above the resonance of suspended
oscillators.. The incoming field quadratures $\hat{u}_{j}$ have uncorrelated
white spectra given by Eq. (36), since we assume the incoming field to be at
vacuum state.

Table 1:  Relevant Experiments  Experiments | $\omega_{m}/(2\pi)$ | $\Omega_{F}/(2\pi)$ | $\Omega_{X}/(2\pi)$ | $\Omega_{q}/(2\pi)$ | $\Omega_{X}/\Omega_{F}$
---|---|---|---|---|---
LIGO | 1 | 10 | 50 | 60 | 5
KAGRA | 1 | 40 | 300 | 80 | 8
VIRGO | 1 | 10 | 40 | 50 | 4
Figure 2: Logarithmic negativity between the mechanical oscillator and the
outgoing light field in the free mass limit as a function of
$\Omega_{X}/\Omega_{F}$ for various $\Omega_{q}/\Omega_{F}$. We plot the
results for Markovian and non-Markovian force and sensing noise with dashed
and plain lines respectively. Note that entanglement does not exist for
$\Omega_{X}/\Omega_{F}\lesssim 1$ with Markovian and for
$\Omega_{X}/\Omega_{F}\lesssim 6.3$ with non-Markoivan noise sources, for
$\omega_{m}/(2\pi)=1$ Hz, $\gamma_{m}/(2\pi)=0.01$ Hz, and a cutoff frequency
of $\Omega_{c}/(2\pi)=0.001$ Hz.

We list some relevant experiments with low mechanical oscillation frequencies
in Table 1. Note that large scale experiments such as aLIGO, KAGRA [52], and
VIRGO [53] are subject to seismic noise in addition to thermal noises, which
usually does not follow structural damping and is disregarded in calculating
$\Omega_{F}$ for this Table. Seismic noise imposes further challenge upon
reaching entanglement, which is explained in detail in Sec. VI.1.

In the limit of large cavity bandwidth, $\gamma\gg\Omega$ the cavity can be
eliminated adiabatically. Then, the equations of motion in (11) are modified
as

$\displaystyle\hat{v}_{1}(\Omega)$ $\displaystyle=\hat{u}_{1}(\Omega)\,,$
(17a) $\displaystyle\hat{v}_{2}(\Omega)$
$\displaystyle=\hat{u}_{2}(\Omega)+\alpha(\hat{x}(\Omega)+\hat{n}_{X}(\Omega))\,,$
(17b)
$\displaystyle\begin{split}-i\Omega\hat{p}(\Omega)&=-\gamma_{m}\hat{p}(\Omega)-M\omega_{m}^{2}\hat{x}(\Omega)\\\
&+\hbar\alpha\hat{u}_{1}(\Omega)+\hat{n}_{F}(\Omega)\,,\end{split}$ (17c)
$\displaystyle-i\Omega\hat{x}(\Omega)$ $\displaystyle=\hat{p}(\Omega)/M\,,$
(17d)

where $\alpha=\Omega_{q}\sqrt{M/\hbar}$ and
$\Omega_{q}=2G\sqrt{\hbar/M\gamma}$ is the characteristic interaction
frequency. We fix $\omega_{m}/(2\pi)=1$ Hz, $\gamma_{m}/(2\pi)=0.01$ Hz,
$\Omega_{F}/(2\pi)=100$ Hz and vary $\Omega_{X}$ and $\Omega_{q}$. Note that
the aim of these choices is to ensure that
$\Omega_{F}\gg\omega_{m}\gg\gamma_{m}$, i.e. that the free-mass limit is
justified and we are in the large mechanical quality factor regime.
Furthermore, $\Omega_{F}\gg\omega_{m}$ implies the large bath temperature
limit where $k_{B}T\gg\hbar\omega_{m}$, $T$ being the temperature of the bath.
Lastly, $\Omega_{q}\gg\omega_{m}$ ensures that the measurement of the system
performed by light is faster than the dynamics of the system.

Since we work in the free-mass limit, the relevant parameters on which the
presence of entanglement depends are $\Omega_{q}$, $\Omega_{F}$, and
$\Omega_{X}$. In this limit, scaling all of them by the same constant factor
does not have an effect on the resulting dynamics. Therefore, we define the
unitless ratios $\Omega_{q}/\Omega_{F}$ and $\Omega_{X}/\Omega_{F}$, which
describe the system fully, in this limit. Note that the amount of classical
noises in the system decreases with increasing $\Omega_{X}/\Omega_{F}$.

We look for entanglement between the oscillator and the outgoing light field
by varying the ratio $\Omega_{X}/\Omega_{F}$ for various $\Omega_{q}$ in Fig.
2 and we plot the results with plain lines. We find that entanglement does not
exist for $\Omega_{X}/\Omega_{F}\lesssim 6.3$ for any value of $\Omega_{q}$.
For $\Omega_{X}/\Omega_{F}\gtrsim 6.3$, the system is entangled for any
(finite) $\Omega_{q}$, and the entanglement increases monotonously with
increasing $\Omega_{X}/\Omega_{F}$. This implies the existence of “universal”
entanglement, meaning that whether the system is entangled or not is
independent from $\Omega_{q}$, the interaction frequency (or, in other words,
how fast the system is measured by light). When the system is entangled, the
amount of entanglement increases when $\Omega_{q}$ is increased. We also note
that the threshold for $\Omega_{X}/\Omega_{F}$ above which the system is
entangled depends on the low-frequency cutoff $\Omega_{c}$ chosen in our model
for the spectra of $n_{F}$ and $n_{X}$. Here, our model is in the form
$1/(\Omega+\Omega_{c})$ with $\Omega_{c}/(2\pi)=0.001$ Hz for both force and
the sensing noises, such that the average noise power at DC is
$1/\Omega_{c}=1000/(2\pi)$ s. We saw that the threshold is inversely
proportional to the cutoff frequency.

In order to see the significance of non-Markovianity on the results, we repeat
the same procedure with a white force and sensing noise. The noise spectra are
given by $S_{n_{F}}(\Omega)=2\hbar M\Omega_{F}^{2}$,
$S_{F}(\Omega)=2\hbar\Omega_{F}^{2}/M\Omega^{4}$, and
$S_{n_{X}}(\Omega)=S_{X}(\Omega)=2\hbar/M\Omega_{X}^{2}$. The results can
again be found in Fig. 2, plotted with dashed lines, where the system is not
entangled for any $\Omega_{q}$ when $\Omega_{X}/\Omega_{F}\lesssim 1$, whereas
for $\Omega_{X}/\Omega_{F}\gtrsim 1$, entanglement exists for all $\Omega_{q}$
and the amount of entanglement increases with increasing $\Omega_{q}$. Since
we see this behavior for both Markovian and non-Markovian noise, we conjecture
that the universality of the entangling-disentangling phase transition is
independent of the power spectral densities of the classical noises, and that
the power spectral densities only determine the threshold above which we have
entanglement for all $\Omega_{q}$.

##  V Noise Model of aLIGO

Figure 3:  aLIGO noise budget obtained from pygwinc. Only the dominant
classical noise sources are plotted, along with the quantum noise. The total
force noise in the system is the sum of the seismic noise and the suspension
thermal noise, which is effective at low frequencies. The coating Brownian
thermal noise is taken as the main constituent of the sensing noise. As can be
seen from the figure, quantum noise dominates over the sensing noise by a
large margin at high frequencies.

The primary noise sources in aLIGO, other than the quantum noise, are the
following [29]: seismic noise and suspension thermal noise are the main
constituents of the force noise, and mirror coating thermal noise constitutes
the sensing noise. The noise spectrum is dominated by seismic and thermal
noise at low frequencies (until 100 Hz), and quantum noise at high
frequencies, cf. Fig. 3. The interferometer noise is stationary and Gaussian
to very good approximation in absence of glitches (i.e. transient noise
artifacts) [54].

Seismic noise occurs because of the ground motion at the interferometer sites.
This motion is $\sim 10^{-9}\text{m}/\sqrt{\text{Hz}}$ at 10 Hz [55]. To
provide isolation from this motion, the mirrors are suspended from quadruple
pendulums [56]. The primary components of thermal noise are suspension thermal
noise and coating Brownian noise. Suspension thermal noise occurs due to loss
in the fused silica fibres used in the final suspension stage [29], whereas
the coating Brownian noise (which is classified as a sensing noise) occurs due
to the mechanical dissipation in the coatings [57]. Other types of sensing
noise comprises of many noise sources that are dominant at high frequencies,
such as thermal fluctuations of the mirror’s shape, optical losses, or
photodetection inefficiency [58]. The noise budget of aLIGO can be found in
Figure 3.

Due to the classification above, we represent the sum of the seismic and the
suspension thermal noise with $\hat{n}_{F}(\Omega)$, and the coating thermal
noise with $\hat{n}_{X}(\Omega)$. We use the aLIGO noise budgets as given by
the Python Gravitational Wave Interferometer Noise Calculator library
(pygwinc) [59] and model them with rational functions of $\Omega^{2}$. The
noise spectra are modelled by

$\displaystyle S^{LIGO}_{F}(\Omega)$
$\displaystyle=\frac{\tau_{F}\alpha_{F_{1}}}{\left(\frac{\Omega}{\omega_{F}}\alpha_{F_{2}}\right)^{14}+1}\,,$
(18a) $\displaystyle S^{LIGO}_{X}(\Omega)$
$\displaystyle=\tau_{X_{1}}\left(\frac{\Omega}{\omega_{X}}\right)^{2}\alpha_{X_{1}}+\tau_{X_{2}}\alpha_{X_{2}}\,,$
(18b)

where $S^{LIGO}_{F}(\Omega)$ is the spectrum of the force noise and
$S^{LIGO}_{X}(\Omega)$ is the spectrum of the sensing noise. We model
$S^{LIGO}_{F}(\Omega)$ to decay as $\Omega^{-14}$ instead of $\Omega^{-16}$
(which is the expected behavior for quadruple suspension systems) since it
performs better at approximating the global behavior. The power spectral
densities are characterized by the time constants $\tau_{F}$, $\tau_{X_{1}}$,
$\tau_{X_{2}}$ and cutoff frequencies $\omega_{F}$, $\omega_{X}$. The values
of these parameters can be found in Appendix B, Table LABEL:tab:noiseParams.
$\alpha_{F_{1}}$, $\alpha_{F_{2}}$, $\alpha_{X_{1}}$, and $\alpha_{X_{2}}$ are
dimensionless constants that will be used to change the noise curves in
Section VI. Their effect on the noise curves can be seen in Fig. 4. If we set
all of them to be unity, we get our model of aLIGO noise curves.

Figure 4:  Force (top row) and sensing noise (bottom row) spectra parametrized
by $\alpha_{F_{1}}$, $\alpha_{F_{2}}$, $\alpha_{X_{1}}$, and $\alpha_{X_{2}}$.
The effect of $\alpha_{F_{1}}$ is to rise and lower the nominal noise strength
below the cutoff frequency and $\alpha_{F_{2}}$ shifts the cutoff frequency.
Similarly, $\alpha_{X_{1}}$ shifts the cutoff frequency where the sensing
noise starts increasing as $\Omega^{2}$, whereas $\alpha_{X_{2}}$ shifts the
nominal noise level.

Even though the parameters of aLIGO are well-known, they need to be recomputed
since we are reducing the antisymmetric mode of the interferometer to a single
cavity. The relation between the parameters of the antisymmetric mode and the
parameters of the reduced cavity have already been computed [32], however, we
choose to work numerically and find the parameters of the reduced cavity $G$,
$\gamma$, $M$, $\gamma_{m}$, $\omega_{m}$, and $L$ by fitting aLIGO’s quantum
noise spectrum found in pygwinc. The modeled classical and quantum noise
curves can be found in Appendix C, Fig. C.1. The fitted parameters can be
found in Appendix B, Table LABEL:tab:ligoParams.

##  VI Entanglement in aLIGO

In aLIGO, the spectrum is dominated by force noise for low frequencies and
quantum noise for high frequencies, therefore we expect the sensing noise to
not affect the entanglement significantly, which we observed with our
numerics. Then, we focus on the effect of the force noise on the entanglement
and use the parameters $\alpha_{F_{1}}$ and $\alpha_{F_{2}}$ to modify the
force noise spectrum, and set $\alpha_{X_{1}}=\alpha_{X_{2}}=1$ throughout all
of the following subsections. We also introduce resonant modes to investigate
the system as accurately as possible. Intuitively, we expect the entanglement
to be destroyed in presence of high classical noise levels.

###  VI.1 Effect of Force Noise

First, we investigate the effect of the force noise spectrum on entanglement
and we calculate the logarithmic negativity $E_{\mathcal{N}}$ as a function of
$\alpha_{F_{1}}$ and $\alpha_{F_{2}}$, for both of the partitions described in
Sec. III. For all pairs $\alpha_{F_{1}}$ and $\alpha_{F_{2}}$ here, we find
larger logarithmic negativity values when we do not trace over the cavity. The
results for the partition where we do not trace over the cavity can be found
in Fig. 5. The amount of entanglement in the system diminishes when the force
noise increases: that is towards the bottom-right of the plot where
$\alpha_{F_{1}}$ increases (proportional to the DC noise power) and
$\alpha_{F_{2}}$ decreases (inversely proportional to the noise bandwidth),
cf. Sec. V and Fig. 4. Our fit of aLIGO’s force noise level is for
$\alpha_{F_{1}}=\alpha_{F_{2}}=1$, hence this plot is for comparatively low
level of force noise. Further to the bottom-right, the numerics become
unstable and do not converge due to the wide range of orders of magnitude
entering the calculation; see Appendix E for a discussion of the numerical
implementation. Therefore, we cannot not give a definite answer about
optomechanical entanglement in aLIGO with our model yet.

Figure 5:  The effect of the force noise spectrum on the logarithmic
negativity when we do not trace over the cavity. Note that the force noise
levels increase toward the bottom-right of the figure and our fit of aLIGO’s
operation regime is for $\alpha_{F_{1}}=\alpha_{F_{2}}=1$.

Next, we look for entanglement in the absence of the seismic noise, since it
is the dominating contribution to the force noise of the system. Not being
subjected to seismic noise is a realistic scenario for space-based
gravitational-wave interferometers, such as the Laser Interferometer Space
Antenna (LISA) [60]. For aLIGO, in absence of seismic noise, suspension
thermal noise dominates the spectrum for low frequencies, which is modelled as

$S_{ST}(\Omega)=\frac{\tau_{ST}}{\left(\frac{\Omega}{\omega_{ST}}\right)^{8}+1}\,,$
(19)

similar to how the total force noise was modelled in Sec. V. The parameters
$\tau_{ST}$ and $\omega_{ST}$ can be found in Table LABEL:tab:noiseParams.
Without changing the other parameters in the system, we find that we can
achieve negativities of 1.52 and 1.72 for the partitions where we do and do
not trace over the cavity respectively. This means that, in the absence of
seismic noise, aLIGO has stationary optomechanical entanglement in its current
operating regime.

###  VI.2 Effect of Low Frequency Resonances

Both the classical and the quantum noise curves contain many resonances, as
can be seen from Fig. 3. The resonances in the seismic and suspension thermal
noise arise from the displacement noises of the rigid-body resonant modes of
the 4-stage suspension system [56]. These modes can be modeled with a sum of
Lorentzians multiplying the force noise spectrum defined in Eq. (18). Then,
the new formula defining the spectrum of the force noise is given as

$\displaystyle\begin{split}S^{R}_{F}(\Omega)&=\frac{\tau_{F}\alpha_{F_{1}}}{\left(\frac{\Omega}{\omega_{F}}\alpha_{F_{2}}\right)^{14}+1}\Biggl{(}1+\\\
&\sum_{i}\frac{A_{v_{i}}^{2}}{(\Omega-\Omega_{v_{i}})^{2}+(\frac{1}{2}\Gamma_{v_{i}})^{2}}\Biggr{)}\quad\text{with}\
\Omega>0\,,\end{split}$ (20a)
$\displaystyle\begin{split}S^{R}_{F}(\Omega)&=S^{R}_{F}(-\Omega)\,,\\\
\end{split}$ (20b)

where the sum is over the resonant modes. The parameters $\Omega_{v_{i}}$,
$\Gamma_{v_{i}}$, and $A_{v_{i}}$ are the mode frequencies, full widths at
half maximum (FWHM), and the amplitudes of the Lorentzians respectively, and
are listed in Appendix B, Table LABEL:tab:violinModeParams for the modes with
the biggest relative amplitudes.

Figure 6: The effect of resonant modes on the logarithmic negativity,
$E_{\mathcal{N}}$, for both partitions. We see that the negativity reduces
with increasing $\alpha_{F_{1}}$ or with the introduction of resonant modes.
Note that $\alpha_{F_{2}}=1$.

We investigate the effect of the resonant modes on the entanglement. In Fig.
6, with low force noise guarantying well-behaved numerics and setting
$\alpha_{F_{2}}=1$, we plot $E_{\mathcal{N}}$ for noise curves with and
without these modes (orange and gray, resp.) and for both partitions (circles
when the cavity is traced out). We see explicitly here that there is more
entanglement in the partition where we do not trace over the cavity. For low
noise (low $\alpha_{F_{1}}$), the negativity remains unchanged, hence the
resonant modes do not affect the entanglement significantly in this regime. As
the level of noise increases, resonant modes cause the logarithmic negativity
to decrease faster than the negativities calculated without the resonant modes
for both partitions. The system becomes separable when resonant modes are
included for $\alpha_{F_{1}}\approx 10^{-15}$ and $\alpha_{F_{1}}\approx
10^{-12}$ when we do and do not trace over the cavity respectively. It seems
reasonable to expect that entanglement will not emerge when force noise
becomes stronger. Hence, extrapolating $\alpha_{F_{1}}\to 1$, this is evidence
(but not a rigorous proof) that aLIGO in its current operation regime (but
without squeezed input) probably contains no optomechanical entanglement.

##  VII Conclusions

In this paper, we developed a framework to determine and quantify bipartite
entanglement in an optomechanical system in the non-sideband-resolved regime,
in presence of non-Markovian Gaussian noises. The main novelty of our work is
to enable the study of non-Markovian noise drives, which are common in devices
with low mechanical frequencies, typically associated with large/macroscopic
masses. Hence, we focused on macroscopic entanglement and used a free mass
with structural damping and coating Brownian noise as an initial example, and
then aLIGO as a more detailed case study.

We tested for bipartite entanglement by looking at the separability between
the mechanical oscillator and (i) the outgoing light field, (ii) the joint
system of the cavity and the outgoing light field. For low levels of classical
noise, we saw that the latter partition is more entangled compared to the
former. However, for high levels of classical noise, we did not see a
significant advantage in using one partition over the other.

In the low mechanical frequency regime, where the free-mass limit
approximation holds, we found that the presence of entanglement is independent
of the coherent optomechanical interaction strength and depends only on the
relative strength of force and sensing noise. This result is similar to that
already found in [4] for white noise drive. Furthermore, by parametrizing the
noise curves of aLIGO, we were able to find a region of noise curves where
entanglement exists, and we showed that there is a trade-off between the
overall noise level and the cutoff frequency. Due to the high level of the
current classical noise in aLIGO, we were not able to reach a definite
conclusion in terms of the existence of entanglement in the system, even
though it’s unlikely to have a significant amount of entanglement based on our
simulations. However, we saw that entanglement exists if we assume a system
without seismic noise, even when the suspension thermal noise is still
present. This is an important result since it shows that classical noises,
even at very low frequencies, are able to demolish entanglement.

We also looked at how resonances in the noise curves of the system affect the
amount of entanglement, and saw that entanglement is more resistant (i.e. it
disappears for higher levels of classical noise) for the partition where we
test for the separability between the mechanical oscillator and the joint
system of the cavity and the outgoing light field. For future work, we plan to
develop better sampling strategies to overcome numerical instabilities.

###### Acknowledgements.

We thank the Chen Quantum Group for helpful discussions. S.D. and Y.C.
acknowledge the support by the Simons Foundation (Award Number 568762). K.W.,
C.G., and M.A. received funding from the European Research Council (ERC) under
the European Union’s Horizon 2020 research and innovation program (grant
agreement No 951234), and from the Research Network Quantum Aspects of
Spacetime (TURIS). K.H. was supported by the Deutsche Forschungsgemeinschaft
(DFG, German Research Foundation) through Project-ID 274200144 – SFB 1227
(projects A06) and Project-ID 390837967 - EXC 2123.

##  Appendix A Transformation of the Optomechanical Hamiltonian

Strictly speaking, the Hamiltonian in Eq. (1) is written in terms of
Schrödinger operators. The symbol $\Omega$ in Eq. (1) is used to label a
spatial mode which has a reduced wavenumber of $\Omega/c$, and a free
oscillation frequency of $\Omega$. More specifically, $\omega_{0}+\Omega$ is
used to indicate a spatial mode whose wave number is $(\omega_{0}+\Omega)/c$,
where $c$ is the speed of light. We can also write the same Hamiltonian in the
spatial domain. Following [33] and setting $c=1$, we define

$\begin{gathered}\hat{c}(z)=\int_{-\infty}^{\infty}\frac{d\Omega}{2\pi}\hat{c}(\omega_{0}+\Omega)e^{i\Omega
z}\,,\end{gathered}$ (21)

which represents a spatial mode of the light field with wave number $\Omega$ —
or a temporal mode of the free light field with frequency $\Omega$, since we
assume no dispersion. The Hamiltonian can then be rewritten as

$\begin{gathered}H=\hbar\omega_{m}\hat{B}^{\dagger}\hat{B}+\hbar\Delta\hat{A}^{\dagger}\hat{A}-\hbar
G\hat{x}(\hat{A}^{\dagger}+\hat{A})\\\
+i\hbar\sqrt{2\gamma}[\hat{A}^{\dagger}\hat{c}(z=0)-\hat{A}\hat{c}^{\dagger}(z=0)]\\\
-i\hbar\int_{-\infty}^{\infty}\hat{c}^{\dagger}(z)\partial_{z}\hat{c}(z)\,dz\,.\end{gathered}$
(22)

The only non-zero commutators for the creation and annihilation operators, as
well as the spatial modes of the light field are

$\displaystyle[\hat{c}(\omega_{0}+\Omega),$
$\displaystyle\hat{c}^{\dagger}(\omega_{0}+\Omega^{\prime})]=2\pi\delta(\Omega-\Omega^{\prime})\,,$
(23a) $\displaystyle[\hat{c}(z),$
$\displaystyle\hat{c}^{\dagger}(z^{\prime})]=\delta(z-z^{\prime})\,,$ (23b)

where $\delta$ is the Dirac delta distribution.

##  Appendix B Parameter Tables

This appendix contains the tables with the numerical values of the parameters
used in modelling the classical and quantum noise curves, including the
resonant modes. The parameters were calculated by minimizing the mean squared
error between the actual noise curves of aLIGO taken from pygwinc and the
theoretical models, characterized by the parameters of interest, sampled
logarithmically in frequency. Table LABEL:tab:noiseParams contains the
parameters for the classical force and sensing noise in aLIGO defined in Eqs.
(18), Table LABEL:tab:ligoParams contains the parameters of aLIGO introduced
in the optomechanical Hamiltonian of Sec. II, and lastly, Table
LABEL:tab:violinModeParams contains the parameters of the resonant modes
present in aLIGO’s noise curves, modeled with Eq. (20).

Table 2: Classical Noise Model Parameters Parameter | Symbol | Value | Units
---|---|---|---
Force Noise Time Constant | $\tau_{F}$ | $1.6\cdot 10^{-20}$ | s
Force Noise Cutoff Frequency | $\omega_{F}$ | $2\pi\cdot 0.25$ | rad/s
Sensing Noise Time Constant 1 | $\tau_{X_{1}}$ | $10^{-50}$ | s
Sensing Noise Time Constant 2 | $\tau_{X_{2}}$ | $10^{-48}$ | s
Sensing Noise Cutoff Frequency | $\omega_{X}$ | $2\pi\cdot 10^{4}$ | rad/s
Suspension Thermal Noise Time Constant | $\tau_{ST}$ | $3.1\cdot 10^{-35}$ | s
Suspension Thermal Noise Cutoff Frequency | $\omega_{ST}$ | $2\pi\cdot 1.9\cdot 10^{3}$ | rad/s
Table 3: aLIGO Parameters Parameter | Symbol | Value | Units
---|---|---|---
Mechanical Resonant Frequency | $\omega_{m}$ | 2$\pi\cdot$ 0.9991 | rad/s
Mirror Mass | $M$ | 9.446 | kg
Cavity Decay Rate | $\gamma$ | 2$\pi\cdot$ 424.6 | rad/s
Arm Length | $L$ | 3.995 | km
Circulating Power | $P_{c}$ | 322.7 | kW
Laser Wavelength | $\lambda$ | 1064 | nm
Mechanical Damping | $\gamma_{m}$ | 2$\pi\cdot 10^{-3}$ | rad/s
Table 4: Resonant Mode Parameters Parameter | Symbol | Value | Units
---|---|---|---
Mode Frequency | $\Omega_{v_{1}}$ | 2$\pi\cdot 0.441$ | rad/s
$\Omega_{v_{2}}$ | 2$\pi\cdot 0.995$ | rad/s
$\Omega_{v_{3}}$ | 2$\pi\cdot 1.98$ | rad/s
$\Omega_{v_{4}}$ | 2$\pi\cdot 2.37$ | rad/s
$\Omega_{v_{5}}$ | 2$\pi\cdot 3.38$ | rad/s
$\Omega_{v_{6}}$ | 2$\pi\cdot 3.81$ | rad/s
$\Omega_{v_{7}}$ | 2$\pi\cdot 9.73$ | rad/s
Full Width at Half Maximum | $\Gamma{v_{1}}$ | 2$\pi\cdot 1.92\cdot 10^{-3}$ | rad/s
$\Gamma{v_{2}}$ | 2$\pi\cdot 5.63\cdot 10^{-5}$ | rad/s
$\Gamma{v_{3}}$ | 2$\pi\cdot 2.11\cdot 10^{-5}$ | rad/s
$\Gamma{v_{4}}$ | 2$\pi\cdot 1.44\cdot 10^{-1}$ | rad/s
$\Gamma{v_{5}}$ | 2$\pi\cdot 1.45\cdot 10^{-4}$ | rad/s
$\Gamma{v_{6}}$ | 2$\pi\cdot 1.65\cdot 10^{-3}$ | rad/s
$\Gamma{v_{7}}$ | 2$\pi\cdot 1.03\cdot 10^{-3}$ | rad/s
Amplitude | $A_{v_{1}}$ | 159 | rad/s
$A_{v_{2}}$ | 93.8 | rad/s
$A_{v_{3}}$ | 538 | rad/s
$A_{v_{4}}$ | 235 | rad/s
$A_{v_{5}}$ | 353 | rad/s
$A_{v_{6}}$ | 27.4 | rad/s
$A_{v_{7}}$ | 78.0 | rad/s

##  Appendix C Noise Model

Our model for the total classical and quantum noise in aLIGO is displayed in
Fig. C.1. Note that the coating Brownian noise, which is the dominating
classical noise source above $\sim 10$ Hz, is modeled in a piecewise manner: a
white noise in the $10-10^{5}$ Hz band and a noise source increasing as
$\Omega^{2}$ for frequencies larger than $10^{5}$ Hz. Even though the sensing
noise in the $10-10^{5}$ Hz band decreases as $\Omega^{0.8}$ in the power
spectral density, we choose to model it with a frequency-independent white
noise for simplified calculations. Our choice is justified since the quantum
noise dominates over the coating Brownian noise in that region.

Figure C.1:  Modeled noise spectrum for aLIGO.

##  Appendix D Structure of the Covariance Matrix

The mirror and the cavity, at $t=0$, constitute two modes, whereas there is an
infinite number of modes in the outgoing light field given by the quadratures
$\hat{v}_{1}(t),\hat{v}_{2}(t)$, $t\in(-\infty,0)$. In continuum coordinates,
we can first write down the commutator matrix,

$\mathbf{K}=\left[\begin{array}[]{ccc}\mathbf{K}^{B}\\\ &\mathbf{K}^{A}\\\
&&\mathbf{K}^{v}\end{array}\right]\,,$ (24)

with

$\displaystyle\mathbf{K}^{A}$
$\displaystyle=\mathbf{K}^{B}=\left[\begin{array}[]{cc}0&i\\\
-i&0\end{array}\right]\,,$ (27) $\displaystyle\mathbf{K}^{v}$
$\displaystyle=\left[\begin{array}[]{cc}0&i\delta(t-t^{\prime})\\\
-i\delta(t-t^{\prime})&0\end{array}\right]\,.$ (30)

Note that $\mathbf{K}^{v}$ is a $2\times 2$ block matrix, but each block is
infinite dimensional, with columns and rows indexed by $t$ and $t^{\prime}$,
respectively. The indices $t$ and $t^{\prime}$ each run through all negative
real numbers, $(-\infty,0]$. In order to represent these covariance matrices
in a less ambiguous and more operational way, we shall adopt an index
notation, in which $j,k,l,m$ are discrete and run through 1 and 2, while $t$
and $t^{\prime}$ are continuous and run through negative real numbers. We can
then write

$K^{A}_{jk}=i\epsilon_{jk}\,,\;K^{B}_{lm}=i\epsilon_{lm}\,,\;K^{v}_{lt,mt^{\prime}}=i\epsilon_{lm}\delta(t-t^{\prime})\,.$
(31)

Note that for $\mathbf{K}^{v}$ we have a two-dimensional row index $(l,t)$, as
well as a two-dimensional column index $(m,t^{\prime})$, to label quadrature
and arrival time.

The covariance matrix $\mathbf{V}$ of the system contains the steady state
correlations between the mechanical mode at $t=0$, the cavity mode at $t=0$,
and the outgoing light field modes for $t<0$. It can be written as

$\mathbf{V}=\left[\begin{array}[]{cc|c}V^{BB}&V^{BA}&V^{Bv}\\\
V^{AB}&V^{AA}&V^{Av}\\\ \hline\cr
V^{vB}&V^{vA}&V^{vv}\end{array}\right]\equiv\left[\begin{array}[]{cc}V^{QQ}&V^{Qv}\\\
V^{vQ}&V^{vv}\end{array}\right]\,.$ (32)

Here each of $A$ and $B$ represent two dimensions, while $v$ represents an
infinite number of dimensions. For computational purposes we can group $A$ and
$B$ together as $Q$, or

$\hat{Q}_{J}=\hat{Q}_{(1,2,3,4)}=(\hat{B}_{1},\hat{B}_{2},\hat{A}_{1},\hat{A}_{2})\,.$
(33)

Here we shall use upper-case Latin indices to run from 1 to 4, to label
quadratures in the mechanical and cavity modes. We can write

$V^{QQ}_{JK}=\frac{1}{2}\langle\hat{Q}_{J}\hat{Q}_{K}+\hat{Q}_{K}\hat{Q}_{J}\rangle=\frac{1}{2}\int_{-\infty}^{+\infty}\frac{d\Omega}{2\pi}S_{Q_{J}Q_{K}}(\Omega)\,,$
(34)

where the (one-sided) cross spectrum is defined via

$\frac{1}{2}S_{XY}(\Omega)\delta(\Omega-\Omega^{\prime})=\frac{1}{2}\langle\hat{X}(\Omega)\hat{Y}^{\dagger}(\Omega^{\prime})+\hat{Y}^{\dagger}(\Omega^{\prime})\hat{X}(\Omega)\rangle\,,$
(35)

and they can be obtained from solutions to the Heisenberg Equations, as well
as the spectra

$S_{u_{i}u_{j}}(\Omega)=\delta_{ij}\,,$ (36)

and the prescriptions we use for the spectra of $S_{n_{X}}(\Omega)$ and
$S_{n_{F}}(\Omega)$. The uncorrelated white spectra between $\hat{u}_{1}$ and
$\hat{u}_{2}$ in Eq. (36) result from the ingoing light field being in its
vacuum state. We shall discuss the magnitude and frequency dependence of
$S_{n_{X}}(\Omega)$ and $S_{n_{F}}(\Omega)$ in depth in the next section.

For elements that involve $v$, we shall still use $lt$ for column indices and
$mt^{\prime}$ for row indices. We then have

$\displaystyle V^{Qv}_{J,mt^{\prime}}$
$\displaystyle=\frac{1}{2}\langle\hat{Q}_{J}\hat{v}_{m}(t^{\prime})+\hat{v}_{m}(t^{\prime})\hat{Q}_{J}\rangle$
$\displaystyle=\frac{1}{2}\int\limits_{-\infty}^{+\infty}\frac{d\Omega}{2\pi}S_{Q_{J}v_{m}}(\Omega)e^{i\Omega
t^{\prime}}\,,$ (37)

and

$\displaystyle V^{vv}_{lt,mt^{\prime}}$
$\displaystyle=\frac{1}{2}\langle\hat{v}_{l}(t)\hat{v}_{m}(t^{\prime})+\hat{v}_{m}(t)\hat{v}_{l}(t^{\prime})\rangle$
$\displaystyle=\frac{1}{2}\int\limits_{-\infty}^{+\infty}\frac{d\Omega}{2\pi}S_{v_{l}v_{m}}(\Omega)e^{-i\Omega(t-t^{\prime})}\,.$
(38)

Note that $V^{Qv}_{J,mt^{\prime}}=V^{vQ}_{mt^{\prime},J}$.

In numerical computations, we will have to convert the continuum of
$t,t^{\prime}\leq 0$ into a finite grid. This means we will sample a finite
duration $T$ with a step size of $\Delta t$. We shall still use lower-case
Latin indices to run from 1 to 2, and upper case Latin indices to run from 1
to 4, while use Greek indices, for example, $\alpha=0,1,2,\ldots,T/\Delta
t\equiv\mathcal{N}-1$ to replace $t$. We shall write

$K_{lt,mt^{\prime}}^{v}\rightarrow
K^{v}_{l\alpha,m\alpha^{\prime}}=i\epsilon_{lm}\delta_{\alpha\alpha^{\prime}}\,.$
(39)

Note that a Kronecker delta now replaces the Dirac delta. For the covariance
matrix, we replace

$\displaystyle V_{J,mt^{\prime}}^{Qv}\rightarrow V_{J,m\alpha^{\prime}}^{Qv}$
$\displaystyle=\frac{\sqrt{\Delta
t}}{2}\langle\hat{Q}_{J}\hat{v}_{m}(t_{\alpha^{\prime}})+\hat{v}_{m}(t_{\alpha^{\prime}})\hat{Q}_{J}\rangle\,,$
(40)

which are $4\times 2\mathcal{N}$ and $2\mathcal{N}\times 4$ dimensional
matrices (with $V_{J,mt^{\prime}}^{vQ}$), and

$\displaystyle V_{lt,mt^{\prime}}^{vv}\rightarrow
V_{l\alpha,m\alpha^{\prime}}^{vv}$ $\displaystyle=\frac{\Delta
t}{2}\langle\hat{v}_{l}(t_{\alpha})\hat{v}_{m}(t_{\alpha^{\prime}})+\hat{v}_{m}(t_{\alpha})\hat{v}_{l}(t_{\alpha^{\prime}})\rangle\,,$
(41)

which is a $2\mathcal{N}\times 2\mathcal{N}$-dimensional matrix. For the
discrete sampling times we have defined

$t_{\alpha}=-\left(\alpha+\frac{1}{2}\right)\Delta t\,,$ (42)

where the additional $\frac{1}{2}\Delta t$ provides a faster convergence in
numerics. The entire covariance matrix is then
$(2\mathcal{N}+4)\times(2\mathcal{N}+4)$ \- dimensional.

Our particular convention of inserting $\Delta t$ at various places of the
matrix is associated with our convention of discretizing vectors. For a
generic variable, in the continuum form, we can always express it as

$\displaystyle X$
$\displaystyle=\alpha^{j}\hat{A}_{j}+\beta^{j}\hat{B}_{j}+\int_{-\infty}^{0}\xi^{j}(t)\hat{v}_{j}(t)dt$
$\displaystyle=\gamma^{J}\hat{Q}_{J}+\int_{-\infty}^{0}\xi^{j}(t)\hat{v}_{j}(t)dt\,,$
(43)

where we have used upper indices for vector components, and have grouped
$\alpha^{j}$ and $\beta^{j}$ into $\gamma^{J}$. The variance of $X$, which is
formally written as $\mathbf{X}^{\dagger}\mathbf{V}\mathbf{X}$, will then be

$\displaystyle\mathbf{X}^{\dagger}\mathbf{V}\mathbf{X}$
$\displaystyle=\frac{1}{2}\gamma^{J}\langle\hat{Q}_{J}\hat{Q}_{K}+\hat{Q}_{K}\hat{Q}_{J}\rangle\gamma^{K}$
$\displaystyle+\int_{-\infty}^{0}\gamma^{J}\langle\hat{Q}_{J}\hat{v}_{m}(t^{\prime})+\hat{v}_{m}(t^{\prime})\hat{Q}_{J}\rangle\xi^{m}(t^{\prime})dt^{\prime}$
$\displaystyle+\frac{1}{2}\iint_{-\infty}^{0}\xi^{l}(t)\langle\hat{v}_{l}(t)\hat{v}_{m}(t^{\prime})+\hat{v}_{m}(t)\hat{v}_{l}(t^{\prime})\rangle\xi^{m}(t^{\prime})dtdt^{\prime}\,.$
(44)

As we convert the integrals in Eq. (D) into summations, $\int$ will become
$\Sigma$, while $dt$ will become $\Delta t$. We shall take

$\xi^{m\alpha}=\xi^{m}(t_{\alpha})\sqrt{\Delta t}.$ (45)

Together with Eq. (40)–(41), the fully discretized version of Eq. (D) will
then be

$\displaystyle\mathbf{X}^{\dagger}\mathbf{V}\mathbf{X}=\gamma^{J}V^{QQ}_{JK}\gamma^{K}+\gamma^{J}V^{Qv}_{J,m\alpha}\xi^{m\alpha}$
$\displaystyle+\xi^{l\beta}V^{vQ}_{l\beta,K}\gamma^{K}$
$\displaystyle+\xi^{l\alpha}V^{vv}_{l\alpha,m\beta}\xi^{m\beta}.$ (46)

In this convention, the usual vector norm for the discretized version of a
function of time coincides with the $L^{2}$-norm of that function. It can also
be checked that discretized matrices in Eqs. (39)–(41), when contracted with
vectors in this convention, lead to the appropriate integrals. Note that if a
$\delta(t_{\alpha}-t_{\alpha}^{\prime})$ shows up in Eq. (41), we will take
$\Delta
t\delta(t_{\alpha}-t_{\alpha}^{\prime})\rightarrow\delta_{\alpha\alpha^{\prime}}$,
as in Eq. (39).

Corresponding to the discussion at the end of Sec. III (also shown in Fig. 1),
here we consider entanglement between: (i) mass at $t=0$ and the out-going
light field that had emerged during $t\leq 0$ and (ii) mass and the joint
system of the cavity mode as well as light that had emerged during $t\leq 0$.
In case (i), we simply throw away elements involving $A$ in both $\mathbf{K}$
and $\mathbf{V}$, while in case (ii) we consider the full matrices. In both
cases, $\mathbf{V}_{\rm pt}$ is obtained by adding a minus sign to the column
involving $\hat{B}_{2}$ and the row involving $\hat{B}_{2}$ – but not the
diagonal element at which they intersect.

##  Appendix E Numerical Implementation for aLIGO’s Noise

In our simulations, we use $dt=0.25$ ms and $T=0.1$ s, which corresponds to
sampling the light field at 4000 Hz and working with the outgoing field
emitted from the cavity between $t=-0.1$ s and $t=0$ s. We achieve numerical
convergence with these parameters. To quantify the amount of entanglement in
the system, we use the logarithmic negativity defined in Eq. (14). However,
this is possible only for low levels of classical noise. For high levels of
classical noise, classical correlations dominate over quantum correlations,
which causes the cross-correlation values in the system to cover a wide-range
of orders of magnitude, mostly due to the 14th power of $\Omega$ in our force
noise model Eq. (18). For aLIGO parameters, the entries of the covariance
matrix span about 20 orders of magnitudes, while we attempt to find a
symplectic eigenvalue of order 1 – this is numerically an extremely
challenging problem.

Figure E.1:  $\lambda_{\mathcal{B}}$ and $\lambda_{\mathcal{N}}$ for
$\alpha_{F_{1}}=10^{-15}$ and $\alpha_{F_{2}}=15$ (top row, Figure a),
$\alpha_{F_{1}}=10^{-8}$, $\alpha_{F_{2}}=10$, and $\alpha_{X_{2}}=10^{3}$
(bottom row, Figure b). We set $T=0.1$ s, and $dt\in[0.1,10]$ ms. Note that
$dt$ decreases to the right of the plots. In Figure a, entanglement exists for
both partitions, and $\lambda_{\mathcal{B}}$($\lambda_{\mathcal{N}}$) is re-
plotted in the inset on the left(right) for $dt\in[0.1,0.5]$ ms. In Figure b,
the state is separable for both partitions, and $\lambda_{\mathcal{B}}$ and
$\lambda_{\mathcal{N}}$ are re-plotted in the inset for $dt\in[0.1,0.5]$ ms.

Numerical errors also arise because of time-binning with an insufficient
resolution. Thus, we lose precision on the numerically determined covariance
matrices. This affects the smallest symplectic eigenvalue $\tilde{\nu}_{min}$
to the point that it cannot be used to measure entanglement with the
logarithmic negativity. Numerical imprecision can lead to covariance matrices
that do not satisfy Heisenberg uncertainty bound Eq. (12), they thus do not
correspond to a bona fide state and we call them _non-physical_. One way to
get around this loss of precision is to use the PPT criterion as a yes/no test
only, renouncing the magnitude information of $\tilde{\nu}_{min}$. The
sampling frequency during time binning should be higher than the Nyquist rate
of the system (i.e. twice the largest frequency in the system), since the
entries of the covariance matrix contain correlations from all frequencies. In
our system, the largest frequency is the cavity decay rate $\gamma=424$ Hz.
Therefore, we choose $dt<1/(2\cdot 424)\approx 1.2$ ms. However, as we
decrease $dt$, we are limited by computational resources, such as the RAM
size, or time. The parameter $dt$ is limited to an optimal range determined by
this trade-off. We thus develop the following strategy: we first quantify the
amount of numerical errors in the system by computing the most negative
eigenvalue (if it exists) of $\mathbf{V}+\frac{1}{2}\mathbf{K}$ before and
after the partial transpose operation, denoted as $\lambda_{\mathcal{B}}$ and
$\lambda_{\mathcal{N}}$ respectively. Then, we decide that entanglement exists
if $\lambda_{\mathcal{B}}>0$ and $\lambda_{\mathcal{N}}<0$, or if
$\lambda_{\mathcal{B}}<0$ and
$|\lambda_{\mathcal{N}}|\gg|\lambda_{\mathcal{B}}|$. Furthermore, we decide
that the system is separable if $\lambda_{\mathcal{B}}>0$ and
$\lambda_{\mathcal{N}}>0$.

Two case studies about this strategy can be found in Fig. E.1 where we fix
$T=0.1$ s, change $dt$, and examine how $\lambda_{\mathcal{N}}$ and
$\lambda_{\mathcal{B}}$ change by plotting $\lambda_{\mathcal{N}}$ and
$\lambda_{\mathcal{B}}$. We decide on entanglement if
$\lambda_{\mathcal{B}}\geq 0$ and $\lambda_{\mathcal{N}}<0$, or
$\lambda_{\mathcal{B}}<0$, $\lambda_{\mathcal{N}}<0$, and
$|\lambda_{\mathcal{N}}|\geq 100|\lambda_{\mathcal{B}}|$. Our criteria for
convergence is a relative change smaller than $5\%$ for both
$\lambda_{\mathcal{N}}$ and $\lambda_{\mathcal{B}}$ as we change $dt$. In Fig.
E.1a, we work with a low level of classical force noise and set
$\alpha_{F_{1}}=10^{-15}$, $\alpha_{F_{2}}=15$ in Eq. (18). We see that
$\lambda_{\mathcal{N}}$ and $\lambda_{\mathcal{B}}$ converge with
$\lambda_{\mathcal{N}}$ changing by $0.053\%$, $0.068\%$, and
$\lambda_{\mathcal{B}}$ changing by $4.9\%$, $0.77\%$ before and after tracing
over the cavity respectively, for $dt=0.1$ ms. The system is entangled for
both partitions since $\lambda_{\mathcal{B}}\geq 0$ and
$\lambda_{\mathcal{N}}<0$. Furthermore, $\lambda_{\mathcal{B}}$ becomes
positive and converges after $dt\sim 1$ ms, or a sampling frequency of 1000
Hz; consistent with the discussion above relating physicality to Nyquist rate
of $\sim 850$ Hz. We also see that $\lambda_{\mathcal{N}}$ converges for
similar values of $dt$ from Fig. E.1a.

Figure E.2: The real and imaginary parts of the eigenvector for the negative
eigenvalue of $\mathbf{V}_{\rm pt}+\frac{1}{2}\mathbf{K}$ for for
$\alpha_{F_{1}}=10^{-15}$, $\alpha_{F_{2}}=15$, denoted by $e_{1}$ and
$e_{2}$, which correspond to the first and second halves of the entire
eigenvector respectively. The reason behind this slicing is the block-matrix
structure of the light-field sector of the covariance matrix. Furthermore,
since the light field modes are continuous in time, $e_{1}$ and $e_{2}$ are
also functions of time.

In Fig. E.2, we plot the light-field section of the eigenvector associated to
the (converged) minimal eigenvalue of $\mathbf{V}_{\rm
pt}+\frac{1}{2}\mathbf{K}$, for the partition where we do not trace over the
cavity and with $\alpha_{F_{1}}=10^{-15}$, $\alpha_{F_{2}}=15$. It corresponds
to a temporal mode of the free electromagnetic field outside the cavity. It is
that particular mode associated to the (sole) negative eigenvalue that is
entangled with the joint system mechanics plus cavity. The four curves
correspond to the real and imaginary parts of the $\hat{v}_{1}(t)$ and
$\hat{v}_{2}(t)$. They have the form of smooth decaying oscillations with the
same frequency and decay rate, but differing by a phase. This form of the mode
functions was predicted for a white force noise in [6]; which gives us
confidence in the correctness of our study. Also, exponentially decaying
demodulation pulses were used to demonstrate optomechanical entanglement [16,
8] and proposed for a demonstration in the stationary regime [23]. We fit
functions in the form of $e^{-\gamma_{*}t}\sin{(\omega_{*}t+\theta)}$ to each
curve, which results in $\omega_{*}/(2\pi)\approx 40$ Hz, and
$\gamma_{*}/(2\pi)\approx 25$ Hz. In the frequency domain, exponentially
decaying harmonic oscillations are Lorentzians, centered at $\pm\omega_{*}$
and with a bandwidth (FWHM) $2\gamma_{*}$. In aLIGO’s noise budget (Fig. 3),
these Lorentzians are on the low frequency side of the low-noise band and
their HWHM to the left crosses the quantum noise, where it is not yet
dominated by suspension thermal and seismic noises – although we saw in Sec.
VI.1 that the latter is probably the main mechanism preventing optomechanical
entanglement. We add that Lorentzians are heavy-tail distributions, being a
possible reason why even lower frequency components matter.

In Fig. E.1b, we set $\alpha_{F_{1}}=10^{-8}$, $\alpha_{F_{2}}=10$, and
$\alpha_{X_{2}}=10^{3}$, causing the sensing noise to dominate over quantum
noise for frequencies in the $30-2000$ Hz band. We again see that
$\lambda_{\mathcal{N}}$ and $\lambda_{\mathcal{B}}$ converge with
$\lambda_{\mathcal{N}}$ changing by $1.3\%$, $1.2\%$, and
$\lambda_{\mathcal{B}}$ changing by $1.2\%$, $1.1\%$ before and after tracing
over the cavity respectively, for $dt=0.1$ ms. Since $\lambda_{\mathcal{B}}$
and $\lambda_{\mathcal{N}}$ are positive for both partitions, we conclude that
there is no entanglement in the system for either partition. When we increase
the classical noise level we see that convergence is much harder to achieve.
Furthermore, $\lambda_{\mathcal{B}}$ and $\lambda_{\mathcal{N}}$ are negative,
and $\lambda_{\mathcal{B}}\sim\lambda_{\mathcal{N}}$ for every value of $dt$
regardless of when we do or do not trace over the cavity. Therefore, we cannot
conclude that there is entanglement for neither of the partitions. These case
studies show that we can use $\lambda_{\mathcal{B}}$ as an indicator of the
“non-physicality” of the covariance matrix $\mathbf{V}$ introduced by finite-
sampling and high levels of classical noise in the system, and that the
negativity of $\lambda_{\mathcal{N}}$ is not enough to decide on entanglement
when we consider the numerics.

Figure E.3:  Contour plot depicting the effect of the force noise spectrum on
our numerical precision for both partitions. The force noise levels increase
toward the bottom-right of the figure. The black and the red dashed lines
separate the regions where numerics converge from the regions where numerics
fail for the partition where we do and do not trace over the cavity
respectively. The region where numerics converge for both partitions is marked
in gray, whereas numerics fail for both partitions toward the bottom right of
the figure, past the red dashed line, in the yellow region.

For our model of aLIGO’s non-Markovian noises, Eqs. (18), we study the
numerical stability of broad noise regimes, parametrized by the pair
$\alpha_{F_{j}}$, $j=1,2$. We set $\alpha_{X_{j}}=1$, $j=1,2$, since we saw
that force noise had a greater impact on numerical stability than sensing
noise in our simulations for aLIGO’s noise. In Fig. E.3, we plot the boundary
between noise regimes where the numerics converge and the computed covariance
matrices are physical (in gray) and those regimes where the numerics fail
(either at converging or at generating physical covariance matrices or both)
with the available computing resources (in yellow). As a matter of fact, in
all the operation regimes in the gray region where the numerics work, the
system is entangled. This means that our model predicts optomechanical
entanglement in aLIGO, if its classical noise is in this gray region. Recall
that the current status of aLIGO corresponds to
$\alpha_{F_{j}}=\alpha_{X_{j}}=1$, $j=1,2$, far to the bottom-right in the
undecidable yellow region.

If we follow the red dashed line in Fig. E.3, we continuously sample the force
noise spectrum over the boundary where our numerics converge, and the system
is entangled, for the partition where we do not trace over the cavity. This
boundary corresponds to a lower limit of the maximum amount of classical noise
allowed in the system in order to have entanglement. We plot some of these
noise curves in Fig. E.4. Note that if we set $\alpha_{F_{1}}$ or
$\alpha_{F_{2}}$ to be unity, the corresponding pairs of $\alpha_{F_{1}}$,
$\alpha_{F_{2}}$ situated on this boundary would be $\alpha_{F_{1}}=1$ and
$\alpha_{F_{2}}\geq 200$, or for $\alpha_{F_{2}}=1$ and $\alpha_{F_{1}}\leq
10^{-10}$. The corresponding noise curves are also plotted on Fig. E.4.

Figure E.4: Examples of force noise spectra along the boundary where we lose
numerical precision (defined with
$|\lambda_{\mathcal{N}}|=100|\lambda_{\mathcal{B}}|$) for the partition where
we do not trace over the cavity, plotted in black. aLIGO noise curves are also
plotted for comparison. The black curves are parametrized with (from left to
right) $\alpha_{F_{1}}=1$ and $\alpha_{F_{2}}=200$, $\alpha_{F_{1}}=10^{-3}$
and $\alpha_{F_{2}}=14$, $\alpha_{F_{1}}=10^{-10}$ and $\alpha_{F_{2}}=1$, and
$\alpha_{F_{1}}=10^{-15}$ and $\alpha_{F_{2}}=0.029$.

## References

  * Genes  _et al._ [2009] C. Genes, A. Mari, D. Vitali, and P. Tombesi, Chapter 2 quantum effects in optomechanical systems, in  _Advances in Atomic Molecular and Optical Physics_, Advances In Atomic, Molecular, and Optical Physics, Vol. 57 (Academic Press, 2009) pp. 33–86.
  * Hofer and Hammerer [2015] S. G. Hofer and K. Hammerer, Entanglement-enhanced time-continuous quantum control in optomechanics, Phys. Rev. A 91, 033822 (2015).
  * Genes  _et al._ [2008] C. Genes, A. Mari, P. Tombesi, and D. Vitali, Robust entanglement of a micromechanical resonator with output optical fields, Phys. Rev. A 78, 032316 (2008).
  * Miao  _et al._ [2010a] H. Miao, S. Danilishin, H. Müller-Ebhardt, and Y. Chen, Achieving ground state and enhancing optomechanical entanglement by recovering information, New Journal of Physics 12, 083032 (2010a).
  * Miao  _et al._ [2010b] H. Miao, S. Danilishin, H. Müller-Ebhardt, H. Rehbein, K. Somiya, and Y. Chen, Probing macroscopic quantum states with a sub-heisenberg accuracy, Phys. Rev. A 81, 012114 (2010b).
  * Miao  _et al._ [2010c] H. Miao, S. Danilishin, and Y. Chen, Universal quantum entanglement between an oscillator and continuous fields, Physical Review A 81, 10.1103/physreva.81.052307 (2010c).
  * Danilishin  _et al._ [2019] S. L. Danilishin, F. Y. Khalili, and H. Miao, Advanced quantum techniques for future gravitational-wave detectors, Living Reviews in Relativity 22, 10.1007/s41114-019-0018-y (2019).
  * Hofer  _et al._ [2011] S. G. Hofer, W. Wieczorek, M. Aspelmeyer, and K. Hammerer, Quantum entanglement and teleportation in pulsed cavity optomechanics, Phys. Rev. A 84, 052327 (2011).
  * Aspelmeyer  _et al._ [2014] M. Aspelmeyer, T. J. Kippenberg, and F. Marquardt, Cavity optomechanics, Reviews of Modern Physics 86, 1391 (2014).
  * Barzanjeh  _et al._ [2019] S. Barzanjeh, E. S. Redchenko, M. Peruzzo, M. Wulf, D. P. Lewis, G. Arnold, and J. M. Fink, Stationary entangled radiation from micromechanical motion, Nature 570, 480 (2019).
  * Chen  _et al._ [2020] J. Chen, M. Rossi, D. Mason, and A. Schliesser, Entanglement of propagating optical modes via a mechanical interface, Nature Communications 11, 943 (2020).
  * Riedinger  _et al._ [2018] R. Riedinger, A. Wallucks, I. Marinković, C. Löschnauer, M. Aspelmeyer, S. Hong, and S. Gröblacher, Remote quantum entanglement between two micromechanical oscillators, Nature 556, 473 (2018).
  * Ockeloen-Korppi  _et al._ [2018] C. F. Ockeloen-Korppi, E. Damskägg, J.-M. Pirkkalainen, M. Asjad, A. A. Clerk, F. Massel, M. J. Woolley, and M. A. Sillanpää, Stabilized entanglement of massive mechanical oscillators, Nature 556, 478 (2018).
  * Kotler  _et al._ [2021] S. Kotler, G. A. Peterson, E. Shojaee, F. Lecocq, K. Cicak, A. Kwiatkowski, S. Geller, S. Glancy, E. Knill, R. W. Simmonds, J. Aumentado, and J. D. Teufel, Direct observation of deterministic macroscopic entanglement, Science 372, 622 (2021).
  * Thomas  _et al._ [2021] R. A. Thomas, M. Parniak, C. Østfeldt, C. B. Møller, C. Bærentsen, Y. Tsaturyan, A. Schliesser, J. Appel, E. Zeuthen, and E. S. Polzik, Entanglement between distant macroscopic mechanical and spin systems, Nature Physics 17, 228 (2021).
  * Palomaki  _et al._ [2013] T. A. Palomaki, J. D. Teufel, R. W. Simmonds, and K. W. Lehnert, Entangling mechanical motion with microwave fields, Science 342, 710 (2013), https://www.science.org/doi/pdf/10.1126/science.1244563 .
  * Romero-Isart [2011] O. Romero-Isart, Quantum superposition of massive objects and collapse models, Phys. Rev. A 84, 052121 (2011).
  * Müller-Ebhardt  _et al._ [2008] H. Müller-Ebhardt, H. Rehbein, R. Schnabel, K. Danzmann, and Y. Chen, Entanglement of macroscopic test masses and the standard quantum limit in laser interferometry, Phys. Rev. Lett. 100, 013601 (2008).
  * Schnabel [2015] R. Schnabel, Einstein-podolsky-rosen–entangled motion of two massive objects, Phys. Rev. A 92, 012126 (2015).
  * Wang  _et al._ [2016] M. Wang, X.-Y. Lü, Y.-D. Wang, J. Q. You, and Y. Wu, Macroscopic quantum entanglement in modulated optomechanics, Phys. Rev. A 94, 053807 (2016).
  * Pirandola  _et al._ [2006] S. Pirandola, D. Vitali, P. Tombesi, and S. Lloyd, Macroscopic entanglement by entanglement swapping, Phys. Rev. Lett. 97, 150403 (2006).
  * Ludwig  _et al._ [2010] M. Ludwig, K. Hammerer, and F. Marquardt, Entanglement of mechanical oscillators coupled to a nonequilibrium environment, Phys. Rev. A 82, 012333 (2010).
  * Gut  _et al._ [2020] C. Gut, K. Winkler, J. Hoelscher-Obermaier, S. G. Hofer, R. M. Nia, N. Walk, A. Steffens, J. Eisert, W. Wieczorek, J. A. Slater, M. Aspelmeyer, and K. Hammerer, Stationary optomechanical entanglement between a mechanical oscillator and its measurement apparatus, Phys. Rev. Res. 2, 033244 (2020).
  * Vitali  _et al._ [2007] D. Vitali, S. Gigan, A. Ferreira, H. R. Böhm, P. Tombesi, A. Guerreiro, V. Vedral, A. Zeilinger, and M. Aspelmeyer, Optomechanical entanglement between a movable mirror and a cavity field, Phys. Rev. Lett. 98, 030405 (2007).
  * Paternostro  _et al._ [2007] M. Paternostro, D. Vitali, S. Gigan, M. S. Kim, C. Brukner, J. Eisert, and M. Aspelmeyer, Creating and probing multipartite macroscopic entanglement with light, Phys. Rev. Lett. 99, 250401 (2007).
  * Giovannetti and Vitali [2001] V. Giovannetti and D. Vitali, Phase-noise measurement in a cavity with a movable mirror undergoing quantum brownian motion, Phys. Rev. A 63, 023812 (2001).
  * Fedorov  _et al._ [2018] S. Fedorov, V. Sudhir, R. Schilling, H. Schütz, D. Wilson, and T. Kippenberg, Evidence for structural damping in a high-stress silicon nitride nanobeam and its implications for quantum optomechanics, Physics Letters A 382, 2251 (2018), special Issue in memory of Professor V.B. Braginsky.
  * Gröblacher  _et al._ [2015] S. Gröblacher, A. Trubarov, N. Prigge, G. D. Cole, M. Aspelmeyer, and J. Eisert, Observation of non-markovian micromechanical brownian motion, Nature Communications 6, 7606 (2015).
  * Aasi  _et al._ [2015] J. Aasi _et al._ , Advanced ligo, Classical and Quantum Gravity 32, 074001 (2015).
  * Yu and members of the LIGO Scientific Collaboration [2020] H. Yu and members of the LIGO Scientific Collaboration, Quantum correlations between light and the kilogram-mass mirrors of ligo, Nature 583, 43 (2020).
  * Note [1] The quantum noise’s beating of the SQL is only inferred by subtracting a classical noise floor that was obtained through calibration.
  * Buonanno and Chen [2003] A. Buonanno and Y. Chen, Scaling law in signal recycled laser-interferometer gravitational-wave detectors, Physical Review D 67, 10.1103/physrevd.67.062002 (2003).
  * Chen [2013] Y. Chen, Macroscopic quantum mechanics: theory and experimental concepts of optomechanics, Journal of Physics B: Atomic, Molecular and Optical Physics 46, 104001 (2013).
  * Note [2] We use the convention $\mathcal{F}\\{f(t)\\}=\intop\ilimits@_{-\infty}^{\infty}f(t)e^{i\omega t}dt$.
  * Caves and Schumaker [1985] C. M. Caves and B. L. Schumaker, New formalism for two-photon quantum optics. i. quadrature phases and squeezed states, Phys. Rev. A 31, 3068 (1985).
  * Schumaker and Caves [1985] B. L. Schumaker and C. M. Caves, New formalism for two-photon quantum optics. ii. mathematical foundation and compact notation, Phys. Rev. A 31, 3093 (1985).
  * Danilishin and Khalili [2012] S. L. Danilishin and F. Y. Khalili, Quantum measurement theory in gravitational-wave detectors, Living Reviews in Relativity 15, 10.12942/lrr-2012-5 (2012).
  * Kubo [1966] R. Kubo, The fluctuation-dissipation theorem, Reports on progress in physics 29, 255 (1966).
  * Kimble  _et al._ [2001] H. J. Kimble, Y. Levin, A. B. Matsko, K. S. Thorne, and S. P. Vyatchanin, Conversion of conventional gravitational-wave interferometers into quantum nondemolition interferometers by modifying their input and/or output optics, Physical Review D 65, 10.1103/physrevd.65.022002 (2001).
  * Peres [1996] A. Peres, Separability criterion for density matrices, Phys. Rev. Lett. 77, 1413 (1996).
  * Adesso and Illuminati [2007] G. Adesso and F. Illuminati, Entanglement in continuous-variable systems: recent advances and current perspectives, Journal of Physics A: Mathematical and Theoretical 40, 7821 (2007).
  * Werner and Wolf [2001] R. F. Werner and M. M. Wolf, Bound entangled gaussian states, Phys. Rev. Lett. 86, 3658 (2001).
  * Simon [2000] R. Simon, Peres-horodecki separability criterion for continuous variable systems, Phys. Rev. Lett. 84, 2726 (2000).
  * Vidal and Werner [2002] G. Vidal and R. F. Werner, Computable measure of entanglement, Phys. Rev. A 65, 032314 (2002).
  * Serafini [2017] A. Serafini,  _Quantum Continuous Variables_ (CRC Press, 2017).
  * Braginsky  _et al._ [1992] V. B. Braginsky, F. Y. Khalili, and K. S. Thorne,  _Quantum Measurement_ (Cambridge University Press, 1992).
  * Saulson [1990] P. R. Saulson, Thermal noise in mechanical experiments, Phys. Rev. D 42, 2437 (1990).
  * González and Saulson [1995] G. I. González and P. R. Saulson, Brownian motion of a torsion pendulum with internal friction, Physics Letters A 201, 12 (1995).
  * Neben  _et al._ [2012] A. R. Neben, T. P. Bodiya, C. Wipf, E. Oelker, T. Corbitt, and N. Mavalvala, Structural thermal noise in gram-scale mirror oscillators, New Journal of Physics 14, 115008 (2012).
  * Note [3] The spectrum of the mechanical displacement per unit Hz, where the transduction of the force and sensing noises by the cavity is effectively factored out.
  * Note [4] To find spectra of the form of Eqs. (16), one assumes the free-mass and high-$Q$ limits where the mechanical susceptibility behaves as $1/\Omega^{2}$; additionally, one must be in frequency ranges where the noise spectra $S_{n_{F}}$ and $S_{n_{X}}$ are constant polynomial roll-off – typical at frequencies above the resonance of suspended oscillators.
  * Somiya  _et al._ [2012] K. Somiya _et al._ , Detector configuration of kagra–the japanese cryogenic gravitational-wave detector, Classical and Quantum Gravity 29, 124007 (2012).
  * Acernese  _et al._ [2014] F. Acernese _et al._ , Advanced virgo: a second-generation interferometric gravitational wave detector, Classical and Quantum Gravity 32, 024001 (2014).
  * Abbott  _et al._ [2020] B. P. Abbott _et al._ , A guide to LIGO–virgo detector noise and extraction of transient gravitational-wave signals, Classical and Quantum Gravity 37, 055002 (2020).
  * Martynov  _et al._ [2016] D. V. Martynov, Hall, _et al._ , Sensitivity of the advanced ligo detectors at the beginning of gravitational wave astronomy, Phys. Rev. D 93, 112004 (2016).
  * Aston  _et al._ [2012] S. M. Aston _et al._ , Update on quadruple suspension design for advanced ligo, Classical and Quantum Gravity 29, 235004 (2012).
  * Liu and Thorne [2000] Y. T. Liu and K. S. Thorne, Thermoelastic noise and homogeneous thermal noise in finite sized gravitational-wave test masses, Phys. Rev. D 62, 122002 (2000).
  * Müller-Ebhardt  _et al._ [2009] H. Müller-Ebhardt, H. Rehbein, C. Li, Y. Mino, K. Somiya, R. Schnabel, K. Danzmann, and Y. Chen, Quantum-state preparation and macroscopic entanglement in gravitational-wave detectors, Physical Review A 80, 10.1103/physreva.80.043802 (2009).
  * Rollins  _et al._ [2020] J. G. Rollins, E. Hall, C. Wipf, and L. McCuller, pygwinc: Gravitational Wave Interferometer Noise Calculator, Astrophysics Source Code Library, record ascl:2007.020 (2020), ascl:2007.020 .
  * Wanner [2019] G. Wanner, Space-based gravitational wave detection and how lisa pathfinder successfully paved the way, Nature Physics 15, 200 (2019).