$\displaystyle Y_{1}$ $\displaystyle=$ $\displaystyle
8\pi^{2/3}\left(\frac{1}{\epsilon}+\log(\mu\,x)+3\log(x)+2\gamma+2\log(\pi)+O\left(\epsilon\right)\right)\,,$
(E.10) $\displaystyle Y_{2}$ $\displaystyle=$
$\displaystyle\frac{1}{\pi^{2/3}}\left(\frac{1}{\epsilon}+\log(\mu\,x)+3\log(x)+2\gamma-\frac{1}{4}+2\log(\pi)+O\left(\epsilon\right)\right)\,,$
(E.11)

where the $d$-dimensional and the 4-dimensional couplings are related through
$g_{d}^{2}=g_{4}^{2}\mu^{\epsilon}$. The $3\log(x)$ term results from the
$\epsilon$-expansion of the tree-level factor, dressing the $1/\epsilon$
piece. It cancels out against the wave function renormalisation of the
operator and hence, it does not contribute to the anomalous dimension. The
resulting eigenvalues of the renormalised matrix of two-point functions are
given by

$Y_{\pm}=\frac{\cancel{x}}{32\pi^{6}(x^{2})^{4}}\left(1\mp
2\sqrt{2}\xi_{4/3}^{2}(\log(\mu\,x)+\text{finite})\right)\,.$ (E.12)

The corresponding one-loop anomalous dimensions read

$\gamma_{\pm}=\pm\sqrt{2}\xi_{4/3}^{2}+O(\xi_{4/3}^{4})\,,$ (E.13)

in agreement with (4.13) at $S={1\over 2}$.

##  Appendix F Embedding Formalism and the Spectrum of ${\rm
tr\,}(F_{\mu\nu}Z)$

Some of the computations of this article are done using the so-called
embedding formalism. We review this formalism here in section F.1, and we use
it in section F.2 to compute the eigenvalues of the graph-building operator
$\widehat{B}_{F}$ of section 2.2.

###  F.1 Embedding Formalism

We introduce here a realisation of the principal series representations of the
conformal group as homogeneous functions on the null cone in six dimensions,
with Lorentzian metric $\eta_{MN}$. The degree of homogeneity is $-\Delta$. We
will denote by $x,\theta,\zeta$, etc. the four-dimensional vectors and by
$Y,\Theta,Z$, etc. the six-dimensional ones. The use of the embedding
formalism is quite standard, see ([40] and references therein. The idea of
working in six dimensions, where conformal transformations are linearly
realised, dating back to [41]. We shall focus here on the scalar and rank-2
antisymmetric tensors representations, (see [42] for more general cases).

For scalar representations, the map between the two realisations sends a
$f:\mathbb{R}^{4}\rightarrow\mathbb{C}$ of dimension $\Delta$ to

$\overline{f}(Y)\equiv\frac{1}{(Y^{+})^{\Delta}}f\left(\frac{Y^{\mu}}{Y^{+}}\right)\,.$
(F.1)

The inverse sends a homogeneous function $F$ on the null cone
$\\{Y\in\mathbb{R}^{6};Y^{M}Y_{M}=0\\}$ to

$\overline{F}(x)=F(\bar{Y})\,,$ (F.2)

where

$\bar{Y}^{\mu}=x^{\mu}\,,\quad\bar{Y}^{+}=1\,,\quad\bar{Y}^{-}=x^{2}\,.$ (F.3)

Clearly, $\overline{\overline{f}}=f$ and $\overline{\overline{F}}=F$, (on the
null cone). It is also straightforward to check that the generators
$q_{MN}^{(\Delta,0,0)}$ given in subsection D.2 are mapped to

$\mathfrak{q}^{(\Delta,0,0)}_{MN}=Y_{M}\partial_{Y^{N}}-Y_{N}\partial_{Y^{M}}\,,$
(F.4)

in embedding space.

For rank-2 antisymmetric representations, there is an additional dependence on
an auxiliary polarisation vector $\theta$ (such that
$\\{\theta^{\mu},\theta^{\nu}\\}=0$) in the usual realisation. In the
embedding formalism, the functions $\Psi$ depend on $Y$ and $\Theta$ and
satisfy the following requirements:

  * •

they are defined on $Y^{M}Y_{M}=Y^{M}\Theta_{M}=0$,

  * •

they are homogeneous of degree $-\Delta$ in $Y$ and of degree $2$ in $\Theta$,

  * •

they are transverse: $\Psi(Y,\Theta+\alpha Y)=\Psi(Y,\Theta)$ for arbitrary
$\alpha$ anticommuting with $\Theta$.

A function $\psi$ is now sent to

$\overline{\psi}(Y,\Theta)\equiv\frac{1}{(Y^{+})^{\Delta}}\psi\left(\frac{Y^{\mu}}{Y^{+}},\Theta^{\mu}-\Theta^{+}\frac{Y^{\mu}}{Y^{+}}\right)\,.$
(F.5)

The inverse map sends $\Psi$ to

$\overline{\Psi}(x,\theta)=\Psi(\bar{Y},\bar{\Theta})\,,$ (F.6)

where

$\bar{\Theta}^{\mu}=\theta^{\mu}\,,\quad\bar{\Theta}^{+}=0\,,\quad\bar{\Theta}^{-}=2\,\theta\cdot
x\,.$ (F.7)

The generators of conformal transformations, given by (D.18)-(D.21) in
physical space, are mapped to

$\mathfrak{q}_{MN}=Y_{M}\partial_{Y^{N}}-Y_{N}\partial_{Y^{M}}+\Theta_{M}\partial_{\Theta^{N}}-\Theta_{N}\partial_{\Theta^{M}}\,,$
(F.8)

in embedding space. Moreover, when $\Delta=2$, the analogue of
$\psi=\theta\cdot\partial_{x}\phi$ is $\Psi=\Theta\cdot\partial_{Y}\,\Phi$
with $\Phi$ homogeneous of degree $-1$ in $Y$.

###  F.2 Spectrum

To compute the spectrum of $\widehat{B}_{F}$ in (2.41), we first translate
(2.42) to embedding space. We start with the transfer matrix. The main
advantage of this new realisation is that the computations are much simpler.
For instance, the Lax matrix (D.5) for a scalar of dimension $1$
(corresponding to the first site of the chain) is easily verified to be

$L_{Y}^{(1,0,0;\mathbf{6})}(u)=\left[u^{2}\eta^{MN}-u\left(Y^{M}\partial_{Y}^{N}-Y^{M}\partial_{Y}^{N}\right)-\frac{1}{2}Y^{M}Y^{N}\Box_{Y}\right]\otimes
e_{MN}\,.$ (F.9)

For rank-2 antisymmetric representations of dimension 2 (second site), one has

$\mathfrak{q}_{MP}\mathfrak{q}\,?\\!^{P}_{N}?-2\,\mathfrak{q}_{MN}=-\Theta\cdot\partial_{Y}\,Y_{M}Y_{N}\,\partial_{Y}\cdot\partial_{\Theta}-\partial_{Y}\cdot\partial_{\Theta}\,Y_{M}Y_{N}\,\Theta\cdot\partial_{Y}\,.$
(F.10)

Consequently, the transfer matrix (D.7) at $u=0$ is

$\widehat{T}_{\mathbf{6}}(0;0,0)=\frac{1}{4}\left[\Theta\cdot\partial_{Y_{2}}\,(Y_{1}\cdot
Y_{2})^{2}\,\partial_{Y_{2}}\cdot\partial_{\Theta}+\partial_{2}\cdot\partial_{\Theta}\,(Y_{1}\cdot
Y_{2})^{2}\,\Theta\cdot\partial_{Y_{2}}\right]\Box_{Y_{1}}\,.$ (F.11)

The eigenvalue equation $\widehat{T}_{\mathbf{6}}(0;0,0)\Psi=\frac{1}{E}\Psi$
for $\Psi=\Theta\cdot\partial_{Y_{2}}\Phi$ becomes

$\Theta\cdot\partial_{2}\left[(Y_{1}\cdot
Y_{2})^{2}\,\partial_{Y_{2}}\cdot\partial_{\Theta}\,\Theta\cdot\partial_{Y_{2}}\Box_{Y_{1}}\Phi-\frac{4}{E}\Phi\right]=0\,.$
(F.12)

As usual, the solutions $\Phi$ should have the form of three-point functions.
One of the operators is a scalar of dimension 1 and another is a vector of
dimension 1. The only possibility for the third operator is either a symmetric
traceless tensor or some mixed symmetry tensor.

##### Symmetric traceless tensors

We use a four-dimensional auxiliary polarisation vector $\zeta$ to
parameterise the tensor structure. Because the tensor is traceless, it
suffices to consider null polarisation vectors $\zeta^{2}=0$. We denote by $Z$
the polarisation vector in embedding space. It is defind such that its
projection to four dimensions flat space, is similar to (F.7), reads

$\bar{Z}^{\mu}=\zeta^{\mu}\,,\quad\bar{Z}^{+}=0\,,\quad\bar{Z}^{-}=2\,\zeta\cdot
x_{0}\,.$ (F.13)

The three-point function $\Phi$ is therefore a function of $Z$ and $Y_{0}$, as
well as $Y_{1}$, $Y_{2}$ and $\Theta$ described above. It is defined on the
space $Z^{2}=Z\cdot Y_{0}=Y_{i}^{2}=0$ and is transverse. Namely,
$\Phi(Z+\alpha Y_{0})=\Phi(Z)$ for arbitrary $\alpha\in\mathbb{R}$. Moreover,
it is an homogeneous function of degree $S$ in $Z$ and of degree $-\Delta$ in
$Y_{0}$.

For $S=0$, the only possibility is

$\Phi\propto\Theta\cdot\partial_{Y_{2}}\left(\frac{1}{(Y_{1}\cdot
Y_{2})^{\frac{1-\Delta}{2}}(Y_{1}\cdot Y_{0})^{\frac{1+\Delta}{2}}(Y_{2}\cdot
Y_{0})^{\frac{\Delta-1}{2}}}\right)\,,$ (F.14)

which trivially satisfies (F.12), and gives $\Psi=0$. We can thus only
consider $S\geqslant 1$. In that case, there are a priori two possible
structures:

$\Phi_{1}=\frac{J_{\Theta,Z}K_{Z}^{S-1}}{(Y_{1}\cdot
Y_{2})^{\frac{1+S-\Delta}{2}}(Y_{1}\cdot
Y_{0})^{\frac{\Delta+S-1}{2}}(Y_{2}\cdot Y_{0})^{\frac{\Delta+S+1}{2}}}$
(F.15)

and

$\Phi_{2}=\frac{K_{\Theta}K_{Z}^{S}}{(Y_{1}\cdot
Y_{2})^{\frac{3+S-\Delta}{2}}(Y_{1}\cdot
Y_{0})^{\frac{\Delta+S+1}{2}}(Y_{2}\cdot Y_{0})^{\frac{\Delta+S+1}{2}}}\,,$
(F.16)

where

$J_{\Theta,Z}=(\Theta\cdot Z)(Y_{2}\cdot Y_{0})-(\Theta\cdot Y_{0})(Z\cdot
Y_{2})\,,\quad J_{\Theta,H}=(\Theta\cdot H)(Y_{2}\cdot Y_{0})-(\Theta\cdot
Y_{0})(H\cdot Y_{2})\,,$ (F.17)

and

$K_{Z}=(Z\cdot Y_{1})(Y_{2}\cdot Y_{0})-(Z\cdot Y_{2})(Y_{1}\cdot
Y_{0})\,,\quad K_{\Theta}=(\Theta\cdot Y_{1})(Y_{2}\cdot Y_{0})-(\Theta\cdot
Y_{0})(Y_{1}\cdot Y_{2})\,.$ (F.18)

However, these two structures are related via

$\Theta\cdot\partial_{Y_{2}}\left(\frac{K_{Z}^{S}}{(Y_{1}\cdot
Y_{2})^{\frac{1+S-\Delta}{2}}(Y_{1}\cdot
Y_{0})^{\frac{1+S+\Delta}{2}}(Y_{2}\cdot
Y_{0})^{\frac{\Delta+S-1}{2}}}\right)=-S\Phi_{1}+\frac{\Delta-1-S}{2}\Phi_{2}\,,$
(F.19)

which means that they give rise to the same $\Psi$, and are therefore
equivalent. We now plug $\Phi_{1}$ in equation (F.12). Firstly,

$\Box_{Y_{1}}\Phi_{1}\approx-\frac{(\Delta-S-1)(\Delta+S-1)}{2}\frac{J_{\Theta,Z}K_{Z}^{S-1}}{(Y_{1}\cdot
Y_{2})^{\frac{3+S-\Delta}{2}}(Y_{1}\cdot
Y_{0})^{\frac{1+S+\Delta}{2}}(Y_{2}\cdot Y_{0})^{\frac{-1+S+\Delta}{2}}}$
(F.20)

where $\approx$ means equal up to terms proportional to
$Y_{2}^{2},Y_{0}^{2},Z\cdot Y_{0}$, etc., which we can consistently ignore as
they project to zero in physical space and are stable under the various
differential operators we apply. Then,

$(Y_{1}\cdot
Y_{2})^{2}\,\partial_{Y_{2}}\cdot\partial_{\Theta}\,\Theta\cdot\partial_{Y_{2}}\Box_{Y_{1}}\Phi_{1}\approx\frac{(\Delta-S-1)(\Delta+S-1)(\Delta-S-3)}{4}\left[(\Delta-3)\Phi_{1}+\frac{\Delta-S-5}{2}\Phi_{2}\right]\,.$
(F.21)

Using (F.19), one verifies that equation (F.12) is satisfied and that the
eigenvalue is

$E(\Delta,S,S)=\frac{16}{(\Delta+S-1)^{2}(\Delta-S-3)^{2}}\,.$ (F.22)

##### Mixed symmetry tensors

In this case, the tensors we are interested in are of the form
$f^{\mu\nu\mu_{1}\cdots\mu_{\ell}}$ with the following symmetry properties:
they are antisymmetric in $(\mu,\nu)$, completely symmetric in
$(\mu_{1},\dots,\mu_{\ell})$, and they also satisfy, for all
$k\in\\{1,\dots,\ell\\}$,

$f^{\mu\nu\mu_{1}\cdots\mu_{\ell}}=f^{\mu_{k}\nu\mu_{1}\cdots\mu\cdots\mu_{\ell}}+f^{\mu\mu_{k}\mu_{1}\cdots\nu\cdots\mu_{\ell}}\,,$
(F.23)

where on the right-hand side, $\mu$ and $\nu$ appear at the $k$’th position.
Additionally, the tensors should be completely traceless. In order to get rid
of the index structure, we consider
$f(\zeta,\eta)=f^{\mu\nu\mu_{1}\cdots\mu_{\ell}}\eta_{\mu}\eta_{\nu}\zeta_{\mu_{1}}\dots\zeta_{\mu_{\ell}}$
for two arbitrary vectors satisfying

$\zeta^{2}=\eta^{2}=\zeta\cdot\eta=0\,,\qquad[\zeta_{\mu},\zeta_{\nu}]=\\{\eta_{\mu},\eta_{\nu}\\}=0\,.$
(F.24)

It is important to notice that not all homogeneous polynomials, of degree $2$
in $\eta$ and $\ell$ in $\zeta$, come from a tensor with the symmetry we want.
Indeed, we have not taken into account the constraints coming from (F.23). In
particular, there should exist a polynomial $\tilde{f}$ such that

$f=\eta\cdot\partial_{\zeta}\tilde{f}\,.$ (F.25)

One could take
$\tilde{f}(\zeta,\eta)=\frac{2}{\ell+2}f^{\mu\nu\mu_{1}\cdots\mu_{\ell}}\zeta_{\mu}\eta_{\nu}\zeta_{\mu_{1}}\dots\zeta_{\mu_{\ell}}$
for instance. In embedding space, the auxiliary vectors are promoted to $Z$
and $H$ satisfying the same properties and, additionally, we restrict
ourselves to $Z\cdot Y_{0}=H\cdot Y_{0}=0$. The three-point function $\Phi$ is
homogeneous of degree $2$ in $H$, of degree $\ell$ in $Z$, and transverse in
both variables

$\Phi(H+\alpha Y_{0},Z)=\Phi(H,Z+Y_{0})=\Phi(H,Z)\,,$ (F.26)

where $\alpha$ anti-commute with $H$. Moreover, there should exist some
function $\widetilde{\Phi}$ such that

$\Phi=H\cdot\partial_{Z}\widetilde{\Phi}\,.$ (F.27)

As is now usual, the polarisation vectors are projected to (F.13) and

$\bar{H}^{\mu}=\eta^{\mu}\,,\qquad\bar{H}^{+}=0\,,\qquad\bar{H}^{-}=2\,\eta\cdot
x_{0}\,.$ (F.28)

Furthermore, one may easily verify that (F.27) implies

$\overline{\Phi}=\eta\cdot\partial_{\zeta}\overline{\widetilde{\Phi}}\,.$
(F.29)

We now examine the possible structures. The only one which does not involve
$\varepsilon^{ABCDEF}$ is

$\Phi^{\prime}_{1}=\frac{K_{H}J_{\Theta,H}K_{Z}^{\ell}}{(Y_{1}\cdot
Y_{2})^{\frac{3+\ell-\Delta}{2}}(Y_{1}\cdot
Y_{0})^{\frac{1+\ell+\Delta}{2}}(Y_{2}\cdot
Y_{0})^{\frac{3+\ell+\Delta}{2}}}\\\
=\frac{H\cdot\partial_{Z}}{\ell+1}\left[\frac{J_{\Theta,H}K_{Z}^{\ell+1}}{(Y_{1}\cdot
Y_{2})^{\frac{3+\ell-\Delta}{2}}(Y_{1}\cdot
Y_{0})^{\frac{1+\ell+\Delta}{2}}(Y_{2}\cdot
Y_{0})^{\frac{3+\ell+\Delta}{2}}}\right]\,,$ (F.30)

where $J_{\Theta,H}$ and $K_{Z}$ are defined in (F.17) and (F.18) above, and

$K_{H}=(H\cdot Y_{1})(Y_{2}\cdot Y_{0})-(H\cdot Y_{2})(Y_{1}\cdot Y_{0})\,.$
(F.31)

Similarly, there is only one structure involving $\varepsilon^{ABCDEF}$

$\Phi^{\prime}_{2}=\frac{\varepsilon(Y_{0},Y_{1},Y_{2},HK_{Z}+\ell
ZK_{H},H,\Theta)K_{Z}^{\ell-1}}{(Y_{1}\cdot
Y_{2})^{\frac{3+\ell-\Delta}{2}}(Y_{1}\cdot
Y_{0})^{\frac{1+\ell+\Delta}{2}}(Y_{2}\cdot
Y_{0})^{\frac{1+\ell+\Delta}{2}}}\\\
=H\cdot\partial_{Z}\left[\frac{\varepsilon(Y_{0},Y_{1},Y_{2},Z,H,\Theta)K_{Z}^{\ell}}{(Y_{1}\cdot
Y_{2})^{\frac{3+\ell-\Delta}{2}}(Y_{1}\cdot
Y_{0})^{\frac{1+\ell+\Delta}{2}}(Y_{2}\cdot
Y_{0})^{\frac{1+\ell+\Delta}{2}}}\right]\,,$ (F.32)

with

$\varepsilon(Y_{0},Y_{1},Y_{2},Z,H,\Theta)=\varepsilon_{ABCDEF}Y_{0}^{A}Y_{1}^{B}Y_{2}^{C}Z^{D}H^{E}\Theta^{F}\,.$
(F.33)

Since $\partial_{2}\cdot\partial_{\Theta}\Phi^{\prime}_{2}=0$ and

$(Y_{1}\cdot
Y_{2})^{2}\Box_{Y_{2}}\Box_{Y_{1}}\Phi^{\prime}_{2}\approx\frac{(\Delta+\ell+1)(\Delta+\ell-1)(\Delta-\ell-3)(\Delta-\ell-5)}{4}\Phi^{\prime}_{2}\,,$
(F.34)

we find that equation (F.12) holds true, with eigenvalue

$E(\Delta,\ell+2,\ell)=\frac{16}{(\Delta+\ell+1)(\Delta+\ell-1)(\Delta-\ell-3)(\Delta-\ell-5)}\,.$
(F.35)

A similar, though more cumbersome, computation, would show that
$\Phi^{\prime}_{1}$ shares the same eigenvalue. Instead of doing this, we will
show that their projections to four dimensional flat space, (2.49) and (2.50),
are related as

$\Phi^{\rho,2}_{x_{0},\zeta,\eta}=\frac{1}{S}\left[\epsilon(\zeta,\eta,\partial_{\zeta},\partial_{\eta})-\frac{1}{2}\epsilon(\eta,\eta,\partial_{\eta},\partial_{\eta})\right]\Phi^{\rho,1}_{x_{0},\zeta,\eta}\,,$
(F.36)

and

$\Phi^{\rho,1}_{x_{0},\zeta,\eta}=\frac{S}{(S+1)^{2}}\left[\epsilon(\zeta,\eta,\partial_{\zeta},\partial_{\eta})-\frac{1}{2}\epsilon(\eta,\eta,\partial_{\eta},\partial_{\eta})\right]\Phi^{\rho,2}_{x_{0},\zeta,\eta}\,.$
(F.37)

We recall that
$\Phi^{\rho,i}_{x_{0},\zeta,\eta}=(\eta\cdot\partial_{\zeta})\widetilde{\Phi}^{\rho,i}_{x_{0},\zeta,\eta}$
with the tilde function being linear in $\eta$. Moreover,

$\left[\epsilon(\zeta,\eta,\partial_{\zeta},\partial_{\eta})-\frac{1}{2}\epsilon(\eta,\eta,\partial_{\eta},\partial_{\eta}),\eta\cdot\partial_{\zeta}\right]=0\,.$
(F.38)

Using this, it is straight forward to check (F.36). Now that we have proven
(F.36), we insert it into the second equation, which is thus equivalent to

$\Phi^{\rho,1}_{x_{0},\zeta,\eta}=\frac{1}{(S+1)^{2}}(\eta\cdot\partial_{\zeta})\left[\epsilon(\zeta,\eta,\partial_{\zeta},\partial_{\eta})\right]^{2}\widetilde{\Phi}^{\rho,1}_{x_{0},\zeta,\eta}\,.$
(F.39)

Since $\widetilde{\Phi}^{\rho,1}_{x_{0},\zeta,\eta}$ is linear in $\eta$, we
can replace
$\left[\epsilon(\zeta,\eta,\partial_{\zeta},\partial_{\eta})\right]^{2}$ with
$\epsilon^{\alpha\beta\gamma\delta}\eta_{\alpha}\zeta_{\beta}\partial_{\zeta^{\gamma}}?\epsilon^{\mu\nu\rho}_{\delta}?\zeta_{\mu}\partial_{\eta^{\nu}}\partial_{\zeta^{\rho}}$
so that the equation becomes

$\displaystyle\Phi^{\rho,1}_{x_{0},\zeta,\eta}$
$\displaystyle=\frac{1}{(S+1)^{2}}(\eta\cdot\partial_{\zeta})\left[\epsilon^{\alpha\beta\gamma\delta}\eta_{\alpha}\zeta_{\beta}\partial_{\zeta^{\gamma}}?\epsilon^{\mu\nu\rho}_{\delta}?\zeta_{\mu}\partial_{\eta^{\nu}}\partial_{\zeta^{\rho}}\right]\widetilde{\Phi}^{\rho,1}_{x_{0},\zeta,\eta}$
$\displaystyle=\frac{1}{S+1}(\eta\cdot\partial_{\zeta})\left[(\eta\cdot\partial_{\eta})(\zeta\cdot\partial_{\zeta})-\zeta^{\beta}(\eta\cdot\partial_{\zeta})\partial_{\eta^{\beta}}\right]\widetilde{\Phi}^{\rho,1}_{x_{0},\zeta,\eta}$
$\displaystyle=(\eta\cdot\partial_{\zeta})\widetilde{\Phi}^{\rho,1}_{x_{0},\zeta,\eta}=\Phi^{\rho,1}_{x_{0},\zeta,\eta}\,,$

where we have used the fact that
$\eta\cdot\zeta=\zeta^{2}=(\eta\cdot\partial_{\zeta})^{2}=0$ and that
$\widetilde{\Phi}^{(\rho^{\prime}),1}_{x_{0},\zeta,\eta}$ is homogeneous of
degree $S$ in $\zeta$.

##  Appendix G Graph-Building Operator for Section 3

In this appendix, we give the explicit expressions for the matrix elements
$\widehat{\mathcal{B}}_{ij}$ of the matrix graph-building operator introduced
in section 3, see (3.3). In order to emphasise that the operators do not all
act on the same space, we will use the notations $\Psi_{\emptyset},\Psi_{A},$
and $\Psi_{X}$ for functions belonging to different spaces. The first one is a
function of $J$ variables $(x_{1},\dots,x_{J})$, but the other two also depend
on $x_{0}$. Moreover, $\Psi_{A}$ is linear in $\theta$, i.e. carries one
index.

There is a small freedom in the way we define the graph-building operator. For
instance, one could replace $\widehat{\mathcal{B}}$ with
$\mathcal{D}\widehat{\mathcal{B}}\mathcal{D}^{-1}$ for an arbitrary constant
diagonal matrix $\mathcal{D}$. In general, we also expect the operator
$\widehat{G}_{0}$ to be promoted to a matrix, and that different matrix
elements of $\frac{1}{1-\widehat{\mathcal{B}}}\widehat{G}_{0}$ correspond to
correlators of different non-local operators. In this paper, we only study the
first few orders in $g$ of the graph-building operator, and it is possible to
take $\widehat{G}_{0}$ to be diagonal. With our convention, we can then write

$\langle{\rm tr\,}(Z(x_{1})\dots Z(x_{J}))\,{\rm tr\,}(Z^{\dagger}(z_{J})\dots
Z^{\dagger}(z_{1}))\rangle\\\
=\sum_{i=1}^{J}\int\frac{\prod_{i=1}^{J}\differential^{4}y_{i}}{\pi^{2J}}\frac{\bra{x_{i},x_{i+1},\dots,x_{i-1}}\left(\frac{1}{1-\widehat{\mathcal{B}}}\right)_{\emptyset\emptyset}\ket{y_{1},\dots,y_{J}}}{(4\pi^{2})^{J}\prod_{j=1}^{J}(y_{j}-z_{j})^{2}}\,,$
(G.1)

where $x_{J+1}=x_{1}$. Similarly,

$\langle{\rm tr\,}(A^{\mu}(x_{0})Z(x_{1})\dots Z(x_{J}))\,{\rm
tr\,}(Z^{\dagger}(z_{J})\dots Z^{\dagger}(z_{1}))\rangle\\\
=-\frac{\operatorname{i}}{2}\int\frac{\prod_{i=1}^{J}\differential^{4}y_{i}}{\pi^{2J}}\frac{\bra{x_{0},x_{1},\dots,x_{J}}\left(\frac{1}{1-\widehat{\mathcal{B}}}\right)^{\mu}_{A\emptyset}\ket{y_{1},\dots,y_{J}}}{(4\pi^{2})^{J}\prod_{i=1}^{J}(y_{i}-z_{i})^{2}}\,,$
(G.2)

and

$\langle{\rm tr\,}(Z(x_{1})\dots Z(x_{J}))\,{\rm
tr\,}(X^{\dagger}(z_{J})Z^{\dagger}(z_{J})\dots
Z^{\dagger}(z_{1})X(z_{1}))\rangle\\\
=g\sum_{i=1}^{J}\int\frac{\prod_{i=0}^{J}\differential^{4}y_{i}}{\pi^{2(J+1)}}\frac{\bra{x_{i},x_{i+1},\dots,x_{i-1}}\left(\frac{1}{1-\widehat{\mathcal{B}}}\right)_{\emptyset
X}\ket{y_{0},y_{1},\dots,y_{J}}}{(4\pi^{2})^{J+2}(y_{0}-z_{1})^{2}(y_{0}-z_{J})^{2}\prod_{j=1}^{J}(y_{j}-z_{j})^{2}}\,,$
(G.3)

for instance. Ultimately, the only important point is to define
$\widehat{\mathcal{B}}$ and $\widehat{G}_{0}$ such that a given matrix element
of $\frac{1}{1-\widehat{\mathcal{B}}}\widehat{G}_{0}$ contains all the Feynman
diagrams appearing in the expansion of a specific correlator, and only those
diagrams.

The fishnet graph-building operator is

$\left[\widehat{\mathcal{B}}_{\emptyset\emptyset}\Psi_{\emptyset}\right](x_{1},\dots,x_{J})=8\pi^{2}\int\frac{\prod_{i=1}^{J}\differential^{4}y_{i}}{\pi^{2J}}\frac{\Psi_{\emptyset}(y_{1},\dots,y_{J})}{\prod_{i=1}^{J}(x_{i}-y_{i})^{2}y_{i,i+1}^{2}}$
(G.4)

with inverse
$\widehat{\mathcal{B}}_{\emptyset\emptyset}^{-1}=(8\pi^{2})^{-1}\prod_{i=1}^{J}x_{i,i+1}^{2}\prod_{i=1}^{J}\frac{\Box_{x_{i}}}{-4}$.
The other operators are

$\displaystyle\left[\widehat{\mathcal{B}}_{\emptyset
A}\Psi_{A}\right](x_{1},\dots,x_{J})$
$\displaystyle=4\pi^{2}\int\frac{\prod_{i=0}^{J}\differential^{4}y_{i}}{\pi^{2(J+1)}}\frac{\Psi_{A}^{\rho}(y_{0},y_{1},\dots,y_{J})}{\prod_{i=1}^{J}(x_{i}-y_{i})^{2}y_{0J}^{2}\prod_{i=0}^{J-1}y_{i,i+1}^{2}}\left(\frac{y_{J0,\rho}}{y_{J0}^{2}}-\frac{y_{10,\rho}}{y_{10}^{2}}\right)\,,$
(G.5) $\displaystyle\left[\widehat{\mathcal{B}}_{\emptyset
X}\Psi_{X}\right](x_{1},\dots,x_{J})$
$\displaystyle=2\pi^{2}\int\frac{\prod_{i=0}^{J}\differential^{4}y_{i}}{\pi^{2(J+1)}}\frac{\Psi_{X}(y_{0},y_{1},\dots,y_{J})}{\prod_{i=1}^{J}(x_{i}-y_{i})^{2}y_{0J}^{2}\prod_{i=0}^{J-1}y_{i,i+1}^{2}}\,,$
(G.6)
$\displaystyle\left[\widehat{\mathcal{B}}_{A\emptyset}\Psi_{\emptyset}\right]^{\mu}(x_{0},\dots,x_{J})$
$\displaystyle=4\pi^{2}\int\frac{\differential^{4}z\prod_{i=1}^{J}\differential^{4}y_{i}}{\pi^{2(J+1)}}\frac{\Delta^{\mu}_{\sigma}(x_{0}-z)\Psi_{\emptyset}(y_{1},\dots,y_{J})}{\prod_{i=1}^{J}y_{i,i+1}^{2}}\sum_{i=1}^{J}\frac{y_{i,i+1}^{2}}{\prod_{j=1}^{J}(x_{j}-y_{i+j})^{2}}$
$\displaystyle\quad\times\frac{1}{(y_{i}-z)^{2}(y_{i+1}-z)^{2}}\left(\frac{(y_{i}-z)^{\sigma}}{(y_{i}-z)^{2}}-\frac{(y_{i+1}-z)^{\sigma}}{(y_{i+1}-z)^{2}}\right)\,,$
(G.7)
$\displaystyle\left[\widehat{\mathcal{B}}_{X\emptyset}\Psi_{\emptyset}\right](x_{0},\dots,x_{J})$
$\displaystyle=2\int\frac{\prod_{i=1}^{J}\differential^{4}y_{i}}{\pi^{2J}}\frac{\Psi_{\emptyset}(y_{1},\dots,y_{J})}{\prod_{i=1}^{J}y_{i,i+1}^{2}}\sum_{i=1}^{J}\frac{y_{i,i+1}^{2}}{\prod_{j=1}^{J}(x_{j}-y_{i+j})^{2}(x_{0}-y_{i})^{2}(x_{0}-y_{i+1})^{2}}\,,$
(G.8)

and

$\left[\widehat{\mathcal{B}}_{AA}\Psi_{A}\right]^{\mu}(x_{0},\dots,x_{J})=4\pi^{2}\int\frac{\differential^{4}z\prod_{i=1}^{J}\differential^{4}y_{i}}{\pi^{2(J+1)}}\left[\frac{\Delta^{\mu}_{\sigma}(x_{0}-z)\Psi_{A}^{\sigma}(z,y_{1},\dots,y_{J})}{\prod_{i=1}^{J}(x_{i}-y_{i})^{2}\prod_{i=1}^{J-1}y_{i,i+1}^{2}(z-y_{1})^{2}(z-y_{J})^{2}}\right.\\\
+\frac{1}{2}\sum_{i=1}^{J-1}\int{\differential^{4}y_{0}\over\pi^{2}}\frac{\Delta^{\mu}_{\sigma}(x_{0}-z)\Psi_{A}^{\rho}(y_{0},\dots,y_{J})}{\prod_{j=1}^{J}(x_{j}-y_{i+j})^{2}y_{0J}^{2}\prod_{i=0}^{J-1}y_{j,j+1}^{2}}\left(\frac{y_{J0,\rho}}{y_{J0}^{2}}-\frac{y_{10,\rho}}{y_{10}^{2}}\right)\\\
\times\frac{y_{i,i+1}^{2}}{(y_{i}-z)^{2}(y_{i+1}-z)^{2}}\left(\frac{(y_{i}-z)^{\sigma}}{(y_{i}-z)^{2}}-\frac{(y_{i+1}-z)^{\sigma}}{(y_{i+1}-z)^{2}}\right)\\\
+\frac{1}{4}\int{\differential^{4}y_{0}\over\pi^{2}}\frac{\Delta^{\mu}_{\sigma}(z-x_{0})\Psi_{A}^{\rho}(y_{0},\dots,y_{J})}{\prod_{j=1}^{J}(x_{j}-y_{j})^{2}\prod_{j=1}^{J-1}y_{j,j+1}^{2}}\Bigg{[}\frac{1}{(y_{J}-z)^{2}}\frac{\partial}{\partial
z_{\sigma}}\left(\frac{y_{10,\rho}}{y_{10}^{4}(z-y_{0})^{2}}-\frac{(z-y_{0})_{\rho}}{(z-y_{0})^{4}y_{10}^{2}}\right)\\\
-\left(\frac{y_{10,\rho}}{y_{10}^{4}(z-y_{0})^{2}}-\frac{(z-y_{0})_{\rho}}{(z-y_{0})^{4}y_{10}^{2}}\right)\frac{\partial}{\partial
z_{\sigma}}\frac{1}{(y_{J}-z)^{2}}+\frac{1}{(y_{1}-z)^{2}}\frac{\partial}{\partial
z_{\sigma}}\left(\frac{y_{J0,\rho}}{y_{J0}^{4}(z-y_{0})^{2}}-\frac{(z-y_{0})_{\rho}}{(z-y_{0})^{4}y_{J0}^{2}}\right)\\\
\left.-\left(\frac{y_{J0,\rho}}{y_{J0}^{4}(z-y_{0})^{2}}-\frac{(z-y_{0})_{\rho}}{(z-y_{0})^{4}y_{J0}^{2}}\right)\frac{\partial}{\partial
z_{\sigma}}\frac{1}{(y_{1}-z)^{2}}\Bigg{]}\right]\,.$ (G.9)

Similarly,

$\left[\widehat{\mathcal{B}}_{AX}\Psi_{X}\right](x_{0},\dots,x_{J})=\pi^{2}\int\frac{\differential^{4}z\prod_{i=0}^{J}\differential^{4}y_{i}}{\pi^{2(J+2)}}\frac{\Delta^{\mu}_{\sigma}(z-x_{0})\Psi_{X}(y_{0},\dots,y_{J})}{y_{0J}^{2}\prod_{i=0}^{J-1}y_{i,i+1}^{2}}\\\
\times\Bigg{[}\sum_{i=0}^{J-1}\frac{y_{i,i+1}^{2}}{\prod_{j=1}^{J}(x_{j}-y_{i+j})^{2}(y_{i}-z)^{2}(y_{i+1}-z)^{2}}\left(\frac{(y_{i}-z)^{\sigma}}{(y_{i}-z)^{2}}-\frac{(y_{i+1}-z)^{\sigma}}{(y_{i+1}-z)^{2}}\right)\\\
+\frac{y_{J0}^{2}}{\prod_{i=1}^{J}(x_{i}-y_{i})^{2}(y_{0}-z)^{2}(y_{J}-z)^{2}}\left(\frac{(y_{J}-z)^{\sigma}}{(y_{J}-z)^{2}}-\frac{(y_{0}-z)^{\sigma}}{(y_{0}-z)^{2}}\right)\Bigg{]}\,,$
(G.10)
$\left[\widehat{\mathcal{B}}_{XA}\Psi_{A}\right]^{\mu}(x_{0},\dots,x_{J})=\int\frac{\prod_{i=0}^{J}\differential^{4}y_{i}}{\pi^{2(J+1)}}\frac{\Psi_{A}^{\rho}(y_{0},\dots,y_{J})}{\prod_{i=1}^{J-1}y_{i,i+1}^{2}}\Bigg{[}\frac{1}{\prod_{i=1}^{J}(x_{i}-y_{i})^{2}(x_{0}-y_{J})^{2}(x_{0}-y_{1})^{2}}\\\
\times\left(\frac{(x_{0}-y_{1})^{2}}{y_{10}^{2}(x_{0}-y_{0})^{2}}\left(\frac{(x_{0}-y_{0})_{\rho}}{(x_{0}-y_{0})^{2}}-\frac{y_{10,\rho}}{y_{10}^{2}}\right)+\frac{(x_{0}-y_{J})^{2}}{y_{J0}^{2}(x_{0}-y_{0})^{2}}\left(\frac{y_{J0,\rho}}{y_{J0}^{2}}-\frac{(x_{0}-y_{0})_{\rho}}{(x_{0}-y_{0})^{2}}\right)\right)\\\
+\frac{1}{y_{10}^{2}y_{J0}^{2}}\left(\frac{y_{J0,\rho}}{y_{J0}^{2}}-\frac{y_{10,\rho}}{y_{10}^{2}}\right)\sum_{i=1}^{J-1}\frac{y_{i,i+1}^{2}}{\prod_{j=1}^{J}(x_{j}-y_{i+j})^{2}(x_{0}-y_{i})^{2}(x_{0}-y_{i+1})^{2}}\Bigg{]}\,,$
(G.11)

and

$\left[\widehat{\mathcal{B}}_{XX}\Psi_{X}\right](x_{0},\dots,x_{J})=\frac{1}{2}\int\frac{\prod_{i=0}^{J}\differential^{4}y_{i}}{\pi^{2(J+1)}}\frac{\Psi_{X}(y_{0},\dots,y_{J})}{y_{0J}^{2}\prod_{i=0}^{J-1}y_{i,i+1}^{2}}\\\
\times\Bigg{[}\frac{1}{\prod_{i=0}^{J}(x_{i}-y_{i})^{2}}\left(\frac{y_{01}^{2}}{(x_{0}-y_{1})^{2}}+\frac{y_{0J}^{2}}{(x_{0}-y_{J})^{2}}\right)+\sum_{i=1}^{J-1}\frac{y_{i,i+1}^{2}}{\prod_{j=1}^{J}(x_{j}-y_{i+j})^{2}(x_{0}-y_{i})^{2}(x_{0}-y_{i+1})^{2}}\Bigg{]}\,,$
(G.12)

Note that the dependence on $g^{2}$ and $\xi^{2}_{1+1/J}$ has been factored
out as in (3.3) and (3.5). Note also that in the matrix elements involving a
gluon, there are also diagrams where the gluon is absorbed or emitted by one
of the $Z$ propagators. These are, however, corrections of order $g^{2}$ to
$\mathcal{B}_{A\emptyset}$ and $\mathcal{B}_{\emptyset A}$. Hence, they do not
contribute at leading order in $g$. As explained in the main text, even if we
had included them here, the subtraction performed in (3.8) would have removed
them.

As explained around (3.7), the spectrum is deducted from the eigenvalues of
the matrix
$\mathfrak{B}_{ij}=\widehat{\mathcal{B}}_{ij}-\widehat{\mathcal{B}}_{i\emptyset}\widehat{\mathcal{B}}_{\emptyset\emptyset}^{-1}\widehat{\mathcal{B}}_{\emptyset
j}$ for $(i,j)\in\\{A,X\\}^{2}$. We shall now compute these matrix elements
and show that they are intertwined with $\mathfrak{B}_{ij}$ for
$(i,j)\in\\{F,X\\}^{2}$ introduced in section 3 in the following way:

$\begin{pmatrix}\theta\cdot\partial_{x_{0}}&0\\\
0&\mathrm{Id}\end{pmatrix}\begin{pmatrix}\mathfrak{B}_{AA}&\mathfrak{B}_{AX}\\\
\mathfrak{B}_{XA}&\mathfrak{B}_{XX}\end{pmatrix}=\begin{pmatrix}\mathfrak{B}_{FF}&\mathfrak{B}_{FX}\\\
\mathfrak{B}_{XF}&\mathfrak{B}_{XX}\end{pmatrix}\begin{pmatrix}\theta\cdot\partial_{x_{0}}&0\\\
0&\mathrm{Id}\end{pmatrix}\,.$ (G.13)

Using the previous expressions, it is clear that $\mathfrak{B}_{XX}$ is given
by equation (3.13). It is also straightforward to check that

$\theta\cdot\partial_{x_{0}}\,\mathfrak{B}_{AX}\Psi_{X}=\mathfrak{B}_{FX}\Psi_{X}\,.$
(G.14)

One then computes

$\left[\mathfrak{B}_{XA}\Psi_{A}\right](x_{0},\dots,x_{J})=\int\frac{\prod_{i=0}^{J}\differential^{4}y_{i}}{\pi^{2(J+1)}}\frac{\Psi_{A}^{\rho}(y_{0},\dots,y_{J})}{\prod_{i=1}^{J}(x_{i}-y_{i})^{2}\prod_{i=1}^{J-1}y_{i,i+1}^{2}(x_{0}-y_{J})^{2}(x_{0}-y_{1})^{2}}\\\
\times\Bigg{[}\frac{y_{1J}^{2}}{y_{10}^{2}y_{J0}^{2}}\left(\frac{y_{10,\rho}}{y_{10}^{2}}-\frac{y_{J0,\rho}}{y_{J0}^{2}}\right)+\frac{(x_{0}-y_{1})^{2}}{y_{10}^{2}(x_{0}-y_{0})^{2}}\left(\frac{(x_{0}-y_{0})_{\rho}}{(x_{0}-y_{0})^{2}}-\frac{y_{10,\rho}}{y_{10}^{2}}\right)+\frac{(x_{0}-y_{J})^{2}}{y_{J0}^{2}(x_{0}-y_{0})^{2}}\left(\frac{y_{J0,\rho}}{y_{J0}^{2}}-\frac{(x_{0}-y_{0})_{\rho}}{(x_{0}-y_{0})^{2}}\right)\Bigg{]}\,.$
(G.15)

If we define
$\Psi_{F}^{\mu\nu}=\frac{1}{2}\left(\partial_{x_{0}}^{\mu}\Psi_{A}^{\nu}-\partial_{x_{0}}^{\nu}\Psi_{A}^{\mu}\right)$,
i.e. $\Psi_{F}=\theta\cdot\partial_{x_{0}}\Psi_{A}$, and use the integral
formula (2.33), the previous equation can be rewritten

$\mathfrak{B}_{XA}\Psi_{A}=\mathfrak{B}_{XF}\Psi_{F}\,.$ (G.16)

Finally,

$\left[\mathfrak{B}_{AA}\Psi_{A}\right]^{\mu}(x_{0},\dots,x_{J})=4\pi^{2}\int\frac{\differential^{4}z\prod_{i=1}^{J}\differential^{4}y_{i}}{\pi^{2(J+1)}}\frac{\Delta^{\mu}_{\sigma}(x_{0}-z)\Psi_{A}^{\sigma}(z,y_{1},\dots,y_{J})}{\prod_{i=1}^{J}(x_{i}-y_{i})^{2}\prod_{i=1}^{J-1}y_{i,i+1}^{2}(z-y_{1})^{2}(z-y_{J})^{2}}\\\
+\pi^{2}\int\frac{\differential^{4}z\prod_{i=0}^{J}\differential^{4}y_{i}}{\pi^{2(J+2)}}\frac{\Delta^{\mu}_{\sigma}(x_{0}-z)\Psi_{A}^{\rho}(y_{0},\dots,y_{J})}{\prod_{j=1}^{J}(x_{j}-y_{j})^{2}\prod_{j=1}^{J-1}y_{j,j+1}^{2}}\Bigg{[}\frac{1}{(y_{J}-z)^{2}}\frac{\partial}{\partial
z_{\sigma}}\left(\frac{y_{10,\rho}}{y_{10}^{4}(z-y_{0})^{2}}-\frac{(z-y_{0})_{\rho}}{(z-y_{0})^{4}y_{10}^{2}}\right)\\\
-\left(\frac{y_{10,\rho}}{y_{10}^{4}(z-y_{0})^{2}}-\frac{(z-y_{0})_{\rho}}{(z-y_{0})^{4}y_{10}^{2}}\right)\frac{\partial}{\partial
z_{\sigma}}\frac{1}{(y_{J}-z)^{2}}+\frac{1}{(y_{1}-z)^{2}}\frac{\partial}{\partial
z_{\sigma}}\left(\frac{y_{J0,\rho}}{y_{J0}^{4}(z-y_{0})^{2}}-\frac{(z-y_{0})_{\rho}}{(z-y_{0})^{4}y_{J0}^{2}}\right)\\\
-\left(\frac{y_{J0,\rho}}{y_{J0}^{4}(z-y_{0})^{2}}-\frac{(z-y_{0})_{\rho}}{(z-y_{0})^{4}y_{J0}^{2}}\right)\frac{\partial}{\partial
z_{\sigma}}\frac{1}{(y_{1}-z)^{2}}\\\
+\frac{2\,y_{1J}^{2}}{(y_{1}-z)^{2}(y_{J}-z)^{2}y_{10}^{2}y_{J0}^{2}}\left(\frac{(y_{1}-z)^{\sigma}}{(y_{1}-z)^{2}}-\frac{(y_{J}-z)^{\sigma}}{(y_{J}-z)^{2}}\right)\left(\frac{y_{J0,\rho}}{y_{J0}^{2}}-\frac{y_{10,\rho}}{y_{10}^{2}}\right)\Bigg{]}\,.$
(G.17)

We first rewrite this expression as

$\left[\mathfrak{B}_{AA}\Psi_{A}\right]^{\mu}(x_{0},\dots,x_{J})=4\pi^{2}\int\frac{\differential^{4}z\prod_{i=1}^{J}\differential^{4}y_{i}}{\pi^{2(J+1)}}\frac{\Delta^{\mu}_{\sigma}(x_{0}-z)\Psi_{A}^{\sigma}(z,y_{1},\dots,y_{J})}{\prod_{i=1}^{J}(x_{i}-y_{i})^{2}\prod_{i=1}^{J-1}y_{i,i+1}^{2}(z-y_{1})^{2}(z-y_{J})^{2}}\\\
+\pi^{2}\int\frac{\Delta^{\mu}_{\sigma}(x_{0}-z)\Psi_{A}^{\rho}(y_{0},\dots,y_{J})}{\prod_{j=1}^{J}(x_{j}-y_{j})^{2}\prod_{j=1}^{J-1}y_{j,j+1}^{2}(z-y_{1})^{2}(z-y_{J})^{2}}\Bigg{[}\frac{\partial}{\partial
z_{\sigma}}\frac{(y_{1}-z)^{2}}{y_{10}^{2}(z-y_{0})^{2}}\left(\frac{y_{10,\rho}}{y_{10}^{2}}-\frac{(z-y_{0})_{\rho}}{(z-y_{0})^{2}}\right)\\\
+\frac{\partial}{\partial
z_{\sigma}}\frac{(y_{J}-z)^{2}}{y_{J0}^{2}(z-y_{0})^{2}}\left(\frac{y_{J0,\rho}}{y_{J0}^{2}}-\frac{(z-y_{0})_{\rho}}{(z-y_{0})^{2}}\right)+2\left(\frac{(y_{1}-z)^{\sigma}}{(y_{1}-z)^{2}}-\frac{(y_{J}-z)^{\sigma}}{(y_{J}-z)^{2}}\right)\bigg{[}\frac{y_{1J}^{2}}{y_{10}^{2}y_{J0}^{2}}\left(\frac{y_{J0,\rho}}{y_{J0}^{2}}-\frac{y_{10,\rho}}{y_{10}^{2}}\right)\\\
+\frac{(z-y_{1})^{2}}{y_{10}^{2}(z-y_{0})^{2}}\left(\frac{y_{10,\rho}}{y_{10}^{2}}-\frac{(z-y_{0})_{\rho}}{(z-y_{0})^{2}}\right)+\frac{(z-y_{J})^{2}}{y_{J0}^{2}(z-y_{0})^{2}}\left(\frac{(z-y_{0})_{\rho}}{(z-y_{0})^{2}}-\frac{y_{J0,\rho}}{y_{J0}^{2}}\right)\bigg{]}\Bigg{]}\frac{\differential^{4}z\prod_{i=0}^{J}\differential^{4}y_{i}}{\pi^{2(J+2)}}\,.$
(G.18)

Applying equations (2.31) and (2.33) allows to write the right-hand side as an
integral over
$\Psi_{F}^{\mu\nu}=\frac{1}{2}\left(\partial_{x_{0}}^{\mu}\Psi_{A}^{\nu}-\partial_{x_{0}}^{\nu}\Psi_{A}^{\mu}\right)$,
the result is

$\left[\mathfrak{B}_{AA}\Psi_{A}\right]^{\mu}(x_{0},\dots,x_{J})=4\pi^{2}\int\frac{\differential^{4}z\prod_{i=0}^{J}\differential^{4}y_{i}}{\pi^{2(J+2)}}\frac{\Delta^{\mu}_{\sigma}(x_{0}-z)\Psi_{F}^{\rho\tau}(y_{0},\dots,y_{J})}{\prod_{i=1}^{J}(x_{i}-y_{i})^{2}\prod_{i=1}^{J-1}y_{i,i+1}^{2}(z-y_{1})^{2}(z-y_{J})^{2}}\\\
\times\Bigg{[}\frac{(z-y_{0})_{\rho}}{(z-y_{0})^{4}}\delta_{\tau}^{\sigma}+\left(\frac{(y_{J}-z)^{\sigma}}{(y_{J}-z)^{2}}-\frac{(y_{1}-z)^{\sigma}}{(y_{1}-z)^{2}}\right)\bigg{(}\frac{(y_{0}-z)_{\rho}y_{0J,\tau}}{(y_{0}-z)^{2}y_{0J}^{2}}+\frac{y_{0J,\rho}y_{01,\tau}}{y_{0J}^{2}y_{01}^{2}}\\\
+\frac{y_{01,\rho}(y_{0}-z)_{\tau}}{y_{01}^{2}(y_{0}-z)^{2}}\bigg{)}+\frac{1}{2}\frac{\partial}{\partial
z_{\sigma}}\frac{(y_{0}-z)_{\rho}}{(y_{0}-z)^{2}}\left(\frac{y_{01,\tau}}{y_{01}^{2}}+\frac{y_{0J,\tau}}{y_{0J}^{2}}\right)\Bigg{]}\,.$
(G.19)

This implies that, if $\Psi_{F}=\theta\cdot\partial_{x_{0}}\Psi_{A}$, then

$\theta\cdot\partial_{x_{0}}\,\mathfrak{B}_{AA}\Psi_{A}=\mathfrak{B}_{FF}\Psi_{F}\,.$
(G.20)

Finally, we state without proof two useful properties:

$\partial_{x_{0}}\cdot\partial_{\theta}\mathfrak{B}_{AX}\Big{|}_{\alpha=1}=\partial_{x_{0}}\cdot\partial_{\theta}\mathfrak{B}_{AA}\Big{|}_{\alpha=1}=0\,,$
(G.21)

where $\alpha$ parameterises the gluon propagator in the $R_{\alpha}$ gauge,
(2.23). These relations are useful when checking that
$\mathfrak{B}^{-1}\mathfrak{B}=\mathrm{Id}$.

## References

  * [1] O. Gürdoğan and V. Kazakov, _New Integrable 4D Quantum Field Theories from Strongly Deformed Planar $\mathcal{N}=$ 4 Supersymmetric Yang-Mills Theory_, _Phys. Rev. Lett._ 117 (2016) 201602 [1512.06704].
  * [2] A. Cavaglia, D. Grabner, N. Gromov and A. Sever, _Colour-twist operators. Part I. Spectrum and wave functions_ , _JHEP_ 06 (2020) 092 [2001.07259].
  * [3] D. Correa, J. Henn, J. Maldacena and A. Sever, _The cusp anomalous dimension at three loops and beyond_ , _JHEP_ 05 (2012) 098 [1203.1019].
  * [4] D. Grabner, N. Gromov, V. Kazakov and G. Korchemsky, _Strongly $\gamma$-Deformed $\mathcal{N}=4$ Supersymmetric Yang-Mills Theory as an Integrable Conformal Field Theory_, _Phys. Rev. Lett._ 120 (2018) 111601 [1711.04786].
  * [5] N. Gromov, V. Kazakov and G. Korchemsky, _Exact Correlation Functions in Conformal Fishnet Theory_ , _JHEP_ 08 (2019) 123 [1808.02688].
  * [6] V. Kazakov, E. Olivucci and M. Preti, _Generalized fishnets and exact four-point correlators in chiral CFT 4_, _JHEP_ 06 (2019) 078 [1901.00011].
  * [7] B. Basso, G. Ferrando, V. Kazakov and D.-l. Zhong, _Thermodynamic Bethe Ansatz for Biscalar Conformal Field Theories in Any Dimension_ , _Phys. Rev. Lett._ 125 (2020) 091601 [1911.10213].
  * [8] J. Caetano, Ö. Gürdoğan and V. Kazakov, _Chiral limit of $\mathcal{N}$ = 4 SYM and ABJM and integrable Feynman graphs_, _JHEP_ 03 (2018) 077 [1612.05895].
  * [9] B. Basso and L. J. Dixon, _Gluing Ladder Feynman Diagrams into Fishnets_ , _Phys. Rev. Lett._ 119 (2017) 071601 [1705.03545].
  * [10] O. Mamroud and G. Torrents, _RG stability of integrable fishnet models_ , _JHEP_ 06 (2017) 012 [1703.04152].
  * [11] B. Basso and D.-l. Zhong, _Continuum limit of fishnet graphs and AdS sigma model_ , _JHEP_ 01 (2019) 002 [1806.04105].
  * [12] G. P. Korchemsky, _Exact scattering amplitudes in conformal fishnet theory_ , _JHEP_ 08 (2019) 028 [1812.06997].
  * [13] B. Basso, J. Caetano and T. Fleury, _Hexagons and correlators in the fishnet theory_ , _JHEP_ 11 (2019) 172 [1812.09794].
  * [14] G. K. Karananas, V. Kazakov and M. Shaposhnikov, _Spontaneous conformal symmetry breaking in fishnet CFT_ , _Phys. Lett. B_ 811 (2020) 135922 [1908.04302].
  * [15] A. C. Ipsen, M. Staudacher and L. Zippelius, _The one-loop spectral problem of strongly twisted $\mathcal{N}$ = 4 Super Yang-Mills theory_, _JHEP_ 04 (2019) 044 [1812.08794].
  * [16] A. Cavaglià, N. Gromov and F. Levkovich-Maslyuk, _Separation of variables in AdS/CFT: functional approach for the fishnet CFT_ , _JHEP_ 06 (2021) 131 [2103.15800].
  * [17] B. Basso, L. J. Dixon, D. A. Kosower, A. Krajenbrink and D.-l. Zhong, _Fishnet four-point integrals: integrable representations and thermodynamic limits_ , _JHEP_ 07 (2021) 168 [2105.10514].
  * [18] O. Shahpo and E. Vescovi, _Correlation functions of determinant operators in conformal fishnet theory_ , _JHEP_ 06 (2022) 070 [2110.09458].
  * [19] E. Olivucci and P. Vieira, _Stampedes I: fishnet OPE and octagon Bootstrap with nonzero bridges_ , _JHEP_ 07 (2022) 017 [2111.12131].
  * [20] N. Gromov, V. Kazakov, G. Korchemsky, S. Negro and G. Sizov, _Integrability of Conformal Fishnet Theory_ , _JHEP_ 01 (2018) 095 [1706.04167].
  * [21] D. Chicherin, V. Kazakov, F. Loebbert, D. Müller and D.-l. Zhong, _Yangian symmetry for bi-scalar loop amplitudes_ , _JHEP_ 05 (2018) 003 [1704.01967].
  * [22] D. Chicherin, V. Kazakov, F. Loebbert, D. Müller and D.-l. Zhong, _Yangian symmetry for fishnet Feynman graphs_ , _Phys. Rev. D_ 96 (2017) 121901 [1708.00007].
  * [23] V. Kazakov, _Quantum Spectral Curve of $\gamma$-Twisted ${\cal N}=4$ SYM Theory and Fishnet CFT_, _Rev. Math. Phys._ 30 (2018) 1840010 [1802.02160].
  * [24] F. Loebbert, D. Müller and H. Münkler, _Yangian Bootstrap for Conformal Feynman Integrals_ , _Phys. Rev. D_ 101 (2020) 066006 [1912.05561].
  * [25] N. Gromov and A. Sever, _Derivation of the Holographic Dual of a Planar Conformal Field Theory in 4D_ , _Phys. Rev. Lett._ 123 (2019) 081602 [1903.10508].
  * [26] N. Gromov and A. Sever, _Quantum fishchain in AdS 5_, _JHEP_ 10 (2019) 085 [1907.01001].
  * [27] N. Gromov and A. Sever, _The holographic dual of strongly $\gamma$-deformed $\mathcal{N}$ = 4 SYM theory: derivation, generalization, integrability and discrete reparametrization symmetry_, _JHEP_ 02 (2020) 035 [1908.10379].
  * [28] N. Gromov, J. Julius and N. Primi, _Open fishchain in N = 4 Supersymmetric Yang-Mills Theory_ , _JHEP_ 07 (2021) 127 [2101.01232].
  * [29] A. B. Zamolodchikov and A. B. Zamolodchikov, _Factorized S-Matrices in Two Dimensions as the Exact Solutions of Certain Relativistic Quantum Field Theory Models_ , _Annals Phys._ 120 (1979) 253.
  * [30] G. F. Cuomo, D. Karateev and P. Kravchuk, _General Bootstrap Equations in 4D CFTs_ , _JHEP_ 01 (2018) 130 [1705.05401].
  * [31] M. Nishida and K. Tamaoka, _Fermions in Geodesic Witten Diagrams_ , _JHEP_ 07 (2018) 149 [1805.00217].
  * [32] K. Tamaoka, _Geodesic Witten diagrams with antisymmetric tensor exchange_ , _Phys. Rev. D_ 96 (2017) 086007 [1707.07934].
  * [33] A. Pittelli and M. Preti, _Integrable fishnet from $\gamma$-deformed $\mathcal{N}=2$ quivers_, _Phys. Lett. B_ 798 (2019) 134971 [1906.03680].
  * [34] J. Fokken, C. Sieg and M. Wilhelm, _Non-conformality of ${{\gamma}_{i}}$-deformed $\mathcal{N}=4$ SYM theory_, _J. Phys. A_ 47 (2014) 455401 [1308.4420].
  * [35] D. Chicherin, S. Derkachov and A. Isaev, _Conformal algebra: R-matrix and star-triangle relation_ , _JHEP_ 04 (2013) 020 [1206.4150].
  * [36] D. Karateev, P. Kravchuk and D. Simmons-Duffin, _Harmonic Analysis and Mean Field Theory_ , _JHEP_ 10 (2019) 217 [1809.05111].
  * [37] V. K. Dobrev, G. Mack, V. B. Petkova, S. G. Petrova and I. T. Todorov, _Harmonic Analysis on the n-Dimensional Lorentz Group and Its Application to Conformal Quantum Field Theory_ , vol. 63 of _Lecture Notes in Physics_. Springer-Verlag Berlin Heidelberg, 1977, 10.1007/BFb0009678.
  * [38] P. P. Kulish, N. Y. Reshetikhin and E. K. Sklyanin, _Yang-Baxter equation and representation theory. I._ , _Lett. Math. Phys._ 5 (1981) 393.
  * [39] P. P. Kulish and E. K. Sklyanin, _Quantum Spectral Transform Method. Recent Developments_ , _Lect. Notes Phys._ 151 (1982) 61.
  * [40] M. S. Costa, J. Penedones, D. Poland and S. Rychkov, _Spinning Conformal Correlators_ , _JHEP_ 11 (2011) 071 [1107.3554].
  * [41] P. A. M. Dirac, _Wave equations in conformal space_ , _Annals Math._ 37 (1936) 429.
  * [42] M. S. Costa and T. Hansen, _Conformal correlators of mixed-symmetry tensors_ , _JHEP_ 02 (2015) 151 [1411.7351].


# Newton’s laws of motion can generate gravity-mediated entanglement

Marta M. Marchese  Naturwissenschaftlich-Technische Fakultät, Universität
Siegen, Walter-Flex-Straße 3, 57068 Siegen, Germany  Martin Plávala
Naturwissenschaftlich-Technische Fakultät, Universität Siegen, Walter-Flex-
Straße 3, 57068 Siegen, Germany  Matthias Kleinmann  Naturwissenschaftlich-
Technische Fakultät, Universität Siegen, Walter-Flex-Straße 3, 57068 Siegen,
Germany  Stefan Nimmrichter  Naturwissenschaftlich-Technische Fakultät,
Universität Siegen, Walter-Flex-Straße 3, 57068 Siegen, Germany

###### Abstract

The interface between quantum theory and gravity represents still uncharted
territory. Recently, some works suggested promising alternative approaches
aimed at witnessing quantum features to test the fundamental nature of gravity
in tabletop experiments: Two masses in an initial superposition of spatially
localized states are allowed to interact only through gravity and it is
measured whether the final state is entangled. Here we show that one can
generate the same amount of entanglement in this setup by using classical time
evolution given by Newton’s laws of motion. We argue that theories of quantum
gravity that can be approximated by the Newtonian potential and classical time
evolution given by Newton’s laws of motion will generate gravity-mediated
entanglement.

## Introduction

Whether gravity is fundamentally quantum or not is still one of the most
debated problems since the lack of experimental evidence does not allow for a
decisive answer in either direction. Gravity-mediated entanglement has
recently been suggested as a testable quantum feature to address this question
[1, 2, 3, 4, 5]. In these schemes, two masses are prepared in a superposition
of spatial states on the two arms of an interferometer, then they are let to
interact through just the gravitational field. If gravity-mediated
entanglement between the two masses is detected at the end of the evolution,
then an argument that local operations and classical communications (LOCC)
cannot generate entanglement [6, 7] is used to conclude that the interaction
cannot be described by classical physics: either gravity is nonlocal, or it is
nonclassical [8]. In this sense the proposed experiments aim to answer the
standing questions whether gravity needs to be quantized or not [9, 10].
Performing the proposed experiment is still outside the realm of possibility
[11], but many variants of the scheme were subsequently proposed [12, 13, 14,
15, 16, 17] and additional theoretical analysis was carried out to further
support this result [18, 19, 20].

Witnessing gravity-mediated entanglement will rule out several alternative
theories such as semi-classical models, that treat gravity as a classical
field acting as a weak continuous measurement on a quantum system [21, 22, 23,
9, 24], and theories as the ones predicting the breaking down of quantum
mechanics at scales where quantum gravity effects should be relevant [25, 26].
Along the same line, some measurement-feedback descriptions have been ruled
out because they are incompatible with experimental data from atom
interferometry experiments [27]. However, there are still classical models
that have not been excluded yet [9, 23]. It has also been argued that
witnessing gravity-induced entanglement does not undoubtedly certify the
nonclassical nature of gravity [28, 29, 30, 31, 32].

Figure 1: Purity marginal as a function of time. $\gamma_{\mathrm{QT}}$ (blue
solid line) is the purity of the marginal of the state evolved according to
the Schrödinger equation; $\gamma_{\mathrm{Swift}}$ (green dotted line) is the
purity of the marginal obtained from the classical time evolution in the phase
space of the Wigner function by using the second-order Taylor expansion of the
potential; $\gamma_{\mathrm{Fit}}$ (red dashed line) is the purity for the
marginal obtained from the classical time evolution with the second-order
polynomial that coincides with the Newtonian potential for
$x_{2}-x_{1}=d\pm\Delta x$ and $x_{2}-x_{1}=d$. Here $x_{1}$ and $x_{2}$ are
the position coordinates of the two particles, $d$ is the distance of centers
of mass of the superpositions and $\Delta x$ is the arm separation. We are
using parameters as proposed in [1], i.e., $m=10^{-14}\text{kg}$, $\Delta
x=250\mu\text{m}$, $d=450\mu\text{m}$, and $\sigma=10\mu\text{m}$. The solid
black line indicates the proposed time $\tau=2.5\text{s}$ [1] over which the
phases $\varphi_{1}=-0.13$ and $\varphi_{2}=0.44$ are accumulated.

Here we show that it is possible to reproduce the gravity-mediated
entanglement between the two simultaneously interfering masses by using a
classical time evolution based on Newton’s laws of motion and a suitable
approximation of the gravitational potential. The amount of entanglement,
quantified using the purity of the marginals, matches the quantum prediction
based on the Schrödinger equation, see Fig. 1. This has two implications:
Firstly, any theory of gravity, for which in the regime of tabletop
experiments Newton’s laws of motion are the leading term of the time evolution
in phase space, is capable of generating entanglement from spatial
superpositions. Secondly, the aforementioned experimental test schemes will
not provide a definitive answer about the quantum nature of gravity, since
they do not rule out a simple classical description.

In the following, we treat the two gravitationally interacting masses in the
nonrelativistic Newtonian limit. In order to calculate the time evolution of a
quantum state according to Newton’s laws of motion, we employ the Wigner
function representation [33] and the Liouville equation. We also use a second-
order polynomial approximation of the interaction potential, which yields an
analytically tractable and physical solution of the masses’ state that
exhibits entanglement without making any assumption on the quantum nature of
gravity. Other approximations of the potential may result in unphysical
dynamics of the Wigner function at large times. For example, a stepwise linear
approximation, subjecting different branches of the superposition to different
gravitational forces, leads to small but noticeable negative probabilities at
$t\gtrsim 1.6\;\text{s}$ for the settings of Fig. 1; see App. B. This could be
remedied by an adequate amount of diffusion noise, which would arise in a
consistent hybrid model of the dynamics of quantum particles interaction with
a classical gravitational field [9, 23].

In the following, we briefly review the interferometric setup and the
wavefunction description of the proposed quantum gravity test schemes and
introduce the purity-based entanglement measure we employ. We then present our
classical phase-space model for gravity-mediated entanglement underlying the
results plotted in Fig. 1 and conclude with an outlook on the future of such
test experiments.

## Experimental scenario

Consider a Mach-Zehnder interferometer setup with two neutral identical masses
$m$ placed at a distance $d$. Each of them is split into a superposition of
spatial states $\ket{L}$ and $\ket{R}$ representing the two arms of the
interferometer with a separation $\Delta x<d$. Each single-arm state is given
by a localized Gaussian wave packet of width $\sigma\ll\Delta x$. Hence, the
two masses are initially in a normalized product state
$\psi(x_{1},x_{2})=\psi_{1}(x_{1})\psi_{2}(x_{2})$ of two spatially
superimposed Gaussian wave packets,

$\psi_{i}(q_{i})=\dfrac{e^{-(q_{i}-\Delta
x/2)^{2}/4\sigma^{2}}+e^{-(q_{i}+\Delta
x/2)^{2}/4\sigma^{2}}}{(2\pi\sigma^{2}N)^{1/4}},$ (1)

where $q_{1}=x_{1}-d/2$ and $q_{2}=x_{2}+d/2$ are the shifted one-dimensional
coordinates of the two masses and $N=4(1+e^{-\Delta x^{2}/8\sigma^{2}})^{2}$
is the normalization constant. We shall omit the other two spatial directions
throughout, assuming that both masses are sufficiently localized in $y,z$ at
all times. We assume that the distance between the two masses is always
sufficiently large so that we can neglect the Casimir-Polder force, and the
only force acting between the two masses is due to their mutual gravitational
attraction. The Hamiltonian for the two masses is then
$H(x_{1},x_{2},p_{1},p_{2})=(p_{1}^{2}+p_{2}^{2})/2m-\kappa/\absolutevalue{x_{1}-x_{2}}$,
with $\kappa=Gm^{2}$. We follow the calculations presented in [1]: The two
Gaussian wave packets now correspond to the left and right arms of an
interferometer represented by the states $\ket{L}$ and $\ket{R}$. Then, using
the full Newtonian potential and neglecting the kinetic part of the
Hamiltonian, the time evolution due purely to gravitational interaction leads
to the state

$\begin{split}\ket{\psi(t)}=\dfrac{e^{i\phi}}{\sqrt{N(t)}}&(\ket{LL}+e^{i\Delta\phi_{LR}t}\ket{LR}\\\
&+e^{i\Delta\phi_{RL}t}\ket{RL}+\ket{RR}),\end{split}$ (2)

where $N(t)=4+4(\cos(\Delta\phi_{LR}t)+\cos(\Delta\phi_{RL}t))e^{-\frac{\Delta
x^{2}}{8\sigma^{2}}})+2\left[1+\cos((\Delta\phi_{LR}-\Delta\phi_{RL})t)\right]e^{-\frac{\Delta
x^{2}}{4\sigma^{2}}}$ is the normalization factor; it depends on the phase
shifts because the states are not perfectly orthogonal. For $\Delta
x\gg\sigma$ we get $N(t)=4$ which we will use from now on. The phases
$\Delta\phi_{LR}=\phi_{LR}-\phi$, $\Delta\phi_{RL}=\phi_{RL}-\phi$ read
$\phi_{LR}=\frac{Gm^{2}}{\hbar(d+\Delta x)}$,
$\phi_{RL}=\frac{Gm^{2}}{\hbar(d-\Delta x)}$, $\phi=\frac{Gm^{2}}{\hbar d}$.
The final state is still pure, since the evolution is unitary. In order to
check if it is still separable or not we look at the purity of marginals. In
particular, the marginal state of particle $2$ is given by tracing out
particle 1, $\rho_{2}(t)=\Tr_{1}[\rho(t)]$, where
$\rho(t)=\outerproduct{\psi(t)}{\psi(t)}$ is the density matrix of the evolved
state from Eq. (2). The reduces states density matrix results in

$\rho_{2}(t)=\dfrac{1}{4}\begin{pmatrix}2&e^{-i\Delta\phi_{LR}t}-e^{i\Delta\phi_{RL}t}\\\
e^{i\Delta\phi_{LR}t}-e^{-i\Delta\phi_{RL}t}&2\end{pmatrix}$ (3)

and it is straightforward to obtain the purity

$\gamma_{\mathrm{QT}}(t)=\Tr[\rho_{2}^{2}(t)]=\dfrac{3+\cos((\Delta\phi_{LR}+\Delta\phi_{RL})t)}{4}.$
(4)

If a quantum state represented by a density matrix $\rho$ is pure, a necessary
and sufficient condition for separability is that all the marginals must also
be pure. Therefore, since the purity of the marginal is found smaller than
$1$, then this implies that $\ket{\psi(t)}$ must be entangled.

There are two ways in which non-classical effects enter the previous
calculations: the initial states are assumed to be superpositions of paths and
we use the Schödinger equation to evolve the state in time. If we changed the
initial states to statistical mixtures instead of superpositions then
entanglement would not be generated, hence the non-classical nature of the
initial state is crucial for obtaining gravity-mediated entanglement. The
question we aim to answer is: what would happen if we changed the quantum time
evolution given by Schödinger equation to a classical time evolution given by
Newton’s laws of motion? This is a relevant question since if one can generate
entanglement using a classical time evolution, then the only necessary non-
classical resource are the initial superpositions.

## Classical time evolution

In phase space, Newton’s laws of motion correspond to the Liouville equation
$\dot{P}=\\{H,P\\}$, where $\\{\cdot,\cdot\\}$ is the Poisson bracket,
$\\{f,g\\}=\frac{\partial f}{\partial q}\frac{\partial g}{\partial
p}-\frac{\partial f}{\partial p}\frac{\partial g}{\partial q}$, $H$ is the
Hamiltonian, and $P$ is the probability distribution of the system in phase
space. In order to apply Newton’s laws of motion to a quantum system, we
simply replace the probability distribution $P$ by its closest quantum
mechanical analog, the Wigner function $W$ [34, 35]. The Wigner function for a
single particle in a pure state is

$W(q,p)=\dfrac{1}{2\pi\hbar}\int
e^{-ips/\hbar}\psi^{*}\left(q-\frac{s}{2}\right)\psi\left(q+\frac{s}{2}\right)\differential
s.$ (5)

and we get the equation

$\dot{W}=\\{H,W\\}$ (6)

where $H$ is still the same Hamiltonian function in phase space. As it is the
case for the Liouville equation, also Eq. (6) is solved by
$W(x_{1},x_{2},p_{1},p_{2},t)=W^{0}(X_{1}(t),X_{2}(t),P_{1}(t),P_{2}(t))$,
where $X_{i}(t)=X_{i}(x_{1},x_{2},p_{1},p_{2},t)$ and
$P_{i}(t)=X_{i}(x_{1},x_{2},p_{1},p_{2},t)$ are the coordinates of the point
in phase space at which pair of classical particles would have to start at
$t=0$ to reach the point $(x_{1},x_{2},p_{1},p_{2})$ at time $t$. It now only
remains to solve the respective Hamilton equations of motion.

At the initial time $t=0$, the total Wigner function for the two-particle
state will be given by the product

$W^{0}(x_{1},x_{2},p_{1},p_{2})=W_{1}^{0}(x_{1},p_{1})W^{0}_{2}(x_{2},p_{2}),$
(7)

where we used Eq. (1) to obtain the explicit expression for each of the two
masses,

$\begin{split}W_{i}^{0}(q_{i},p_{i})=\;&\dfrac{e^{\frac{-2\sigma^{2}p_{i}^{2}}{\hbar^{2}}}}{\pi\hbar\sqrt{N}}\Bigg{[}e^{-\frac{(q_{i}+\frac{\Delta
x}{2})^{2}}{2\sigma^{2}}}+e^{-\frac{(q_{i}-\frac{\Delta
x}{2})^{2}}{2\sigma^{2}}}\\\
&+2e^{-\frac{q_{i}}{2\sigma^{2}}}\cos(\frac{\Delta
x\;p_{i}}{\hbar})\Bigg{]}.\end{split}$ (8)

see Appendix A for detailed calculations.

In order to make the problem tractable, we will approximate the Newtonian
gravitational potential by a second-order polynomial. This ensures that the
time-evolution yields a pure quantum state at all times [33]. To do so, we
first apply a symplectic transformation that preserves both the Poisson
bracket and the commutator. The new coordinates are

$\begin{split}x_{\mathrm{avg}}=\dfrac{1}{\sqrt{2}}(x_{1}+x_{2}),\qquad&p_{\mathrm{avg}}=\dfrac{1}{\sqrt{2}}(p_{1}+p_{2}),\\\
x_{\mathrm{rel}}=\dfrac{1}{\sqrt{2}}(x_{2}-x_{1}-d),\qquad&p_{\mathrm{rel}}=\dfrac{1}{\sqrt{2}}(p_{2}-p_{1}),\end{split}$
(9)

and the Hamiltonian is

$H(x_{\mathrm{avg}},x_{\mathrm{rel}},p_{\mathrm{avg}},p_{\mathrm{rel}})=\dfrac{p_{\mathrm{avg}}^{2}+p_{\mathrm{rel}}^{2}}{2m}-\dfrac{\kappa}{d+\sqrt{2}x_{\mathrm{rel}}}.$
(10)

Here, we assumed that during the whole experiment, the Wigner function
approximately vanishes at $x_{1}\geq 0$ or at $x_{2}\leq 0$. Therefore,
$d+\sqrt{2}x_{\mathrm{rel}}\geq 0$ and so we dropped the absolute value.

The Taylor expansion of the potential up to second order in $x_{\mathrm{rel}}$
is

$V_{\mathrm{Swift}}(x_{\mathrm{rel}})=-\dfrac{\kappa}{d}+\dfrac{\sqrt{2}\kappa}{d^{2}}x_{\mathrm{rel}}-\dfrac{2\kappa}{d^{3}}x_{\mathrm{rel}}^{2},$
(11)

Alternatively, we will consider a quadratic fit of the Newtonian potential
such that the fitting polynomial coincides with the potential at the most
relevant points $x_{\mathrm{rel}}=\pm\Delta x/\sqrt{2}$ and
$x_{\mathrm{rel}}=0$,

$V_{\mathrm{Fit}}(x_{\mathrm{rel}})=-\dfrac{\kappa}{d}+\dfrac{\sqrt{2}\kappa}{d^{2}-\Delta
x^{2}}x_{\mathrm{rel}}-\dfrac{2\kappa}{d(d^{2}-\Delta
x^{2})}x_{\mathrm{rel}}^{2}.$ (12)

Both of the quadratic potentials $V_{\mathrm{Swift}}$ and $V_{\mathrm{Fit}}$
give reasonable approximations to the Newtonian potential and the respective
force, see Appendix C. In both cases, the respective Hamilton equations of
motion can be easily solved, see Appendix D. The two approximations will in
general lead to different time evolutions, but note that for $d\gg\Delta x$ we
have $V_{\mathrm{Swift}}(x_{\mathrm{rel}})\approx
V_{\mathrm{Fit}}(x_{\mathrm{rel}})$. We will use $W^{\mathrm{Swift}}$ and
$W^{\mathrm{Fit}}$ to denote the Wigner functions evolved according to
Hamiltonians with the approximate potentials $V_{\mathrm{Swift}}$ and
$V_{\mathrm{Fit}}$, respectively.

We will again use purity of the marginals to quantify entanglement of the
states; we use this method since it is easy to compute from the Wigner
function [36]. We want to compute the quantities

$\displaystyle\gamma_{\mathrm{Swift}}(t)$
$\displaystyle=2\pi\hbar\int(W^{\mathrm{Swift}}_{2}(x_{2},p_{2},t))^{2}\differential
x_{2}\differential p_{2}$ (13) $\displaystyle\gamma_{\mathrm{Fit}}(t)$
$\displaystyle=2\pi\hbar\int(W^{\mathrm{Fit}}_{2}(x_{2},p_{2},t))^{2}\differential
x_{2}\differential p_{2}$ (14)

where

$\displaystyle W^{\mathrm{Swift}}_{2}(x_{2},p_{2},t)$ $\displaystyle=\int
W^{\mathrm{Swift}}(x_{1},x_{2},p_{1},p_{2},t)\differential x_{1}\differential
p_{1}$ (15) $\displaystyle W^{\mathrm{Fit}}_{2}(x_{2},p_{2},t)$
$\displaystyle=\int W^{\mathrm{Fit}}(x_{1},x_{2},p_{1},p_{2},t)\differential
x_{1}\differential p_{1}$ (16)

are the respective marginals. Since we have already solved the time evolution,
we can compute the purities numerically.

In Fig. 1 we show the purity of the reduced state as a function of time for
the three different time evolutions: quantum dynamics given by Schrödinger
equation and classical dynamics given by Newton’s laws of motion and the
Taylor approximated potential or the second-order polynomial fit of the
potential. The state is pure at the initial time $t=0$, since the total state
was just the product of the two particle states. After a certain time, all the
curves reach the minimum value of purity $\frac{1}{2}$ meaning that at this
time the entanglement between the two subsystems is maximal. The classical
time evolution with the Taylor expansion of the potential has a slower rate of
entanglement, since the minimum is reached between $30\unit{s}$ and
$40\unit{s}$, compared to the quantum evolution and classical evolution with
the polynomial fit of the potential, where the minimum is reached for
$t=\frac{\pi}{\Delta\phi_{LR}+\Delta\phi_{RL}}\approx 25\unit{s}$. The
difference between $\gamma_{\mathrm{QT}}$ and $\gamma_{\mathrm{Fit}}$ for
$t\geq 40\text{s}$ is likely due to an addition approximation underlying the
quantum model in Eq. (2): The deflection and dispersion of the wavepackets is
neglected, contrary to the classical dynamics that we employ.

## Conclusions

In this paper we have shown that gravity-mediated entanglement can be
generated by using Newton’s laws of motion, here represented via Liouville
equation using a suitable approximation of the Newtonian potential. This shows
that observing gravity-mediated entanglement does not certify the non-
classicality of the gravitational interaction, in the sense that initial
superpositions and classical time evolution are already sufficien. New
experiments should be designed that will certify non-classicality or
quantumness of the gravitational interaction by ruling out simple classical
models, such as the one we presented. For example, one could consider
different arm separations $\Delta x$ for the two particles and quantitatively
analyze the entanglement rate in the setup of Ref. [1]. A even more
conclusive, but much more challenging, version of the experiment could use the
gravitational interaction to create the initial superpositions.

In order to tackle the problem more systematically, an approach based on
operational theories [35, 37] similar to the approach to violations of Bell
inequalities can be adopted: This approach was recently used to prove that
quantum theory based on real numbers can be experimentally ruled out [38], to
analyze quantumness of uncertainty relations [39] and interference phenomena
[40, 41, 42]. A similar approach to quantum gravity may be able to rule out
large classes of classical description and settle the question of whether
gravity is non-classical. We have already made first step in this direction in
this work by showing that simple classical models can describe current
experimental setups and offer different explanations of the experimental
results.

A potential criticism of our approach is that we are using the Newtonian
potential to model the gravitational interaction which is by definition
nonlocal. This is the standard approach in description of the proposed
experiments [1, 2, 3, 4, 5] and one can easily make the interaction non-
signaling by considering retarded potentials, although retardation will not
play any significant role at the distances of the order $10^{-12}$ light-
seconds considered within current experimental proposals. Moreover, using
essentially the contrapositive of the LOCC argument, the observation of
gravity-mediated entanglement should prove that the gravitational field was at
some point entangled with the particles, but this is impossible since the
state of the particles is always pure even if one describes the experiment
using a first-principles approach and path integrals [19]. This is a clear
contradiction since monogamy of entanglement states that a pure state cannot
be entangled with any other state, even if one would consider the
gravitational field to be described by an operational theory potentially
different from quantum theory [43]. A possible resolution of this conundrum is
that, in the field theoretic approach in the Lorentz gauge where all
interactions are local, the entanglement is mediated via auxiliary,
unobservable mediators that do not correspond to physical degrees of freedom
themselves [30]. Nevertheless, as also pointed out in Ref. [30], one can
always describe the experiment in the Poisson gauge using the non-local
potentials, at which point, as we showed, one can obtain gravity-mediated
entanglement even using only the classical time evolution given by Newton’s
laws of motion.

## Acknowledgments

We acknowledge support from the Deutsche Forschungsgemeinschaft (DFG, German
Research Foundation, project numbers 447948357 and 440958198), the Sino-German
Center for Research Promotion (Project M-0294), the German Ministry of
Education and Research (Project QuKuK, BMBF Grant No. 16KIS1618K), the DAAD.
MP acknowledges support from the Alexander von Humboldt Foundation. MMM
acknowledges support from the Walter Benjamin Programme.

## References

  * Bose  _et al._ [2017] S. Bose, A. Mazumdar, G. W. Morley, H. Ulbricht, M. Toroš, M. Paternostro, A. A. Geraci, P. F. Barker, M. Kim, and G. Milburn, Spin entanglement witness for quantum gravity, Physical review letters 119, 240401 (2017).
  * Marletto and Vedral [2017] 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).
  * Al Balushi  _et al._ [2018] A. Al Balushi, W. Cong, and R. B. Mann, Optomechanical quantum cavendish experiment, Physical Review A 98, 043811 (2018).
  * Krisnanda  _et al._ [2020] T. Krisnanda, G. Y. Tham, M. Paternostro, and T. Paterek, Observable quantum entanglement due to gravity, npj Quantum Information 6, 12 (2020).
  * Weiss  _et al._ [2021] T. Weiss, M. Roda-Llordes, E. Torrontegui, M. Aspelmeyer, and O. Romero-Isart, Large quantum delocalization of a levitated nanoparticle using optimal control: Applications for force sensing and entangling via weak forces, Physical Review Letters 127, 023601 (2021).
  * Horodecki  _et al._ [2009] R. Horodecki, P. Horodecki, M. Horodecki, and K. Horodecki, Quantum entanglement, Reviews of modern physics 81, 865 (2009).
  * Gühne and Tóth [2009] O. Gühne and G. Tóth, Entanglement detection, Physics Reports 474, 1 (2009).
  * Galley  _et al._ [2022] T. D. Galley, F. Giacomini, and J. H. Selby, A no-go theorem on the nature of the gravitational field beyond quantum theory, Quantum 6, 779 (2022).
  * Tilloy [2019] A. Tilloy, Does gravity have to be quantized? lessons from non-relativistic toy models, in  _Journal of Physics: Conference Series_ , Vol. 1275 (IOP Publishing, 2019) p. 012006.
  * Oppenheim [2023a] J. Oppenheim, Is it time to rethink quantum gravity? (2023a).
  * Aspelmeyer [2022] M. Aspelmeyer, When zeh meets feynman: How to avoid the appearance of a classical world in gravity experiments, in  _From Quantum to Classical: Essays in Honour of H.-Dieter Zeh_ (Springer, 2022) pp. 85–95.
  * Nguyen and Bernards [2020] H. C. Nguyen and F. Bernards, Entanglement dynamics of two mesoscopic objects with gravitational interaction, The European Physical Journal D 74, 1 (2020).
  * Lami  _et al._ [2023] L. Lami, J. S. Pedernales, and M. B. Plenio, Testing the quantumness of gravity without entanglement (2023).
  * Howl  _et al._ [2023] R. Howl, N. Cooper, and L. Hackermüller, Gravitationally-induced entanglement in cold atoms (2023).
  * Hanif  _et al._ [2023] F. Hanif, D. Das, J. Halliwell, D. Home, A. Mazumdar, H. Ulbricht, and S. Bose, Testing whether gravity acts as a quantum entity when measured (2023).
  * Pesci and Pieri [2023] A. Pesci and P. Pieri, Testing the nonclassicality of gravity with the field of a single delocalized mass, Physical Review A 108, 062801 (2023).
  * Pedernales  _et al._ [2022] 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).
  * Danielson  _et al._ [2022] D. L. Danielson, G. Satishchandran, and R. M. Wald, Gravitationally mediated entanglement: Newtonian field versus gravitons, Physical Review D 105, 086001 (2022).
  * Christodoulou  _et al._ [2023] M. Christodoulou, A. Di Biagio, M. Aspelmeyer, Č. Brukner, C. Rovelli, and R. Howl, Locally mediated entanglement in linearized quantum gravity, Physical Review Letters 130, 100202 (2023).
  * Feng  _et al._ [2023] T. Feng, C. Marletto, and V. Vedral, Conservation laws reveal the quantumness of gravity (2023).
  * Kafri  _et al._ [2014] D. Kafri, J. Taylor, and G. Milburn, A classical channel model for gravitational decoherence, New Journal of Physics 16, 065020 (2014).
  * Khosla and Nimmrichter [2018] K. E. Khosla and S. Nimmrichter, Classical Channel Gravity in the Newtonian Limit, arXiv: 1812.03118 10.48550/arXiv.1812.03118 (2018).
  * Oppenheim [2023b] J. Oppenheim, A postquantum theory of classical gravity?, Physical Review X 13, 041040 (2023b).
  * Carney and Taylor [2023] D. Carney and J. M. Taylor, Strongly incoherent gravity, arXiv: 2301.08378 10.48550/arXiv.2301.08378 (2023), arXiv:2301.08378 [quant-ph] .
  * Diósi [1989] L. Diósi, Models for universal reduction of macroscopic quantum fluctuations, Physical Review A 40, 1165 (1989).
  * Tilloy and Diósi [2016] A. Tilloy and L. Diósi, Sourcing semiclassical gravity from spontaneously localized quantum matter, Physical Review D 93, 024026 (2016).
  * Altamirano  _et al._ [2018] N. Altamirano, P. Corona-Ugalde, R. B. Mann, and M. Zych, Gravity is not a pairwise local classical channel, Classical and Quantum Gravity 35, 145005 (2018).
  * Hall and Reginatto [2018] M. J. Hall and M. Reginatto, On two recent proposals for witnessing nonclassical gravity, Journal of Physics A: Mathematical and Theoretical 51, 085303 (2018).
  * Anastopoulos  _et al._ [2021] C. Anastopoulos, M. Lagouvardos, and K. Savvidou, Gravitational effects in macroscopic quantum systems: a first-principles analysis, Classical and Quantum Gravity 38, 155012 (2021).
  * Fragkos  _et al._ [2022] V. Fragkos, M. Kopp, and I. Pikovski, On inference of quantization from gravitationally induced entanglement, AVS Quantum Science 4 (2022).
  * Döner and Großardt [2022] M. K. Döner and A. Großardt, Is gravitational entanglement evidence for the quantization of spacetime?, Foundations of Physics 52, 101 (2022).
  * Spaventa  _et al._ [2023] G. Spaventa, L. Lami, and M. B. Plenio, On tests of the quantum nature of gravitational interactions in presence of non-linear corrections to quantum mechanics, arXiv preprint arXiv:2302.00365 (2023).
  * Case [2008] W. B. Case, Wigner functions and weyl transforms for pedestrians, American Journal of Physics 76, 937 (2008).
  * Schleich [2011] W. P. Schleich,  _Quantum optics in phase space_ (John Wiley & Sons, 2011).
  * Plávala and Kleinmann [2022] M. Plávala and M. Kleinmann, Operational theories in phase space: Toy model for the harmonic oscillator, Physical Review Letters 128, 040405 (2022).
  * Almeida [2009] A. O. d. Almeida, Entanglement in phase space, in  _Entanglement and Decoherence: Foundations and Modern Trends_ (Springer, 2009) pp. 157–219.
  * Plávala and Kleinmann [2023] M. Plávala and M. Kleinmann, Generalized dynamical theories in phase space and the hydrogen atom, Physical Review A 108, 052212 (2023).
  * Renou  _et al._ [2021] M.-O. Renou, D. Trillo, M. Weilenmann, T. P. Le, A. Tavakoli, N. Gisin, A. Acín, and M. Navascués, Quantum theory based on real numbers can be experimentally falsified, Nature 600, 625 (2021).
  * Catani  _et al._ [2022] L. Catani, M. Leifer, G. Scala, D. Schmid, and R. W. Spekkens, What is nonclassical about uncertainty relations?, Physical Review Letters 129, 240401 (2022).
  * Blasiak [2015] P. Blasiak, Local model of a qubit in the interferometric setup, New Journal of Physics 17, 113043 (2015).
  * Catani  _et al._ [2023a] L. Catani, M. Leifer, G. Scala, D. Schmid, and R. W. Spekkens, Aspects of the phenomenology of interference that are genuinely nonclassical, Physical Review A 108, 022207 (2023a).
  * Catani  _et al._ [2023b] L. Catani, M. Leifer, D. Schmid, and R. W. Spekkens, Why interference phenomena do not capture the essence of quantum theory, Quantum 7, 1119 (2023b).
  * Plávala [2023] M. Plávala, General probabilistic theories: An introduction, Physics Reports 1033, 1 (2023).

##  Appendix A Wigner function for the initial state

The Wigner function for a single-particle quantum state described by a
wavefunction $\psi(q)$ is

$W(q,p)=\dfrac{1}{2\pi\hbar}\int
e^{-i\frac{ps}{\hbar}}\psi^{*}(q-\frac{s}{2})\psi(q+\frac{s}{2})\differential
s.$ (17)

Then if we have Gaussian wavepackets of the form

$\psi(q)=\left(\dfrac{1}{2\pi\sigma^{2}N}\right)^{1/4}\left(e^{-\frac{1}{4\sigma^{2}}\left(q+\frac{\Delta
x}{2}\right)^{2}}+e^{-\frac{1}{4\sigma^{2}}\left(q-\frac{\Delta
x}{2}\right)^{2}}\right),$ (18)

we obtain

$\begin{split}W(q,p)=\dfrac{1}{2\pi\hbar}\left(\dfrac{1}{2\pi\sigma^{2}N}\right)^{1/2}\int
e^{-i\frac{ps}{\hbar}}&\Bigg{[}e^{-\frac{1}{4\sigma^{2}}\left(q+\frac{\Delta
x-s}{2}\right)^{2}-\frac{1}{4\sigma^{2}}\left(q+\frac{\Delta
x+s}{2}\right)^{2}}\\\ &+e^{-\frac{1}{4\sigma^{2}}\left(q+\frac{\Delta
x-s}{2}\right)^{2}-\frac{1}{4\sigma^{2}}\left(q-\frac{\Delta
x-s}{2}\right)^{2}}\\\ &+e^{-\frac{1}{4\sigma^{2}}\left(q-\frac{\Delta
x-s}{2}\right)^{2}-\frac{1}{4\sigma^{2}}\left(q+\frac{\Delta
x+s}{2}\right)^{2}}\\\ &+e^{-\frac{1}{4\sigma^{2}}\left(q-\frac{\Delta
x+s}{2}\right)^{2}-\frac{1}{4\sigma^{2}}\left(q-\frac{\Delta
x-s}{2}\right)^{2}}\Bigg{]}\differential s.\end{split}$ (19)

The integral will yield a real function in the end, such that the result for
one particle is

$W(q,p)=\dfrac{1}{\pi\hbar\sqrt{N}}e^{\frac{-2\sigma^{2}}{\hbar^{2}}p^{2}}\left[e^{-\frac{(q+\Delta
x/2)^{2}}{2\sigma^{2}}}+e^{-\frac{(q-\Delta
x/2)^{2}}{2\sigma^{2}}}+2e^{-\frac{q^{2}}{2\sigma^{2}}}\cos\left(\frac{\Delta
x\;p}{\hbar}\right)\right].$ (20)

For our two-particle state, the total Wigner function is just the product of
single-particle Wigner functions,
$W^{0}(x_{1},x_{2},p_{1},p_{2})=W^{0}_{1}(x_{1},p_{1})W^{0}_{2}(x_{2},p_{2})$.
The final expression is obtained by substituting in Eq. (20), $q=x_{1}+d/2$
for particle 1 and $q=x_{2}-d/2$ for particle 2, respectively,

$\begin{split}W^{0}(x_{1},x_{2},p_{1},p_{2})=\left(\dfrac{1}{\pi\hbar\sqrt{N}}\right)^{2}e^{-\frac{2\sigma^{2}}{\hbar^{2}}(p_{1}^{2}+p_{2}^{2})^{2}}&\left[e^{-\frac{(x_{1}+\frac{d+\Delta
x}{2})^{2}}{2\sigma^{2}}}+e^{-\frac{(x_{1}+\frac{d-\Delta
x}{2})^{2}}{2\sigma^{2}}}+2e^{-\frac{(x_{1}+\frac{d}{2})^{2}}{2\sigma^{2}}}\cos\left(\frac{\Delta
xp}{\hbar}\right)\right]\\\ &\left[e^{-\frac{(x_{2}-\frac{d-\Delta
x}{2})^{2}}{2\sigma^{2}}}+e^{-\frac{(x_{2}-\frac{d+\Delta
x}{2})^{2}}{2\sigma^{2}}}+2e^{-\frac{(x_{2}-\frac{d}{2})^{2}}{2\sigma^{2}}}\cos\left(\frac{\Delta
xp}{\hbar}\right)\right].\end{split}$ (21)

The expansion of this product will result in nine contributions that will
depend on different distances,

$W^{0}_{\text{TOT}}=\sum_{j=1}^{9}w_{j}(x_{1},x_{2},p_{1},p_{2}),$ (22)

omitting the arguments $(x_{1},x_{2},p_{1},p_{2})$ for readability. The
individual terms are

$\begin{split}&w_{1}=\dfrac{1}{N}\left(\dfrac{1}{\pi\hbar}\right)^{2}e^{-\frac{2\sigma^{2}}{\hbar^{2}}(p_{1}^{2}+p_{2}^{2})}e^{-\frac{1}{2\sigma^{2}}(x_{1}+\frac{d+\Delta
x}{2})^{2}}e^{-\frac{1}{2\sigma^{2}}(x_{2}-\frac{d-\Delta x}{2})^{2}}\\\
&w_{2}=\dfrac{1}{N}\left(\dfrac{1}{\pi\hbar}\right)^{2}e^{-\frac{2\sigma^{2}}{\hbar^{2}}(p_{1}^{2}+p_{2}^{2})}e^{-\frac{1}{2\sigma^{2}}(x_{1}+\frac{d+\Delta
x}{2})^{2}}e^{-\frac{1}{2\sigma^{2}}(x_{2}-\frac{d+\Delta x}{2})^{2}}\\\
&w_{3}=\dfrac{1}{N}\left(\dfrac{1}{\pi\hbar}\right)^{2}e^{-\frac{2\sigma^{2}}{\hbar^{2}}(p_{1}^{2}+p_{2}^{2})}e^{-\frac{1}{2\sigma^{2}}(x_{1}+\frac{d+\Delta
x}{2})^{2}}e^{-\frac{1}{2\sigma^{2}}(x_{2}-\frac{d}{2})^{2}}2\text{cos}\left(\frac{p_{2}\Delta
x}{\hbar}\right)\\\
&w_{4}=\dfrac{1}{N}\left(\dfrac{1}{\pi\hbar}\right)^{2}e^{-\frac{2\sigma^{2}}{\hbar^{2}}(p_{1}^{2}+p_{2}^{2})}e^{-\frac{1}{2\sigma^{2}}(x_{1}+\frac{d-\Delta
x}{2})^{2}}e^{-\frac{1}{2\sigma^{2}}(x_{2}-\frac{d-\Delta x}{2})^{2}}\\\
&w_{5}=\dfrac{1}{N}\left(\dfrac{1}{\pi\hbar}\right)^{2}e^{-\frac{2\sigma^{2}}{\hbar^{2}}(p_{1}^{2}+p_{2}^{2})}e^{-\frac{1}{2\sigma^{2}}(x_{1}+\frac{d-\Delta
x}{2})^{2}}e^{-\frac{1}{2\sigma^{2}}(x_{2}-\frac{d+\Delta x}{2})^{2}}\\\
&w_{6}=\dfrac{1}{N}\left(\dfrac{1}{\pi\hbar}\right)^{2}e^{-\frac{2\sigma^{2}}{\hbar^{2}}(p_{1}^{2}+p_{2}^{2})}e^{-\frac{1}{2\sigma^{2}}(x_{1}+\frac{d-\Delta
x}{2})^{2}}e^{-\frac{1}{2\sigma^{2}}(x_{2}-\frac{d}{2})^{2}}2\text{cos}\left(\frac{p_{2}\Delta
x}{\hbar}\right)\\\
&w_{7}=\dfrac{1}{N}\left(\dfrac{1}{\pi\hbar}\right)^{2}e^{-\frac{2\sigma^{2}}{\hbar^{2}}(p_{1}^{2}+p_{2}^{2})}e^{-\frac{1}{2\sigma^{2}}(x_{1}+\frac{d}{2})^{2}}e^{-\frac{1}{2\sigma^{2}}(x_{2}-\frac{d-\Delta
x}{2})^{2}}2\text{cos}\left(\frac{p_{1}\Delta x}{\hbar}\right)\\\
&w_{8}=\dfrac{1}{N}\left(\dfrac{1}{\pi\hbar}\right)^{2}e^{-\frac{2\sigma^{2}}{\hbar^{2}}(p_{1}^{2}+p_{2}^{2})}e^{-\frac{1}{2\sigma^{2}}(x_{1}+\frac{d}{2})^{2}}e^{-\frac{1}{2\sigma^{2}}(x_{2}-\frac{d+\Delta
x}{2})^{2}}2\text{cos}\left(\frac{p_{1}\Delta x}{\hbar}\right)\\\
&w_{9}=\dfrac{1}{N}\left(\dfrac{1}{\pi\hbar}\right)^{2}e^{-\frac{2\sigma^{2}}{\hbar^{2}}(p_{1}^{2}+p_{2}^{2})}e^{-\frac{1}{2\sigma^{2}}(x_{1}+\frac{d}{2})^{2}}e^{-\frac{1}{2\sigma^{2}}(x_{2}-\frac{d}{2})^{2}}2\text{cos}\left(\frac{p_{1}\Delta
x}{\hbar}\right)2\text{cos}\left(\frac{p_{2}\Delta x}{\hbar}\right)\\\
\end{split}$ (23)

All these contributions can be formally expressed as

$w_{j}(x_{1},x_{2},p_{1},p_{2})=\frac{4}{N(\pi\hbar)^{2}}e^{-\beta(x_{1}-x_{0,j})^{2}-\beta(x_{2}-y_{0,j})^{2}-\alpha(p_{1}^{2}+p_{2}^{2})}\;\frac{\cos(\gamma_{1,j}\,p_{1})}{D_{1,j}}\;\frac{\cos(\gamma_{2,j}\,p_{2})}{D_{2,j}},$
(24)

abbreviating $\alpha=2\sigma^{2}/\hbar^{2}$, $\beta=1/2\sigma^{2}$, and
introducing the the parameter arrays

$\begin{split}[\gamma_{1,j}]_{j=1}^{9}&=\left[0,0,0,0,0,0,\frac{\Delta
x}{\hbar},\frac{\Delta x}{\hbar},\frac{\Delta
x}{\hbar}\right],\qquad[\gamma_{2,j}]_{j=1}^{9}=\left[0,0,\frac{\Delta
x}{\hbar},0,0,\frac{\Delta x}{\hbar},0,0,\frac{\Delta x}{\hbar}\right];\\\
[D_{1,j}]_{j=1}^{9}&=[2,2,2,2,2,2,1,1,1],\qquad[D_{2,j}]_{j=1}^{9}=[2,2,1,2,2,1,2,2,1];\\\
[x_{0,j}]_{j=1}^{9}&=\left[-\frac{d+\Delta x}{2},-\frac{d+\Delta
x}{2},-\frac{d+\Delta x}{2},-\frac{d-\Delta x}{2},-\frac{d-\Delta
x}{2},-\frac{d-\Delta x}{2},-\frac{d}{2},-\frac{d}{2},-\frac{d}{2}\right]\\\
[y_{0,j}]_{j=1}^{9}&=\left[\frac{d-\Delta x}{2},\frac{d+\Delta
x}{2},\frac{d}{2},\frac{d-\Delta x}{2},\frac{d+\Delta
x}{2},\frac{d}{2},\frac{d-\Delta x}{2},\frac{d+\Delta
x}{2},\frac{d}{2}\right]\end{split}$ (25)

##  Appendix B Stepwise linear potential

To model the classical evolution of the Wigner function in phase space under
the gravitational interaction, one could use a stepwise linear potential.
Assuming that the mean distance $d$ between the two interfering particles is
greater than the arm separation $\Delta x$ and that the Gaussian wave packets
are well localized, $\sigma\ll\Delta x\ll d$, the Wigner function of the two-
particle state consists of the non-overlapping contributions (23). Each
contribution $w_{j}$ corresponds to the two particles being localized around
an average distance given by
$[\bar{x}_{j}]_{j=1}^{9}=\absolutevalue{x_{0,j}-y_{0,j}}$. Hence we can safely
approximate the gravitational force by opposite constant values
$F_{j}=m^{2}G/\bar{x}_{j}^{2}$ and $-F_{j}$ acting on particle 1 and 2,
respectively, as resulting from a stepwise linear interaction potential.

The classical time evolution of a particle under a constant force $F_{j}$
amounts to $x(t)=x(0)+p(0)t/m+F_{j}t^{2}/2m$ and $p(t)=p(0)+F_{j}t$. If we
neglect the deflection of the position coordinate and only account for the
momentum change—a reasonable approximation for the short experimental time
scales considered in Ref. [1]—, the classical time evolution of each Wigner
function contribution reads as

$w_{j}(x_{1},x_{2},p_{1},p_{2},t)=w_{j}(x_{1},x_{2},p_{1}-F_{j}t,p_{2}+F_{j}t).$
(26)

The problem with this solution of the classical Liouville equation based on
stepwise constant forces $F_{j}$ is that it will generally no longer
correspond to a valid state transformation. Indeed, we find that the momentum
marginal of the two-particle Wigner function, $W(p_{1},p_{2},t)=\int
W(x_{1},x_{2},p_{1},p_{2},t)\differential x_{1}\differential
x_{2}=\sum_{j}w_{j}(p_{1},p_{2},t)$, ceases to be a valid probability
distribution of observable momenta as it takes on negative values with growing
$t$. Specifically, we get

$w_{j}(p_{1},p_{2},t)=\frac{4}{N\pi\beta\hbar^{2}}e^{-\alpha(p_{1}-F_{j}t)^{2}-\alpha(p_{2}+F_{j}t)^{2}}\;\frac{\cos[\gamma_{1,j}(p_{1}-F_{j}t)]}{D_{1,j}}\;\frac{\cos[\gamma_{2,j}(p_{2}+F_{j}t)]}{D_{2,j}}.$
(27)

The negative regions of these functions get displaced with respect to one
another so that they no longer cancel under the $j$-sum. Our numerical
simulation predicts that a small but noticeable negativity of
$[\int|W(p_{1},p_{2},t)|\differential p_{1}\differential p_{2}-1]/2\approx
0.1\,\%$ arises at $t\geq 1.6\;\text{s}$ for the parameters used in the main
text.

##  Appendix C Approximations of the potential

In order to find the classical time evolution of the Wigner function we need
to solve the Hamilton’s equation of motions with the respective approximation
of the potential. From the Hamiltonian in the new coordinates given by Eq. (9)
we can expand the Newtonian potential $V_{\mathrm{N}}$ in Taylor series in
$x_{\mathrm{rel}}$

$V_{\mathrm{N}}(x_{\mathrm{rel}})=-\dfrac{\kappa}{d+\sqrt{2}x_{\mathrm{rel}}}=\kappa\sum_{n=0}^{\infty}\dfrac{(-1)^{n+1}2^{\frac{n}{2}}x_{\mathrm{rel}}^{n}}{d^{n+1}}=-\dfrac{\kappa}{d}\sum_{n=0}^{\infty}\left(\dfrac{-\sqrt{2}x_{\mathrm{rel}}}{d}\right)^{n}$
(28)

where the expansion converges for
$\absolutevalue{x_{\mathrm{rel}}}<\frac{d}{\sqrt{2}}$. After dropping the
terms of the order $\order{x_{\mathrm{rel}}^{3}}$ and higher we get the
potential

$V_{\mathrm{Swift}}(x_{\mathrm{rel}})=-\dfrac{\kappa}{d}+\dfrac{\sqrt{2}\kappa}{d^{2}}x_{\mathrm{rel}}-\dfrac{2\kappa}{d^{3}}x_{\mathrm{rel}}^{2},$
(29)

Between two out of the four possible combinations of paths we have
$x_{\mathrm{rel}}=0$ and so the expansion gives the exact result, in the other
two cases we have $x_{\mathrm{rel}}=\pm\frac{\Delta x}{\sqrt{2}}$ and we get
$\frac{\Delta x}{d}\approx 0.556$, where we have used the previously proposed
values $\Delta x=250\mu\text{m}$ and $d=450\mu\text{m}$. It follows that while
the Taylor expansion converges, we should not expect to obtain exact results
by replacing the Newtonian potential $V_{\mathrm{N}}$ with
$V_{\mathrm{Swift}}$.

The potential $V_{\mathrm{Fit}}$, given as,

$V_{\mathrm{Fit}}(x_{\mathrm{rel}})=-\dfrac{\kappa}{d}+\dfrac{\sqrt{2}\kappa}{d^{2}-\Delta
x^{2}}x_{\mathrm{rel}}-\dfrac{2\kappa}{d(d^{2}-\Delta
x^{2})}x_{\mathrm{rel}}^{2},$ (30)

was computed such that it coincides with Newtonian potential for
$x_{\mathrm{rel}}=0$ and $x_{\mathrm{rel}}=\pm\frac{\Delta x}{\sqrt{2}}$.
Nevertheless one must note that it is not the potential, but the derivative of
the potential that appears in the Liouville equation. Let
$V^{\prime}_{\mathrm{Swift}}$ and $V^{\prime}_{\mathrm{Fit}}$ denote the
derivatives of the respective potential, then one should ask how well
$V^{\prime}_{\mathrm{Swift}}$ and $V^{\prime}_{\mathrm{Fit}}$ approximate the
derivative of the Newtonian potential $V^{\prime}_{\mathrm{N}}$. The
potentials are compared in Fig. 2a and their derivatives are compared in Fig.
2b, one can see that the approximations are reasonably close and that
$V^{\prime}_{\mathrm{Fit}}$ is a better approximation of
$V^{\prime}_{\mathrm{N}}$ than $V^{\prime}_{\mathrm{Swift}}$.

(a) Comparison of the Newtonian potential $V_{\mathrm{N}}$ with the Taylor
expansion approximation $V_{\mathrm{Swift}}$ and with the approximation
$V_{\mathrm{Fit}}$. (b) Comparison of the derivative of the Newtonian
potential $V^{\prime}_{\mathrm{N}}$ with the derivatives of the Taylor
expansion approximation $V^{\prime}_{\mathrm{Swift}}$ and with the
approximation $V^{\prime}_{\mathrm{Fit}}$.

Figure 2: Comparison of the potentials and their derivatives. Previously
proposed values $\Delta x=250\mu\text{m}$ and $d=450\mu\text{m}$ were used;
the horizontal axis is in the units of $\Delta x$ while the vertical axis is
in the units of $\kappa$.

##  Appendix D Classical time evolution

Using the potential $V_{\mathrm{Swift}}(x_{\mathrm{rel}})$ we get the
Hamiltonian

$H(x_{\mathrm{avg}},x_{\mathrm{rel}},p_{\mathrm{avg}},p_{\mathrm{rel}})=\dfrac{p_{\mathrm{avg}}^{2}+p_{\mathrm{rel}}^{2}}{2m}-\dfrac{\kappa}{d}+\dfrac{\sqrt{2}\kappa}{d^{2}}x_{\mathrm{rel}}-\dfrac{2\kappa}{d^{3}}x_{\mathrm{rel}}^{2}$
(31)

which yields the equations of motion

$\displaystyle\dot{x}_{\mathrm{avg}}=\dfrac{p_{\mathrm{avg}}}{m},$
$\displaystyle\dot{p}_{\mathrm{avg}}=0,$ (32)
$\displaystyle\dot{x}_{\mathrm{rel}}=\dfrac{p_{\mathrm{rel}}}{m},$
$\displaystyle\dot{p}_{\mathrm{rel}}=-\dfrac{\sqrt{2}\kappa}{d^{2}}+\dfrac{4\kappa}{d^{3}}x_{\mathrm{rel}}.$
(33)

We immediately have

$\displaystyle
x_{\mathrm{avg}}(t)=\dfrac{p_{\mathrm{avg}}(0)}{m}t+x_{\mathrm{avg}}(0),$
$\displaystyle p_{\mathrm{avg}}(t)=p_{\mathrm{avg}}(0).$ (34)

To solve the other two equations we compute the second time derivative of
$p_{\mathrm{rel}}$, we get

$\ddot{p}_{\mathrm{rel}}=\dfrac{4\kappa}{d^{3}}\dot{x}_{\mathrm{rel}}=\dfrac{4\kappa}{md^{3}}p_{\mathrm{rel}}$
(35)

which is solved by the hyperbolic functions

$p_{\mathrm{rel}}(t)=A\sinh\left(2\sqrt{\dfrac{\kappa}{md^{3}}}t\right)+p_{\mathrm{rel}}(0)\cosh\left(2\sqrt{\dfrac{\kappa}{md^{3}}}t\right)$
(36)

where $A$ is a constant to be determined from the initial conditions. We also
get

$x_{\mathrm{rel}}(t)=\dfrac{A}{2}\sqrt{\dfrac{d^{3}}{\kappa
m}}\cosh\left(2\sqrt{\dfrac{\kappa}{md^{3}}}t\right)+\dfrac{p_{\mathrm{rel}}(0)}{2}\sqrt{\dfrac{d^{3}}{\kappa
m}}\sinh\left(2\sqrt{\dfrac{\kappa}{md^{3}}}t\right)+B,$ (37)

where $B$ is an integration constant. First of all we need to enforce that
(33) holds, we have

$\dot{p}_{\mathrm{rel}}(t)=2A\sqrt{\dfrac{\kappa}{md^{3}}}\cosh\left(2\sqrt{\dfrac{\kappa}{md^{3}}}t\right)+2p_{\mathrm{rel}}(0)\sqrt{\dfrac{\kappa}{md^{3}}}\cosh\left(2\sqrt{\dfrac{\kappa}{md^{3}}}t\right)$
(38)

and

$-\dfrac{\sqrt{2}\kappa}{d^{2}}+\dfrac{4\kappa}{d^{3}}x_{\mathrm{rel}}=-\dfrac{\sqrt{2}\kappa}{d^{2}}+2A\sqrt{\dfrac{\kappa}{md^{3}}}\cosh\left(2\sqrt{\dfrac{\kappa}{md^{3}}}t\right)+2p_{\mathrm{rel}}(0)\sqrt{\dfrac{\kappa}{md^{3}}}\sinh\left(2\sqrt{\dfrac{\kappa}{md^{3}}}t\right)+\dfrac{4\kappa}{d^{3}}B,$
(39)

By comparing the two expressions we get

$\dfrac{4\kappa}{d^{3}}B=\dfrac{\sqrt{2}\kappa}{d^{2}},\qquad\text{i.e.,}\quad
B=\dfrac{d}{2\sqrt{2}}.$ (40)

Finally we obtain $A$ from the initial condition

$x_{\mathrm{rel}}(0)=\dfrac{A}{2}\sqrt{\dfrac{d^{3}}{\kappa
m}}+\dfrac{d}{2\sqrt{2}},\qquad\text{which gives}\quad A=2\sqrt{\dfrac{\kappa
m}{d^{3}}}\left(x_{\mathrm{rel}}(0)-\dfrac{d}{2\sqrt{2}}\right).$ (41)

This leaves us with

$\displaystyle x_{\mathrm{rel}}(t)$
$\displaystyle=\left(x_{\mathrm{rel}}(0)-\dfrac{d}{2\sqrt{2}}\right)\cosh\left(2\sqrt{\dfrac{\kappa}{md^{3}}}t\right)+\dfrac{p_{\mathrm{rel}}(0)}{2}\sqrt{\dfrac{d^{3}}{\kappa
m}}\sinh\left(2\sqrt{\dfrac{\kappa}{md^{3}}}t\right)+\dfrac{d}{2\sqrt{2}},$
(42) $\displaystyle p_{\mathrm{rel}}(t)$ $\displaystyle=2\sqrt{\dfrac{\kappa
m}{d^{3}}}\left(x_{\mathrm{rel}}(0)-\dfrac{d}{2\sqrt{2}}\right)\sinh\left(2\sqrt{\dfrac{\kappa}{md^{3}}}t\right)+p_{\mathrm{rel}}(0)\cosh\left(2\sqrt{\dfrac{\kappa}{md^{3}}}t\right).$
(43)

We obtain the solution in terms of the original coordinates by using the
inverse transformation

$\displaystyle
x_{1}=\dfrac{1}{\sqrt{2}}(x_{\mathrm{avg}}-x_{\mathrm{rel}}-\dfrac{d}{\sqrt{2}}),$
$\displaystyle p_{1}=\dfrac{1}{\sqrt{2}}(p_{\mathrm{avg}}-p_{\mathrm{rel}}),$
(44) $\displaystyle
x_{2}=\dfrac{1}{\sqrt{2}}(x_{\mathrm{avg}}+x_{\mathrm{rel}}+\dfrac{d}{\sqrt{2}}),$
$\displaystyle p_{2}=\dfrac{1}{\sqrt{2}}(p_{\mathrm{avg}}+p_{\mathrm{rel}}).$
(45)

We arrive at

$\displaystyle x_{1}(t)$
$\displaystyle=\dfrac{1}{\sqrt{2}}\left(\dfrac{p_{\mathrm{avg}}(0)}{m}t+x_{\mathrm{avg}}(0)-\left(x_{\mathrm{rel}}(0)-\dfrac{d}{2\sqrt{2}}\right)\cosh\left(2\sqrt{\dfrac{\kappa}{md^{3}}}t\right)-\dfrac{p_{\mathrm{rel}}(0)}{2}\sqrt{\dfrac{d^{3}}{\kappa
m}}\sinh\left(2\sqrt{\dfrac{\kappa}{md^{3}}}t\right)-\dfrac{3d}{2\sqrt{2}}\right),$
(46) $\displaystyle x_{2}(t)$
$\displaystyle=\dfrac{1}{\sqrt{2}}\left(\dfrac{p_{\mathrm{avg}}(0)}{m}t+x_{\mathrm{avg}}(0)+\left(x_{\mathrm{rel}}(0)-\dfrac{d}{2\sqrt{2}}\right)\cosh\left(2\sqrt{\dfrac{\kappa}{md^{3}}}t\right)+\dfrac{p_{\mathrm{rel}}(0)}{2}\sqrt{\dfrac{d^{3}}{\kappa
m}}\sinh\left(2\sqrt{\dfrac{\kappa}{md^{3}}}t\right)+\dfrac{3d}{2\sqrt{2}}\right),$
(47) $\displaystyle p_{1}(t)$
$\displaystyle=\dfrac{1}{\sqrt{2}}\left(p_{\mathrm{avg}}(0)-2\sqrt{\dfrac{\kappa
m}{d^{3}}}\left(x_{\mathrm{rel}}(0)-\dfrac{d}{2\sqrt{2}}\right)\sinh\left(2\sqrt{\dfrac{\kappa}{md^{3}}}t\right)-p_{\mathrm{rel}}(0)\cosh\left(2\sqrt{\dfrac{\kappa}{md^{3}}}t\right)\right),$
(48) $\displaystyle p_{2}(t)$
$\displaystyle=\dfrac{1}{\sqrt{2}}\left(p_{\mathrm{avg}}(0)+2\sqrt{\dfrac{\kappa
m}{d^{3}}}\left(x_{\mathrm{rel}}(0)-\dfrac{d}{2\sqrt{2}}\right)\sinh\left(2\sqrt{\dfrac{\kappa}{md^{3}}}t\right)+p_{\mathrm{rel}}(0)\cosh\left(2\sqrt{\dfrac{\kappa}{md^{3}}}t\right)\right),$
(49)

which simplifies to

$\displaystyle x_{1}(t)$
$\displaystyle=\dfrac{p_{1}(0)+p_{2}(0)}{2m}t+\dfrac{x_{1}(0)+x_{2}(0)}{2}-\dfrac{x_{2}(0)-x_{1}(0)-\frac{3d}{2}}{2}\cosh\left(2\sqrt{\dfrac{\kappa}{md^{3}}}t\right)$
$\displaystyle-\dfrac{p_{2}(0)-p_{1}(0)}{4}\sqrt{\dfrac{d^{3}}{\kappa
m}}\sinh\left(2\sqrt{\dfrac{\kappa}{md^{3}}}t\right)-\dfrac{3d}{4},$ (50)
$\displaystyle x_{2}(t)$
$\displaystyle=\dfrac{p_{1}(0)+p_{2}(0)}{2m}t+\dfrac{x_{1}(0)+x_{2}(0)}{2}+\dfrac{x_{2}(0)-x_{1}(0)-\frac{3d}{2}}{2}\cosh\left(2\sqrt{\dfrac{\kappa}{md^{3}}}t\right)$
$\displaystyle+\dfrac{p_{2}(0)-p_{1}(0)}{4}\sqrt{\dfrac{d^{3}}{\kappa
m}}\sinh\left(2\sqrt{\dfrac{\kappa}{md^{3}}}t\right)+\dfrac{3d}{4},$ (51)
$\displaystyle p_{1}(t)$
$\displaystyle=\dfrac{p_{1}(0)+p_{2}(0)}{2}-2\sqrt{\dfrac{\kappa
m}{d^{3}}}\dfrac{x_{2}(0)-x_{1}(0)-\frac{3d}{2}}{2}\sinh\left(2\sqrt{\dfrac{\kappa}{md^{3}}}t\right)-\dfrac{p_{2}(0)-p_{1}(0)}{2}\cosh\left(2\sqrt{\dfrac{\kappa}{md^{3}}}t\right),$
(52) $\displaystyle p_{2}(t)$
$\displaystyle=\dfrac{p_{1}(0)+p_{2}(0)}{2}+2\sqrt{\dfrac{\kappa
m}{d^{3}}}\dfrac{x_{2}(0)-x_{1}(0)-\frac{3d}{2}}{2}\sinh\left(2\sqrt{\dfrac{\kappa}{md^{3}}}t\right)+\dfrac{p_{2}(0)-p_{1}(0)}{2}\cosh\left(2\sqrt{\dfrac{\kappa}{md^{3}}}t\right).$
(53)

where we have also used

$x_{\mathrm{rel}}(0)-\dfrac{d}{2\sqrt{2}}=\dfrac{1}{\sqrt{2}}\left(x_{2}(0)-x_{1}(0)-\dfrac{3d}{2}\right).$
(54)

In an analogous way, it is possible to derive the equation of motions under
the second-order polynomial approximation of the potential $V_{\mathrm{Fit}}$.
By using the approximated potential from Eq. (12) we obtain the Hamiltonian as

$H(x_{\mathrm{avg}},x_{\mathrm{rel}},p_{\mathrm{avg}},p_{\mathrm{rel}})=\dfrac{p_{\mathrm{avg}}^{2}+p_{\mathrm{rel}}^{2}}{2m}-\dfrac{\kappa}{d}+\dfrac{\sqrt{2}\kappa}{d^{2}-\Delta
x^{2}}x_{\mathrm{rel}}-\dfrac{2\kappa}{d(d^{2}-\Delta
x^{2})}x_{\mathrm{rel}}^{2}.$ (55)

Using analogical calculations as in the previous case, we get the solutions to
the Hamilton’s equations

$\displaystyle x_{1}(t)$
$\displaystyle=\dfrac{p_{1}(0)+p_{2}(0)}{2m}t+\dfrac{x_{1}(0)+x_{2}(0)}{2}-\dfrac{x_{2}(0)-x_{1}(0)-\frac{3d}{2}}{2}\cosh\left(2\sqrt{\dfrac{\kappa}{d(d^{2}-\Delta
x^{2})m}}t\right)$
$\displaystyle-\dfrac{p_{2}(0)-p_{1}(0)}{4}\sqrt{\dfrac{d(d^{2}-\Delta
x^{2})}{\kappa m}}\sinh\left(2\sqrt{\dfrac{\kappa}{d(d^{2}-\Delta
x^{2})m}}t\right)-\dfrac{3d}{4}$ (56) $\displaystyle x_{2}(t)$
$\displaystyle=\dfrac{p_{1}(0)+p_{2}(0)}{4}t+\dfrac{x_{1}(0_{+}x_{2}(0)}{2}+\dfrac{x_{2}(0)-x_{1}(0)-\frac{3d}{2}}{2}\cosh\left(2\sqrt{\dfrac{\kappa}{d(d^{2}-\Delta
x^{2})m}}t\right)$
$\displaystyle+\dfrac{p_{2}(0)-p_{1}(0)}{4}\sqrt{\dfrac{d(d^{2}-\Delta
x^{2})}{\kappa m}}\sinh\left(2\sqrt{\dfrac{\kappa}{d(d^{2}-\Delta
x^{2})m}}t\right)+\dfrac{3d}{4}$ (57) $\displaystyle p_{1}(t)$
$\displaystyle=\dfrac{p_{1}(0)+p_{2}(0)}{2}-2\sqrt{\dfrac{\kappa
m}{d(d^{2}-\Delta
x^{2})}}\dfrac{x_{2}(0)-x_{1}(0)-\frac{3d}{2}}{2}\sinh\left(2\sqrt{\dfrac{\kappa}{d(d^{2}-\Delta
x^{2})m}}t\right)$
$\displaystyle-\dfrac{p_{2}(0)-p_{1}(0)}{2}\cosh\left(2\sqrt{\dfrac{\kappa}{d(d^{2}-\Delta
x^{2})m}}t\right)$ (58) $\displaystyle p_{2}(t)$
$\displaystyle=\dfrac{p_{1}(0)+p_{2}(0)}{2}+2\sqrt{\dfrac{\kappa
m}{d(d^{2}-\Delta
x^{2})}}\dfrac{x_{2}(0)-x_{1}(0)-\frac{3d}{2}}{2}\sinh\left(2\sqrt{\dfrac{\kappa}{d(d^{2}-\Delta
x^{2})m}}t\right)$
$\displaystyle+\dfrac{p_{2}(0)-p_{1}(0)}{2}\cosh\left(2\sqrt{\dfrac{\kappa}{d(d^{2}-\Delta
x^{2})m}}t\right)$ (59)