Continuous-Variable Entanglement through Central Forces: Application to Gravity between Quantum Masses

We describe a complete method for a precise study of gravitational interaction between two nearby quantum masses. Since the displacements of these masses are much smaller than the initial separation between their centers, the displacement-to-separation ratio is a natural parameter in which the gravitational potential can be expanded. We show that entanglement in such experiments is sensitive to initial relative momentum only when the system evolves into non-Gaussian states, i.e., when the potential is expanded at least up to the cubic term. A pivotal role of force gradient as the dominant contributor to position-momentum correlations is demonstrated. We establish a closed-form expression for the entanglement gain, which shows that the contribution from the cubic term is proportional to momentum and from the quartic term is proportional to momentum squared. From a quantum information perspective, the results find applications as a momentum witness of non-Gaussian entanglement. Our methods are versatile and apply to any number of central interactions expanded to any order.


Introduction
Due to the weakness of gravitational coupling, all quantum experiments up to date in which gravity plays a role utilized the field of the Earth, see Refs. [1][2][3][4][5][6] for milestone examples. Since this field undergoes practically undetectable back-action from quantum particles, it effectively admits a classical description either in terms of a fixed background Newtonian field [3][4][5] or as a fixed background spacetime [1,2,6]. This argument strongly motivates theoretical and experimental research towards a demonstration of gravitation between two quantum masses, as this is one of the simplest scenarios where quantum properties Ankit Kumar: kumar.ankit.vyas@gmail.com of gravity could be observed. Along this line, several proposals studied the possibility of the generation of quantum entanglement between two massive objects [7][8][9][10][11][12][13][14][15][16]. Our aim here is to build upon these ideas and provide a simple precision test of gravitational coupling between two nearby quantum particles. The methods we develop are generic and apply to arbitrary central interactions, including situations where many of them are present at the same time.
The experiment we have in mind could be realized within the field of optomechanics [17], which already succeeded in cooling individual massive particles near their motional ground state [18][19][20], and in entangling cantilevers to light and themselves [21][22][23]. In such a setup the particles are separated much more than their displacements. For example, two Osmium spheres (the densest natural material) each of mass 100 µg (radius 0.1 mm) with an initial inter-surface distance of 0.1 mm move by less than a nanometer within 1 second of evolution [11], vividly illustrating the weakness of gravity. Since the situation is nonrelativistic, the relevant interaction is characterised by quantum Newtonian potential. Given that the displacements are small compared to the initial separation between the two spheres, a natural parameter in which the potential can be expanded is the displacement-to-separation ratio [11,12,[24][25][26]. We Pr ob e Pr ob e propose to identify phenomena that can only occur if the potential is expanded to a particular order, thus witnessing the relevance of at least this order in experiments. From this perspective, the gravitational entanglement proposals, in addition to providing clues about the quantum nature of gravity, also supply tests of the form of gravitational interaction. For example, entangling two initially disentangled masses requires at least the second-order term [7,8,11]. Here we show a method that witnesses the third-and fourth-order terms and has an advantage of a simple modification of the entanglement scheme with confined particles. Hence, an experiment designed to probe gravitational entanglement can also be used to witness even weaker gravitational coupling.
Our basic idea is to push the particles towards each other as it is intuitively expected that such obtained stronger gravity will lead to higher accumulated entanglement. Yet, we demonstrate that the quantum entanglement generated by gravitational potential truncated at the second order is insensitive to any relative motion of the two masses. This is shown explicitly with an exact analytical solution for the time evolution of the corresponding covariance matrix [27][28][29]. Our intuition is only recovered with the potential containing at least the third-order term, i.e., when the system evolves into non-Gaussian states. Closed-form expressions for the amount of entanglement are established for the potential expanded to any order. The introduced methods apply to any central interaction, even when many are present side by side. They also show remarkable robustness, e.g., even the impact of the fourth-order term on the non-Gaussianity quantifier and the amount of entanglement can be captured numerically despite an astonishingly weak gravitational interaction. Moreover, the derived closed forms can be extrapolated for expansions to arbitrary order. They are in remarkable quantitative agreement with numerical simulations, which show that the contribution from the cubic term is proportional to momentum. In contrast, the contribution from the quartic term is proportional to momentum squared. Accordingly, the cubic correction decreases the entanglement gain when the two particles are moving away, and the quartic correction increases the entanglement irrespective of whether they are moving towards or away from each other.

Experimental setup
Consider the setup schematically represented in Fig. 1, where we introduce our notation. The initial wave function is assumed to describe two independent masses, each in a natural Gaussian state with position spread σ: Note that without loss of generality we chose the momenta to be opposite and equal. The Hamiltonian in the non-relativistic regime is given bŷ Since this is a two-body problem, it is well-known that the center-of-mass (COM) motion separates from the relative movement. Accordingly, we introduce the usual change of variables from the LAB frame to the COM frame: R = (x A +x B )/2 and r = x B −x A , where R and r denote the respective displacements of the COM and the reduced mass from their initial average positions (see Appendix A for details). As a result, the initial wave function separates as The wave functions φ and ψ describe the motion of the COM and the reduced mass, respectively. Compared to the original wave packets, the COM wave packet admits a smaller width of σ/ √ 2, and the reduced mass wave packet has a larger width of σ √ 2. The corresponding relations are illustrated in Fig. 2. In this frame the Hamiltonian decouples aŝ whereP = −i ∂/∂R andp = −i ∂/∂r are the momentum operators for the COM and the reduced mass, respectively. A separable Hamiltonian implies that the two-body wave function retains its product form at all times, i.e., Ψ(x A , x B , t) = φ(R, t) ψ(r, t). Furthermore, the COM wave packet evolves like a free particle, i.e., its Gaussianity is preserved [30][31][32]. The Figure 2: From LAB frame to COM frame. Gaussianity of the initial state is preserved as well as the product form. The widths, however, are different in different frames as marked.
first two statistical moments characterize the quantum state fully, and for completeness, they are given in Appendix B.
The state ψ evolves in the gravitational potential, which we now expand in the powers of the displacement-to-separation ratio r/L: where N is the order of approximation, and we defined ω 2 = 4Gm/L 3 for later convenience. In Appendix B, we derive exact analytical expressions for the statistical moments of ψ by solving the related Ehrenfest equations in the case of N = 2. Together with the statistical moments for the COM, these determine the covariance matrix in an exact closed-form [see Appendix C for the methodology].
With the inclusion of higher-order terms in the potential, i.e., N > 2, the corresponding Ehrenfest's equations cannot be solved analytically due to the emergence of an infinite set of coupled differential equations involving ever-increasing statistical moments. We therefore resort to numerical methods and calculate the time evolution of ψ by implementing Cayley's form of evolution operator [33]. The numerical evaluations for ψ are combined with analytical solutions for the COM to construct the covariance matrix and the two-body wave function. In order to deal with weak gravitational interaction, we improve the accuracy of Cayley's method by utilizing the five-point stencil and discretise onto a pentadiagonal Crank-Nicolson scheme, which is further solved by implementing the LU factorization techniques. The code is publicly available at Zenodo [34], with the corresponding documentation in Ref. [35] where we demonstrate our superior accuracy as compared to the standard tridiagonal solutions. We also implemented a dynamic grid allocation, as described in Ref. [36], to avoid any reflections from numerical boundaries.

Results
The methodology just described returns an analytical form of the covariance matrix at time t for potentials truncated at N = 2 and a numerical form of the two-body wave function for all N . These are thereafter used for computing the entanglement between two masses [see Appendix C for the methodology]. In particular, we use logarithmic negativity and entropy of entanglement as entanglement quantifiers. While logarithmic negativity is known to be a faithful entanglement quantifier for Gaussian states [27][28][29], we will also discuss non-Gaussian pure states and hence the inclusion of the entropy of entanglement. We first give the results for N = 2, emphasizing the independence of relative momentum and its origin. Then we move to N = 3 and demonstrate that entanglement is linearly dependent on the initial momentum. We also analyze an indicator of non-Gaussianity (skewness) and demonstrate the precision of our methods by calculating the marginal impacts of the fourth-order term in the potential expansion. A methodology to obtain closed-form expressions for the entanglement gain through potentials expanded to arbitrary order is presented. Quantitative comparisons are made with numerical simulations for the fourth-order potential, which show that the contribution from the quartic term is proportional to relative momentum squared.

Quadratic interactions
Consider first the gravitational potential truncated at the second order. We obtained exact analytical forms for the independent elements of the covariance matrix σ. The solutions simplify if they are written in terms of already introduced ω and in terms of ω 0 = /2mσ 2 , which is the frequency of harmonic trap for which the initial state is the ground state: The logarithmic negativity for p 0 = 0, in the regime of ω ω 0 and ωt 1, was already approximated to [11] E(σ) ≈ − log 2 1 + 2Ω 6 t 6 − 2Ω 3 t 3 1 + Ω 6 t 6 , (9) where Ω 3 = ω 0 ω 2 /6 ≡ G/3σ 2 L 3 . We verified that this formula indeed matches our results and emphasize that the solutions obtained in this work are exact. Hence, they can quantify entanglement outside the demanding constraints that led to Eq. (9). An example is presented below.
The most striking feature of the covariance matrix is its insensitivity to the initial momentum p 0 . Accordingly, all quantities derived from the covariance matrix, say entanglement or squeezing [24,37], are independent of the initial momentum. In this approximation, the two initially moving masses accumulate the same amount of entanglement as when they start from rest. Furthermore, the amount of entanglement is the same irrespective of whether the masses are moving toward or away from each other. This is confirmed by the numerical simulations presented in Fig. 3. Not only there is no momentum dependence in the dynamics of logarithmic negativity and entropy of entanglement, they also perfectly overlap with analytical results showing that our methods are reliable and consistent. We emphasize that the configurations considered here are non-relativistic. Field theory calculations imply momentum-dependent relativistic corrections to the Newtonian potential [38,39], and accordingly, we verify that second-order quantum entanglement generated by relativistic particles is, in principle, momentum dependent. However, for the parameters in Fig. 3, such corrections to the Newtonian potential energy are sixteen orders of magnitude smaller, hence not discussed in this work. We also note that Eq. (9) is not applicable to the configuration considered in Fig. 3 because ω 0 ≈ 25ω.

Relevance of force gradient
We now move to explanations of the observed results. Mathematically, it is clear that a non-zero force gradient across the size of the wave packet is a necessary condition for entanglement. Without it the potential is effectively truncated at N = 1, and the total Hamiltonian is the sum of local terms. Physically, entanglement is caused by correlations in complementary variables, here between positions and momenta. Due to a force gradient, the parts of the wave packets that are closer are gravitationally attracted more than the parts which are further apart. Hence a moment later, different momentum is developed across different positions within the wave packets, leading to quantum entanglement.
Furthermore, assuming that the force gradient is the dominant contributor to entanglement gain explains the insensitivity to the initial momentum.
Since the potential is truncated at N = 2, the force gradient is constant in space. Therefore, it is irrelevant if the particle moves to a different location in the meantime; hence the initial momentum does not play a role in entanglement dynamics. Quantitatively, the force gradient is F 2 = mω 2 /2, and therefore it fully describes entanglement in Eq. (9) since now In the following section we provide further evidence for the pivotal role of force gradient in the entanglement dynamics due to higher-order interactions.

Higher-order interactions
Let us continue with the working hypothesis that the force gradient is the dominant factor in entanglement dynamics. For the cubic potential, N = 3, the gra-  Figure 4: Accumulation of entanglement with gravitational potential truncated at the cubic term (N = 3). The same physical configuration as in Fig. 3. E denotes the logarithmic negativity, and S is the entanglement entropy. Panels (a) on the left show the dependence of entanglement on the initial momentum. Analytical results are calculated from the closed-form of the covariance matrix for N = 2, and coincide with the numerical results for N = 3 and p0 = 0. Panels (b) on the right compare entanglement accumulated with non-zero momentum to entanglement gained from rest. The ratios show a linear dependence on the momentum and is very well approximated in the regime of positive p0 (masses moving towards each other) with Eqs. (11) and (12). The values of p0/mL in the legends are written in multiples of 6.18082292 × 10 −3 s −1 . dient is given by F 3 (r) = (1 − 3r/L)mω 2 /2 and importantly it admits a position dependence. Accordingly, entanglement should be sensitive to the initial momentum as the gradients differ at different locations. This is indeed observed in Fig. 4a for gravitational potential truncated at the cubic term. Furthermore, when the two masses move towards each other, p 0 > 0 and r < 0, the gradient increases, matching the growing entanglement. Conversely, when the masses move away, p 0 < 0 and r > 0, the force gradient decreases, matching the slower entanglement gain. Quantitative statements can also be achieved. Fig. 4b shows experimentally friendly plots of the ratio of entanglement accumulated within time t with non-zero initial momentum to entanglement gained from rest. The numerically calculated linear dependence (solid lines) can be explained with closed expressions (dotted lines) that we now explain. The force gradients for the quadratic and the cubic interactions are related by the following factor:F 3 (r) = (1 − 3r/L)F 2 . The average factor therefore reads where we utilized the fact that p 0 is much larger than the momenta generated by gravity and the wave packet, on average, practically follows a free particle trajectory: r ≈ r cl = −2p 0 t/m. Fig. 4a shows that for vanishing initial momentum, p 0 = 0, the entanglement obtained with cubic and quadratic potentials is practically the same. We therefore extrapolate that entanglement for non-zero initial momentum is related to entanglement from rest by a simple function of the conversion factor. The plots of Fig. 4b are fitted with Note that the factor of 1/2 next to 3 in the logarithmic negativity is causing a departure from the exact conversion factor between the force gradients. These formulae are in remarkable agreement with the computational results in the regime of positive initial momentum (masses moving towards each other, the regime of experimental interest) and also work quite well for negative initial momenta. This again affirms that the force gradient is the primary driver of gravitational entanglement. Furthermore, these closed forms can now be used in many configurations to estimate the amplification of entanglement for a non-zero initial momentum given entanglement from rest. For the potentials expanded to even higher-order terms, their contribution can also be incorporated with an appropriate conversion factor between the force gradients. Note that the entanglement entropy in Eq. (11) is amplified in the same way as the force gradient in Eq. (10). Assuming that this holds for higher-order terms, the entropy amplification factor can be written as where the corrections for gravity-like interactions (inverse-square forces) arising due to the n th term in the potential expansion is Since the gravitational force between two quantum masses is weak, for the estimation of r n we approximate the reduced mass wave packet to be a Gaussian, with the average position following classical trajectory and the width following the free evolution: where ∆r 2 0 = 2σ 2 1 + ω 2 0 t 2 . With this approximation one obtains the correction terms for n ≥ 3 as where Γ is the gamma function, and the summation is only over even m. Note that the gravitational interaction is already included in F 2 , and the present estimation is for the ratio of the force gradients of different orders only, F N /F 2 . Since n ∝ 1/L n−2 , each consecutive term is diminished by a factor of L. Hence, a cubic order correction should be sufficient for practical applications in the near future. Nevertheless, one can see that the fourth-order correction to entanglement entropy is given by Unlike the third-order term, which was sensitive to the direction of momentum, the fourth-order one depends on the momentum squared, leading to a positive correction in both the scenarios of masses moving towards and away from each other. This prediction is confirmed in Fig. 5, where we show the entanglement accumulated with the gravitational potential expanded up to the fourth order. The derived formulae exactly recover the entanglement gain in the regime of positive momentum, and they work quite well in the case of negative momentum. Note that 4 also depends on the position spread, hence it might be important even for stationary configurations where the wave packet undergoes a fast expansion. Compared to the quadratic case, the cubic term lowers entanglement between particles that move away from each other. Compared to the cubic case, the quartic term adds a positive correction irrespective of the particles moving towards or away from each other. The values of p0/mL in the legends are written in multiples of 6.18082292 × 10 −3 s −1 .

Galilean relativity and a drifting COM
We made a change of reference frames to dissect the bipartite evolution into two independent singleparticle dynamics. The first one is the free evolution of the COM, and the second one is the evolution of reduced mass in gravitational potential. Under the assumption that the two spheres are imparted with equal and opposite momentum, the COM is stationary on average. While this simplifies our theoretical framework substantially, such a configuration may be difficult to achieve in experiments. It is much easier to push one of the masses while the other is at rest. In such a case the COM moves rectilinearly with a constant velocity.
The Galilean principle of relativity demands that the laws of non-relativistic physics must be invariant in all inertial frames of reference. Consequently, the centered moments of the moving COM should evolve in the same way as for the stationary COM [40]. This is readily cross-checked as we get exactly the same correlation and variances after incorporating a nonzero momentum in the initial conditions for solving COM Ehrenfest's equations. In conclusion, a uniformly moving COM has no role in generating quantum (or, for that matter, classical) correlations. Only the relative momentum matters, and as long as it remains the same, the individual momenta can be tweaked as per convenience.

Discussion
The results presented so far could also be perceived as a simple momentum-based witness of non-Gaussianity in a quantum state. Indeed the cubic term is responsible for non-Gaussian evolution that we now quantify in more detail. Fig. 6 presents the skewnessμ 3 in the evolution of the reduced mass wave function ψ. Whileμ 3 vanishes for N = 2, as it should, it rises steeply for N = 3. The physical reason is clear from Fig. 2 describing the change of variables between LAB and COM frames. The left end of wave function ψ is attracted towards the COM much more than the right end. Over time this makes the probability density function negatively skewed, which is indicated byμ 3 < 0. Just like Fig. 5, Fig. 6 also demonstrates the precision of our numerical methods, which even capture marginal contributions of the fourth-order term.
Let us summarise the interplay between non-Gaussianity, initial momentum dependence of entanglement, and the amount of entanglement. For the quadratic Hamiltonian the skewness of course vanishes, and we derived an exact analytical solution for the covariance matrix that shows no dependence on the initial momentum. For the cubic Hamiltonian the skewness rises but its value is small, and within the first few seconds of the evolution the wave function is very close to a Gaussian. Nevertheless, this is already sufficient for non-zero force gradient, and we obtained closed-form equations for the amount of entanglement that show a linear dependence on the initial momentum and do not feature skewness, see Eqs. (10) and (12). Therefore, while non-zero skewness enables entanglement dependence on the initial momentum, it plays a marginal role in the amount of entanglement. Quantitative estimation of this small contribution of skewness to the amount of entanglement is left as an open problem.
All these mutual dependencies can be seen in our data. Fig. 6 shows that skewness is practically the same for the two considered relative momenta. For the same momenta, Fig. 4a demonstrates that entanglement entropy accumulated after 5 seconds is different by 30% and linear in initial momentum. The amount of entanglement is therefore determined by the initial momentum only. Similarly, skewness is non-zero for initially stationary particles but entanglement dynamics with and without non-Gaussianity look practically the same. To give quantitative values, consider the stationary configuration of two Osmium spheres with m = 1 pg separated by a distance of L = 2.1 times their radius. After an evolution for 5 seconds, the entanglement gain with cubic potential is larger than the entanglement accumulated with quadratic potential by only ≈ 0.001, 0.002, and 0.003%, for an initial spread of σ = 5.00, 0.50, and 0.05 nm, respectively. We emphasize that the force gradient plays a pivotal role in entanglement dynamics.
We would also like to address the question of whether a simpler method for detecting the thirdorder coupling exists than based on measurements of entanglement. Indeed, note that solely the mean relative momentum signal could be used for such purposes. One verifies that by truncating the potential in Eq. (7) at N = 2, the relative momentum satisfies the condition¨ p / p = ω 2 , i.e., it is time-independent. Any time dependence of this ratio reveals third-order coupling. In cases where the center of mass is stationary, instead of the relative momentum, the local momentum of any particle could be used.
Finally, a word on decoherence effects is in place. The common decoherence mechanisms, due to thermal photons and air molecules, have already been studied in the considered setup [7,8,11,13,24]. The experiment was found to be challenging, but the required coherence times are in principle realisable, e.g., for freely-falling particles in a high vacuum. The calculations presented here only relax these requirements as the entanglement is improved when the two masses are pushed toward each other. For example, in the configuration considered in this work, the entanglement gain of E ≈ 1.75 × 10 −4 is relaxed from 5 seconds to 4 seconds with an initial momentum of p 0 /mL ≈ +0.022 s −1 . Note that an entanglement detection scheme achieving a precision of E ∼ 10 −4 has recently been put forward in Ref. [41].

Versatility
While our discussions were mainly focused on gravityinduced entanglement, the methods we have presented are applicable more generally. First of all, they hold for arbitrary central interactions. We need to expand the potential in a binomial series similar to  Eq. (7) and the entanglement is characterised by the parameters ω and 3 . As we derived, for identical masses coupled via Newtonian gravity, The general rule is quite simple. Once the potential is expanded in a binomial series of the relative displacement, the coefficient of r 2 is to be compared with −mω 2 /4, and 3 is to be calculated by comparing the force gradients as F 3 /F 2 = 1 + 3 (t). One can then verify that for the Coulomb potential between charges q 1 and q 2 embedded into the masses, we obtain where α is the fine structure constant and e is the electronic charge. For the Casimir interaction between their surfaces (under proximity force approximation [42]) we arrive at where R 0 is the radius of each sphere.
For an arbitrary central interaction with a potential of the form we obtain In some situations the force is known, but solving for the potential is difficult or uncertain due to nonunique boundary conditions. In those cases the general rule would be to expand the force in a binomial series and compare the coefficient of r with mω 2 /2. For example, if the force is of the form one arrives at Note that the functional forms of 3 in this section are valid only for weak interactions. For stronger potentials one has to take a step back and use 3 = −3 r /L, where r has to be approximated either analytically or numerically. Furthermore, the methods established also work for multiple central forces acting simultaneously. If we write the interaction as a sum (25) the total ω characterising the Gaussian covariance matrix is given by a Pythagoras-like theorem, and the equivalent 3 governing the entanglement amplification due to the cubic-order term is calculated to be a weighted sum: where ω k and k3 (t) characterise the individual interactions. This is particularly useful from an experimental point of view as, in practice, it might be difficult to screen all interactions except gravity. For example, the gravitational and the Casimir interaction will likely act side by side.

Conclusions
We have shown that experiments aimed at observing gravitational entanglement can also be used as precision tests of gravitational coupling. In particular, entanglement dependence on the relative momentum of interacting particles indicates a coupling which is higher than quadratic. Furthermore, for the potential expanded to the cubic term, the amount of entanglement accumulated in a fixed time interval grows linearly with the relative momentum when the particles are pushed toward each other. We presented a closed expression for the amount of entanglement as a function of relative momentum based on the derived exact covariance matrix for Gaussian dynamics extended to higher-order couplings. The methods introduced apply to arbitrary central interactions, even when many are present side by side.

A Coordinate transformations
The time-dependent Schrödinger equation (TDSE) for two identical particles A and B is given by (27) wherep where R and r are the displacements of the COM (mass 2m) and the reduced mass (mass m/2), respectively. One can take time derivatives to arrive at their respective momenta as Accordingly, the inverse transformations are The displacements x A and x B are functions of R and r, and the rules of differentiation imply which means that the kinetic energy is equivalent to , and hence the TDSE is transformed to

t). (33)
In this work the initial wave function is Ψ(x A , x B , t = 0) = φ(R, t = 0) ψ(r, t = 0), and the separation of variables in Eq. (33) ensures that the product form is maintained at all times. Accordingly, the problem decouples into two independent single-particle TDSEs: whereP = −i ∂/∂R andp = −i ∂/∂r can now be identified as the momentum operators for the COM and the reduced mass, respectively. The COM is a free particle (but not a plane wave), and the reduced mass undergoes an evolution under the influence of the interaction. The bipartite wave function is simply the product of the two wave functions

B The Ehrenfest dynamics
Ehrenfest's theorem relates the time derivative of the expectation value of an operatorÂ to the expectation of its commutator with the HamiltonianĤ: In this work we focus on the noiseless dynamics, i.e., when the system is isolated from the environment. Accordingly, the COM and the reduced mass evolve exclusively into pure states, with their covariance matrices satisfying Det(σ R ) = Det(σ r ) = ( /2) 2 [43].

B.1 Evolution of the COM
The COM evolution corresponds to the free expansion of a Gaussian wave packet: The initial is characterised by R = 0, P = 0, R ,P = 0, R 2 = σ 2 /2, and P 2 = 2 /2σ 2 [see Fig. 2 or Eq. (4)]. The solution to the corresponding Ehrenfest's differential equations for the first two statistical moments, imply Alternatively, one arrives at the same results by utilizing the functional form of the wave function [40]:

C Quantification of entanglement
We have employed the formalism based on the covariance matrix to quantify entanglement gain via logarithmic negativity and additionally used the density matrix to compute the entropy of entanglement.

C.1 Covariance matrix
The covariance matrix formalism is based on the first two statistical moments of a quantum state. Given a bipartite system AB withû = (x A ,p A ,x B ,p B ) T , the covariance matrix is defined as [27][28][29]: In the block form we can write where α(β) contains the local mode correlation for A(B), and γ describes the intermodal correlation. In our setting the local modes are identical, i.e., α = β, and a coordinate change to the COM frame implies σ 01 (σ 03 ) = 1 2 Cov(R, P ) +(−) 1 2 Cov(r, p).

C.3 Entropy of entanglement
For a pure bipartite system described by a density matrix ρ AB , the entanglement entropy is defined as the von Neumann entropy for any one of the subsystems, e.g., S(ρ A ) = − Tr [ρ A log 2 (ρ A )], where ρ A = Tr B (ρ AB ) is the reduced density matrix for subsystem A. In order to calculate S(ρ A ) we start with the two-body wave function of Eq. (36). The COM wave function φ(R, t) is derived analytically in Eq. (40), and we calculate ψ(r, t) numerically by implementing the improved Cayley's propagator [34]. Once this is available at a given time t, we perform a singular value decomposition [45,46]: where χ (A) j and χ (B) j are orthonormal states in subsystems A and B, respectively, and λ j are the Schmidt coefficients. A numerical implementation utilizes the algorithms in Google TensorNetwork [47,48] with the help of open source libraries hosted at GitHub [49]. Note that the total number of Schmidt coefficients is not fixed. At any given time, the number is dynamically increased until the norm bipartite wave function is recovered up to seven decimal places, i.e., until 1 − |Ψ(x A , x B , t)| 2 = 1 − j λ j (t) 10 −7 . With this decomposition, the entanglement entropy reduces to S(ρ A ) = − j λ j log 2 (λ j ). (51) In the case of Gaussian evolution, the entanglement entropy is calculable using the covariance matrix [43]: where f (x) = x+ 1 2 log 2 x+

D Numerical details
Numerical calculations are performed in natural units of c = 1, hence the conversion constant c = 197.3269804 keV pm. The density of Osmium is 22.5872 g/cm 3 . An error analysis implies that, in the numerical time evolution of the reduced mass wave function, a grid size of 0.2 pm with a time step of 10 µs is required to maintain accuracy in the extreme cases of the largest momentum considered in this work. Accordingly, we set a grid size of 0.1 pm and a time step of 5 µs throughout this work. Note that the first term in the gravitational potential of Eq. (7) is just a constant energy offset, which only contributes to an irrelevant global phase in the quantum dynamics. Moreover, this term is the most prominent of all magnitude-wise. We have therefore ignored it in numerical simulations to maximally utilize the computer precision for the (relevant) higher-order terms.