Observable quantum entanglement due to gravity

No experiment to date has provided evidence for quantum features of the gravitational interaction. Recently proposed tests suggest looking for the generation of quantum entanglement between massive objects as a possible route towards the observation of such features. Motivated by advances in optical cooling of mirrors, here we provide a systematic study of entanglement between two masses that are coupled gravitationally. We first consider the masses trapped at all times in harmonic potentials (optomechanics) and then the masses released from the traps. This leads to the estimate of the experimental parameters required for the observation of gravitationally induced entanglement. The optomechanical setup demands LIGO-like mirrors and squeezing or long coherence times, but the released masses can be light and accumulate detectable entanglement in a timescale shorter than their coherence times. No macroscopic quantum superposition develops during the evolution. We discuss the implications from such thought experiments regarding the nature of the gravitational coupling.


INTRODUCTION
The successful unification of electromagnetic, weak and strong interactions within the quantum framework strongly suggests that gravity should also be quantised. Up to date, however, there is no experimental evidence of quantum features of gravity. In numerous experiments, gravity is key to the interpretation of the observed data, but it is sufficient to use Newtonian theory (quantum particle moving in a background classical field) or general relativity (quantum particle moving in a fixed spacetime) to gather a meaningful understanding of such data. Milestone experiments described within Newtonian framework include gravity-induced quantum phase shift in a vertical neutron interferometer, 1 precise measurement of gravitational acceleration by dropping atoms, 2 or quantum bound states of neutrons in a confining potential created by the gravitational field and a horizontal mirror. 3 Quantum experiments that require general relativity include gravitational redshift of electromagnetic radiation 4 or time dilation of atomic clocks at different heights. 5 A number of theoretical proposals discussed scenarios capable of revealing quantumness of gravity. For example, refs. [6][7][8][9][10][11][12][13][14] proposed the observation of a probe mass interacting with the gravitational field generated by another mass. More recent proposals put gravity in a role of mediator of quantum correlations and are based on the fact that quantum entanglement between otherwise non-interacting objects can only increase via a quantum mediator. [15][16][17] Motivated by these proposals and by advances in optomechanics, 18 in particular the cooling of massive mechanical (macroscopic) oscillators close to their quantum ground state [19][20][21] and the measurement of quantum entanglement of a two-mode system, [22][23][24] we study two nearby cooled masses interacting gravitationally.
We propose two scenarios capable of increasing gravitational entanglement between the masses. In the first scenario, we consider the masses trapped at all times in 1D harmonic potentials (optomechanics). In the second one, the masses are released from the optical traps. For both settings, we derive an analytic figure of merit characterising the amount of gravitationally induced entanglement and the time it takes to observe it. The derivation includes various initial states and shows that the objects have to be cooled down very close to their ground states and that squeezing of their initial state significantly enhances the amount of generated entanglement. We then formulate a numerical approach, which accounts for all the relevant sources of noise affecting the settings that we propose, to identify a set of parameters required for the observation of such entanglement. Finally, we discuss the conclusions that can be drawn from this experiment with emphasis on the need for independent laboratory verification that the gravitational interaction between nearby objects is indeed mediated.

Proposed setup
Consider two particles, separated by a distance L, as depicted in Fig. 1. In what follows, we study the setting where the massive particles are either held or released from unidimensional harmonic traps. In the former case, one can treat the particles as identical harmonic oscillators, with the same shape, mass m, and vibrational frequency ω. The two oscillators and the gravitational interaction between them give rise to the total Hamiltonian H = H 0 + H g , where and H g describes the gravitational term. If the harmonic traps are removed, the corresponding Hamiltonian simplifies to . Before we proceed with detailed calculations, we shall discuss generic features of the gravitational term and the conditions required for the creation of entanglement.
In general, the gravitational term H g depends on the geometry of the objects. Various configurations have been analysed in the Supplementary Information accompanying this paper. The results of such analysis suggest that spherical masses give rise to the 1 highest amount of generated entanglement. The Newtonian gravitational energy of this setting is the same as if the two objects were point-like masses, that is where L is the distance between the objects at equilibrium and x A (x B ) is the displacement of mass A (B) from equilibrium. By expanding the energy in the limit x A − x B ≪ L, which is well justified for oscillators that are cooled down close to their ground state, one gets The first term is a rigid energy offset, while the second is a bi-local term and cannot thus give rise to quantum entanglement. The third term, which is proportional to ðx A À x B Þ 2 , is the first that couples the masses. When written in second quantisation, it becomes apparent that this term includes contributions responsible for the correlated creation of excitations in both oscillators. In the quantum optics language, this is commonly referred to as a "two-mode squeezing" operation, which can in principle entangle the masses provided a sufficient strength of their mutual coupling. Based on this observation, we provide an intuitive argument setting the scales of experimentally relevant parameters, which will then be proven rigorously.

Calculations of entanglement: oscillators
In order to achieve considerable entanglement, we should ensure that the coupling (third term) in Eq. (2) is comparable to the energy hω of each oscillator, that is Gm 2 ðx A À x B Þ 2 =L 3 $ _ω. As we assume that the oscillators are near their ground state, we estimate their displacements by the ground state extension, ðx A À x B Þ 2 $ 2_=mω. We thus introduce the (dimensionless) figure of merit We should have η~1 in order for the oscillators to be significantly entangled. This sets the requested values of the experimentally relevant parameters m, ω, and L.
In what follows, we will demonstrate the following results, which embody the key findings of our investigation: (i) Starting from the ground state of each oscillator and assuming (for the sake of argument) only negligible environmental noise, the maximum entanglement (as quantified by the logarithmic negativity 25,26 ) generated during the dynamics is given by E max th % η=ln2. Moreover, the time taken for entanglement to reach such maximum value is t max th ¼ π=2ð1 À ηÞω. (ii) Single-mode squeezing of the initial ground state of each oscillator substantially enhances the gravity-induced entanglement. The corresponding maximum entanglement becomes E max sq % js A þ s B j=ln2, where s j (j = A, B) is the degree of squeezing of the jth oscillator, and we assume η ≪ s A , s B . In this case, the maximum entanglement is reached in a time t max sq ¼ π=2ηω. (iii) Weaker entanglement is generated with increasing temperature of the masses or coupling to the environment.
As the third term in Eq. (2) is already very small under usual experimental conditions, we neglect all terms of order higher than the second in the displacement from equilibrium. Note that the ratio between any two consecutive terms in Eq. (2) is given by For instance, taking m = 100 μg, ω = 100 kHz, and L = 0.1 mm gives this ratio~10 −12 , and for macroscopic values m = 1 kg, ω = 0.1 Hz, and L = 1 cm the ratio is~10 −15 . We note ref. 27 for similar treatment of linearised central-potential interactions. By taking the total Hamiltonian with a suitably truncated gravitational term H g , one gets a set of Langevin equations in Heisenberg picture where we have introduced the constant frequency ν ¼ . These equations incorporate Brownian-like noise-described by the noise operators ξ j -and damping (at rate γ) affecting the dynamics of the mechanical oscillators, due to their interactions with their respective environment. We assume the (high mechanical quality) conditions Q ¼ ω=γ ) 1, as it is the case experimentally, so that the Brownian noise operators can de facto be treated as uncoloured noise and we can write À1 is the thermal phonon number with β ¼ hω=k B T and T the temperature of the environment with which the oscillators are in contact. The linearity of Eqs. (4) and the Gaussian nature of the noise make the theory of continuous variable Gaussian systems very well suited to the description of the dynamics and properties of the oscillators under scrutiny. In this respect, the key tool to use is embodied by the covariance matrix V(t) associated with the state of the system, whose elements V ij (t) = 〈u i (t)u j (t) + u j (t)u i (t)〉∕2 − 〈u i (t)〉〈u j (t)〉 encompass the variances and correlations of the elements of the quadrature vector uðtÞ ¼ ðX A ðtÞ; P A ðtÞ; X B ðtÞ; P B ðtÞÞ T . The temporal behaviour of physically relevant quantities for our system of mechanical oscillators can be drawn from V(t) by making use of the approach for the solution of the dynamics that is illustrated in Methods.
Due to weakness of the gravitational coupling, we have η ≪ 1 in practically any realistic experimental situation, and we thus assume such condition throughout. In the case of no damping (i.e., γ = 0) and assuming an initial (uncorrelated) thermal state of the oscillators, a tedious but otherwise straightforward analytical derivation shows that the entanglement between the mechanical systems, as quantified by the logarithmic negativity, oscillates in time with an amplitude of η=ln2 À log 2 ð2n þ 1Þ. At low operating temperature, a condition achieved through a combination of passive and radiation-pressure cooling, 18 n % 0 and the maximum entanglement between the oscillators is E max th % η=ln2, a value reached at a time t max th such that ωt max th ¼ π=2ð1 À ηÞ. An analytic solution is also possible for the case of mechanical systems initially prepared in squeezed thermal states, a situation that can be arranged by suitable optical driving. 30,31 Each mass is prepared in a state Sρ th S † , where ρ th is a thermal state and S ¼ expðÀi sðX 2 À P 2 Þ=2Þ is the squeezing operator with strength s. This operator corresponds to anti-squeezing (squeezing) the position quadrature for s > 0 (s < 0). By writing individualoscillator squeezing as s j and assuming s j ≫ η, the entanglement is again observed to oscillate, but with amplitude js A þ s B j=ln2 À log 2 ð2n þ 1Þ. Note that it is irrelevant whether Fig. 1 Proposed experimental setup. Two masses, placed at a distance L, are either trapped with harmonic potentials at all times or released after cooling has been achieved. The particles are assumed to be cooled down near the ground state of their trapping potentials. We study entanglement generated in both scenarios and note that it can be probed with weak light fields. Our model includes gravitational coupling (dominant), noise, damping, decoherence and Casimir forces.
T. Krisnanda et al. the quadratures of both masses are squeezed or anti-squeezed. We provide an explanation in the Supplementary Information. Therefore, only the degree of pre-available single-oscillator squeezing and the environmental temperature set a limit to the amount of entanglement that can be generated between the mechanical systems through the gravitational interaction. In the low temperature limit, where E max sq % js A þ s B j=ln2, which is in principle arbitrarily larger than the case without squeezing, a time t max sq ¼ π=ð2ηωÞ ) t max th would be required for such entanglement to accumulate. Needless to say, long accumulation times are far from the possibilities offered by state-of-the-art optomechamical experiments, which prompts an assessment that includes ab initio the effects of environmental interactions.
In the case of noisy dynamics, however, an analytical solution is no longer available and we have to resort to a numerical analysis. Let us therefore consider the figure of merit η in order to set the parameters for numerical investigation. We consider two oscillators of spherical shape with uniform density ρ and radius R, which are separated by a distance L = 2.1R. This might be a situation matching current experiments in levitated optomechanics, 32,33 which are rapidly evolving towards the possibility of trapping multiple dielectric nano-spheres in common optical traps and controlling their relative positions. However, low-frequency oscillators, which are favourable for the figure of merit and typically associated with large masses, are unsuited to such platforms and would require a different arrangement, such as LIGO-like ones. 21 In terms of the density ρ, we have η = 8πGρ / 3(2.1) 3 ω 2 , which does not depend on the dimensions of the oscillators nor their mass. As the density of materials currently available for such experiments varies within a range of only two orders of magnitude, the linear dependence on ρ sets a considerable restriction on the values that η can take. The densest naturally available material is Osmium, which has ρ = 22.59 g/cm 3 and, in order to provide an upper bound to the generated entanglement that would be attainable using other materials, we shall use this density in our numerical simulations. Accordingly, η = 1.36 × 10 −6 / ω 2 , where ω is in Hertz. Figure 2 shows exemplary entanglement dynamics for different values of the thermal phonon number n and mechanical quality factor Q. The frequency has been fixed to ω = 0.1 Hz (cf. Discussion section). As expected, higher damping (lower Q) results in the decay of entanglement, and the higher the temperature of the mirror (higher n), the higher the mechanical quality factor needed to maintain entanglement. The setup allows for high entanglement, even with low coupling strength η~10 −4 . However, this comes at the expense of the time for which the dynamics of the oscillators should be kept coherent. It is also evident that cooling down the masses close to their ground state, n % 0, is crucial for the reduction of the required coherence time. The oscillations of entanglement for unsqueezed initial state are still present in this dynamics, showing repeating pattern with a period of π / [(1 − η)ω] ≈ 31 s.
Calculations of entanglement: released masses As seen, the experimental parameters required for detectable gravitational entanglement of masses in harmonic traps are demanding. We therefore study one more feasible system, where the traps are switched off after cooling the masses. Similar to the treatment of two oscillators, one starts with the total Hamiltonian for free masses and truncated gravitational term, and obtains the following equations of motion: Note that ω here just sets the conversion between x j , p j and their dimensionless counterparts X j , P j . In what follows, we will consider starting the dynamics with thermal state for each mass. For example, the ground state is a Gaussian state with width Δxð0Þ ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi _=2mω p . This way, one can think of ω as a parameter characterising the initial spread of the wave function.
One can obtain the covariance matrix V(t) from Eqs. (5) and consequently derive the entanglement dynamics using the approach discussed in Methods section. After imposing the limits η ≪ 1 and ffiffiffi η p ωt ( 1, which apply in typical experimental situations, one obtains an analytical expression for the entanglement dynamics as follows: E th ðtÞ ¼ maxf0; E gnd ðtÞ À log 2 ð2n þ 1Þg; where E gnd ðtÞ is the entanglement with initial ground state for each mass and σ(t) = 4G 2 m 2 ω 2 t 6 / 9L 6 . Since entanglement is an increasing function of σ(t), the latter is a figure of merit for entanglement gain relevant in the case of released masses. We present exemplary entanglement dynamics in Fig. 3 for which entanglement~10 −2 is achieved within seconds. The parameters used here are m = 100 μg, ω = 100 kHz, and L = 3R. We will show later that with these values gravity is the dominant interaction and coherence times are much longer than 1 s. Note that this setup does not require any squeezing. These improvements over the scheme with trapped masses are the result of unlimited expansion of the wave functions. For example, for initial ground state, the evolution of the width of each sphere closely follows ΔxðtÞ % ffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi _=2mω p ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 1 þ ω 2 t 2 p , which is an exact solution to a free non-interacting mass. The effect of gravity is stronger attraction of parts of the spatial superposition that are closer, and hence generation of position and momentum correlations, leading to the growth of quantum entanglement; see inset in Fig. 3.
In order to understand the effect of squeezing in this setup, let us suppose, for simplicity, the squeezing strengths s A,B = s. It is as if one initially prepared each mass in a Gaussian state with a new initial spread Δx 0 ð0Þ ¼ Δxð0Þ expðsÞ. One can then calculate the entanglement dynamics using Eq. (6) with a new frequency ω 0 ¼ ω expðÀ2sÞ. This means that anti-squeezing the initial position quadrature (s > 0) would decrease entanglement gain, a situation opposite to the oscillators setup. This is because a Gaussian state with smaller Δx(0) spreads faster, such that during the majority of the evolution, the width is larger than that if one started with larger Δx(0). In principle, one obtains higher entanglement gain by squeezing the position quadrature (s < 0). However, this will result in higher final width, making it more susceptible to decoherence by environmental particles (see Discussion).
From numerical simulations, one confirms that, within t = [0, 10] s, the displacements of the two masses follow x A − x B ≪ L. Furthermore, the trajectories coincide for both quantum treatment with truncated gravitational energy and classical treatment with full H g (see Supplementary Information). This justifies the approximations used.

Discussion
We have shown that two nearby masses-both trapped and released-can become entangled via gravitational interaction. Let us now discuss the conditions required to observe this entanglement in light of recent experimental achievements.
Logarithmic negativity in the order 10 −2 has already been observed between mechanical motion and microwave cavity field. 22 Extrapolating the same entanglement resolution to the case of two massive oscillators sets the required frequency to ω1 0 −2 Hz, see Eq. (3) and its expression in terms of ω. Interestingly, kilogram-scale mirrors of similar frequency (ω~10 −1 Hz) were recently cooled down near their quantum ground state. 21 Furthermore, recent experiments on squeezed light have reported high squeezing strength 34,35 (see also a review in this context 36 ), up to 15 dB, which corresponds to s ≈ 1.73. Advances in the state transfer between light and optomechanical mirrors 18 make this high squeezing promising also for mechanical systems.
For released masses, the experimental requirements are more relaxed. Their mass can be considerably smaller while the frequency for initial trapping considerably higher, which is close to common experimental parameters used for optomechanical system. [18][19][20] Note that higher frequency (lower Δx(0)) improves entanglement gain, unlike in the oscillators case where small ω is preferable. However, one has to be cautious of decoherence mechanisms as a result of faster spreading rate of the wavefunctions. For future experiments, an improvement in the sensitivity of entanglement detection will also be beneficial.
In light of their proximity, apart from gravitational interaction, the two masses can also interact via Casimir force. It has been shown that the Casimir energy between two nearby spheres is given by a fraction of the "proximity force approximation" E ¼ Àf 0 ðπ 3 =1440Þ_cR=ðL À 2R þ x B À x A Þ 2 , with the factor 0 ≤ f 0 ≤ 1. 37,38 As typically x A − x B ≪ L − 2R, we expand such expression to find a quadratic term in x A − x B that can produce entanglement between the masses. The strength of this term, however, is much weaker than the strength of the corresponding entangling term of gravitational origin: for Osmium oscillators with mass~1 kg separated by L = 2.1R, the ratio between the Casimir and gravitational term is r cg~1 0 −12 . Similar calculations made for released masses of the same material with m = 100 μg and L = 3R, give r cg~1 0 −2 . It is thus legitimate to ignore Casimir interaction in both schemes.
Let us also discuss common decoherence mechanisms, i.e., due to interactions with thermal photons and air molecules; 39 see Supplementary Information for details. We take the average width of the wave function as an estimate for the superposition that is subjected to decoherence. All the situations we consider follow the limit Δx ≪ λ, i.e., the "size" of superposition is much smaller than the wavelength of the particles causing the decoherence. For oscillators made of Osmium, we use m~1 kg and frequency ω~0.1 Hz. Taking L = 2.1R and starting with ground state give Δx ≈ 8 × 10 −17 m. From interactions with thermal photons at environmental temperature of 4 K (liquid Helium), the coherence time for the oscillators is τ ph~5 s. The coherence time due to collisions with air molecules can be improved by evacuating the chamber with the oscillators-for about 10 12 molecules/m 3 (ultrahigh vacuum), the coherence time is τ am~5 s. One could also consider performing these experiments in space. Taking the temperature as 2.7 K (cosmic microwave background) and assuming 10 7 molecules/m 3 (space pressure~10 −15 Pa) we obtain τ ph~1 70 s and τ am~1 0 6 s.
By making similar calculations for released masses, with parameters considered in Fig. 3, one obtains τ ph~1 0 5 s and τ am 10 −4 s for the experiment on Earth with liquid Helium temperature and ultrahigh vacuum. For the space experiment, the coherence times improve to τ ph~1 0 7 s and τ am~4 1 s, respectively.
Other schemes have been proposed for gravitationally induced entanglement. 16,17 They are based on Newtonian interaction between spatially superposed microspheres with embedded spins. In those proposals, entanglement is reached faster and the small size of the experiment is the main advantage. However, in order to separate gravitational and Casimir contributions in that setup, each diamond sphere with mass m = 10 −14 kg has to be superposed across 250 μm. Decoherence due to scattering of molecules then becomes the main limiting factor. The schemes we discussed here are complementary in a sense that vibrations of each oscillator are minute (no macroscopic superposition) but larger mass, 100 μg, is required for observable entanglement.
We conclude with analysis of implications on the nature of the gravitational coupling that one can draw from such experiments and future research directions. In laboratory, we deal with two nearby masses which are experimentally shown to become entangled. These setups can be theoretically treated in different ways depending on the assumptions one makes about gravity. In a "conservative approach", the two masses are coupled via Newtonian potential. As seen from our calculations and those in refs. 16,17 this indeed leads to gravitationally induced entanglement. In this picture, gravity is a direct interaction and hence it is difficult to draw conclusions about the form of quantumness of the gravitational field. We note that even in this conservative approach such an experiment has considerable value as it would Fig. 3 Entanglement dynamics between two released Osmium spheres. Each mass has m = 100 μg and is initially prepared in a Gaussian thermal state from a harmonic trap with ω = 100 kHz. The two masses are separated by L = 3R ≈ 0.3 mm. With these parameters, gravity dominates Casimir interactions and observable entanglement is generated in seconds, which is shorter than the coherence times. In particular, entanglement in the order 0.01 is achieved within 0.8 s with initial ground state, and in 4.5 (7.5) s when starting with thermal states of n ¼ 1 (5).
show the necessity of at least the quadratic term in the expansion of the Newtonian potential for generating entanglement.
The objection to the conservative approach is instantaneity of gravitational interaction: Newtonian potential directly couples masses independently of their separation. On the other hand, it has been shown that gravitational waves travel with finite speed. 40 For nearby masses, this retardation is hardly measurable and Newtonian potential is dominant and expected to correctly describe the amount of generated entanglement. A more consistent option in our opinion, motivated by quantum formalism and comparison with other fundamental interactions, is to treat gravitational field as a separate physical object. In this picture, the masses are not directly coupled, but each of them individually interacts with the field. It has been argued within this mass-field-mass setting that entanglement gain between the masses implies non-classical features of the field. [15][16][17] This discussion shows that it would be useful to provide methods for independent verification of the presence or absence of a physical object mediating the interaction. We finish with a toy example of a condition capable of revealing that there was no mediator. To this end, we consider two scenarios: (i) evolving a bipartite system described at time t by a density matrix ρ 12 (t); (ii) two objects interacting via a mediator M, i.e., with Hamiltonian H 1M + H M2 , described by a tripartite state ρ 1M2 (t). We ask whether there exists bipartite quantum dynamics ρ 12 (t) that cannot be obtained by tracing out the mediator in scenario (ii). Indeed, if ρ 12 (t) is a pure state at all times and entanglement increases, the dynamics could not have been mediated. The purity assumption requires the mediator to be uncorrelated from ρ 12 (t), and uncorrelated mediator is not capable of entangling the principal system, composed of particles 1 and 2. 15 It would be valuable to generalise this argument to mixed states measured at finite number of time instances.
For example, in an experimental situation, the state of particles 1 and 2 might only be close to a pure state (with purity 1 − ε, where ε is a small positive parameter) and therefore they could be weakly entangled to the mediator M (with entanglement $ OðεÞ). One would naturally expect that entanglement gain between particles 1 and 2 is bounded, e.g., as a function of ε. An observation of higher entanglement gain therefore excludes the possibility of mediated dynamics.

Langevin equations and covariance matrix
Let us first consider the setup with oscillators. One can rewrite the equations in (4) as a single matrix equation _ uðtÞ ¼ KuðtÞ þ lðtÞ, with the vector uðtÞ ¼ ðX A ðtÞ; P A ðtÞ; X B ðtÞ; P B ðtÞÞ T and a drift matrix We split the last term in the matrix equation into two parts, representing the noise and constant term, respectively, i.e., l(t) = υ(t) + κ, where υðtÞ ¼ ð0; ξ A ðtÞ; 0; ξ B ðtÞÞ T and the constant κ ¼ νð0; 1; 0; À1Þ T with ν ¼ Gm 2 = ffiffiffiffiffiffiffiffiffiffiffiffiffiffi _mωL 4 p . The solution to the Langevin equations is given by where W ± ðtÞ ¼ expð±KtÞ. This allows one to calculate the expectation value of the ith quadrature 〈u i (t)〉 numerically, which is given by the ith element of where we have used the fact that the noises have zero mean, i.e., 〈υ i (t)〉 = 0 and that 〈κ〉 = tr(κρ) = κ. From Eq. (8), one can also calculate other important quantities via the covariance matrix as shown below. Covariance matrix of our system is defined as This means that κ does not contribute to Δu i (t) (and hence the covariance matrix) since 〈κ〉 = κ. We can then construct the covariance matrix at time t from Eq. (8) without considering κ as follows V ij ðtÞ ¼ hU i ðtÞu j ðtÞ þ u j ðtÞu i ðtÞi=2 À hu i ðtÞihu j ðtÞi where D ¼ Diag ½0; γð2n þ 1Þ; 0; γð2n þ 1Þ and we have assumed that the initial quadratures are not correlated with the noise quadratures such that the mean value of the cross terms are zero. A more explicit solution of the covariance matrix, after integration in Eq. (10), is given by which is linear and can be solved numerically. Consider a special case, in which the damping term γ is negligible, giving D = 0. In this case, Eq. (11) simplifies to In this regime, we have obtained analytical results. For free masses, one can follow similar treatments as above, keeping in mind γ = 0 and υ(t) = (0, 0, 0, 0) T such that the solution to quadrature dynamics and covariance matrix is given by Eqs. (9) and (12), respectively, with a new drift matrix Entanglement from covariance matrix The covariance matrix V(t) describing our two-mode system can be written in a block form where the component I A (I B ) is a 2 × 2 matrix describing local mode correlation for A (B) while L is a 2 × 2 matrix characterising the intermodal correlation. A two-mode covariance matrix has two symplectic eigenvalues fν 1 ; ν 2 g. A physical system has ν 1 ; ν 2 ! 1=2. 41 For entangled modes, the covariance matrix will not be physical after partial transposition with respect to mode B (this is equivalent to flipping the sign of the oscillator's momentum operator P B in V(t)). This unphysical VðtÞ TB is shown by the minimum symplectic eigenvalueν min <1=2. The explicit expression is given byν min ¼ ðΣ À ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi Σ 2 À 4 det V p Þ 1=2 = ffiffi ffi 2 p , where Σ = det I A + det I B − 2 det L. Entanglement between modes A and B is then quantified by logarithmic negativity as E ¼ maxf0; Àlog 2 ð2ν min Þg. 25,26 Note that the separability condition, when VðtÞ TB hasν min ! 1=2, is sufficient and necessary for two-mode systems. 42

DATA AVAILABILITY
All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Information.

CODE AVAILABILITY
Codes are available upon request from the authors.