[1]Department of Physics, Indian Institute of Technology Roorkee, Roorkee
247667, India

[2]School of Physical and Mathematical Sciences, Nanyang Technological
University, Singapore 637371, Singapore

[3]Centre for Quantum Technologies, National University of Singapore,
Singapore 117543, Singapore

[4]Centre for Photonics and Quantum Communication Technology, Indian Institute
of Technology Roorkee, Roorkee 247667, India

[5]Institute of Theoretical Physics and Astrophysics, University of Gdańsk,
80-308 Gdańsk, Poland

[6]School of Mathematics and Physics, Xiamen University Malaysia, 43900
Sepang, Malaysia

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

Ankit Kumar<EMAIL_ADDRESS>Tanjung Krisnanda  P. Arumugam  Tomasz
Paterek

###### Abstract

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.

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

Figure 1: Setup under consideration. Two spheres of mass $m$ are released from
the ground state of identical harmonic traps with an equal and opposite
momentum $p_{0}$ along the line joining their centers. The centers are
initially separated by a distance $L$, and displacements from them are denoted
by $x_{A}$ and $x_{B}$. After time $t$ entanglement is estimated with the help
of the probing lasers.

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 $\mu$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 non-relativistic, 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 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.

##  2 Experimental setup

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.

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 $\sigma$:
$\Psi(x_{A},x_{B},t=0)=\psi_{A}(x_{A})\ \psi_{B}(x_{B})$, where

$\displaystyle\psi_{A}(x_{A})=\quantity(\frac{1}{2\pi\sigma^{2}})^{1/4}\exp(-\frac{x_{A}^{2}}{4\sigma^{2}}+i\frac{p_{0}}{\hbar}x_{A}),$
(1)
$\displaystyle\psi_{B}(x_{B})=\quantity(\frac{1}{2\pi\sigma^{2}})^{1/4}\exp(-\frac{x_{B}^{2}}{4\sigma^{2}}-i\frac{p_{0}}{\hbar}x_{B}).$
(2)

Note that without loss of generality we chose the momenta to be opposite and
equal. The Hamiltonian in the non-relativistic regime is given by

$\hat{H}=\frac{\hat{p}_{A}^{2}}{2m}+\frac{\hat{p}_{B}^{2}}{2m}-\frac{Gm^{2}}{L+(\hat{x}_{B}-\hat{x}_{A})}.$
(3)

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 $\Psi(x_{A},x_{B},t=0)=\phi(R,t=0)\ \psi(r,t=0)$, where

$\displaystyle\phi(R,t=0)=\quantity(\frac{1}{\pi\sigma^{2}})^{1/4}\exp(-\frac{R^{2}}{2\sigma^{2}}),$
(4)
$\displaystyle\psi(r,t=0)=\quantity(\frac{1}{4\pi\sigma^{2}})^{1/4}\exp(-\frac{r^{2}}{8\sigma^{2}}-i\frac{p_{0}}{\hbar}r).\
$ (5)

The wave functions $\phi$ and $\psi$ 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 $\sigma/\sqrt{2}$, and the reduced mass
wave packet has a larger width of $\sigma\sqrt{2}$. The corresponding
relations are illustrated in Fig. 2. In this frame the Hamiltonian decouples
as

$\hat{H}=\hat{H}_{R}+\hat{H}_{r}=\quantity(\frac{\hat{P}^{2}}{4m})+\quantity(\frac{\hat{p}^{2}}{m}-\frac{Gm^{2}}{L+\hat{r}}),$
(6)

where $\hat{P}=-i\hbar\partial/\partial R$ and
$\hat{p}=-i\hbar\partial/\partial 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.,
$\Psi(x_{A},x_{B},t)=\phi(R,t)\ \psi(r,t)$. Furthermore, the COM wave packet
evolves like a free particle, i.e., its Gaussianity is preserved [30, 31, 32].
The first two statistical moments characterize the quantum state fully, and
for completeness, they are given in Appendix B.

The state $\psi$ evolves in the gravitational potential, which we now expand
in the powers of the displacement-to-separation ratio $r/L$:

$V(\hat{r})=-\frac{Gm^{2}}{L+\hat{r}}\approx-\frac{1}{4}m\omega^{2}\sum_{n=0}^{N}\frac{(-1)^{n}}{L^{n-2}}\hat{r}^{n},$
(7)

where $N$ is the order of approximation, and we defined $\omega^{2}=4Gm/L^{3}$
for later convenience. In Appendix B, we derive exact analytical expressions
for the statistical moments of $\psi$ 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 $\psi$ by implementing Cayley’s form of
evolution operator [33]. The numerical evaluations for $\psi$ 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.

##  3 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.

###  3.1 Quadratic interactions

Figure 3: Accumulation of entanglement with the gravitational potential
truncated at the quadratic term ($N=2$). The configuration consists of
identical Osmium spheres with $m=0.25$ pg, $L=2.5$ times their radius, and
$\sigma=2.5$ nm. $p_{0}$ is the initial momentum. Analytical results are
calculated with the closed form of the covariance matrix in Eq. (8). $E$
denotes the logarithmic negativity, and $S$ is the entanglement entropy. The
values of $p_{0}/mL$ in the legends are written in multiples of
$6.18082292\times 10^{-3}$ s-1.

Consider first the gravitational potential truncated at the second order. We
obtained exact analytical forms for the independent elements of the covariance
matrix $\bm{\sigma}$. The solutions simplify if they are written in terms of
already introduced $\omega$ and in terms of $\omega_{0}=\hbar/2m\sigma^{2}$,
which is the frequency of harmonic trap for which the initial state is the
ground state:

$\displaystyle\bm{\sigma}_{00}=\frac{\hbar}{4m\omega_{0}}\quantity[2+\omega_{0}^{2}t^{2}+\quantity(1+\frac{\omega_{0}^{2}}{\omega^{2}})\sinh^{2}(\omega
t)],$
$\displaystyle\bm{\sigma}_{02}=\frac{\hbar}{4m\omega_{0}}\quantity[\omega_{0}^{2}t^{2}-\quantity(1+\frac{\omega_{0}^{2}}{\omega^{2}})\sinh^{2}(\omega
t)],$
$\displaystyle\bm{\sigma}_{11}=\frac{m\hbar\omega_{0}}{4}\quantity[2+\quantity(1+\frac{\omega^{2}}{\omega_{0}^{2}})\sinh^{2}(\omega
t)],$
$\displaystyle\bm{\sigma}_{13}=-\frac{m\hbar\omega_{0}}{4}\quantity(1+\frac{\omega^{2}}{\omega_{0}^{2}})\sinh^{2}(\omega
t),$
$\displaystyle\bm{\sigma}_{01}=\frac{\hbar}{8}\quantity[2\omega_{0}t+\quantity(\frac{\omega_{0}}{\omega}+\frac{\omega}{\omega_{0}})\sinh(2\omega
t)],$
$\displaystyle\bm{\sigma}_{03}=\frac{\hbar}{8}\quantity[2\omega_{0}t-\quantity(\frac{\omega_{0}}{\omega}+\frac{\omega}{\omega_{0}})\sinh(2\omega
t)].$ (8)

The logarithmic negativity for $p_{0}=0$, in the regime of
$\omega\ll\omega_{0}$ and $\omega t\ll 1$, was already approximated to [11]

$E(\bm{\sigma})\approx-\log_{2}\sqrt{1+2\Omega^{6}t^{6}-2\Omega^{3}t^{3}\sqrt{1+\Omega^{6}t^{6}}},$
(9)

where $\Omega^{3}=\omega_{0}\omega^{2}/6\equiv\hbar G/3\sigma^{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 $\omega_{0}\approx 25\omega$.

###  3.2 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}^{\prime}=m\omega^{2}/2$, and therefore it fully describes entanglement
in Eq. (9) since now
$\Omega^{3}=\omega_{0}\omega^{2}/6\equiv(\omega_{0}/3m)F_{2}^{\prime}$. In the
following section we provide further evidence for the pivotal role of force
gradient in the entanglement dynamics due to higher-order interactions.

(a)

(b)

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 $p_{0}=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 $p_{0}$ (masses moving towards each other) with Eqs. (11) and (12).
The values of $p_{0}/mL$ in the legends are written in multiples of
$6.18082292\times 10^{-3}$ s-1.

###  3.3 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
gradient is given by $F_{3}^{\prime}({\hat{r}})=(1-3\hat{r}/L)m\omega^{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 $\expectationvalue{\hat{r}}<0$, the gradient
increases, matching the growing entanglement. Conversely, when the masses move
away, $p_{0}<0$ and $\expectationvalue{\hat{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:
$\hat{F}_{3}^{\prime}(\hat{r})=(1-3\hat{r}/L)F_{2}^{\prime}$. The average
factor therefore reads

$\frac{\expectationvalue{F_{3}^{\prime}}}{F_{2}^{\prime}}=1-\frac{3}{L}\expectationvalue{\hat{r}}\approx
1+\frac{6p_{0}t}{mL}\equiv 1+\epsilon_{3}(t),$ (10)

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: $\expectationvalue{\hat{r}}\approx
r_{\text{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

$\displaystyle S(\rho_{A})=\Big{[}1+\epsilon_{3}(t)\Big{]}\
S(\rho_{A},p_{0}=0),$ (11) $\displaystyle
E(\bm{\sigma})=\Bigg{[}1+\frac{1}{2}\epsilon_{3}(t)\Bigg{]}\
E(\bm{\sigma},p_{0}=0).$ (12)

Note that the factor of $1/2$ next to $\epsilon_{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

$\frac{S(\rho_{A})}{S(\rho_{A},p_{0}=0)}=\frac{\expectationvalue{F_{N}^{\prime}}}{F_{2}^{\prime}}=1+\sum_{n=3}^{N}\epsilon_{n}(t),$
(13)

where the corrections for gravity-like interactions (inverse-square forces)
arising due to the $n^{\text{th}}$ term in the potential expansion is

$\epsilon_{n}(t)=\frac{(-1)^{n}}{2L^{n-2}}\
n(n-1)\expectationvalue{\hat{r}^{n-2}}.$ (14)

Since the gravitational force between two quantum masses is weak, for the
estimation of $\expectationvalue{\hat{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:

$\absolutevalue{\psi_{0}(r,t)}^{2}\approx\frac{1}{\bm{\Delta}r_{0}\sqrt{2\pi}}\exp\quantity(-\frac{\quantity(r-r_{\text{cl}})^{2}}{2\bm{\Delta}r_{0}^{2}}),$
(15)

where $\bm{\Delta}r_{0}^{2}=2\sigma^{2}\quantity(1+\omega_{0}^{2}t^{2})$. With
this approximation one obtains the correction terms for $n\geq 3$ as

$\displaystyle\epsilon_{n}(t)=\frac{(-1)^{n}}{2\sqrt{\pi}L^{n-2}}\
n(n-1)\hskip 71.13188pt$ (16)
$\displaystyle\times\sum_{m=0,2,}^{n-2}{n-2\choose m}r_{\text{cl}}^{n-m-2}\
\quantity(\sqrt{2}\bm{\Delta}r_{0})^{m}\ \Gamma\quantity(\frac{m+1}{2}),$

where $\Gamma$ is the gamma function, and the summation is only over even $m$.
Note that the gravitational interaction is already included in
$F_{2}^{\prime}$, and the present estimation is for the ratio of the force
gradients of different orders only,
$\expectationvalue{F_{N}^{\prime}}/F_{2}^{\prime}$. Since $\epsilon_{n}\propto
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

$\epsilon_{4}(t)=24\frac{p_{0}^{2}t^{2}}{m^{2}L^{2}}+12\frac{\sigma^{2}}{L^{2}}\quantity(1+\omega_{0}^{2}t^{2}).$
(17)

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 $\epsilon_{4}$ also depends on the position spread, hence it might
be important even for stationary configurations where the wave packet
undergoes a fast expansion.

Figure 5: Comparison of entanglement accumulated with the gravitational
potential expanded up to quadratic ($N=2$), cubic ($N=3$) and quartic ($N=4$)
term, respectively. Solid lines show the results for positive momentum (masses
moving toward each other), and dashed lines are for negative momentum. The
dots represent the entanglement entropy ($S$) computed with the closed
formulae derived in this work. 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 $p_{0}/mL$
in the legends are written in multiples of $6.18082292\times 10^{-3}$ s-1.

###  3.4 Galilean relativity and a drifting COM

We made a change of reference frames to dissect the bipartite evolution into
two independent single-particle 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 non-zero 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.

##  4 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 $\tilde{\mu}_{3}$ in the evolution of the reduced
mass wave function $\psi$. While $\tilde{\mu}_{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 $\psi$ 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 $\tilde{\mu}_{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 $\approx 0.001,\ 0.002$, and $0.003\%$, for an initial
spread of $\sigma=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 third-order 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
$\ddot{\expectationvalue{p}}/\expectationvalue{p}=\omega^{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\approx
1.75\times 10^{-4}$ is relaxed from 5 seconds to 4 seconds with an initial
momentum of $p_{0}/mL\approx+0.022$ s-1. Note that an entanglement detection
scheme achieving a precision of $E\sim 10^{-4}$ has recently been put forward
in Ref. [41].

Figure 6: Non-Gaussian reduced mass dynamics. The physical situation is as in
Fig. 3. Skewness ($\tilde{\mu}_{3}$) is computed for the position space
distribution. $N$ denotes the order of approximation, see Eq. (7). The left
panel is for particles initially at rest, and the right panel is for masses
moving toward each other. The values of $p_{0}/mL$ in the legends are written
in multiples of $6.18082292\times 10^{-3}$ s-1.

##  5 Versatility

While our discussions were mainly focused on gravity-induced 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 $\omega$ and $\epsilon_{3}$. As we derived, for identical
masses coupled via Newtonian gravity,

$\omega^{2}=\frac{4Gm}{L^{3}},\hskip
14.22636pt\epsilon_{3}(t)=\frac{6p_{0}t}{mL}.$ (18)

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\omega^{2}/4$, and $\epsilon_{3}$ is to be calculated by
comparing the force gradients as
$\expectationvalue{F_{3}^{\prime}}/F_{2}^{\prime}=1+\epsilon_{3}(t)$. One can
then verify that for the Coulomb potential between charges $q_{1}$ and $q_{2}$
embedded into the masses, we obtain

$\omega^{2}=\frac{4q_{1}q_{2}\alpha\hbar c}{e^{2}mL^{3}},\hskip
14.22636pt\epsilon_{3}(t)=\frac{6p_{0}t}{mL},$ (19)

where $\alpha$ 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

$\omega^{2}=\frac{\pi^{3}\hbar cR_{0}}{120m(L-2R_{0})^{4}},\hskip
5.69046pt\epsilon_{3}(t)=\frac{8p_{0}t}{m(L-2R_{0})},$ (20)

where $R_{0}$ is the radius of each sphere.

For an arbitrary central interaction with a potential of the form

$V(x_{A},x_{B})=-\frac{C}{(X+x_{B}-x_{A})^{j}},$ (21)

we obtain

$\omega^{2}=\frac{2j(j+1)C}{mX^{j+2}},\hskip
14.22636pt\epsilon_{3}(t)=\frac{2(j+2)p_{0}t}{mX}.$ (22)

In some situations the force is known, but solving for the potential is
difficult or uncertain due to non-unique 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\omega^{2}/2$. For example, if the force is of
the form

$F(x_{A},x_{B})=-\frac{C}{(X+x_{B}-x_{A})^{j}},$ (23)

one arrives at

$\omega^{2}=\frac{2jC\ }{mX^{j+1}},\hskip
14.22636pt\epsilon_{3}(t)=\frac{2(j+1)p_{0}t}{mX}.$ (24)

Note that the functional forms of $\epsilon_{3}$ in this section are valid
only for weak interactions. For stronger potentials one has to take a step
back and use $\epsilon_{3}=-3\expectationvalue{\hat{r}}/L$, where
$\expectationvalue{\hat{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

$V=\sum_{k}V_{k}(x_{B}-x_{A}),\hskip 7.11317ptF=\sum_{k}F_{k}(x_{B}-x_{A}),$
(25)

the total $\omega$ characterising the Gaussian covariance matrix is given by a
Pythagoras-like theorem, and the equivalent $\epsilon_{3}$ governing the
entanglement amplification due to the cubic-order term is calculated to be a
weighted sum:

$\omega^{2}=\sum_{k}\omega^{2}_{k},\hskip
14.22636pt\epsilon_{3}(t)=\frac{1}{\omega^{2}}\sum_{k}\omega_{k}^{2}\epsilon_{k3}(t).$
(26)

where $\omega_{k}$ and $\epsilon_{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.

##  6 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.

## Acknowledgements

This work is jointly supported by (i) NAWA, Poland, via project
PPN/PPO/2018/1/00007/U/00001, (ii) XMUM, Malaysia, via project
XMUMRF/2022-C10/IPHY/0002, and (iii) DORA office of IIT Roorkee, India. A.K.
thanks the IIT Roorkee Heritage Foundation, USA, for the ‘Pledge a Dream’
grant. Discussions with Andy Chia (CQT-NUS, Singapore) are acknowledged. T.K.
thanks Timothy Liew for their hospitality at NTU, Singapore, and Yvonne Gao
for their hospitality at NUS, Singapore. We acknowledge the National
Supercomputing Mission (NSM) for providing computing resources of ‘PARAM
Ganga’ at IIT Roorkee, India, which is implemented by C-DAC and supported by
MeitY and DST, Govt. of India. The authors gratefully thank the reviewers for
their comments and suggestions which have significantly improved the
presentation of this article.

##  Appendix A Coordinate transformations

The time-dependent Schrödinger equation (TDSE) for two identical particles $A$
and $B$ is given by

$\displaystyle\quantity(-\frac{\hbar^{2}}{2m}\partialderivative[2]{x_{A}}-\frac{\hbar^{2}}{2m}\partialderivative[2]{x_{B}}+V(\hat{x}_{A},\hat{x}_{B}))\Psi(x_{A},x_{B},t)$
$\displaystyle=i\hbar\partialderivative{t}\Psi(x_{A},x_{B},t),\hskip
14.22636pt$ (27)

where $\hat{p}_{A}=-i\hbar\partial/\partial x_{A}$ and
$\hat{p}_{B}=-i\hbar\partial/\partial x_{B}$ are their respective momentum
operators. Coordinate transformations between the LAB the COM frame of
reference are defined by

$R=\frac{x_{A}+x_{B}}{2},\hskip 14.22636ptr=x_{B}-x_{A},$ (28)

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

$P=p_{A}+p_{B},\hskip 14.22636ptp=\frac{p_{B}-p_{A}}{2}.$ (29)

Accordingly, the inverse transformations are

$x_{A}(x_{B})=R-\\!(\\!+\\!)\ \frac{r}{2},\hskip
14.22636ptp_{A}(p_{B})=\frac{P}{2}-\\!(\\!+\\!)\ p.$ (30)

The displacements $x_{A}$ and $x_{B}$ are functions of $R$ and $r$, and the
rules of differentiation imply

$\partialderivative[2]{x_{A}}\quantity(\partialderivative[2]{x_{B}})=\frac{1}{4}\partialderivative[2]{R}+\partialderivative[2]{r}-\\!(\\!+\\!)\partialderivative{}{R}{r},$
(31)

which means that the kinetic energy is equivalent to

$-\frac{\hbar^{2}}{2m}\partialderivative[2]{x_{A}}-\frac{\hbar^{2}}{2m}\partialderivative[2]{x_{B}}=-\frac{\hbar^{2}}{4m}\partialderivative[2]{R}-\frac{\hbar^{2}}{m}\partialderivative[2]{r}.$
(32)

For central potentials $V(x_{A},x_{B})=V(x_{B}-x_{A})\equiv V(r)$, and hence
the TDSE is transformed to

$\displaystyle\quantity(-\frac{\hbar^{2}}{4m}\partialderivative[2]{R}-\frac{\hbar^{2}}{m}\partialderivative[2]{r}+V(\hat{r}))\Psi(x_{A},x_{B},t)$
$\displaystyle=i\hbar\partialderivative{t}\Psi(x_{A},x_{B},t).$ (33)

In this work the initial wave function is $\Psi(x_{A},x_{B},t=0)=\phi(R,t=0)\
\psi(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:

$\displaystyle-\frac{\hbar^{2}}{4m}\partialderivative[2]{R}\phi(R,t)=i\hbar\partialderivative{t}\phi(R,t),$
(34)
$\displaystyle\quantity(-\frac{\hbar^{2}}{m}\partialderivative[2]{r}+V(\hat{r}))\psi(r,t)=i\hbar\partialderivative{t}\psi(r,t),$
(35)

where $\hat{P}=-i\hbar\partial/\partial R$ and
$\hat{p}=-i\hbar\partial/\partial 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

$\Psi(x_{A},x_{B},t)=\phi\quantity(\frac{x_{A}+x_{B}}{2},t)\
\psi\quantity(x_{B}-x_{A},t).$ (36)

##  Appendix B The Ehrenfest dynamics

Ehrenfest’s theorem relates the time derivative of the expectation value of an
operator $\hat{A}$ to the expectation of its commutator with the Hamiltonian
$\hat{H}$:

$\derivative{t}\expectationvalue{\hat{A}}=\frac{1}{i\hbar}\expectationvalue{\commutator{\hat{A}}{\hat{H}}}+\expectationvalue{\partialderivative{\hat{A}}{t}}.$
(37)

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
$\text{Det}(\bm{\sigma}_{R})=\text{Det}(\bm{\sigma}_{r})=(\hbar/2)^{2}$ [43].

###  B.1 Evolution of the COM

The COM evolution corresponds to the free expansion of a Gaussian wave packet:
$\hat{H}_{R}=\hat{P}^{2}/4m$. The initial is characterised by
$\expectationvalue{\hat{R}}=0$, $\expectationvalue{\hat{P}}=0$,
$\expectationvalue{\anticommutator{\hat{R}}{\hat{P}}}=0$,
$\expectationvalue{\hat{R}^{2}}=\sigma^{2}/2$, and
$\expectationvalue{\hat{P}^{2}}=\hbar^{2}/2\sigma^{2}$ [see Fig. 2 or Eq.
(4)]. The solution to the corresponding Ehrenfest’s differential equations
for the first two statistical moments,

$\displaystyle\derivative{t}\expectationvalue{\hat{R}}=\frac{1}{2m}\expectationvalue{\hat{P}},$
$\displaystyle\derivative{t}\expectationvalue{\hat{P}}=0,$
$\displaystyle\derivative{t}\expectationvalue{\anticommutator{\hat{R}}{\hat{P}}}=\frac{1}{m}\expectationvalue{P^{2}},$
$\displaystyle\derivative{t}\expectationvalue{\hat{R}^{2}}=\frac{1}{2m}\expectationvalue{\anticommutator{\hat{R}}{\hat{P}}},$
$\displaystyle\derivative{t}\expectationvalue{\hat{P}^{2}}=0,$ (38)

imply

$\displaystyle\bm{\Delta}R^{2}=\frac{1}{2}\sigma^{2}(1+\omega_{0}^{2}t^{2}),$
$\displaystyle\bm{\Delta}P^{2}=\frac{\hbar^{2}}{2\sigma^{2}},$
$\displaystyle\textbf{Cov}({R},{P})=\frac{1}{2}\hbar\omega_{0}t.$ (39)

Alternatively, one arrives at the same results by utilizing the functional
form of the wave function [40]:

$\phi(R,t)=\frac{1}{\sqrt{\sigma(1+i\omega_{0}t)\sqrt{\pi}}}\exp\quantity(-\frac{R^{2}}{2\sigma^{2}(1+i\omega_{0}t)}).$
(40)

###  B.2 Evolution of the reduced mass

The reduced mass Hamiltonian can be represented by a binomial series as in Eq.
(7). The initial state is characterised by $\expectationvalue{\hat{r}}=0$,
$\expectationvalue{\hat{p}}=p_{0}$,
$\expectationvalue{\anticommutator{\hat{r}}{\hat{p}}}=0$,
$\expectationvalue{\hat{r}^{2}}=2\sigma^{2}$, and
$\expectationvalue{\hat{p}^{2}}=p_{0}^{2}+\hbar^{2}/8\sigma^{2}$ [see Fig. 2
or Eq. (5)]. For $N=2$, i.e., a quadratic Hamiltonian,

$\hat{H}_{r}=\frac{\hat{p}^{2}}{m}-\frac{1}{4}m\omega^{2}\quantity(L^{2}-Lr+r^{2}),$
(41)

the relevant Ehrenfest’s coupled differential equations for the first two
statistical moments are given by

$\displaystyle\derivative{t}\expectationvalue{\hat{r}}=\frac{2}{m}\expectationvalue{\hat{p}},$
$\displaystyle\derivative{t}\expectationvalue{\hat{p}}=-\frac{1}{4}m\omega^{2}L+\frac{1}{2}m\omega^{2}\expectationvalue{\hat{r}},$
$\displaystyle\derivative{t}\expectationvalue{\anticommutator{\hat{r}}{\hat{p}}}=\frac{4}{m}\expectationvalue{\hat{p}^{2}}-\frac{1}{2}m\omega^{2}L\expectationvalue{\hat{r}}+m\omega^{2}\expectationvalue{\hat{r}^{2}},$
$\displaystyle\derivative{t}\expectationvalue{r^{2}}=\frac{2}{m}\expectationvalue{\anticommutator{\hat{r}}{\hat{p}}},$
$\displaystyle\derivative{t}\expectationvalue{\hat{p}^{2}}=\frac{1}{2}m\omega^{2}\expectationvalue{\anticommutator{\hat{r}}{\hat{p}}}-\frac{1}{2}m\omega^{2}L\expectationvalue{\hat{p}},$
(42)

Their exact analytical solution reads

$\displaystyle\expectationvalue{\hat{r}}$ $\displaystyle=$
$\displaystyle\frac{1}{2}L\Big{(}1-\cosh(\omega
t)\Big{)}-\frac{2p_{0}}{m\omega}\sinh(\omega t),$
$\displaystyle\expectationvalue{\hat{p}}$ $\displaystyle=$ $\displaystyle-
p_{0}\cosh(\omega t)-\frac{1}{4}m\omega L\sinh(\omega t),$
$\displaystyle\expectationvalue{\anticommutator{\hat{r}}{\hat{p}}}$
$\displaystyle=$ $\displaystyle Lp_{0}\Big{(}\cosh(2\omega t)-\cosh(\omega
t)\Big{)}$
$\displaystyle+\frac{2}{m\omega}\quantity(p_{0}^{2}+\frac{\hbar^{2}}{8\sigma^{2}}+\frac{1}{2}m^{2}\omega^{2}\sigma^{2})\sinh(2\omega
t)$ $\displaystyle+\frac{1}{8}m\omega L^{2}\Big{(}\sinh(2\omega
t)-2\sinh(\omega t)\Big{)},$ $\displaystyle\expectationvalue{\hat{r}^{2}}$
$\displaystyle=$ $\displaystyle 2\sigma^{2}\Big{(}1+\sinh[2](\omega t)\Big{)}$
$\displaystyle+\frac{1}{8}L^{2}\Big{(}3+\cosh(2\omega t)-4\cosh(\omega
t)\Big{)}$ $\displaystyle+\frac{Lp_{0}}{m\omega}\Big{(}\sinh(2\omega
t)-2\sinh(\omega t)\Big{)}$
$\displaystyle+\frac{4}{m^{2}\omega^{2}}\quantity(p_{0}^{2}+\frac{\hbar^{2}}{8\sigma^{2}})\sinh[2](\omega
t),$ $\displaystyle\expectationvalue{\hat{p}^{2}}$ $\displaystyle=$
$\displaystyle\quantity(p_{0}^{2}+\frac{\hbar^{2}}{8\sigma^{2}})\Big{(}1+\sinh[2](\omega
t)\Big{)}$ (43)
$\displaystyle+\frac{1}{4}m^{2}\omega^{2}\quantity(2\sigma^{2}+\frac{1}{4}L^{2})\sinh[2](\omega
t)$ $\displaystyle+\frac{1}{4}m\omega Lp_{0}\sinh(2\omega t),$

which implies that variances and the correlation are

$\displaystyle\bm{\Delta}r^{2}=2\sigma^{2}\quantity(\cosh[2](\omega
t)+\frac{\omega_{0}^{2}}{\omega^{2}}\sinh[2](\omega t)),$
$\displaystyle\bm{\Delta}p^{2}=\frac{\hbar^{2}}{8\sigma^{2}}\quantity(\cosh[2](\omega
t)+\frac{\omega^{2}}{\omega_{0}^{2}}\sinh[2](\omega t)),$
$\displaystyle\textbf{Cov}({r},{p})=\frac{\hbar}{4}\quantity(\frac{\omega_{0}}{\omega}+\frac{\omega}{\omega_{0}})\sinh(2\omega
t).$ (44)

##  Appendix 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
$\hat{u}=(\hat{x}_{A},\hat{p}_{A},\hat{x}_{B},\hat{p}_{B})^{T}$, the
covariance matrix is defined as [27, 28, 29]:

$\bm{\sigma}_{jk}=\frac{1}{2}\expectationvalue{\anticommutator{\hat{u}_{j}}{\hat{u}_{k}}}-\expectationvalue{\hat{u}_{j}}\expectationvalue{\hat{u}_{k}}.$
(45)

In the block form we can write

$\bm{\sigma}\equiv\matrixquantity(\bm{\alpha}&\bm{\gamma}\\\
\bm{\gamma}^{T}&\bm{\beta}),$ (46)

where $\bm{\alpha}(\bm{\beta}$) contains the local mode correlation for
$A(B)$, and $\bm{\gamma}$ describes the intermodal correlation. In our setting
the local modes are identical, i.e., $\bm{\alpha}=\bm{\beta}$, and a
coordinate change to the COM frame implies

$\displaystyle\bm{\sigma}_{00}(\bm{\sigma}_{02})=\bm{\Delta}R^{2}+\\!(\\!-\\!)\
\frac{1}{4}\bm{\Delta}r^{2},$
$\displaystyle\bm{\sigma}_{11}(\bm{\sigma}_{13})=\frac{1}{4}\bm{\Delta}P^{2}+\\!(\\!-\\!)\
\bm{\Delta}p^{2},$
$\displaystyle\bm{\sigma}_{01}(\bm{\sigma}_{03})=\frac{1}{2}\textbf{Cov}({R},{P})+\\!(\\!-\\!)\
\frac{1}{2}\textbf{Cov}({r},{p}).$ (47)

Given the symmetry of the problem we have $\bm{\sigma}_{22}=\bm{\sigma}_{00}$,
$\bm{\sigma}_{33}=\bm{\sigma}_{11}$, $\bm{\sigma}_{23}=\bm{\sigma}_{01}$,
$\bm{\sigma}_{12}=\bm{\sigma}_{03}$. The rest of the elements are constrained
due to the symmetry property $\bm{\sigma}_{jk}=\bm{\sigma}_{kj}$.

###  C.2 Logarithmic negativity

The negativity of partially transposed density matrix is a necessary and
sufficient condition for entanglement in two–mode Gaussian states [44]. As a
result of partial transposition, the covariance matrix is transformed to
$\tilde{\bm{\sigma}}$, which differs from $\bm{\sigma}$ by a sign-flip of
$\text{Det}(\bm{\gamma})$ [29]. The symplectic eigenvalues of the covariance
matrix, $\tilde{\nu}_{\pm}(\bm{\sigma})$, are given by [28, 27]

$2\tilde{\nu}^{2}_{\pm}(\bm{\sigma})=\tilde{\Sigma}(\bm{\sigma})\pm\sqrt{\tilde{\Sigma}^{2}(\bm{\sigma})-4\
\text{Det}\quantity(\bm{\sigma})},$ (48)

where
$\tilde{\Sigma}(\bm{\sigma})=\text{Det}(\bm{\alpha})+\text{Det}(\bm{\beta})-2\
\text{Det}(\bm{\gamma})$. For the symmetric problem as considered in this work
the local modes $\bm{\alpha}$ and $\bm{\beta}$ are identical, and hence
$\tilde{\Sigma}(\bm{\sigma})=2\
\quantity[\text{Det}(\bm{\alpha})-\text{Det}(\bm{\gamma})]$. Entanglement is
quantified by the minimum symplectic eigenvalue via logarithmic negativity:

$E(\bm{\sigma})=\max\Bigg{[}0,\
-\log_{2}\quantity(\frac{\tilde{\nu}_{-}(\bm{\sigma})}{\hbar/2})\Bigg{]}.$
(49)

###  C.3 Entropy of entanglement

For a pure bipartite system described by a density matrix $\rho_{AB}$, the
entanglement entropy is defined as the von Neumann entropy for any one of the
subsystems, e.g., $S(\rho_{A})=-\Tr\quantity[\rho_{A}\log_{2}(\rho_{A})]$,
where $\rho_{A}=\Tr_{B}\quantity(\rho_{AB})$ is the reduced density matrix for
subsystem $A$. In order to calculate $S(\rho_{A})$ we start with the two-body
wave function of Eq. (36). The COM wave function $\phi(R,t)$ is derived
analytically in Eq. (40), and we calculate $\psi(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]:

$\Psi(x_{A},x_{B},t)=\sum_{j}\sqrt{\lambda_{j}(t)}\ \chi^{(A)}_{j}(x_{A},t)\
\chi^{(B)}_{j}(x_{B},t),$ (50)

where $\quantity{\chi^{(A)}_{j}}$ and $\quantity{\chi^{(B)}_{j}}$ are
orthonormal states in subsystems $A$ and $B$, respectively, and $\lambda_{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-\absolutevalue{\Psi(x_{A},x_{B},t)}^{2}=1-\sum_{j}\lambda_{j}(t)\lesssim
10^{-7}$. With this decomposition, the entanglement entropy reduces to

$S(\rho_{A})=-\sum_{j}\lambda_{j}\log_{2}(\lambda_{j}).$ (51)

In the case of Gaussian evolution, the entanglement entropy is calculable
using the covariance matrix [43]:

$S(\bm{\alpha})=f\quantity(\frac{1}{\hbar}\sqrt{\text{Det}(\bm{\alpha})})\equiv
S(\rho_{A}),$ (52)

where

$f(x)=\quantity(x\\!+\frac{1}{2})\log_{2}\quantity(x\\!+\frac{1}{2})-\quantity(x\\!-\frac{1}{2})\log_{2}\quantity(x\\!-\frac{1}{2}).$
(53)

##  Appendix D Numerical details

Numerical calculations are performed in natural units of $c=1$, hence the
conversion constant $\hbar c=197.3269804$ keV pm. The density of Osmium is
$22.5872$ g/cm3. An error analysis implies that, in the numerical time
evolution of the reduced mass wave function, a grid size of $\lesssim 0.2$ pm
with a time step of $\lesssim 10\ \mu$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\ \mu$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.

## References

  * [1] R. V. Pound and G. A. Rebka.  “Apparent Weight of Photons”.  Physical Review Letters 4, 337 (1960).
  * [2] J. C. Hafele and R. E. Keating.  “Around-the-World Atomic Clocks: Observed Relativistic Time Gains”.  Science 177, 168 (1972).
  * [3] R. Colella, A. W. Overhauser, and S. A. Werner.  “Observation of Gravitationally Induced Quantum Interference”.  Physical Review Letters 34, 1472 (1975).
  * [4] A. Peters, K. Y. Chung, and S. Chu.  “Measurement of gravitational acceleration by dropping atoms”.  Nature 400, 849 (1999).
  * [5] V. V. Nesvizhevsky, H. G. Börner, A. K. Petukhov, H. Abele, S. Baeßler, F. J. Rueß, T. Stöferle, A. Westphal, A. M. Gagarski, G. A. Petrov, and A. V. Strelkov.  “Quantum states of neutrons in the Earth’s gravitational field”.  Nature 415, 297 (2002).
  * [6] P. Asenbaum, C. Overstreet, T. Kovachy, D. D. Brown, J. M. Hogan, and M. A. Kasevich.  “Phase Shift in an Atom Interferometer due to Spacetime Curvature across its Wave Function”.  Physical Review Letters 118, 183602 (2017).
  * [7] S. Bose, A. Mazumdar, G. W. Morley, H. Ulbricht, M. Toroš, M. Paternostro, A. A. Geraci, P. F. Barker, M. S. Kim, and G. Milburn.  “Spin Entanglement Witness for Quantum Gravity”.  Physical Review Letters 119, 240401 (2017).
  * [8] C. Marletto and V. Vedral.  “Gravitationally Induced Entanglement between Two Massive Particles is Sufficient Evidence of Quantum Effects in Gravity”.  Physical Review Letters 119, 240402 (2017).
  * [9] A. Al Balushi, W. Cong, and R. B. Mann.  “Optomechanical quantum Cavendish experiment”.  Physical Review A 98, 043811 (2018).
  * [10] T. Krisnanda, M. Zuppardo, M. Paternostro, and T. Paterek.  “Revealing Nonclassicality of Inaccessible Objects”.  Physical Review Letters 119, 120402 (2017).
  * [11] T. Krisnanda, G. Y. Tham, M. Paternostro, and T. Paterek.  “Observable quantum entanglement due to gravity”.  npj Quantum Information 6, 12 (2020).
  * [12] S. Qvarfort, S. Bose, and A. Serafini.  “Mesoscopic entanglement through central–potential interactions”.  Journal of Physics B: Atomic, Molecular and Optical Physics 53, 235501 (2020).
  * [13] S. Rijavec, M. Carlesso, A. Bassi, V. Vedral, and C. Marletto.  “Decoherence effects in non-classicality tests of gravity”.  New Journal of Physics 23, 043040 (2021).
  * [14] K. Kustura, C. Gonzalez-Ballestero, Andrés de los Ríos Sommer, N. Meyer, R. Quidant, and O. Romero-Isart.  “Mechanical Squeezing via Unstable Dynamics in a Microcavity”.  Physical Review Letters 128, 143601 (2022).
  * [15] T. W. van de Kamp, R. J. Marshman, S. Bose, and A. Mazumdar.  “Quantum gravity witness via entanglement of masses: Casimir screening”.  Physical Review A 102, 062807 (2020).
  * [16] R. J. Marshman, A. Mazumdar, R. Folman, and S. Bose.  “Constructing nano-object quantum superpositions with a Stern-Gerlach interferometer”.  Physical Review Research 4, 023087 (2022).
  * [17] M. Aspelmeyer, T. J. Kippenberg, and F. Marquardt.  “Cavity optomechanics”.  Reviews of Modern Physics 86, 1391 (2014).
  * [18] B. Abbott et al.  “Observation of a kilogram-scale oscillator near its quantum ground state”.  New Journal of Physics 11, 073032 (2009).
  * [19] C. Whittle et al.  “Approaching the motional ground state of a 10-kg object”.  Science 372, 1333 (2021).
  * [20] J. Chan, T. P. M. Alegre, A. H. Safavi-Naeini, J. T. Hill, A. Krause, S. Gröblacher, M. Aspelmeyer, and O. Painter.  “Laser cooling of a nanomechanical oscillator into its quantum ground state”.  Nature 478, 89 (2011).
  * [21] T. A. Palomaki, J. D. Teufel, R. W. Simmonds, and K. W. Lehnert.  “Entangling Mechanical Motion with Microwave Fields”.  Science 342, 710 (2013).
  * [22] N. Fiaschi, B. Hensen, A. Wallucks, R. Benevides, J. Li, T. P. M. Alegre, and S. Gröblacher.  “Optomechanical quantum teleportation”.  Nature Photonics 15, 817 (2021).
  * [23] R. Riedinger, A. Wallucks, Igor Marinković, C. Löschnauer, M. Aspelmeyer, S. Hong, and S. Gröblacher.  “Remote quantum entanglement between two micromechanical oscillators”.  Nature 556, 473 (2018).
  * [24] A. Datta and H. Miao.  “Signatures of the quantum nature of gravity in the differential motion of two masses”.  Quantum Science and Technology 6, 045014 (2021).
  * [25] J. S. Pedernales, K. Streltsov, and M. B. Plenio.  “Enhancing Gravitational Interaction between Quantum Systems by a Massive Mediator”.  Physical Review Letters 128, 110401 (2022).
  * [26] F. Roccati, B. Militello, E. Fiordilino, R. Iaria, L. Burderi, T. Di Salvo, and F. Ciccarello.  “Quantum correlations beyond entanglement in a classical-channel model of gravity”.  Scientific Reports 12, 17641 (2022).
  * [27] G. Vidal and R. F. Werner.  “Computable measure of entanglement”.  Physical Review A 65, 032314 (2002).
  * [28] G. Adesso, A. Serafini, and F. Illuminati.  “Extremal entanglement and mixedness in continuous variable systems”.  Physical Review A 70, 022318 (2004).
  * [29] G. Adesso and F. Illuminati.  “Gaussian measures of entanglement versus negativities: Ordering of two-mode gaussian states”.  Physical Review A 72, 032334 (2005).
  * [30] B. M. Garraway.  “Extended Gaussian wavepacket dynamics”.  Journal of Physics B: Atomic, Molecular and Optical Physics 33, 4447 (2000).
  * [31] C. Weedbrook, S. Pirandola, R. García-Patrón, N. J. Cerf, T. C. Ralph, J. H. Shapiro, and S. Lloyd.  “Gaussian quantum information”.  Reviews of Modern Physics 84, 621 (2012).
  * [32] M. A. Schlosshauer.  “Decoherence and the Quantum-To-Classical Transition”.  Springer-Verlag, Germany. (2007).  1st edition.
  * [33] W. H. Press, S. A. Teukolsky, H. A. Bethe, W. T. Vetterling, and B. P. Flannery.  “Numerical Recipes”.  Cambridge University Press, USA. (2007).  3rd edition.
  * [34] A. Kumar.  “An accurate pentadiagonal matrix solution for the time-dependent Schrödinger equation”.  Zenodo/7275668 (2022).
  * [35] A. Kumar and P. Arumugam.  “An accurate pentadiagonal matrix solution for the time-dependent Schrödinger equation”.  arXiv:2205.13467 (2022).
  * [36] A. Kumar, T. Krisnanda, P. Arumugam, and T. Paterek.  “Nonclassical trajectories in head-on collisions”.  Quantum 5, 506 (2021).
  * [37] H. Miao, D. Martynov, H. Yang, and A. Datta.  “Quantum correlations of light mediated by gravity”.  Physical Review A 101, 063804 (2020).
  * [38] S. Bose, A. Mazumdar, M. Schut, and M. Toroš.  “Mechanism for the quantum natured gravitons to entangle masses”.  Physical Review D 105, 106028 (2022).
  * [39] R. J. Marshman, A. Mazumdar, and S. Bose.  “Locality and entanglement in table-top testing of the quantum nature of linearized gravity”.  Physical Review A 101, 052110 (2020).
  * [40] S. M. Blinder.  “Evolution of a Gaussian Wavepacket”.  American Journal of Physics 36, 525 (1968).
  * [41] T. Krisnanda, T. Paterek, M. Paternostro, and Timothy C. H. Liew.  “Quantum neuromorphic approach to efficient sensing of gravity-induced entanglement”.  Physical Review D 107, 086014 (2023).
  * [42] T. Emig, N. Graham, R. L. Jaffe, and M. Kardar.  “Casimir Forces between Arbitrary Compact Objects”.  Physical Review Letters 99, 170403 (2007).
  * [43] A. Serafini, F. Illuminati, and S. D. Siena.  “Symplectic invariants, entropic measures and correlations of Gaussian states”.  Journal of Physics B: Atomic, Molecular and Optical Physics 37, L21 (2003).
  * [44] R. Simon.  “Peres-Horodecki Separability Criterion for Continuous Variable Systems”.  Physical Review Letters 84, 2726 (2000).
  * [45] A. Ekert and P. L. Knight.  “Entangled quantum systems and the Schmidt decomposition”.  American Journal of Physics 63, 415 (1995).
  * [46] Ho Kwan Yuet.  “Quantum Entanglement in Continuous System”.  PhD thesis.  The Chinese University of Hong Kong, China.  (2004).
  * [47] C. Roberts, A. Milsted, M. Ganahl, A. Zalcman, B. Fontaine, Y. Zou, J. Hidary, G. Vidal, and S. Leichenauer.  “TensorNetwork: A Library for Physics and Machine Learning”.  arXiv:1905.01330 (2019).
  * [48] Google.  “GitHub/google/TensorNetwork”.
  * [49] Kwan-Yuet “Stephen” Ho.  “GitHub/stephenhky/pyqentangle”.