# Measurement-free error correction with coherent ancillas

Vickram N. Premakumar  Department of Physics, University of Wisconsin-Madison,
1150 University Avenue, Madison, WI, USA  M. Saffman  Department of Physics,
University of Wisconsin-Madison, 1150 University Avenue, Madison, WI, USA
Robert Joynt  Department of Physics, University of Wisconsin-Madison, 1150
University Avenue, Madison, WI, USA  Kavli Institute for Theoretical Sciences,
University of Chinese Academy of Sciences, Beijing 100190, China

###### Abstract

We show that coherent error correction with redundant syndrome extraction is
feasible even when realistic error models for the ancilla qubits are used.
This is done by means of a fault-tolerant circuit that has a favorable error
threshold for a conservative error model. It is optimal when Toffoli gates are
available, as they are in many implementations. The possibility of coherent
error correction is important for those quantum computing platforms in which
measurements in the course of the computation are undesirable.

##  I Introduction

Future large scale quantum computers based on existing qubits will require
quantum error correction. The actual error correction protocol depends
strongly on the platform used. For some implementations, measurement during
the course of a computation is not advisable. For example, in neutral atom
quantum computation the measurement time is much longer than other gate times
and crosstalk due to light scattering in the course of the measurement is
difficult to control Beterov and Saffman (2015). For spin qubits in
semiconductors, long measurement times also make traditional quantum error
correction less appealing Petta  _et al._ (2005); Thorgrimsson  _et al._
(2017). There is a known solution to this problem, which is coherent error
correction (CEC) Paz-Silva  _et al._ (2010); Li  _et al._ (2011); Crow  _et
al._ (2016) . CEC protocols do not require measurements to detect errors.
Instead the error information is transferred to ancilla qubits during a
syndrome detection cycle, and this quantum information is then used for
correction using conditional gates. The ancillas are then reset, which flushes
the entropy from the computer. Generally speaking, the overhead involved in
fault-tolerant CEC has been thought to be greater than in fault-tolerant
measurement-based error correction (MEC) Paz-Silva  _et al._ (2017). A recent
calculation suggested that this may not be the case Crow  _et al._ (2016),
which would make neutral-atom and electron-spin implementations more
attractive.

However, the calculation in Ref. Crow  _et al._ (2016) along these lines used
an error model that is not appropriate for the envisioned applications. It did
not allow for phase-flip errors on the ancilla qubits. Such phase errors will
certainly occur and should be given finite probability. In fact, we show below
that the circuit used in Ref. Crow  _et al._ (2016) is actually not fault-
tolerant when these errors are included, i.e., there is no error threshold. It
is also shown, however, that a modified circuit is fault-tolerant when all
errors are given reasonable probabilities. This circuit retains the innovative
feature of redundant syndrome extraction. We calculate the threshold and show
that it is somewhat smaller than in Ref. Crow  _et al._ (2016), but still
favorable.

##  II Circuit design

The error-correction circuit from Ref. Crow  _et al._ (2016) is shown in Fig.
1. Data qubits are in the top half and ancillas on the bottom. The circuit is
based on the Steane 7-qubit error correction code Steane (1996) but there are
two major differences. The first is that ancillas are used to detect the
error, as is standard in CEC. The second is that there is redundancy in the
set of stabilizer operators associated with the code. That is, instead of the
usual 6 operators that are usually measured for syndrome extraction, 14
diagnostic results are encoded in the ancillas. This redundancy builds fault
tolerance into the circuit, obviating the need for repetition of the syndrome
extraction.

If phase errors on the ancillas are included then we can get back-action on
the data qubits by a well-known mechanism Preskill (1998). This is shown in
Fig. 1, where a phase error on the 1st ancilla (the blue Z) produces data
errors in the data qubits that ultimately result in a logical error at the
end, indicated by the 3 red X’s at the top right of the figure. Thus the phase
errors spread through the computer, rendering it faulty: the probability of a
logical error is linearly proportional to the probability of a phase-flip
error on the ancillas. We calculate this probability below.

Figure 1: (color online) The coherent error correction circuit with redundant
stabilizer extraction proposed in Ref. Crow  _et al._ (2016). Moving from left
to right, the first quarter of the circuit extracts syndromes for bit flip
errors, the second quarter corrects them, while the third and fourth do the
same for phase errors. This circuit is not fault tolerant when phase errors on
the ancilla qubits are included. An example of a single phase error (blue Z)
in the ancilla block propagating (red X’s) to a logical error at the end is
shown. $R$ gates reset the ancilla qubits to $\ket{0}$.

To identify the single qubit errors that destroy the fault tolerance of the
previously proposed circuit we must consider the correction procedure for both
bit- and phase-flip type errors. While simulations of CSS codes typically take
advantage of the natural division of the stabilizer group into multiqubit
Pauli $X$ and $Z$ operators, doing so would blind us to the crucial fault of
this implementation. When a phase error happens on an ancilla qubit halfway
through a stabilizer measurement, the entangling gates propagate two phase
errors into the data register. These are rotated into bit-flip errors by the
time the circuit attempts to correct the errors. As the stabilizers form a
group, these two errors together are interpreted as a third: the bitwise sum
of the two syndrome strings of the propagated phase errors is equal to the
syndrome of another, and therefore the circuit erroneously corrects a qubit
which had no error, and we have a fatal logical error. After identifying the
dangerous errors, we can proceed to reconfigure the circuit to avoid them.

To do this we must insist that there is no way for an error to happen during
the syndrome extraction for a single stabilizer element. By using a textbook
circuit identity, we can re-order the gates and flip the role of target and
control qubits during extraction to ensure this. The problem, though not the
solution, is somewhat similar to “hook” errors Knill (2005).

In Fig. 2, we introduce a new circuit that does not allow phase errors to
propagate and produce logical errors. It retains the important feature of
redundant syndrome extraction. The crucial operational change is to use the
ancillas as the control qubits for CkNOT (Toffoli) gates rather than the data
qubits. With this circuit design, any errors that happen in the syndrome
extraction do not spread because the syndrome information for each stabilizer
is extracted in a single time step. This step does introduce errors in the
data qubits and joint errors in a data qubit and the appropriate ancilla. The
point, however, is that these errors are controlled by the redundancy -
multiqubit errors are still possible, but all of them lie within the
stabilizer group. Hence they do not destroy fault tolerance.

##  III Results

We test these ideas with a sequence of 3 simulations and calculations. With
the exception of the inclusion of the ancilla phase errors, the details of the
error model used in the simulations are the same as in Ref. Crow  _et al._
(2016). In particular, we take the gate error rate to be the same as the
memory or ”idle” error rate, both denoted by $p$. The results are presented in
Fig. 3. For each, the logical error probability $p_{log}$ is plotted against
the physical error probability $p$.

Figure 2: (color online) Revised circuit that avoids logical errors due to
back-action by extracting each stabilizer to an ancilla in a single time step.
The main change from the circuit of Fig.1 is the reversal of the roles of data
qubits and ancilla qubits in the extraction process. This circuit is fault
tolerant even in the presence of phase errors on the ancilla qubits.

1\. For comparison to earlier results, we simulate the circuit in Fig. 1
allowing only bit flip errors on the ancillas. This means treating them as
essentially classical, as was done in Ref. Crow  _et al._ (2016). We find a
threshold value of about $6\times 10^{-5}$ but we stress once more that this
is not a physically realistic model. This is the green curve, shown together
with a fit, in Fig. 3. The threshold occurs when the curve crosses the dashed
line that represents the equation $p=p_{log}$ .

2\. We simulate the circuit in Fig. 1 allowing bit flip and phase flip errors
and we calculate the logical error rate by counting the error sites between
the first and second pair of CNOT gates on each ancilla qubit. A fraction of
these errors involve a phase flip. Two of the three single qubit errors and
eight of the 15 two-qubit errors result in logical errors. This gives a
prediction for the logical error rate as a function of the physical error
rate. If the dangerous memory and gate error sites are $N_{m}$ and $N_{g}$ we
can estimate that

$p_{log}(p)=\left(N_{m}\times\frac{2}{3}+N_{g}\times\frac{8}{15}\right)p+\mathcal{O}(p^{2})\approx
11.73\;p.$

This form fits the simulations well. The data and the fit are the red curve in
Fig.3. Note that it does not intersect the dashed line, indicating that there
is no threshold: the circuit is not fault tolerant.

3\. We simulate the revised circuit of Fig.2 allowing for both flip and phase
errors on the ancillas. The results and a fit are shown as the blue curve in
Fig. 3. (It lies very close to the green curve.) There is a threshold, which
we compute to be about $5\times 10^{-5}$. This means that the revised circuit
is indeed fault tolerant with a realistic error model.

Figure 3: (color online) Data and fits for three simulations. The dashed line
represents $p_{log}=p$. The green curve and points represent the results for
the circuit in Fig. 1 without allowing for phase errors on ancilla qubits.
There is a threshold of about $6\times 10^{-5}$. The circuit is fault tolerant
but the error model is not realistic. The red curve and points represent the
results for the circuit in Fig. 1 allowing for phase errors on ancilla qubits.
There is no threshold, since the curve does not intersect the dashed line. The
circuit is not fault tolerant with this more realistic error model. The blue
curve and points represent the results for the circuit in Fig. 2 allowing for
phase errors on ancilla qubits. There is a threshold of about $5\times
10^{-5}$. The circuit is fault tolerant with a realistic error model.

This work shows that CEC can be performed with thresholds that are comparable
to those of MEC. It is important to use redundancy in the syndrome extraction.
This increases the threshold. We have assumed in this work that there is a
native CkNOT gate available. This is likely to be the case in, for example,
neutral atom set-ups, where such gates have been demonstrated Levine  _et al._
(2019). If only strictly 2-qubit gates are used, then the CkNOTs used in the
syndrome extraction can be replaced by CNOT gates at the cost of $k-1$
additional ancillas per extraction Isenhower  _et al._ (2011). This
modification would not destroy the fault-tolerant nature of the circuit.
However, it would involve a considerable number of additional qubits.

###### Acknowledgements.

We would like to thank D. Crow for useful discussions. This research was
sponsored in part by the Army Research Office (ARO) under Grant Number
W911NF-17-1-0274. The views and conclusions contained in this document are
those of the authors and should not be interpreted as representing the
official policies, either expressed or implied, of the Army Research Office
(ARO), or the U.S. Government. The U.S. Government is authorized to reproduce
and distribute reprints for Government purposes notwithstanding any copyright
notation herein. MS was supported by NSF PHY-1720220.

## References

  * Beterov and Saffman (2015) I. I. Beterov and M. Saffman, “Rydberg blockade, Förster resonances, and quantum state measurements with different atomic species,” Phys. Rev. A 92, 042710 (2015).
  * Petta  _et al._ (2005) J. R. Petta, A. C. Johnson, J. M. Taylor, E. A. Laird, A. Yacoby, M. D. Lukin, C. M. Marcus, M. P. Hanson, and A. C. Gossard, “Coherent manipulation of coupled electron spins in semiconductor quantum dots,” Science 309, 2180–2184 (2005).
  * Thorgrimsson  _et al._ (2017) B. Thorgrimsson, D. Kim, Y.-C. Yang, L. W. Smith, C. B. Simmons, D. R. Ward, R. H. Foote, J. Corrigan, D. E. Savage, M. G. Lagally, M. Friesen, S. N. Coppersmith, and M. A. Eriksson, “Extending the coherence of a quantum dot hybrid qubit,” npj Qu. Infor. 3, 32 (2017).
  * Paz-Silva  _et al._ (2010) G. A. Paz-Silva, G. K. Brennen, and J. Twamley, “Fault tolerance with noisy and slow measurements and preparation,” Phys. Rev. Lett. 105, 100501 (2010).
  * Li  _et al._ (2011) C.-K. Li, M. Nakahara, Y. Poon, N.-S. Sze, and H. Tomita, “Recovery in quantum error correction for general noise without measurement,” Qu. Inf. & Comp. 12, 0149 (2011).
  * Crow  _et al._ (2016) D. Crow, R. Joynt, and M. Saffman, “Improved error thresholds for measurement-free error correction,” Phys. Rev. Lett. 117, 130503 (2016).
  * Paz-Silva  _et al._ (2017) G. A. Paz-Silva, L. M. Norris, and L. Viola, “Multiqubit spectroscopy of gaussian quantum noise,” Phys. Rev. A 95, 022121 (2017).
  * Steane (1996) A. Steane, “Multiple-particle interference and quantum error correction,” Proc. R. Soc. Lond. A 452, 2551 (1996).
  * Preskill (1998) J. Preskill, “Fault-tolerant quantum computation,” in  _Introduction to Quantum Computation_ , edited by H.K. Lo, S. Popescu, and T.P. Spiller (World Scientific, Singapore, 1998) pp. 213–269.
  * Knill (2005) E. Knill, “Scalable quantum computing in the presence of large detected-error rates,” Phys. Rev. A 71, 042322 (2005).
  * Levine  _et al._ (2019) H. Levine, A. Keesling, G. Semeghini, A. Omran, T. T. Wang, S. Ebadi, H. Bernien, M. Greiner, V. Vuletić, H. Pichler, and M. D. Lukin, “Parallel implementation of high-fidelity multiqubit gates with neutral atoms,” Phys. Rev. Lett. 123, 170503 (2019).
  * Isenhower  _et al._ (2011) L. Isenhower, M. Saffman, and K. Mølmer, “Multibit CkNOT quantum gates via Rydberg blockade,” Quant. Inf. Proc. 10, 755 (2011).


#  Relativistic Dips in Entangling Power of Gravity

Marko Toroš  University College London, Gower Street, WC1E 6BT London, United
Kingdom  Martine Schut  Van Swinderen Institute, University of Groningen, 9747
AG, The Netherlands  Patrick Andriolo  Physics Institute, University of São
Paulo, Rua do Matão, 1371, São Paulo, Brazil  Sougato Bose  University College
London, Gower Street, WC1E 6BT London, United Kingdom  Anupam Mazumdar  Van
Swinderen Institute, University of Groningen, 9747 AG, The Netherlands

###### Abstract

The salient feature of both classical and quantum gravity is its universal and
attractive character. However, less is known about the behaviour and build-up
of quantum correlations when quantum systems interact via graviton exchange.
In this work, we show that quantum correlations can remain strongly suppressed
for certain choices of parameters even when considering two adjacent quantum
systems in delocalized states. Using the framework of linearized quantum
gravity with post-Newtonian contributions, we find that there are special
values of delocalization where gravitationally induced entanglement drops to
negligible values, albeit non-vanishing. We find a pronounced cancellation
point far from the Planck scale, where the system tends towards
classicalization. In addition, we show that quantum correlations begin to
reemerge for large and tiny delocalizations due to Heisenberg’s uncertainty
principle and the universal coupling of gravity to the energy-momentum tensor,
forming a valley of gravitational entanglement.

##  I Introduction

The coupling of classical gravity to the stress-energy tensor has been probed
in numerous experiments and has withstood the test of time in all astronomical
observations [1]. One of its most distinguishing consequences is the universal
and attractive character of the induced gravitational matter-matter
interaction. Such behaviour is manifest by looking at Newton’s 1/r potential,
but persists also when including post-Newtonian (PN) corrections depending on
the particle momenta [2, 3, 4, 5, 6].

The same universal and attractive behaviour is also a feature of an effective
field theory of quantum gravity [4], where the gravitational field is not
real-valued but rather operator-valued [7, 8, 9, 10]. Such a quantum
interaction can generate non-classical correlations between quantum systems
with no classical analogue, making it ideal for testing genuinely quantum
aspects. This observation was critical in conceiving a protocol to test the
quantum nature of gravity using two massive particles [11, 12] 111The results
of Ref. [11] were already known earlier, see [13].. This protocol, known as
the quantum gravity-induced entanglement of masses (QGEM), is in this regard
akin to Bell’s original idea of testing quantum correlations between two
spatially separated systems [14, 15].

According to the Local Operations and Classical Communication (LOCC) principle
[16], entanglement can only be generated by ostensibly quantum interactions
between the test particles. Hence, only if gravity is a quantum entity will it
generate an entangled state of the two masses [11, 17]. Within the context of
an effective field theory of quantum gravity the gravitational interaction is
being mediated by the massless spin-2 graviton, see [18, 19, 4, 17, 20, 21,
22, 23, 24, 7, 8, 9, 10], for a textbook, see [25].

Theoretical works and feasibility studies about the QGEM proposal have mainly
focused on the static regime where the momenta of the particles (i.e., the PN
corrections) are neglected in the gravitational interaction. In the static
Newtonian limit, the interaction is in position, and the spatial
delocalizations of the quantum states control entanglement generation.
Generally speaking, increasing the spatial delocalization $\Delta x$ will
increase the overall generated entanglement. Based on the intuition from the
static Newtonian limit, decreasing the spatial delocalization $\Delta x$ would
suggest that the generated entanglement is bound to decrease. We will show
that this is not the case and that this naive picture breaks down when the PN
corrections are included in the analysis [26].

In this paper, we will show that the generated entanglement in general
increases for both very small and very large spatial delocalizations $\Delta
x$. This observation is, in a way, simple, and it follows directly from
Heisenberg’s uncertainty principle [27] and the universal coupling of gravity
to energy [28]. The generated entanglement entropy, quantifying the degree of
entanglement, scales as a function of $\Delta x$ for large spatial
delocalizations and as a function of $\Delta p\propto\hbar/\Delta x$ for small
spatial delocalizations. We will illustrate this result using a toy setup of
two harmonic oscillators interacting gravitationally for the initial state of
the product of ground states.

There is also an unexpected twist in the story. We find small pockets in the
parameter space of $\Delta x$ – far from the Planck scale – where the
entanglement is strongly suppressed. This decrease is pronounced when the
dominant 0PN and 1PN terms entering in the generation of quantum correlations
cancel each other’s contributions. The explanation lies in the opposite sign
of the two-mode squeezing parameter induced by the 0PN and 1PN couplings. We
show that this behaviour persists, and can be controlled, by varying the
degree of squeezing of the initial states. We conclude by arguing that
analogous cancellations of quantum correlations, where the system tends
towards classicalization, should appear also in any other quantum mechanical
theory of gravity.

##  II Quantum gravitational potential and harmonic oscillators

We consider the simple 1D toy model of two identical harmonic oscillators, A
and B, characterized by the mass $m$ and angular frequency
$\omega_{\text{m}}$, oscillating along the x-axis. We will assume that the
centres of the two harmonic traps are separated by a distance $d$, and the two
particles interact only gravitationally. Using Gupta’s framework of linearized
quantum gravity [18, 19] and perturbation theory, we can then obtain the
gravitational matter-matter potential up to order 2PN [7] (see also [8, 10,
9]):

$\displaystyle\hat{H}_{\text{grav}}=$
$\displaystyle-\frac{Gm^{2}}{|\hat{r}_{A}-\hat{r}_{B}|}-\frac{G(3\hat{p}_{A}^{2}-8\hat{p}_{A}\hat{p}_{B}+3\hat{p}_{B}^{2})}{2c^{2}|\hat{r}_{A}-\hat{r}_{B}|}$
$\displaystyle+\frac{G(5\hat{p}_{A}^{4}-18\hat{p}_{A}^{2}\hat{p}_{B}^{2}+5p_{B}^{4})}{8c^{4}m^{2}|\hat{r}_{A}-\hat{r}_{B}|}\,,$
(1)

where $c$ ($G$) denotes the speed of light (the gravitational constant). Here,
we will implicitly assume that the momenta are sufficiently small such that
higher order terms, i.e., the terms
$\propto\hat{p}_{j}^{n}\hat{p}_{k}^{n^{\prime}}$ (with $j,k=A,B$ and
$n+n^{\prime}>4$), can be neglected. Higher order terms would only modify the
quantitative results for relativistic momenta without affecting the features
in the regime where the velocities are small compared to the speed of light
(i.e., we neglect terms beyond $\mathcal{O}(c^{-4})$). To keep the expressions
short, we also implicitly assume the convention that unsymmetrized expressions
(e.g., $\hat{x}\hat{p}$) are to be interpreted in the symmetrized ordering
(e.g., $(\hat{x}\hat{p}+\hat{p}\hat{x})/2$).

Post-Newtonian corrections have, of course, been analyzed extensively in
previous works in the center of momentum frame [2, 3, 5, 6, 29]. If we set
$\hat{p}\equiv\hat{p}_{A}=-\hat{p}_{B}$, and denote
$\hat{r}\equiv|{\hat{r}}_{A}-{\hat{r}}_{B}|$, we recover from Eq. (1) the
known result in the literature:

$\hat{H}_{\text{grav}}=-\frac{Gm^{2}}{\hat{r}}-7\frac{G\hat{p}^{2}}{c^{2}\hat{r}}-\frac{G\hat{p}^{4}}{c^{4}m^{2}\hat{r}}.$
(2)

In other words, Eq. (1) can be seen as the adaptation of Eq. (2) to the
specific case of two quantum harmonic oscillators. Generalizations to more
particles, e.g., four harmonic oscillators, forming pairs of particles and
detectors [30], modified gravity scenarios such as in the case of a massive
graviton [10], fat graviton with nonlocal interaction [9, 31], or a dilaton-
graviton combination [32], can also be analyzed using similar methods.

We now suppose that the two trap centres are located at $\pm d/2$ and write
$\hat{r}_{A}=-d/2+{\hat{x}}_{A}$, and $\hat{r}_{B}=d/2+{\hat{x}}_{B}$. The
operators ${\hat{x}}_{A}$ and $\hat{x}_{B}$ denote small displacements from
the equilibrium position, while the corresponding conjugate momenta are given
by $\hat{p}_{A}$ and $\hat{p}_{B}$, respectively. The Hamiltonian of the two
harmonic oscillators is given by

$\hat{H}_{\text{matter}}=\frac{\hat{p}_{A}^{2}}{2m}+\frac{\hat{p}_{B}^{2}}{2m}+\frac{m\omega_{\text{m}}^{2}}{2}\hat{x}_{A}^{2}+\frac{m\omega_{\text{m}}^{2}}{2}\hat{x}_{B}^{2},$
(3)

where $\omega_{\text{m}}$ and $m$ denote the harmonic frequency and mass,
respectively (assumed for simplicity to be the same for the two harmonic
oscillators). For later convenience, we introduce the mode decompositions

$\displaystyle\hat{x}_{\text{A}}$ $\displaystyle=\delta
x(\hat{a}+\hat{a}^{\dagger}),\,\,\,\,$ $\displaystyle\hat{x}_{\text{B}}$
$\displaystyle=\delta x(\hat{b}+\hat{b}^{\dagger}),$ (4)
$\displaystyle\hat{p}_{\text{A}}$ $\displaystyle=i\delta
p(\hat{a}^{\dagger}-\hat{a}),\,\,\,\,$ $\displaystyle\hat{p}_{\text{B}}$
$\displaystyle=i\delta p(\hat{b}^{\dagger}-\hat{b}),$ (5)

where $\hat{a},\hat{b}$ ($\hat{a}^{\dagger},\hat{b}^{\dagger}$) denote the
annihilation (creation) operators, and $\delta
x=\sqrt{\frac{\hbar}{2m\omega_{\text{m}}}}$, $\delta p=\sqrt{\frac{\hbar
m\omega_{\text{m}}}{2}}$ are the position, momentum zero-point-motions,
respectively.

Here, we will be interested in the phenomenology of the gravitational
potential in Eq. (1) up to the quartic order in the operators
$\propto\hat{O}_{i}\hat{O}_{j}\hat{O}_{k}\hat{O}_{l}$ (with
$\hat{O}_{i,j,k,l}=\hat{x}_{A},\hat{x}_{B},\hat{p}_{A},\hat{p}_{B}$)
corresponding to small position and momentum fluctuations. Taylor expanding
Eq. (1) around the equilibrium positions and considering small fluctuations,
we find the following Hamiltonian:

$\hat{H}_{\text{grav}}=-\frac{Gm^{2}}{d}+\hat{H}_{\text{A}}+\hat{H}_{\text{B}}+\hat{H}_{\text{AB}},$
(6)

where the first term only produces a global phase, and $\hat{H}_{A}$ and
$\hat{H}_{B}$ depend only on the operators of particle A and B, respectively.
The leading order cross-coupling terms between the two particles are given by

$\hat{H}_{\text{AB}}=\frac{Gm^{2}}{d}\left[\frac{2\hat{x}_{\text{A}}\hat{x}_{\text{B}}}{d^{2}}+\frac{4\hat{p}_{\text{A}}\hat{p}_{\text{B}}}{m^{2}c^{2}}-\frac{9\hat{p}_{A}^{2}\hat{p}_{B}^{2}}{4m^{4}c^{4}}\right],$
(7)

where the 0PN, 1PN, and 2PN contributions appear from left to right,
respectively.

While the leading order gravitational force arises from the terms
$\hat{H}_{\text{A}}$ and $\hat{H}_{\text{B}}$ (as these terms contain the
uniform gravitational fields affecting the motion of the individual particle),
the leading order contribution for entanglement generation arises from the
cross-couplings in $\hat{H}_{\text{AB}}$ [33]. This is already hinting that
the generation of gravitationally induced entanglement might be hiding some
surprises.

##  III Gravitational entanglement entropy

In this section we illustrate how the gravitational entanglement entropy
depends on the delocalization using simple first order perturbation theory. We
assume that the initial state is a product of the two ground states

$|\psi_{\text{i}}\rangle=|0\rangle_{\text{A}}|0\rangle_{\text{B}},$ (8)

which is perturbed by the gravitational interaction $\hat{H}_{\text{grav}}$ in
Eq. (6). We can always decompose the perturbed state vector in the number
basis:

$|\psi_{AB}\rangle=\frac{1}{\sqrt{\mathcal{N}}}\sum_{n,N}C_{nN}|n\rangle|N\rangle,$
(9)

where $C_{nN}$ denote the coefficients (with $C_{00}=1$),
$\mathcal{N}=\sum_{n,N}|C_{nN}|^{2}$ is the normalization, and $|n\rangle$,
$|N\rangle$ denote the number states of the two harmonic oscillators. In
particular, the coefficients appearing in Eq. (9) can be computed using

$C_{nN}=\frac{\langle n|\langle
N|\hat{H}_{AB}|0\rangle|0\rangle}{2E_{0}-E_{n}-E_{N}},$ (10)

where $E_{0}$ and $E_{n}$, $E_{N}$ denote the energy of the ground state
$|0\rangle$ and of the excited states $|n\rangle$, $|N\rangle$, respectively
(we recall that for a harmonic oscillator, we have $E_{j}=E_{0}+\hbar\omega j$
with $j$ denoting the occupation number).

Figure 1: Entanglement entropy $S$ as a function of the spatial delocalization
$\Delta x$ (bottom) or the momentum delocalization $\Delta p=\hbar/\Delta x$
(top). We have expressed the spatial (momentum) superposition in units of the
Planck length (momentum) given by $x_{\text{P}}\equiv\sqrt{\hbar G/c^{3}}$
($p_{\text{P}}\equiv\hbar/x_{\text{P}}$). The behaviour on the right is
determined by the static contribution
$\propto\hat{x}_{\text{A}}\hat{x}_{\text{B}}$. The plateau in the middle is
determined by the non-static 1PN contribution
$\propto\hat{p}_{\text{A}}\hat{p}_{\text{B}}$. The dominant static and the
non-static couplings cancel each other’s contribution to the entanglement
entropy at $\sqrt{2}d\omega_{\text{m}}\sim c$, and we observe that the
entanglement entropy $S$ drops to negligible values (the three pronounced dips
indicated by the arrow). We note that these dips of quantum correlations
indicate a _classicalization_ of the system and that they occur far from the
Planck length $x_{p}$. Finally, as we approach relativistic velocities
(corresponding to large momentum superposition sizes on the left side of the
figure), we find that higher-order momentum contributions become important.
Here we have plotted the contribution from the term
$\propto\hat{p}_{\text{A}}^{2}\hat{p}_{\text{B}}^{2}$ arising at order
$\mathcal{O}(1/c^{4})$. We plot the curves for different values of the trap
distance $d$. Note that the qualitative behaviour remains the same in all
three cases, with the location of the dip shifting to the left (right) for
smaller (larger) values of $d$ as expected (i.e., larger distances weaken the
position couplings, hence requiring larger spatial superposition sizes). The
entanglement entropy is expressed in units of the entanglement entropy when we
set the spatial superposition to that of the Planck length, i.e.,
$S_{p}=S(\Delta
x=x_{\text{P}})=-(\frac{Gm}{c^{2}d})^{2}\log((\frac{Gm}{c^{2}d})^{2})$ which
is achieved on the plateau.

We will quantify the degree of entanglement using the von Neumann entanglement
entropy given by $S=-\text{tr}[\rho_{A}\text{ln}\rho_{A}]$, where
$\rho_{A}=\text{tr}_{B}[\rho_{AB}]$ is the reduced density matrix of subsystem
$A$, and $\rho_{AB}=|\psi_{AB}\rangle\langle\psi_{AB}|$ is the total density
matrix of the system (see Appendix A for a short review on the entanglement
entropy). Using Eqs. (7)-(10) find a simple formula for the steady-state
entanglement entropy:

$\displaystyle S\approx$
$\displaystyle-\left(\frac{Gm}{c^{2}d}-\frac{Gm}{2d^{3}\text{$\omega_{\text{m}}^{2}$}}\right)^{2}\log\left(\left(\frac{Gm}{c^{2}d}-\frac{Gm}{2d^{3}\omega_{\text{m}}^{2}}\right)^{2}\right)$
$\displaystyle-\frac{81G^{2}\omega_{\text{m}}^{2}\hbar^{2}}{1024c^{8}d^{2}}\log\left(\frac{81G^{2}\omega_{\text{m}}^{2}\hbar^{2}}{1024c^{8}d^{2}}\right),$
(11)

which is plotted in Fig. 1 (see Appendix B for the detailed derivation). The
first line of Eq. (11) captures the right-most part of the plot (from the 0PN
$\hat{x}_{A}\hat{x}_{B}$ coupling) as well as the plateau (from the 1PN
$\hat{p}_{A}\hat{p}_{B}$ coupling). The right-most part of the plot in Fig. 1
is captured by the 2PN contribution in the second line of Eq. (11) arising
from the coupling $\propto\hat{p}_{A}^{2}\hat{p}_{B}^{2}$.

We find that the entanglement entropy grows for very large and very small
positions delocalizations $\Delta x$ given by the zero-point-motion. For large
$\Delta x$ the entanglement entropy grows as a consequence of the 0PN static
position couplings arising from the familiar 1/r potential, while for small
$\Delta x$ the momentum delocalization $\Delta p\sim\hbar/\Delta x$ becomes
large, and the post-Newtonian momentum couplings start increasing the
entanglement entropy. In other words, the landscape of entanglement entropy as
a function of delocalizations $\Delta x$ or $\Delta p$ forms a valley. This
can be seen as a consequence of Heisenberg’s minimum uncertainty relation
$\Delta p\Delta x=\hbar/2$ and the universal coupling of gravity to all forms
of energy.

In addition, Fig. 1 reveals unexpected dips in entanglement entropy, which
appear far the Planck length. These dips arise because of the cancellation
between the 0PN and 1PN terms in the first line of Eq. (11) corresponding to
the couplings $\propto\hat{x}_{A}\hat{x}_{B}$ and
$\propto\hat{p}_{A}\hat{p}_{B}$, respectively. The first line vanishes when
the product of the harmonic frequency, $\omega_{\text{m}}$, and of the
distance between the two traps, $d$, becomes comparable to the speed of light,
i.e., $\sqrt{2}\,\omega_{\text{m}}d=c$. In this case, the position and
momentum delocalizations are given by $\Delta x=\sqrt{\frac{\hbar
d}{mc}}/\sqrt[4]{2}$ and $\Delta p=\sqrt{\frac{\hbar mc}{d}}/\sqrt[4]{2}$,
respectively. At this point, the two adimensional parameters governing the 0PN
and 1PN match, i.e., $\Delta x/d=\sqrt{2}\Delta p/(mc)$.

The adimensional parameters however only explain that the 0PN and 1PN
contributions are of equal magnitude without revealing the origin of the
cancellation. As we will see in the next section, this cancellation arises
from the opposite sign of the two-mode squeezing (TMS) character of the
position and momentum couplings.

##  IV Dips of gravitational entanglement

Figure 2: The maximum entanglement value achieved during the time evolution
starting with a product of two single-mode-squeezed-vacuum (SMSV) states as a
function of the spatial delocalization $\Delta x=\delta xe^{-r}$ (bottom axis)
and momentum delocalization $\Delta p=\delta pe^{r}$ (top axis). By setting
$\omega_{\text{m}}=\omega_{0}\equiv c/(\sqrt{2}d)$ we find the case of equal
couplings $g_{x}=g_{p}=g_{0}\equiv\sqrt{2}Gm/(c^{2}d)$ (green line). We also
consider the case $\omega_{\text{m}}=4\omega_{0}$ producing the couplings
$4g_{x}=g_{p}/4=g_{0}$ (orange line), and the case
$\omega_{\text{m}}=\omega_{0}/4$ producing the coupling $g_{x}/4=4g_{p}=g_{0}$
(black line). The entanglement entropy is normalized to the maximum value of
the green line to ease the comparison. We note that both the horizontal
location and the depth can be changed by tuning the couplings $g_{x}$, $g_{p}$
using the mechanical frequency $\omega_{\text{m}}$. At the location of the
dips the system tends towards classicalization as there is a strong
suppression of quantum correlations.

In this section we look closer at the unexpected dips of gravitational
entanglement uncovered in the previous section. We restrict the analysis only
to the leading order terms identified in the previous section as its origin,
i.e.,

$\hat{H}_{\text{AB}}=\frac{2Gm^{2}}{d^{3}}\hat{x}_{\text{A}}\hat{x}_{\text{B}}+\frac{4G}{c^{2}d}\hat{p}_{\text{A}}\hat{p}_{\text{B}}.$
(12)

Using Eqs. (4) and (5) in Eq. (12) we then find

$\hat{H}_{\text{AB}}=\hbar
g_{-}(\underbrace{ab+a^{\dagger}b^{\dagger}}_{\text{TMS}})+\hbar
g_{+}(\underbrace{ab^{\dagger}+a^{\dagger}b}_{\text{BS}}),$ (13)

where the coupling rates are

$\displaystyle g_{-}=g_{x}-g_{p}$ $\displaystyle,\qquad g_{+}=g_{x}+g_{p},$
(14) $\displaystyle g_{x}=\frac{Gm}{d^{3}\omega_{\text{m}}}$
$\displaystyle,\qquad g_{p}=\frac{2Gm\omega_{\text{m}}}{c^{2}d},$ (15)

and TMS (BS) labels the two-mode squeezing (beam-splitter) contribution (for
an introduction on quantum optics transformations see for example [34, 35]).

For an initial product of ground states $|0\rangle|0\rangle$ the BS part of
the Hamiltonian has no effect, while the TMS-induced change depends on the PN
order: at 0PN we have the coupling $+g_{x}$, while at 1PN we have the coupling
$-g_{p}$. The explanation for the dips in Fig. 1 thus lies in the opposite
sign of the TMS transformation arising from the position and momentum
couplings of gravity. When $g_{x}=g_{p}$ the 0PN and 1PN contributions to TMS
cancel, i.e., $g_{-}=0$, and we find a dip.

The observed cancellation point persists even if we consider the initial state
to be the product of two single-mode-squeezed-vacuum (SMSV) states which have
enhanced position or momentum delocalizations given by

$\Delta x=\delta xe^{-r},\qquad\Delta p=\delta pe^{r},$ (16)

respectively (with $\delta x$ and $\delta p$ denoting the zero-point-motions
defined below Eq. (5)). Such states are not energy eigenstates, and hence, we
have to take into account their time-evolution in the harmonic traps (see
derivation in Appendix C). In Fig. 2 we show how the maximum achieved
entanglement entropy changes by varying the degree of squeezing $r$ of the
initial state. We note that by tuning the mechanical frequency
$\omega_{\text{m}}$ (and hence the couplings $g_{x},g_{p}$ in Eq. (15)) we can
also shift to profile horizontally and increase its depth. The location of the
dip is now determined by the condition $g_{x}e^{-2r}=g_{p}e^{2r}$, which
depends on the squeezing parameter $r$ as expected.

##  V Discussion

In this work, we have discussed entanglement generation within the context of
a linearized quantum gravity with post-Newtonian (PN) momentum contributions.
We have considered a simple toy model of two nearby harmonic oscillators and
computed the entanglement entropy arising from their gravitational
interaction. We uncovered a surprising interplay between the Heisenberg’s
uncertainty relation, the squeezing character of gravity, and the generation
of non-classical correlations when 1PN and 2PN contributions are taken into
account.

We found that the entanglement entropy drops to negligible values when the 0PN
contribution and the 1PN are of the same magnitude, i.e., when $\Delta
x/d\sim\Delta p/(mc)$, where $m$ is the mass of each harmonic oscillator, and
$d$ denotes the distance between the two trap centres. Although the classical
gravitational force is universally attractive, the generation of
gravitationally induced entanglement can be suppressed to negligible values
for specific states. As this will remain true within any theory that recovers
the predictions of perturbative quantum gravity, such states could. provide a
method to distinguish perturbative quantum gravity from other classes of
gravitational theories, such as a scalar-tensor theory with different PN
contributions. While this work does not answer questions of experimental
feasibility, it nonetheless uncovers experimentally defining features of
infrared quantum gravity. The entanglement dips provide a distinct signature
of classicalization, showing how quantum correlations can become suppressed in
a fully quantum mechanical framework.

###### Acknowledgements.

MT would like to acknowledge funding from ST/W006227/1. MS is supported by the
Fundamentals of the Universe research program at the University of Groningen.
PA acknowledges funding from Brazilian Coordination for the Improvement of
Higher Education Personnel (CAPES), Erasmus+ Programme and the Strategic
Partnership Framework at University of Groningen. SB thanks EPSRC grants
EP/R029075/1, EP/X009467/1, and ST/W006227/1. SB and AM’s research is
supported by the Sloan, and Betty and Moore foundations.

## References

  * [1] Clifford M Will.  The confrontation between general relativity and experiment.  Living reviews in relativity, 17(1):1–117, 2014.
  * [2] Yoichi Iwasaki.  Quantum theory of gravitation vs. classical theory: fourth-order potential.  Progress of Theoretical Physics, 46(5):1587–1609, 1971.
  * [3] Kichiro Hiida and Hiroshi Okamura.  Gauge transformation and gravitational potentials.  Progress of Theoretical Physics, 47(5):1743–1757, 1972.
  * [4] John F. Donoghue.  General relativity as an effective field theory: The leading quantum corrections.  Phys. Rev. D, 50:3874–3888, 1994.
  * [5] Andrea Cristofoli, NEJ Bjerrum-Bohr, Poul H Damgaard, and Pierre Vanhove.  Post-minkowskian hamiltonians in general relativity.  Physical Review D, 100(8):084040, 2019.
  * [6] Gianluca Grignani, Troels Harmark, Marta Orselli, and Andrea Placidi.  Fixing the non-relativistic expansion of the 1pm potential.  Journal of High Energy Physics, 2020(12):1–18, 2020.
  * [7] Sougato Bose, Anupam Mazumdar, Martine Schut, and Marko Toroš.  Mechanism for the quantum natured gravitons to entangle masses.  Phys. Rev. D, 105:106028, 2022.
  * [8] Dripto Biswas, Sougato Bose, Anupam Mazumdar, and Marko Toroš.  Gravitational optomechanics: Photon-matter entanglement via graviton exchange.  Phys. Rev. D, 108(6):064023, 2023.
  * [9] Ulrich K. Beckering Vinckers, Álvaro de la Cruz-Dombriz, and Anupam Mazumdar.  Quantum entanglement of masses with nonlocal gravitational interaction.  Phys. Rev. D, 107(12):124036, 2023.
  * [10] Shafaq Gulzar Elahi and Anupam Mazumdar.  Probing massless and massive gravitons via entanglement in a warped extra dimension.  Phys. Rev. D, 108(3):035018, 2023.
  * [11] Sougato Bose, Anupam Mazumdar, Gavin W. Morley, Hendrik Ulbricht, Marko Toroš, Mauro Paternostro, Andrew Geraci, Peter Barker, M. S. Kim, and Gerard Milburn.  Spin Entanglement Witness for Quantum Gravity.  Phys. Rev. Lett., 119:240401, 2017.
  * [12] Chiara Marletto and Vlatko Vedral.  Gravitationally induced entanglement between two massive particles is sufficient evidence of quantum effects in gravity.  Phys. Rev. Lett., 119:240402, 2017.
  * [13] Sougato Bose, 2016.  https://www.youtube.com/watch?v=0Fv-0k13s_k, Accessed 1/11/22, on behalf of QGEM collaboration.
  * [14] John S Bell.  On the einstein podolsky rosen paradox.  Physics Physique Fizika, 1(3):195, 1964.
  * [15] Nicolas Brunner, Daniel Cavalcanti, Stefano Pironio, Valerio Scarani, and Stephanie Wehner.  Bell nonlocality.  Reviews of modern physics, 86(2):419, 2014.
  * [16] Charles H. Bennett et al.  Mixed-state entanglement and quantum error correction.  Phys. Rev. A, 54:3824, 1996.
  * [17] Ryan J. Marshman, Anupam Mazumdar, and Sougato Bose.  Locality and entanglement in table-top testing of the quantum nature of linearized gravity.  Phys. Rev. A, 101:052110, 2020.
  * [18] Suraj N Gupta.  Quantization of einstein’s gravitational field: linear approximation.  Proceedings of the Physical Society. Section A, 65(3):161, 1952\.
  * [19] Suraj N Gupta.  Quantization of einstein’s gravitational field: General treatment.  Proceedings of the Physical Society. Section A, 65(8):608, 1952\.
  * [20] Daniel Carney, Philip C E Stamp, and Jacob M Taylor.  Tabletop experiments for quantum gravity: a user’s manual.  Class. Quant. Grav., 36:034001, 2019.
  * [21] Marios Christodoulou, Andrea Di Biagio, Markus Aspelmeyer, Časlav Brukner, Carlo Rovelli, and Richard Howl.  Locally mediated entanglement in linearized quantum gravity.  Phys. Rev. Lett., 130:100202, Mar 2023.
  * [22] Marios Christodoulou and Carlo Rovelli.  On the possibility of laboratory evidence for quantum superposition of geometries.  Physics Letters B, 792:64–68, 2019.
  * [23] Alessio Belenchia et al.  Quantum Superposition of Massive Objects and the Quantization of Gravity.  Phys. Rev. D, 98:126009, 2018.
  * [24] Daine L. Danielson, Gautam Satishchandran, and Robert M. Wald.  Gravitationally mediated entanglement: Newtonian field versus gravitons.  Phys. Rev. D, 105:086001, 2022.
  * [25] M. D. Scadron.  Advanced quantum theory.  World Scientific Publishing Company, 2007.
  * [26] Eric Poisson and Clifford M Will.  Gravity: Newtonian, post-newtonian, relativistic.  Cambridge University Press, 2014.
  * [27] Paul Busch, Teiko Heinonen, and Pekka Lahti.  Heisenberg’s uncertainty principle.  Physics reports, 452(6):155–176, 2007.
  * [28] Charles W Misner, Kip S Thorne, and John Archibald Wheeler.  Gravitation.  Macmillan, 1973.
  * [29] Luc Blanchet.  Gravitational radiation from post-newtonian sources and inspiralling compact binaries.  Living Reviews in Relativity, 17(1), February 2014.
  * [30] Fabian Gunnink, Anupam Mazumdar, Martine Schut, and Marko Toroš.  Gravitational decoherence by the apparatus in the quantum-gravity-induced entanglement of masses.  Classical and Quantum Gravity, 40(23):235006, 2023.
  * [31] Ulrich K Beckering Vinckers, Álvaro de la Cruz-Dombriz, and Anupam Mazumdar.  Smearing out contact terms in ghost-free infinite derivative quantum gravity.  arXiv preprint arXiv:2402.18694, 2024.
  * [32] Sumanta Chakraborty, Anupam Mazumdar, and Ritapriya Pradhan.  Distinguishing Jordan and Einstein frames in gravity through entanglement.  Phys. Rev. D, 108(12):L121505, 2023.
  * [33] Vijay Balasubramanian, Michael B. McDermott, and Mark Van Raamsdonk.  Momentum-space entanglement and renormalization in quantum field theory.  Phys. Rev. D, 86:045014, 2012.
  * [34] Christopher C Gerry and Peter L Knight.  Introductory quantum optics.  Cambridge university press, 2023.
  * [35] Ulf Leonhardt.  Essential quantum optics: from quantum measurements to black holes.  Cambridge University Press, 2010.
  * [36] Michael A Nielsen and Isaac L Chuang.  Quantum computation and quantum information.  Cambridge university press, 2010.
  * [37] Koenraad Audenaert, Jens Eisert, Martin B Plenio, and Reinhard F Werner.  Entanglement properties of the harmonic chain.  Physical Review A, 66(4):042327, 2002.
  * [38] Gerardo Adesso, Alessio Serafini, and Fabrizio Illuminati.  Extremal entanglement and mixedness in continuous variable systems.  Phys. Rev. A, 70:022318, Aug 2004.
  * [39] Alessio Serafini, Gerardo Adesso, and Fabrizio Illuminati.  Unitarily localizable entanglement of gaussian states.  Phys. Rev. A, 71:032349, Mar 2005.
  * [40] Christian Weedbrook, Stefano Pirandola, Raúl García-Patrón, Nicolas J. Cerf, Timothy C. Ralph, Jeffrey H. Shapiro, and Seth Lloyd.  Gaussian quantum information.  Rev. Mod. Phys., 84:621–669, May 2012.
  * [41] Samuel L. Braunstein and Peter van Loock.  Quantum information with continuous variables.  Rev. Mod. Phys., 77:513–577, Jun 2005.
  * [42] N Canosa, Swapan Mandal, and R Rossignoli.  Exact dynamics and squeezing in two harmonic modes coupled through angular momentum.  Journal of Physics B: Atomic, Molecular and Optical Physics, 48(16):165501, 2015.

##  Appendix A Brief overview of entanglement entropy formulae

The von Neumann entanglement entropy is given by

$S=-\text{tr}[\rho_{A}\text{ln}\rho_{A}],$ (17)

where $\rho_{A}=\text{tr}_{B}[\rho_{AB}]$ is the reduced density matrix of
subsystem $A$, $\rho_{AB}=|\psi_{AB}\rangle\langle\psi_{AB}|$ is the total
density matrix of the system. The same form of the expression would be
obtained by tracing over system A by formally exchanging the subsystems A and
B, i.e., $S=-\text{tr}[\rho_{B}\text{ln}\rho_{B}]$ and
$\rho_{B}=\text{tr}_{A}[\rho_{AB}]$. The general expression in Eq. (17) can be
however rewritten in a more convenient way for the two cases we analyse in
Secs. III and IV.

In Sec. III, we will suppose that the initial state is the product of the
ground states of two harmonic oscillators, i.e., $|0\rangle|0\rangle$ and use
time-independent perturbation theory. Using the Schmidt decomposition of the
state in Eq. (9) it can be shown that the entanglement entropy from Eq. (17)
reduces to

$S=-\sum_{j}|\alpha_{j}|^{2}\text{log}\left[|\alpha_{j}|^{2}\right],$ (18)

where $\alpha_{j}$ are the Schmidt coefficients. The Schmidt decomposition
can, in general, be performed using the linear algebra technique of singular
value decomposition (SVD) [36]. However, we will be primarily interested in
the states of the form $\sum_{j}\alpha_{j}|j\rangle|j\rangle$ (with
$|j\rangle$ denoting the number states), where one can readily read-off the
Schmidt coefficients $\alpha_{j}$.

In Sec. IV, we will be interested in Gaussian states and their time-evolution.
The formula for the entanglement entropy simplifies to [37, 38, 39, 40, 41]:

$S=-f\text{ln}f+(1+f)(1+\text{ln}f),$ (19)

where $f(t)$ is the symplectic eigenvalue of the single-mode covariance matrix
of subsystem A given by

$f(t)=\frac{1}{\hbar}\sqrt{4\langle\hat{x}_{A}^{2}\rangle\langle\hat{p}_{A}^{2}\rangle-\langle\hat{x}_{A}\hat{p}_{A}+\hat{p}_{A}\hat{x}_{A}\rangle^{2}}-\frac{1}{2}.$
(20)

Also, the formula in Eq. (19) remains valid if we formally replace the
quantities related to subsystem A with the ones for subsystem B in Eq. (20)
(i.e., we compute the symplectic eigenvalue of the single mode covariance
matrix of subsystem B).

##  Appendix B List of expansion coefficients up to quartic order in the
operators and up to 2PN

order | coefficient | coupling | Non-zero coefficients $C_{nN}$ from first-order perturbation theory
---|---|---|---
$C_{11}$ | $C_{21}$ | $C_{12}$ | $C_{31}$ | $C_{13}$ | $C_{22}$
0PN | $\frac{2Gm^{2}}{d^{3}}$ | $\mathbf{\hat{x}_{\text{A}}\hat{x}_{\text{B}}}$ | $\boldsymbol{-\frac{Gm}{2d^{3}\omega_{\text{m}}^{2}}}$ |  |  |  |  |
$\frac{3Gm^{2}}{d^{4}}$ | $\hat{x}_{\text{A}}^{2}\hat{x}_{\text{B}}$ |  | $-\frac{Gm\delta x}{\sqrt{2}d^{4}\omega_{\text{m}}^{2}}$ |  |  |  |
$-\frac{3Gm^{2}}{d^{4}}$ | $\hat{x}_{\text{A}}\hat{x}_{\text{B}}^{2}$ |  |  | $\frac{Gm\delta x}{\sqrt{2}d^{4}\omega_{\text{m}}^{2}}$ |  |  |
$\frac{4Gm^{2}}{d^{5}}$ | $\hat{x}_{\text{A}}^{3}\hat{x}_{\text{B}}$ | $-\frac{3G\hbar}{2d^{5}\omega_{\text{m}}^{3}}$ |  |  | $-\frac{\sqrt{6}G\hbar}{4d^{5}\omega_{\text{m}}^{3}}$ |  |
$-\frac{6Gm^{2}}{d^{5}}$ | $\hat{x}_{\text{A}}^{2}\hat{x}_{\text{B}}^{2}$ |  |  |  |  |  | $\frac{3G\hbar}{4d^{5}\omega_{\text{m}}^{3}}$
$\frac{4Gm^{2}}{d^{5}}$ | $\hat{x}_{\text{A}}\hat{x}_{\text{B}}^{3}$ | $-\frac{3G\hbar}{2d^{5}\omega_{\text{m}}^{3}}$ |  |  |  | $-\frac{\sqrt{6}G\hbar}{4d^{5}\omega_{\text{m}}^{3}}$ |
1PN | $\frac{4G}{c^{2}d}$ | $\mathbf{\hat{p}_{\text{A}}\hat{p}_{\text{B}}}$ | $\boldsymbol{\frac{Gm}{c^{2}d}}$ |  |  |  |  |
$\frac{4G}{c^{2}d^{2}}$ | $\hat{p}_{\text{A}}\hat{p}_{\text{B}}\hat{x}_{\text{A}}$ |  | $\frac{2\sqrt{2}Gm\delta x}{3c^{2}d^{2}}$ |  |  |  |
$-\frac{4G}{c^{2}d^{2}}$ | $\hat{p}_{\text{A}}\hat{p}_{\text{B}}\hat{x}_{\text{B}}$ |  |  | $-\frac{2\sqrt{2}Gm\delta x}{3c^{2}d^{2}}$ |  |  |
$\frac{4G}{c^{2}d^{3}}$ | $\hat{p}_{\text{A}}\hat{p}_{\text{B}}\hat{x}_{\text{A}}^{2}$ | $\frac{G\hbar}{2c^{2}d^{3}\omega_{\text{m}}}$ |  |  | $\frac{\sqrt{6}G\hbar}{4c^{2}d^{3}\omega_{\text{m}}}$ |  |
$-\frac{8G}{c^{2}d^{3}}$ | $\hat{p}_{\text{A}}\hat{p}_{\text{B}}\hat{x}_{\text{A}}\hat{x}_{B}$ |  |  |  |  |  | $-\frac{G\hbar}{c^{2}d^{3}\omega_{\text{m}}}$
$\frac{4G}{c^{2}d^{3}}$ | $\hat{p}_{\text{A}}\hat{p}_{\text{B}}\hat{x}_{B}^{2}$ | $\frac{G\hbar}{2c^{2}d^{3}\omega_{\text{m}}}$ |  |  |  | $\frac{\sqrt{6}G\hbar}{4c^{2}d^{3}\omega_{\text{m}}}$ |
$\frac{3G}{2c^{2}d^{2}}$ | $\hat{p}_{\text{A}}^{2}\hat{x}_{\text{B}}$ |  | $\frac{\sqrt{2}Gm\delta x}{4c^{2}d^{2}}$ |  |  |  |
$-\frac{3G}{2c^{2}d^{3}}$ | $\hat{p}_{\text{A}}^{2}\hat{x}_{\text{B}}^{2}$ |  |  |  |  |  | $-\frac{3G\hbar}{16c^{2}d^{3}\omega_{\text{m}}}$
$\frac{3G}{c^{2}d^{3}}$ | $\hat{p}_{\text{A}}^{2}\hat{x}_{\text{A}}\hat{x}_{\text{B}}$ | $-\frac{3G\hbar}{8c^{2}d^{3}\omega_{\text{m}}}$ |  |  | $\frac{3\sqrt{6}G\hbar}{16c^{2}d^{2}\omega_{\text{m}}}$ |  |
$\frac{3G}{2c^{2}d^{2}}$ | $\hat{p}_{\text{B}}^{2}\hat{x}_{\text{A}}$ |  |  | $\frac{\sqrt{2}Gm\delta x}{4c^{2}d^{2}}$ |  |  |
$-\frac{3G}{2c^{2}d^{2}}$ | $\hat{p}_{\text{B}}^{2}\hat{x}_{\text{A}}^{2}$ |  |  |  |  |  | $-\frac{3G\hbar}{16c^{2}d^{3}\omega_{\text{m}}}$
$\frac{3G}{c^{2}d^{3}}$ | $\hat{p}_{\text{B}}^{2}\hat{x}_{\text{B}}\hat{x}_{\text{A}}$ | $-\frac{3G\hbar}{8c^{2}d^{3}\omega_{\text{m}}}$ |  |  |  | $\frac{3\sqrt{6}G\hbar}{16c^{2}d^{2}\omega_{\text{m}}}$ |
2PN | $-\frac{9G}{4c^{4}m^{2}d}$ | $\mathbf{\hat{p}_{A}^{2}\hat{p}_{B}^{2}}$ |  |  |  |  |  | $\boldsymbol{\frac{9G\boldsymbol{\hbar}\omega_{\text{m}}}{32c^{4}d}}$
Table 1: List of all gravitational couplings up to quartic order in operators
and up to order 2PN (first three columns). Using first-order perturbation
theory in Eq. (10) we obtain the 25 non-zero coefficients $C_{nN}$ in columns
four to nine (the coefficients $C_{0N}$ and $C_{n0}$ will not contribute to
the entanglement at first order in perturbation theory and are left out of the
table [33]). The _local dip feature_ in the generation of entanglement occurs
as a result of the cancellation arising from the terms
$\propto\hat{x}_{\text{A}}\hat{x}_{\text{B}}$ and
$\propto\hat{p}_{\text{A}}\hat{p}_{\text{B}}$. As discussed in the main text,
this is a result of the different signs of the induced two-mode-squeezing
(TMS) transformation from the leading order 0PN and 1PN terms with the hint in
the different sign of the corresponding $C_{11}$ coefficients. We note that
additional cancellations occur in every column; we see that the $C_{nN}$
coefficients in all columns have both positive and negative values. The _broad
valley feature_ can be understood directly from the couplings. On the one
hand, for large-position delocalization $\Delta x$ the $0$PN term
$\propto\hat{x}_{\text{A}}\hat{x}_{\text{B}}$ generates a rapid increase of
the entanglement entropy with increasing $\Delta x$, and, on the other hand,
for large momentum delocalization $\Delta p$ (i.e., tiny spatial
delocalization) the term $\propto\hat{p}_{A}^{2}\hat{p}_{B}^{2}$ also produces
a fast growth of the entanglement entropy with increasing $\Delta p$. In other
words, entanglement entropy as function of position delocalization is loosely
speaking “U” shaped, forming a valley of entanglement.

From Eq. (1) we find that the cross-coupling terms between the two particles
are given by

$\displaystyle\Delta$
$\displaystyle\hat{H}_{\text{AB}}=\frac{Gm^{2}}{d^{3}}\bigg{(}2\hat{x}_{\text{A}}\hat{x}_{\text{B}}+3\frac{(\hat{x}_{\text{A}}^{2}\hat{x}_{\text{B}}-\hat{x}_{\text{A}}\hat{x}_{\text{B}}^{2})}{d}$
$\displaystyle\quad+\frac{4\hat{x}_{\text{A}}^{3}\hat{x}_{\text{B}}-6\hat{x}_{\text{A}}^{2}\hat{x}_{\text{B}}^{2}+4\hat{x}_{\text{A}}\hat{x}_{\text{B}}^{3}}{d^{2}}\bigg{)}$
$\displaystyle+\frac{G}{c^{2}d}\bigg{(}4\hat{p}_{\text{A}}\hat{p}_{\text{B}}(1+\frac{\hat{x}_{\text{A}}-\hat{x}_{\text{B}}}{d}+\frac{(\hat{x}_{\text{A}}-\hat{x}_{\text{B}})^{2}}{d^{2}})$
$\displaystyle\quad+\frac{3\hat{p}_{\text{A}}^{2}}{2}(\frac{\hat{x}_{\text{B}}}{d}+\frac{2\hat{x}_{\text{A}}\hat{x}_{\text{B}}-\hat{x}_{\text{B}}^{2}}{d^{2}})$
$\displaystyle\quad+\frac{3\hat{p}_{\text{B}}^{2}}{2}(\frac{\hat{x}_{A}}{d}+\frac{2\hat{x}_{\text{A}}\hat{x}_{\text{B}}-\hat{x}_{\text{A}}^{2}}{d^{2}})\bigg{)}-\frac{9G}{4c^{4}m^{2}d}\hat{p}_{A}^{2}\hat{p}_{B}^{2},$
(21)

where the first two lines contain the static limit 0PN contribution, lines
three to five contain the 1PN contribution and the last term corresponds to
the 2PN contribution. While the leading order gravitational force arises from
the terms $\hat{H}_{\text{A}}$ and $\hat{H}_{\text{B}}$ (as these terms
contain the uniform gravitational fields affecting the motion of the
individual particle), the leading order contribution for entanglement
generation arises from the cross-couplings in $\hat{H}_{\text{AB}}$.

In Table. 1, we have listed the couplings from Eq. (7) and applied Eq. (10) to
obtain the non-zero expansion coefficient $C_{nN}$ for $n,N>0$ (the
coefficients $C_{0N}$ and $C_{n0}$ will not contribute to the entanglement at
first order in perturbation theory and are left out of the computation). The
resulting entanglement entropy computed using Eq. (18) as a function of
position delocalization $\Delta x\equiv\delta x$ (bottom axis) and momentum
delocalization $\Delta p\equiv\delta p=\hbar/\Delta x$ (top axis) is plotted
in Fig. 1.

We can, however, capture both qualitatively and quantitatively the
entanglement entropy shown in Fig. 1 by considering only three couplings from
Table. 1. In the right-most part of the figure, the momentum couplings (1PN
and 2PN effects) are negligible as the momentum delocalization $\Delta p$ is
tiny compared to $mc$, and hence dominated by the position 0PN couplings.
Furthermore, if the size of the delocalization $\Delta x$ is also small
compared to the distance between the traps $d$, then we can further neglect
the cubic and quartic couplings, leaving us with the coupling
$\propto\hat{x}_{A}\hat{x}_{B}$. The intermediate plateau is captured by the
1PN coupling $\propto\hat{p}_{A}\hat{p}_{B}$ as the ratio $\Delta x/d$ becomes
tiny, and $\Delta p$ becomes non-negligible compared to $mc$. Finally, the
left-most part of the figure is dominated by the 2PN coupling
$\propto\hat{p}_{A}^{2}\hat{p}_{B}^{2}$, and eventually by higher order PN
corrections as we would further increase $\Delta p$. By making such
simplifications, the perturbed state can be written as

$|\psi_{\text{AB}}\rangle\approx\frac{1}{\mathcal{N}}[|0\rangle|0\rangle-
C_{11}|1\rangle|1\rangle-C_{22}|2\rangle|2\rangle],$ (22)

where $|0\rangle$, $|1\rangle$, $|2\rangle$ denote the number states, and
$\mathcal{N}$ denotes the overall normalization. To compute the entanglement
entropy in Eq. (11), we can now readily use Eq. (18), where we can make the
further approximation $\mathcal{N}\approx 1$ as we have $C_{00}\approx 1$ and
$C_{11},C_{22}\ll 1$ (while the terms $C_{11}^{2},C_{22}^{2}$ appearing in
$\mathcal{N}$ would only contribute higher order corrections).

Let us briefly comment how to see the entanglement dip from Table. 1. We first
recall that the TMS operator can be written in the form
$\hat{S(\xi)}=\exp(-i(\xi^{*}\hat{a}\hat{b}+\xi\hat{a}^{\dagger}\hat{b}^{\dagger}))$,
where the TMS generator is
$\propto\xi^{*}\hat{a}\hat{b}+\xi\hat{a}^{\dagger}\hat{b}^{\dagger}$. From Eq.
(10) we immediately find

$C_{nN}\propto\langle n|\langle
N|\hat{H}_{AB}|0\rangle|0\rangle\propto\pm\langle n|\langle
N|\hat{a}^{\dagger}\hat{b}^{\dagger}+\hat{a}\hat{b}|0\rangle|0\rangle,$ (23)

with the plus (minus) sign corresponding to the 0PN coupling
$\hat{x}_{A}\hat{x}_{B}$ (1PN coupling $\hat{p}_{A}\hat{p}_{B}$). In other
words, the 0PN position coupling would like to squeeze with TMS parameter
$\xi=+1$ while the 1PN momentum coupling would like to squeeze in the opposite
direction with TMS parameter $\xi=-1$. The 0PN and 1PN two-mode squeezing
contributions cancel when $\Delta x/d=\sqrt{2}\Delta p/(mc)$. Summing the two
contributions, we find a total squeezing parameter $\xi=0$ and a suppression
of the gravitationally induced entanglement. The reason for the entanglement
suppression thus lies in the opposite sign of the two-mode squeezing (TMS)
parameter at the leading order 0PN and 1PN gravitational interaction.

##  Appendix C Derivation of the time-dependent entanglement entropy

Figure 3: Entanglement entropy $S$ as a function of time $t$ for different
values of the couplings $g_{x}$, $g_{p}$ defined in Eq. (15) (which scale as a
function of the harmonic frequency $\omega_{\text{m}}$). The SMSV squeezing
parameter is set to $r=-3$ corresponding to initial position squeezing $\Delta
x=\delta xe^{-r}$ and momentum delocalization $\Delta p=\delta pe^{r}$. By
setting $\omega_{\text{m}}=\omega_{0}\equiv c/(\sqrt{2}d)$ we find the case of
equal couplings $g_{x}=g_{p}=g_{0}\equiv\sqrt{2}Gm/(c^{2}d)$ (green line). We
also consider the case $\omega_{\text{m}}=4\omega_{0}$ producing the couplings
$4g_{x}=g_{p}/4=g_{0}$ (orange line), and the case
$\omega_{\text{m}}=\omega_{0}/4$ producing the coupling $g_{x}/4=4g_{p}=g_{0}$
(black line). In all cases, the maximum entanglement is generated at
$\omega_{\text{m}}t=\pi/2$. The entanglement entropy is normalized to the
maximum value $S_{\text{max}}$ of the green curve to ease the comparison with
Fig. 2.

Here, we further explore the dip in entanglement generation by considering the
initial state to be the product of two single-mode-squeezed-vacuum (SMSV)
states:

$|\psi_{\text{i}}\rangle=|r\rangle_{\text{A}}|r\rangle_{\text{B}},$ (24)

where $r\in\text{Re}$ is the SMSV squeezing parameter. The single mode
squeezed state is given by

$|r\rangle=\frac{1}{\sqrt{\text{cosh}r}}\sum_{n=0}^{\infty}(\text{tanh}r)^{n}\frac{\sqrt{(2n)!}}{2^{n}n!}|2n\rangle,$
(25)

where $|n\rangle$ denotes the number state of the considered harmonic
oscillator. The state in Eq. (25) has enhanced position or momentum
delocalization is given by

$\Delta x=\delta xe^{-r},\qquad\Delta p=\delta pe^{r},$ (26)

respectively (with $\delta x$ and $\delta p$ denoting the zero-point-motions
defined below Eq. (5)). Such states are not energy eigenstates, and hence, we
have to take into account their time-evolution in the harmonic traps. However,
as the initial state in Eq. (24) is Gaussian, and the interaction in Eq. (12)
is quadratic in the operators, the state will remain Gaussian also at any
later time.

The Heisenberg equations of motion for the modes of the harmonic oscillators
evolve as [42]:

$\displaystyle\hat{a}(t)$
$\displaystyle=c_{0}(t)\hat{a}+c_{+}(t)\hat{b}+c_{-}(t)\hat{b}^{\dagger},$
(27) $\displaystyle\hat{b}(t)$
$\displaystyle=c_{0}(t)\hat{b}-c_{+}(t)\hat{a}-c_{-}(t)\hat{a}^{\dagger},$
(28)

where $\hat{a}\equiv\hat{a}(0)$, $\hat{b}\equiv\hat{b}(0)$. The time-dependent
coefficients are given by: [42]

$\displaystyle c_{0}(t)$
$\displaystyle=\text{cos}(\omega_{\text{e}}t)-i\frac{\omega_{\text{m}}}{\omega_{\text{e}}}\text{sin}(\omega_{\text{e}}t),$
(29) $\displaystyle c_{\pm}(t)$
$\displaystyle=g_{\pm}\frac{\omega_{\text{m}}}{\omega_{\text{e}}}\text{sin}(\omega_{\text{e}}t),$
(30)

where we have defined the effective frequency
$\omega_{\text{e}}=\sqrt{\omega_{\text{m}}^{2}+g_{+}^{2}-g_{-}^{2}}$. We can
now readily compute the time dependency of the entanglement entropy. Inserting
Eq. (27) in Eq.(20) we find that the problem reduces to evaluating the
expectation values of the initial state in Eq. (24). In particular, to
complete the analysis, we use the following expectation values [34, 35]:

$\displaystyle\langle\hat{a}\hat{a}\rangle$
$\displaystyle=\langle\hat{b}\hat{b}\rangle=\text{-sinh}(r)\text{cosh}(r),$
(31) $\displaystyle\langle\hat{a}{}^{\dagger}\hat{a}{}^{\dagger}\rangle$
$\displaystyle=\langle\hat{b}^{\dagger}\hat{b}^{\dagger}\rangle=\text{-sinh}(r)\text{cosh}(r),$
(32) $\displaystyle\langle\hat{a}\hat{a}{}^{\dagger}\rangle$
$\displaystyle=\langle\hat{b}\hat{b}^{\dagger}\rangle=\text{cosh}^{2}(r),$
(33) $\displaystyle\langle\hat{a}{}^{\dagger}\hat{a}\rangle$
$\displaystyle=\langle\hat{b}^{\dagger}\hat{b}\rangle=\text{sinh}^{2}(r),$
(34)

where the expectation values are computed with respect to the SMSV state in
Eq. (25).

Inserting Eq. (20) in Eq. (19), and neglecting higher order terms in
$g_{x,p}/\omega_{\text{m}}$, we eventually find the formula for the time-
dependent entanglement entropy $S(t)$:

$\displaystyle S(t)\approx$
$\displaystyle-\frac{A(t)}{2\text{$\omega_{\text{m}}^{2}$}}\sin^{2}(\text{$\omega_{\text{m}}$}t)\left(\text{ln}\left[\frac{A(t)}{2\omega_{\text{m}}^{2}}\sin^{2}(\text{$\omega_{\text{m}}$}t)\right]-1\right).$
(35)

where

$\displaystyle A(t)=$ $\displaystyle
g_{p}^{2}-4g_{p}\text{$g_{x}$}+g_{x}^{2}+\left(g_{p}^{2}+g_{x}^{2}\right)\cos(2\text{$\omega_{\text{m}}$}t)$
$\displaystyle+2\left(g_{p}^{2}e^{4r}+g_{x}^{2}e^{-4r}\right)\sin^{2}(\text{$\omega_{\text{m}}$}t).$
(36)

We analyze the temporal behaviour of Eq. (35) in Fig. 3. We first find that
the behaviour remains qualitatively similar as we change the frequency
$\omega_{\text{m}}$. The entanglement entropy $S(t)$ has the maximum at
$\omega_{\text{m}}t=\pi/2$ in all cases. Setting $\omega_{\text{m}}t=\pi/2$ we
find that the maximum entanglement increases both for $r<0$ (i.e., spatial
delocalization $\Delta x>\delta x$) as well as for $r>0$ (momentum
delocalization $\Delta p>\delta p$). This is not surprising as a squeezed
state breaks the symmetry of the ground state, leaving it more exposed to TMS
with either positive or negative values.

The generated entanglement becomes negligible when the condition $A=0$ is met
(at $t=\pi/(2\omega_{\text{m}})$). In this case we are in an entanglement dip,
such that $S(t)=0\,\,\forall t$, as can be noted by computing the limit of Eq.
(35):

$\lim_{A\rightarrow 0}S(t)=0.$ (37)

In particular, if we set $t=\pi/(2\omega_{\text{m}})$ in Eq. (36), and impose
$A=0$, we find the simple condition for the location of the dip:

$g_{x}e^{-2r}=g_{p}e^{2r}.$ (38)

If we set $r=0$ in Eq. (38) we recover the condition $g_{x}=g_{p}$ emerging at
the level of the Hamiltonian in Eqs. (13)-(15) as highlighted in the main
text. In other words, the squeezing effectively changes the quadratic position
and momentum couplings resulting in the modified condition for the dip. We
observe the shifted dip location according to Eq. (38) in Fig. 2.

To summarize, the 0PN position coupling $\hat{x}_{\text{A}}\hat{x}_{\text{B}}$
and the 1PN momentum coupling $\hat{p}_{\text{A}}\hat{p}_{\text{B}}$ induce
two-mode squeezing (TMS) with an opposite sign of the squeezing parameter
$\xi$. The 0PN and 1PN coupling are a source of TMS with squeezing parameter
$\xi=+1$ and $\xi=-1$, which would individually generate TMS entangled states.
However, when the two effects combine, they cancel, resulting in a strong
suppression of gravitationally induced entanglement.

As a safety check, we consider the case where the position and momentum
coupling match such that $g_{-}=g_{x}-g_{p}$ vanishes, and we would thus
expect a strong suppression of entanglement generation. We define

$\omega_{0}\equiv c/(\sqrt{2}d),\qquad g_{0}\equiv\sqrt{2}Gm/(c^{2}d),$ (39)

and find that when $\omega_{\text{m}}=\omega_{0}$ we have equal couplings
$g_{x}=g_{p}=g_{0}$. Using this symmetric coupling regime, we find from Eq.
(36) a simplified expression.

$A(t)\approx
8g_{0}^{2}\sinh^{2}(2r)\sin^{2}(\text{{\hbox{\omega_{\text{m}}}}}t).$ (40)

Setting $\omega_{\text{m}}t=\pi/2$ and taking the limit $r\rightarrow 0$ in
Eq. (35) with $A(t)$ from Eq. (40) we recover that the entanglement entropy
vanishes (case corresponding to the dip found in Fig. 1 for the ground state).