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

Probing quantum entanglement with macroscopic objects allows us 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 and discuss the parameter regimes where entanglement exists, even in the presence of quantum and classical noises.


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 the 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 mechanical decoherence imposed by coupling to a possibly highly populated thermal bath [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, Refs.[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 * sdirekci@caltech.edumany proposals in the literature to further study macroscopic quantum phenomena in optomechanical systems [17][18][19][20][21] and entanglement between coupled oscillators in the 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 [4,6,[23][24][25], whereas to the best of our knowledge an experimental demonstration has 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 ′ ≤ 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 ≫ ℏω m , T is the temperature of the bath and ω 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-noise processes, such as structural damping [27][28][29][30], are relevant.We extend the methods developed in [6], incorporating a cavity, and, more importantly, non-Markovian noise processes.The technique consists in computing the minimal symplectic 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.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.We work in this limit to study the general behavior of heavy suspended oscillators affected by environmental decoherence, however the algorithm does not require this assumption.Subsequently, we stop working in the free-mass limit and focus our attention on the Advanced Laser Interferometer Gravitational-Wave Observatory (aLIGO), [31] using 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 [32].
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 [32], we see that the level of classical noise is not yet below the SQL [33].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 ω 0 close to one of the resonant frequencies of the cavity, ω 0 + ∆ [34].The quantity ∆ 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 [35] where B and B † are the annihilation and creation operators of the mechanical mode (center of mass motion of the mirror), x is the position of the center of mass of the mirror, ω m is the mechanical resonance frequency, Â and Â † are the annihilation and creation operators of the cavity mode, ĉ(ω 0 + Ω) and ĉ † (ω 0 + Ω) are the annihilation and creation operators of the external vacuum light field at frequency ω 0 +Ω, G is the linear optomechanical coupling constant, and γ 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 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, 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 [34].From the scaling-law, the coupling constant G is related to the parameters of the interferometer by, 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 † ) 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 û and v to label the field right before entering and right after exiting the cavity, respectively.
We include two classes of "classical" noises [39] in our system: a force noise nF and a sensing noise nX , originally arising from a quantum treatment of the interaction of the system with its environment.We label them as "classical" noises because we assume that the details -and possible quantum limit -of a microscopic model of these noises are irrelevant (typically because they arise from thermal baths in the high temperature limit) unlike the noise arising from vacuum fluctuations, û1 (Ω) and û2 (Ω).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 γ m .γ m and n F are associated with the heat bath(s) the mass is coupled to, with the value of γ m and the spectrum of n F related by the fluctuationdissipation theorem [40].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 We refer to Eqs. (11) as the Heisenberg equations for the rest of the article.It is straightforward to solve them to obtain ( x, p, Â1,2 , v1,2 ) in terms of the input fields, (û 1,2 , nX , nF ), referred as the input-output relations of the system.More specifically, quantum fluctuations in the in-going quadratures û1,2 (Ω) drive the system's quantum noise [41].From Eq. (11e) and (11f), reading the out-going field quadratures are subject to noises in û1 and û2 , giving rise to the shot noise (SN) for that readout strategy [41].On the other hand, from Eq. (11a), we see that û1 drives Â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 V + 1 2 K is positive semidefinite, where V is the covariance matrix with V i j = ⟨{ Xi − ⟨ Xi ⟩, X j − ⟨ X j ⟩}⟩/2 and K i j = [ Xi , X j ] is the commutator matrix of the quadratures in the system.This relation can be stated as Here for an N-partite system containing N harmonic oscillators, the matrices V and K are 2N-dimensional.
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 [42][43][44].
To use the PPT criterion in this context, one obtains the partial-transposed covariance matrix V 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) [45].The PPT criterion for separability is expressed as The amount of entanglement is quantified by the logarithmic negativity, E N [46].For a Gaussian state of N modes, it is defined as where ν j , j = 1, . . .N are the symplectic eigenvalues of the partially transposed covariance matrix, V pt , which are given by the absolute values of the eigenvalues of K −1 V pt .For 1 vs. N − 1 mode partitions, only one of the symplectic eigenvalues of V pt can have a magnitude smaller than 1 [47], therefore there can be at most one negative eigenvalue of V pt + 1 2 K.We label the corresponding symplectic eigenvalue as ν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 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 V for both types of partitions, as well as the discretization of V, can be found in Appendix D.

IV. ENTANGLEMENT IN THE PRESENCE OF NON-MARKOVIAN NOISES
Due to the numerical nature of the algorithm, we can tackle any noise spectral density associated with û1 (Ω), û2 (Ω), nF (Ω) and nX (Ω) 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.Examples of such systems are aLIGO, KAGRA [48], and VIRGO [49], but also smaller devices such as those in [50,51].For such systems, the mechanical resonance frequency, ω m , is much smaller than the other frequencies of the system, which is referred to as the free-mass limit.In this setting, ω m essentially does not affect the dynamics.Furthermore, in this limit where Ω ≫ ω m , γ m , the tradeoff between shot noise and QRPN gives rise to the SQL [52], given by In the context of suspended oscillators, nF (Ω) is the force that gives rise to the suspension thermal noise, whereas nX (Ω) 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 nF (Ω) and nX (Ω) decreases as 1/Ω above internal resonances, which is referred to as structural damping [27,28,53].Evidently, structural damping gives rise to non-Markovian noises, and the position-referred noise spectral densities of nF (Ω) and nX (Ω) are given, in the free-mass limit, by where Ω F and Ω 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 nF (Ω), the position-referred spectrum is related to the noise spectral density, labeled as , whereas for nX (Ω), the position-referred spectrum S X (Ω) is also the noise spectral density S n X (Ω) [54].The incoming field quadratures û j have uncorrelated white spectra given by Eq. (D8), since we assume the incoming field to be at vacuum state.
In the limit of large cavity bandwidth, γ ≫ Ω the cavity can be eliminated adiabatically.Then, the equations of motion in (11) are modified as where α = Ω q √ M/ℏ and Ω q = 2G ℏ/Mγ is the characteristic interaction frequency.In the context of structural damping, there is no velocity damping, Instead, the damping arises from a complex spring constant associated with the mechanical oscillator.Accordingly, Eq. (17b) is modified as where ϕ(Ω) is referred to as the loss angle.When structural damping is present, ϕ(Ω) is constant for a large band of frequencies and goes to zero as Ω → 0, however the dependence of ϕ(Ω) to Ω depends on the properties of the material [27].Numerically, we choose to model this with FIG. 2. Logarithmic negativity between the mechanical oscillator and the outgoing light field in the free mass limit as a function of Ω X /Ω F for various Ω q /Ω 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 Ω X /Ω F ≲ 1 with Markovian and for Ω X /Ω F ≲ 2.6 with non-Markoivan noise sources, for ω m /(2π) = 1 Hz, γ m /(2π) = 0.01 Hz, ϕ = 0.05, and a cutoff frequency of Ω c /(2π) = 0.05 Hz.
ϕ(Ω) = ϕ • Ω/(Ω + Ω c ), so that ϕ(Ω) ≈ ϕ for Ω ≫ Ω c and ϕ(0) = 0 for some cutoff frequency Ω c .Then, Ω c determines the noise power of nF and nX at 0 Hz.We fix ω m /(2π) = 1 Hz, Ω F /(2π) = 100 Hz, Ω c /(2π) = 0.05 Hz, ϕ = 0.05, and vary Ω X and Ω q .Note that the aim of these choices is to ensure that Ω F ≫ ω m , i.e. that the freemass limit is justified.This choice also implies the large bath temperature limit where k B T ≫ ℏω m , T being the temperature of the bath.Furthermore, Ω q ≫ ω m ensures that the measurement of the system performed by light is faster than the dynamics of the system.Lastly, we perform the simulations by sampling the covariance matrix for τ = 0.1 s, from t = −0.1 s to t=0 s, where t indicates the time (for the sampling of the covariance matrix, see Appendix D).Hence, we cannot resolve the very low frequency regime of Ω c .Physically, this corresponds to the finite detection frequency resolution that renders the behavior of ϕ at the cutoff inaccessible.In this sense, one might expect Ω c to be irrelevant.However, we see that it does affect the final characterisation of entanglement, due to its contribution to the total variance of the quadratures of the mechanical oscillator.
After specifying the low-frequency behavior of ϕ with Ω c , the relevant parameters on which the presence of entanglement depends are Ω q , Ω F , and Ω X .Then, working in the free-mass limit enables us to examine the general behavior of suspended oscillators with low mechanical resonance frequencies, classified by their coating materials and suspension systems (i.e. the low-frequency behavior of ϕ).Therefore, we look for entanglement between the oscillator and the outgoing light field by varying the ratio Ω X /Ω F for various Ω q in Fig. 2, and we plot the results with plain lines.We find that for Ω c /(2π) = 0.05 Hz, entanglement does not exist for Ω X /Ω F ≲ 2.6 for any value of Ω q .For Ω X /Ω F ≳ 2.6, the system is entangled for any (finite) Ω q , and the entanglement increases monotonously with increasing Ω X /Ω F .This im-FIG.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.
plies the existence of "universal" entanglement, meaning that whether the system is entangled or not is independent of Ω 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 Ω q is increased.We also note that the threshold for Ω X /Ω F above which the system is entangled depends on the low-frequency cutoff Ω c chosen in our model for the spectra of nF and nX .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 The results can again be found in Fig. 2, plotted with dashed lines, where the system is not entangled for any Ω q when Ω X /Ω F ≲ 1, whereas for Ω X /Ω F ≳ 1, entanglement exists for all Ω q and the amount of entanglement increases with increasing Ω q .Since we see this behavior for both Markovian and non-Markovian noise, we prove 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 Ω q in a manuscript that is currently in preparation.

V. NOISE MODEL OF aLIGO
The primary noise sources in aLIGO, other than the quantum noise, are the following [31]: 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 the absence of glitches (i.e.transient noise artifacts) [55].
Seismic noise occurs because of the ground motion at the interferometer sites.This motion is ∼ 10 −9 m/ √ Hz at 10 Hz [56].To provide isolation from this motion, the mirrors are suspended from quadruple pendulums [57].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 fibers used in the final suspension stage [31], whereas the coating Brownian noise (which is classified as a sensing noise) occurs due to the mechanical dissipation in the coatings [58].Other types of sensing noise comprise many noise sources that are dominant at high frequencies, such as thermal fluctuations of the mirror's shape, optical losses, or photodetection inefficiency [59].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 nF (Ω), and the coating thermal noise with nX (Ω).We use the aLIGO noise budgets as given by the Python Gravitational Wave Interferometer Noise Calculator library (pygwinc) [60] and we model them with rational functions of Ω 2 .The noise spectra are modelled by where S LIGO F (Ω) is the spectrum of the force noise, and S LIGO X (Ω) is the spectrum of the sensing noise.We model S LIGO F (Ω) to decay as Ω −14 instead of Ω −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 τ F , τ X 1 , τ X 2 and cutoff frequencies ω F , ω X .The values of these parameters can be found in Appendix B, Table I and α 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.
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 [34], however, we choose to work numerically and find the parameters of the reduced cavity G, γ, M, γ m , ω 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. 7.The fitted parameters can be found in Appendix B, Table II.

VI. ENTANGLEMENT IN aLIGO
In aLIGO, the spectrum is dominated by force noise for low frequencies and quantum noise for high frequencies, therefore FIG. 4. Force (top row) and sensing noise (bottom row) spectra parametrized by α F 1 , α F 2 , α X 1 , and α X 2 .The effect of α F 1 is to rise and lower the nominal noise strength below the cutoff frequency and α F 2 shifts the cutoff frequency.Similarly, α X 1 shifts the cutoff frequency where the sensing noise starts increasing as Ω 2 , whereas α X 2 shifts the nominal noise level.
we expect the sensing noise not to 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 α F 1 and α F 2 to modify the force noise spectrum, and set α X 1 = α 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 the presence of high classical noise levels.

A. Effect of Force Noise
First, we investigate the effect of the force noise spectrum on entanglement and we calculate the logarithmic negativity E N as a function of α F 1 and α F 2 , for both of the partitions described in Sec.III.For all pairs α F 1 and α F 2 here, we find larger logarithmic negativity values when we do not trace over FIG. 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 α 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 α F 1 increases (proportional to the DC noise power) and α 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 α F 1 = α F 2 = 1, hence this plot is for a 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 give a definite answer about optomechanical entanglement in aLIGO with our model yet.
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) [61].For aLIGO, in the absence of seismic noise, suspension thermal noise dominates the spectrum for low frequencies, which is modelled as similar to how the total force noise was modelled in Sec.V.The parameters τ S T and ω S T can be found in Table I.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.

B. 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 FIG. 6.The effect of resonant modes on the logarithmic negativity, E N , for both partitions.We see that the negativity reduces with increasing α F 1 or with the introduction of resonant modes.Note that in the seismic and suspension thermal noise arise from the displacement noises of the rigid-body resonant modes of the 4-stage suspension system [57].These modes can be modeled with a sum of Lorentzians multiplying the force noise spectrum defined in Eq. ( 19).Then, the new formula defining the spectrum of the force noise is given as with Ω > 0 , where the sum is over the resonant modes.The parameters Ω v i , Γ v i , and A v i are the mode frequencies, full widths at half maximum (FWHM), and the amplitudes of the Lorentzians, respectively, and they are listed in Appendix B, Table III for the modes with the biggest relative amplitudes.We investigate the effect of the resonant modes on the entanglement.In Fig. 6, with low force noise guaranteeing wellbehaved numerics and setting α F 2 = 1, we plot E 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 α 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 α F 1 ≈ 10 −15 and α F 1 ≈ 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 α F 1 → 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. DISCUSSION AND 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 the 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.However, non-Markovian noise can also be found in optomechanical systems with higher frequencies [29,30,62,63], and our framework can be applied to them in the same way.Besides improving the quality and accuracy of predictions, modeling non-Markovian dynamics is rich in interesting physics and possibly useful phenomena: for example, in Ref. [63], they find that squashing/squeezing the mechanical state is less demanding in the presence of structural damping compared to Markovian viscous damping.
We tested for bipartite entanglement by looking at the separability between the mechanical oscillator and (i) the outgoing light field, and (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 freemass 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.Beyond the free-mass limit approximation, 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 we 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.Our model for the total classical and quantum noise in aLIGO is displayed in Fig. 7.Note that the coating Brownian noise, which is the dominating classical noise source above ∼ 10 Hz, is modeled in a piecewise manner: a white noise in the 10 − 10 5 Hz band and a noise source increasing as Ω 2 for frequencies larger than 10 5 Hz.Even though the sensing noise in the 10−10 5 Hz band decreases as Ω 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.

Appendix D: Structure of the Covariance Matrix
The mirror and the cavity, at t = 0, constitute two modes, whereas there are an infinite number of modes in the outgoing light field given by the quadratures v1 (t), v2 (t), t ∈ (−∞, 0).In continuum coordinates, we can first write down the commutator matrix, with Note that K v is a 2 × 2 block matrix, but each block is infinite dimensional, with columns and rows indexed by t and t ′ , respectively.The indices t and t ′ each run through all negative real numbers, (−∞, 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 ′ are continuous and run through negative real numbers.We can then write Note that for K v we have a two-dimensional row index (l, t), as well as a two-dimensional column index (m, t ′ ), to label quadrature and arrival time.
The covariance matrix 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 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 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 where the (one-sided) cross spectrum is defined via and they can be obtained from solutions to the Heisenberg Equations, as well as the spectra and the prescriptions we use for the spectra of S n X (Ω) and S n F (Ω).The uncorrelated white spectra between û1 and û2 in Eq. (D8) result from the ingoing light field being in its vacuum state.We shall discuss the magnitude and frequency dependence of S n X (Ω) and S n F (Ω) in depth in the next appendix.
For elements that involve v, we shall still use lt for column indices and mt ′ for row indices.We then have and Note that V Qv J,mt ′ = V vQ mt ′ ,J .In numerical computations, we will have to convert the continuum of t, t ′ ≤ 0 into a finite grid.This means we will sample a finite duration T with a step size of ∆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 we use Greek indices, for example, α = 0, 1, 2, . . ., T/∆t ≡ N − 1, to replace t.We shall write Note that a Kronecker delta now replaces the Dirac delta.For the covariance matrix, we replace which are 4 × 2N and 2N × 4 dimensional matrices (with V vQ J,mt ′ ), and which is a 2N × 2N-dimensional matrix.For the discrete sampling times we have defined where the additional 1 2 ∆t provides a faster convergence in numerics.The entire covariance matrix is then (2N+4)×(2N+4) -dimensional.
Our particular convention of inserting ∆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 where we have used upper indices for vector components, and we have grouped α j and β j into γ J .The variance of X, which is formally written as X † VX, will then be As we convert the integrals in Eq. (D16) into summations, will become Σ, while dt will become ∆t.We shall take Together with Eq. (D12)-(D13), the fully discretized version of Eq. (D16) will then be 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.(D11)-(D13), when contracted with vectors in this convention, lead to the appropriate integrals.Note that if a δ(t α − t ′ α ) shows up in Eq. (D13), we will take ∆tδ(t α − t ′ α ) → δ αα ′ , as in Eq. (D11).
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 ≤ 0 and (ii) mass and the joint system of the cavity mode as well as light that had emerged during t ≤ 0. In case (i), we simply throw away elements involving A in both K and V, while in case (ii) we consider the full matrices.In both cases, V pt is obtained by adding a minus sign to the column involving B2 and the row involving B2 -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 Ω in our force noise model Eq.(19).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.
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 ν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 nonphysical.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 ν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 γ = 424 Hz.Therefore, we choose dt < 1/(2 • 424) ≈ 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 V + 1 2 K before and after the partial transpose operation, denoted as λ B and λ N , respectively.Then, we decide that entanglement exists if λ B > 0 and λ N < 0, or if λ B < 0 and |λ N | ≫ |λ B |. Furthermore, we decide that the system is separable if λ B > 0 and λ N > 0.
Two case studies about this strategy can be found in Fig. 8 where we fix T = 0.1 s, change dt, and examine how λ N and λ B change by plotting λ N and λ B .We decide on entanglement if λ B ≥ 0 and λ N < 0, or λ B < 0, λ N < 0, and |λ N | ≥ 100|λ B |. Our criteria for convergence is a relative change smaller than 5% for both λ N and λ B as we change dt.In Fig. 8a, we work with a low level of classical force noise and set α F 1 = 10 −15 , α F 2 = 15 in Eq. (19).We see that λ N and λ B converge with λ N changing by 0.053%, 0.068%, and λ 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 λ B ≥ 0 and λ N < 0. Furthermore, λ B becomes positive and converges after dt ∼ 1 ms, or a sampling frequency of 1000 Hz; consistent with the discussion above relating physicality to Nyquist rate of ∼ 850 Hz.We also see that λ N con-FIG.9.The real and imaginary parts of the eigenvector for the negative eigenvalue of V pt + 1 2 K for α F 1 = 10 −15 , α 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.
verges for similar values of dt from Fig. 8a.
In Fig. 9, we plot the light-field section of the eigenvector associated with the (converged) minimal eigenvalue of V pt + 1 2 K, for the partition where we do not trace over the cavity and with α F 1 = 10 −15 , α F 2 = 15.It corresponds to a temporal mode of the free electromagnetic field outside the cavity.It is that particular mode associated with 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 v1 (t) and v2 (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 dicted 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 [8,16] and proposed for a demonstration in the stationary regime [23].We fit functions in the form of e −γ * t sin (ω * t + θ) to each curve, which results in ω * /(2π) ≈ 40 Hz, and γ * /(2π) ≈ 25 Hz.In the frequency domain, exponentially decaying harmonic oscillations are Lorentzians, centered at ±ω * and with a bandwidth (FWHM) 2γ * .In aLIGO's noise budget (Fig. 3), these Lorentzians are on the low frequency side of the low-noise band and their halfwidth at halfmaximum to the left crosses the quantum noise, where it is not yet dominated by suspension thermal and seismic noisesalthough we saw in Sec.VI A 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. 8b, we set α F 1 = 10 −8 , α F 2 = 10, and α 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 λ N and λ B converge with λ N changing by 1.3%, 1.2%, and λ B changing by 1.2%, 1.1% before and after tracing over the FIG.10.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.cavity, respectively, for dt = 0.1 ms.Since λ B and λ 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, λ B and λ N are negative, and λ B ∼ λ 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 either of the partitions.These case studies show that we can use λ B as an indicator of the "non-physicality" of the covariance matrix V introduced by finite-sampling and high levels of classical noise in the sys-tem, and that the negativity of λ N is not enough to decide on entanglement when we consider the numerics.
For our model of aLIGO's non-Markovian noises, Eqs. ( 19), we study the numerical stability of broad noise regimes, parametrized by the pair α F j , j = 1, 2. We set α 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. 10, 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 α F j = α 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. 10, 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. 11.Note that if we set α F 1 or α F 2 to be unity, the corresponding pairs of α F 1 , α F 2 situated on this boundary would be α F 1 = 1 and α F 2 ≥ 200, or for α F 2 = 1 and α F 1 ≤ 10 −10 .The corresponding noise curves are also plotted in Fig. 11.

FIG. 1 .
FIG. 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.

TABLE I
, sampled logarithmically in frequency.Table I contains the parameters for the classical force and sensing noise in aLIGO defined in Eqs.(19), Table II contains the parameters of aLIGO introduced in the optomechanical Hamiltonian of Sec.II, and lastly, Table III contains the parameters of the resonant modes present in aLIGO's noise curves, modeled with Eq. (21). interest