# Large Spin Stern-Gerlach Interferometry for Gravitational Entanglement

Lorenzo Braccini<EMAIL_ADDRESS>Department of Physics and
Astronomy, University College London, Gower Street, WC1E 6BT London, United
Kingdom  Martine Schut  Van Swinderen Institute for Particle Physics and
Gravity, University of Groningen, 9747 AG Groningen, Netherlands  Alessio
Serafini  Department of Physics and Astronomy, University College London,
Gower Street, WC1E 6BT London, United Kingdom  Anupam Mazumdar  Van Swinderen
Institute for Particle Physics and Gravity, University of Groningen, 9747 AG
Groningen, Netherlands  Sougato Bose  Department of Physics and Astronomy,
University College London, Gower Street, WC1E 6BT London, United Kingdom

###### Abstract

Recently, there has been a proposal to test the quantum nature of gravity in
the laboratory by witnessing the growth of entanglement between two masses in
spatial quantum superpositions. The required superpositions can be created via
Stern-Gerlach interferometers, which couple an embedded spin qubit quantum
state to the spatial dynamics of each mass. The masses would entangle only if
gravity is quantum in nature. Here, we generalise the experiment to an
arbitrary spin $j$ or equivalently to an ensemble of uniformly coupled spins.
We first exemplify how to create a generalized Stern-Gerlach interferometer,
which splits the mass into $2j+1$ trajectories. This shows that a controlled
protocol can be formulated to encode the amplitudes of any spin state to a
spatial superposition. Secondly, two masses in spatial superpositions of the
above form are left to interact via gravity, and the entanglement is computed.
Different families of initial spin states are varied to find the optimal spin
state that maximizes the entanglement. We conclude that larger spins can offer
a modest advantage in enhancing gravity-induced entanglement.

††preprint: APS/123-QED

##  I Introduction

In the last decades, the search for a unifying theory of general relativity
and quantum mechanics has been central in discussions of theoretical physics.
Historically, a crucial difficulty faced by the Quantum Gravity (QG) community
has been the absence of related phenomena detectable in a laboratory setting.
In the low energy regime, where most of the experiments are conceived, it is
tough to test directly the carrier of the gravitational force, known as the
massless spin-2 graviton [1, 2]. In the high energy regime ($M_{p}\sim
10^{19}$ GeV), it is inconceivable to probe the time scale ($t_{p}\sim
10^{-44}$ s) and length scale ($\ell_{p}\sim 10^{-35}$ m) of spacetime where
it is believed that actual non-perturbative QG effects would be relevant.

Recently, a promising proposal has emerged which will enable us to witness the
quantum coherent nature of gravitational superpositions while still staying
within the weak field regime. The proposed experiment uses a quantum spatial
superposition of mesoscopic objects (say, via a Stern-Gerlach interferometer
[3]). The state describing the whole system is assumed to be initially in a
product state, and only the gravitational interaction is present. The proposed
experiment is based on the fact that local operation and classical
communication cannot yield entanglement [4, 5]. Thus, if gravity - the sole
channel between the masses - is inherently quantum, it would entangle two test
masses in a quantum superposition, while a classical natured gravity cannot.
Measurement of entanglement would genuinely witness the quantum nature of
gravity. This scheme can be dubbed as Quantum Gravity-induced Entanglement of
Masses (QGEM) [6, 7, 8].

The lowest order interaction between the masses is Newtonian, which in a
perturbative quantum theory of gravity results from the interaction between
the masses through virtual graviton exchanges. That this process can provide
the quantum communication necessary to entangle the two objects as highlighted
in Ref. [9, 10]. In this regard witnessing the entanglement between the two
test masses due to their Newtonian interaction is similar in spirit to Bell’s
test of quantum mechanics [11]. Another way of justifying the protocol is the
observation that entanglement generation requires an operator valued
interaction between the masses which is not possible with a classical
mediator, as was shown explicitly by quantizing the spin-2 graviton around a
Minkowski background [10]. The idea that entanglement detection would witness
QG is supported by several independent works, namely that it is an evidence of
the quantum superposition of geometries [12], by the relation between
different degrees of freedom of the gravitational field [13, 14, 15], path
integral calculations [16], generalized probabilistic theories [17], quantum
reference frames [18], and canonical quantization of gravity [19, 20, 21].

This class of proposed experiments have many other potentialities to study
fundamental questions: spin-2 character of the exchanged graviton via light-
matter entanglement (a quantum counterpart of the bending of light) [22],
discreteness of time [23, 24], the weak equivalence principle [25], fifth
force [26], massive graviton [27], and non-local gravitational interactions
arising from string theory [20]. Alternative apparatus, such as optomechanical
systems [28, 29, 30] and quantum gasses [31], have been studied in the
literature, showing both potentialities and difficulties for detecting
entanglement. Non-Gaussianity of Bose-Einstein condensates [32] may also be
proof of QG. Instead of witnessing gravitational entanglement, a recent
proposal has been developed that emphasizes testing the quantum nature of
gravity via the measurement-induced disturbance on a quantum system [33].

The aim of our paper is to investigate whether large spins embedded in
nanocrystals could increase the entanglement generated by the gravitational
interaction. This is motivated by several considerations. First of all, in the
quantum sensing community, it has been noted that large spins (and spins
ensembles) are excellent sensors [34]. Squeezed state and superposition of
coherent states of the spin system reach the fundamental Heisenberg limit [35,
36]. Thus, large spins are a promising system for sensing the weak
entanglement generated by the gravitational interaction between two mesoscopic
masses. Secondly, large spins in ferromagnetic nanoparticles could be used to
create large superpositions [37]. Thirdly, so far, the spin qubit based
gravitational entanglement experiments have been formulated in terms of a
single Nitrogen Vacancy (NV) centre in a mesoscopic diamond crystal. However,
it is an experimental challenge to create a single NV-center in a nano-
diamond, as multiple NVs (spin qubits) may be present. In cases when the
simplifying assumption that the coupling with the magnetic gradient is the
same for each spin qubit is valid, and when we can neglect the spin-spin
coupling, the system can be treated as a large effective spin.

We consider a large spin undergoing a spatial splitting: we will call this a
“Generalized Stern-Gerlach” (GSG) interferometer. As the first part of the
paper, we describe the dynamics of this process and find that the
superposition distance spreads proportionally to total angular momentum $j$,
while the splitting time remains constant. The mass is superposed in $2j+1$
coherent states, i.e. semi-classical paths. The protocol exemplifies a way to
encode any spin state to a spatial superposition, fully capturing the ability
of large spins in spatial-dependent protocols. As the second part, two such
states, each with $2j+1$ semi-classical paths generated due to spin-motion
interaction (for example, as generated in adjacent GSGs, but could be
otherwise), is let interact via gravity. The gravitational entanglement
generated is computed and compared for various values of spins.

The paper is organized as follows. In Sec. II, the GSG interferometer is
presented, deriving the splitting and recombining dynamics that create the
one-loop interferometer. We find that GSG is a map from any spin state to a
spatial superposition. Sec. III introduces the spins states that are used in
this study (coherent spin states, a superposition of coherence spin states and
squeezed spin states). In Sec. IV, two different two-interferometer geometries
are presented. The dynamics of these set-ups under gravitational interaction
is described in Sec. V, and the optimal spin states that maximize the
gravitational entanglement are numerically computed in Sec. VI. A suitable
witness and the influence of decoherence is discussed in Sec. VII and Sec.
VIII, respectively. The concluding remarks are given in Sec. IX.

##  II Generalized Stern-Gerlach Interferometer with Large Spins

Let us consider a large spin $j$ embedded in a mass $M$. The spin Hilbert
space $\mathbb{C}^{2j+1}$ is spanned in the Dicke basis $\\{\ket{j,m}\\}$ with
$m\in\left[-j,j\right]$, which are the eigenvectors of the $J_{z}$ operator,
where $z$ is the quantization axis. The $j$-dependency of the quantum states
will not be carried forward, as $j$ is assumed to be constant over the
experiment. Such a system can physically arise either from a single large spin
$j$ or an ensemble of $n$ spin–1/2 particles, with $n=j/2$, embedded in a
nanocrystal. The latter requires that the spin ensemble is permutationally
invariant, i.e. each spin couples uniformly with the external forces, and no
spin-spin interaction is present [38, 39]. While this may not be always an
accurate assumption, it is sometimes used for spins ensembles in a nano-
crystal [40], and we are going to make this simplifying assumption.

We consider a protocol that can be described in $6$ stages as follows:

(1) System Initiation: The system is initially prepared in the ground state
for both spin and position degrees of freedom, i.e.
$\ket{\psi}=\ket{\psi_{S}}\otimes\ket{\psi_{x}}=\ket{m=-j}\otimes\ket{0}$. For
the position degree of freedom, we assume that the ground state can be
achieved by usual feedback cooling of the mass in a trap [41]. The spin state
can be created by applying a strong magnetic field in the direction of the
quantization axis. Subsequently, both the trapping mechanism and the magnetic
field are turned off.

(2) Spin State Preparation: A desired spin state can be created by applying
secondary magnetic fields or by using microwave and radio-frequency field
pulses. Generally, the spin state will take the form:

$\ket{\psi_{S}}=\sum_{m=-j}^{j}c_{m}\ket{m}\,,$ (1)

where $\\{c_{m}\\}$ are the coefficients in the Dicke basis. A review of the
spin states considered in this work is given in Sec III, while the following
dynamics is kept in full generality.

(3) Splitting Process: A uniform magnetic gradient
$\left(\partial_{x}B\right)$ is applied to the mass. We will assume that there
is no external rotation of the mass [42, 43]. The splitting Hamiltonian is
then given by:

$H=\hbar\omega_{M}a^{\dagger}a-\hbar gJ_{z}\left(a+a^{\dagger}\right)\,.$ (2)

The first term is that of the quantum harmonic oscillator in which the mass is
trapped, with frequency $\omega_{M}$ and ladder operator $a$. The Stern-
Gerlach spin-position coupling is the second term, where $J_{z}$ is the
$z$-spin operator and the coupling constant for the mass $M$ is given by

$g=\tilde{g}\mu_{B}\sqrt{\frac{1}{2\hbar
M\omega_{M}}}\left(\partial_{x}B\right)\,,$ (3)

where $\mu_{B}$ is the Bohr magneton and $\tilde{g}$ is the Landé g-factor.

It has to be noted that this work is the first theoretical description of GSG
interferometry, and does not aim to resolve well-known experimental challenges
of creating a mesoscopic superposition with the use of Stern-Gerlach
interferometery. As pointed out in Ref. [44, 45, 46], the diamagnetic response
of the mass – a nanodiamond with a single NV center – placed in the magnetic
gradient significantly increases the splitting time, representing one of the
greatest experimental issues111However, once the initial superposition is
created of decent splitting, the diamagnetic induced potential can also aid to
create a larger superposition with the help of wires, see [47, 48].. We find
that this challenge is also present for large spins, or multiple NV centers
making an effective large spin; in Appendix A, the Hamiltonian in Eq. (2) is
derived by considering the diamagnetic response of the mesoscopic mass and a
linear magnetic gradient. Here we proceed by assuming that $\omega_{M}$ can be
engineered to be independent of $\left(\partial_{x}B\right)$. If this
assumption is not made, all our results will still be valid except occurring
over much longer time scales as in the spin-1/2 case. In essence, the
assumption of the independence of $\omega_{M}$ to $\left(\partial_{x}B\right)$
is difficult to fulfill, without being able to nullify the effect of
diamagnetism. Hence, this assumption is equivalent to be able to engineer a
composite material with extremely low magnetic susceptibility. We are making
this simplifying assumption herewith in order to make progress to be able to
compare the efficacy of large vs small spins in equal footing in generating
gravitationally induced entanglement.

The unitary evolution and the time evolution of the state resulting from the
Hamiltonian in Eq. (2) are derived in Appendix B. It follows that the quantum
state at time $t$ is given by:

$\Ket{\psi(t)}=\sum_{m=-j}^{j}c_{m}e^{i\frac{g^{2}}{\omega_{M}^{2}}m^{2}\left(\omega_{M}t-\sin(\omega_{M}t)\right)}\ket{m}\otimes\Ket{\alpha_{m}(t)}\,,$
(4)

where the position coherent states222Mathematically, position coherent states
are defined by the application of a displacement operator to the vacuum state
$\ket{\alpha}=\mathcal{D}(\alpha)\ket{0}=e^{\alpha^{\dagger}a}\ket{\alpha}$.
They are eigenstates of the annihilation operator
($a\ket{\alpha}=\alpha\ket{\alpha}$) and Gaussian states in both position and
momentum representations. They form an over-complete basis, i.e.
$\int_{\mathcal{C}}\ket{\alpha}\bra{\alpha}d^{2}\alpha$ and they are states of
minimum uncertainty. have parameters:

$\alpha_{m}(t)=m\frac{g}{\omega_{M}}\left(1-e^{-i\omega_{M}t}\right)\,.$ (5)

Under the evolution, the spin and position degrees of freedom entangle –
represented by the $m$-dependency in $\Ket{\alpha_{m}(t)}$ state of Eq. (4) –,
and the mass’ position degree of freedom branches out over $2j+1$ paths, i.e.
coherent state trajectories. Each path has minimum uncertainty; its position
and momentum expectation values are given by:

$\displaystyle\braket{x_{m}(t)}$
$\displaystyle=m\frac{g}{\omega_{M}}\sqrt{\frac{2\hbar}{M\omega_{M}}}\left(1-\cos(\omega_{M}t)\right),$
(6) $\displaystyle\braket{p_{m}(t)}$ $\displaystyle=mg\sqrt{\frac{2\hbar
M}{\omega_{M}}}\sin(\omega_{M}t)\,.$ (7)

Let us define the spin probability distribution
$P_{s}(m):=\left|\braket{\psi_{S}}{m}\right|^{2}=|c_{m}|^{2}$. We denote a
normal distribution with expectation value $\mu$ and standard deviation
$\sigma$ as $\mathcal{N}\left(\mu,\sigma^{2}\right)$. The probability of
finding the mass at $x$ at time $t$ in the position representation is:

$P\left(x,t\right)=\sum_{m=-j}^{j}P_{s}(m)\mathcal{N}\left(\braket{x_{m}(t)},\sigma_{x_{m}}^{2}=\frac{\hbar}{M\omega_{M}}\right)\,,$
(8)

which represents a sum of Gaussians, weighted by the discrete spin probability
distribution $P_{s}(m)$. An example of $P\left(x,t\right)$ is given in Fig. 1.

Figure 1: Probability distribution $P(x,t)$ as a function of time (in
frequency unit $\omega_{M}$) in the position representation (in $\Delta x$
units) of a Coherent Spin State $\Ket{j=5;\phi=0,\theta=\pi/2}$ in a GSG
Interferomenter. The mass undergoes the splitting and recombining process in
$11$ semi-classical paths, weighted by the probability distribution of the
$m$-states.

(4) Maximum Displacement and Interaction: The maximum displacement between
adjacent paths is achieved at the splitting time $t_{s}:=\pi/\omega_{M}$, and
takes the value:

$\Delta
x:=\braket{x_{m+1}\left(t_{s}\right)}-\braket{x_{m}\left(t_{s}\right)}=2\sqrt{\frac{2\hbar}{M\omega_{M}}}\frac{g}{\omega_{M}}\,,$
(9)

This splitting is independent on which adjacent paths $m$ and $m+1$ are taken,
as the distance between each adjacent path is the same.

The spatial superposition spreads over the total range

$\Delta
D:=\braket{x_{j}\left(t_{s}\right)}-\braket{x_{-j}\left(t_{s}\right)}=2j\Delta
x\,,$ (10)

i.e. it scales linearly with the total angular momentum, while the splitting
time is independent of $j$, promising a larger total splitting within the same
$t_{s}$ for higher $j$.

At the time $t_{s}$, the magnetic gradient is switched off and the system is
described by the quantum state:

$\Ket{\psi(t_{s})}=\sum_{m=-j}^{j}c_{m}e^{i\pi\frac{g^{2}}{\omega_{M}^{2}}m^{2}}\ket{m}\otimes\Ket{\alpha_{m}=m\Delta
x}\,.$ (11)

During this stage, the system can interact with external forces, such as
gravity, and, as we will see, generate entanglement.

(5) Recombination Process: Assuming that during the interaction the masses
were not displaced, the state can be recombined by applying a magnetic
gradient with the same profile as during the splitting process stage. After a
time $t_{s}$ the mass is brought back to the center of the trap and the
quantum state is give by:

$\Ket{\psi(2t_{s})}=\left(\sum_{m=-j}^{j}c_{m}e^{i2\pi\frac{g^{2}}{\omega_{M}^{2}}m^{2}}\ket{m}\right)\otimes\Ket{0}\,,$
(12)

where we assumed, for now, that no interaction was present during the previous
stage. The recombining magnetic field has to be, with extreme precision,
identical to the splitting one. Small differences between the two fields would
lead to the Humpty-Dumpty effect [49], as the recombination process would be
required to be a perfect time reversal of the splitting. Recently, a full-loop
has been achieved in experiments [3, 50, 51], in the context of atom
interferometry [52].

(6) Measurements: One can note that there is no spatial dependency Eq. (12):
the spatial degrees of freedom are disentangled from the spin ones and the
mass is placed back to the center of the trap. the position component of the
wavefunction can be trivially traced out, and, thus, general spin measurements
can be carried out on the spin embedded in the mass.

##  III  Families of Spin States

(a) Spin Coherent State

(b) Superposition of Two CSS

(c) Squeezed Spin State

Figure 2: Examples of Husimi-Q quasiprobability distributions of (a) a Spin
Coherent State, (b) a superposition of two CSS, and (c) a Spin Squeezed State.
The axis are the spin axis with spin m from -j to j.

The GSG interferometer represents a map from a spin state to an entangled
superposition of spin and position coherent states. The spins states
$\ket{\psi_{S}}$ used in this study are: Coherent Spin States (CSS), a
superposition of CSS, and Squeezed Spin State (SSS). In this Section a brief
overview of the aforementioned states is presented.

###  III.1 Coherent Spin States

A position coherent state of a quantum harmonic oscillator is a semi-classical
state describing a translation of the vacuum state. Similarly, a CSS is
defined as the state arising from an arbitrary rotation of the spin ground
state, $\ket{-j}$333The CSS can be written in terms of Euler angles
$\ket{\alpha,\beta,\gamma}:=\mathcal{R}\left(\alpha,\beta,\gamma\right)\ket{m=j}$,
where we used the Euler’s angles $\alpha,\beta,\gamma$ and the arbitrary
rotation $\mathcal{R}\left(\alpha,\beta,\gamma\right)\coloneqq e^{-i\alpha
J_{z}}e^{-i\beta J_{y}}e^{-i\gamma J_{z}}$. . This can be experimentally
prepared by applying a magnetic field on an arbitrary axis or alternatively by
applying a cavity field in resonance with the spin frequency. CSS states can
be expanded on the Dicke basis in terms of the Wigner D-Matrix. By defining
the quantity $\mu=e^{i\phi}\tan{\theta/2}$, where $\theta$ is the azimuth
angle with respect to the $z$-axis and $\phi$ is the $xy$-plane phase, a
general CSS is:

$\displaystyle\begin{split}\ket{\phi,\theta}&:=\mathcal{N}e^{\mu
J_{-}}\ket{m=-j}\\\
&=\mathcal{N}\sum_{m=-j}^{j}\mu^{j+m}\sqrt{\frac{2j!}{(j+m)!(j-m)!}}\ket{m}\,,\end{split}$
(13)

where $J_{-}$ is the spin lowering operator and the normalization is given by
$P_{s}\approx\mathcal{N}=\left(1-\left|\mu\right|^{2}\right)^{-j}$ [53]. This
parametrization has been chosen for its geometrical interpretation: the
created CSS is a semi-classical spin state pointing in the $\phi,\theta$
direction, i.e. the state with the maximum eigenvalue of the $\phi\theta-$spin
projection operator $J_{\phi,\theta}\ket{\phi,\theta}=j\ket{\phi,\theta}$. The
probability distribution in the $m$ representation is binomial and, thus, for
large $j$, it approaches a Gaussian distribution
$\mathcal{N}\left(\braket{\theta}=j\cos\theta,\sigma_{s}^{2}=j/2\sin^{2}\theta\right)$[54].
Hence, in a GSG interferometer, the embedded CSS allows to manipulate the mass
probability distribution. By varying the rotation angle $\theta$, the mass is
centered at $\Delta x\braket{\theta}$ with spread $\Delta x\sigma_{s}$444Given
three orthogonal axis ($\alpha$, $\beta$, and $\gamma$), the uncertainty
relation for a spin state is given by:  $\left(\Delta
J_{\alpha}\right)^{2}\left(\Delta
J_{\beta}\right)^{2}\geq\frac{\left|\braket{J_{\gamma}}\right|^{2}}{4}\,.$
Analogous to position coherent states, CSS are of minimum uncertainty $(\Delta
J_{\tilde{\alpha}})^{2}=(\Delta J_{\tilde{\beta}})^{2}=\hbar j/2$, where
$\tilde{\alpha}$ and $\tilde{\beta}$ are orthogonal to the CSS axis (pointing
in the $\theta$-$\phi$ direction) [55].. In the example given in Fig. 1, the
mass is centered and of maximum spread, as $\theta=\pi/2$.

###  III.2 Superpositions of Coherent Spin States

A superposition of CSS can be considered, and we will limit ourselves to the
simple case of a superposition of two states:

$\ket{\phi_{1},\theta_{1};\phi_{2},\theta_{2}}\coloneqq\mathcal{\tilde{N}}\left(\ket{\phi_{1},\theta_{1}}+\ket{\phi_{2},\theta_{2}}\right),$

where $\mathcal{\tilde{N}}$ is a normalization constant. Recently, such states
has been of interest in the quantum-sensing community [36] and the
NOON555$\ket{NOON}\equiv\frac{1}{\sqrt{2}}(\ket{m=-j}+\ket{m=j})=\frac{1}{\sqrt{2}}(\ket{\theta=0}+\ket{\theta=\pi})$
with $\phi=0$. state employed to create large mass superposition [37].

###  III.3 Squeezed Spin States

Finally, the last spin states that will be considered are SSS: states with
$\left(\Delta
J_{\alpha}\right)^{2}<\frac{\left|\braket{J_{\gamma}}\right|}{2}$, with
($\alpha\neq\gamma\in(x,y,z)$. Contrary to the quantum harmonic oscillator
case, for spin states there is not a unique notion of the squeezing parameter,
in literature several have been proposed [35]. Being agnostic on such an
issue, we define these states from a state-preparation perspective. A SSS is
created by the evolution of a CSS under the unitary evolution arising from
non-linear Hamiltonians in the spin operators: the one-axis twisting
Hamiltonian $H_{\text{OA}}\propto J_{y}^{2}$, which is experimentally easier
to implement, and the two-axis twisting Hamiltonian $H_{\text{TA}}\propto
J_{+}^{2}-J_{-}^{2}$[55]. For instance, we can define a one-axis twisted SSS
as:

$\ket{\chi_{OAT};\phi,\theta}\coloneqq
e^{-i\chi_{OAT}J_{y}^{2}}\ket{\phi,\theta}\,,$ (14)

and similarly for two-axis twisted SSS.

The probability distribution of these states in the $m$ basis is sub-binomial
[35]. Thus, with the use of a GSG, these states represents a controlled way to
increase and decrease the spread of the weights of the paths.

In order to picture the aforementioned spin states $\ket{\psi_{S}}$, we
compute the quasiprobability distribution known as the Husimi-Q function,
defined according to
$Q(\theta,\phi)\coloneqq\left|\braket{\theta,\phi}{\psi_{S}}\right|^{2}$. The
Husimi-Q functions of each of the aforementioned spin states is presented in
Fig. 2. To conclude the section, CSS represents a 2-parameter family of states
(with parameters $\theta$ and $\phi$), superposition of two CSSs are
4-parameter family ($\theta_{i}$ and $\phi_{i}$, with $i\in\\{1,2\\}$), and
the SSS is a 3-parameter family ($\chi$, $\theta$ and $\phi$). These
parameters can be varied to find the optimal states to maximize the detectable
entanglement, for instance, generated by gravity.

##  IV Geometries for the Intereferometers

In the original experimental protocol [7], two masses of
$M\sim\mathcal{O}(10^{-14}){\rm Kg}$, initially placed at distance $\Delta
s\sim{\cal O}(250){\rm\mu m}$, are individually put in a spatial superposition
of $\Delta x\sim{\cal O}(200){\rm\mu m}$. In order to witness QG, two GSG
interferometers have to be placed at distance $\Delta s$, and let interacting
via gravity. In this section, two geometries of the set-up will be presented
(represented in Fig. 3 and 4).

The photon-induced Casimir-Polder (CP) interaction between two spheres [56,
57] is the main source of unwanted entanglement between the neutral masses.
The CP potential is described by

$V_{\text{CP}}(r)\approx-\frac{23\hbar
c}{4\pi}\frac{R^{6}}{r^{7}}\left(\frac{\epsilon-1}{\epsilon+2}\right)^{2}\;,$

where $R$ is the radius of the mass, $\epsilon$ its dielectric constant, $r$
the distance between the two masses666The Casimir-Polder potential is given in
Ref. [56]. We assume that the separation between the plate and the mass ($r$)
is much larger than the radius of the mass ($R$), i.e. $r\gg R$. We have
furthermore assumed that the polarizability can be expressed as
$R^{3}(\epsilon-1)/(\epsilon+2)$, which is known to be the case for e.g.
diamond in a high vacuum [58]. The equation in this paper is found based on
these assumptions, see also Ref. [59]..

Since the CP potential scales with $r^{-7}$, while the Newtonian potential
scales with $r^{-1}$, this electromagnetic phenomenon restricted the
aforementioned experimental parameters to a large distance between the two
interferometers ($\Delta s\sim{\cal O}(200){\rm\mu m}$), such that the
gravitational-induced entanglement dominates at least ten times over the
electromagnetic one. However, it was later realized that separating the two
interferometers by a conducting plate ameliorates the parameter space
tremendously [59, 60], reducing the separation distance within $\Delta
s\sim{\cal O}(50){\rm\mu m}$. We will consider screening in our proposal,
which will allow for a higher gravitational entanglement.

Furthermore, following the insights [61, 62, 63], two interferometric
geometries will be used, namely, ‘linear’ and ‘parallel’ set-ups. They are
illustrated in Fig. 3 and Fig. 4, respectively.

Figure 3: Linear set-up: two spins ($J=5/2$) at distance $d$ are spatially
superposed and recombined on the $x$-axis via two GSG interferometers, with
superposition distance $\Delta x$. The minimum distance between the two
closest paths is $\Delta s=d-j\Delta x$. A Casimir shield (CS) is placed
between the masses and aligned with the $y$-axis. Figure 4: Parallel set-up:
two spins ($J=5/2$) at distance $\Delta s$ are spatially superposed and
recombined on the $y$-axis via two GSG interferometers, with superposition
distance $\Delta y$. The minimum distance between the two closest paths is
$\Delta s=d/2$. A Casimir shield (CS) is placed between the masses and aligned
with the $y$-axis.

##  V Gravitational Entanglement

As exemplified in Sec. II, it is possible to create an arbitrary spin-position
quantum state of the form in Eq. (11). We assume that, after the splitting,
the two masses are in tensor product of these quantum states777Hence, we are
neglecting entanglement generated during the splitting process. This was done
to not commit to a specific trajectory (in our case a sinusoidal function),
which may change with different splitting processes in experimental setting,
and ensure to give a comparative analysis between different spins.. The
quantum state of the composite system is
$\ket{\Psi_{AB}(t_{s})}=\Ket{\psi_{A}(t_{s})}\otimes\Ket{\psi_{B}(t_{s})}$,
where:

$\Ket{\psi_{I}(t_{s})}=\sum_{m=-j}^{j}c_{m}^{(I)}e^{i\pi\frac{g^{2}}{\omega_{M}^{2}}m^{2}}\ket{m}\otimes\Ket{x_{m}^{(I)}(t_{s})}\,,$

whit $x_{m}^{(I)}(t_{s})=m\;\Delta x+x_{0}^{{(I)}}$, $I\in\\{A,B\\}$. The
position $x_{0}^{A}$ and $x_{0}^{B}$ are the centers of the traps of each GSG
interferometer: in the linear set-up $x_{0}^{A}-x_{0}^{B}=d$, while in the
case of parallel $x_{0}^{A}-x_{0}^{B}=\Delta s$.

The system is let to interact via gravity for a time $\tau$. In the low energy
regime of the experiment, the gravitational interaction can be approximated by
the operator:

$\hat{H}_{G}=\frac{GM_{A}M_{B}}{\left|\hat{x}_{A}-\hat{x}_{B}\right|}\,,$

where $\hat{x}_{I}$ ($I=A,B$) are the position operators of each mass. We are
going to use the approximation that the position coherent state representing
each of the semi-classical paths is an eigenstate of the corresponding
position operator, i.e. $\hat{x}_{I}\ket{x_{m}^{(I)}}\approx
x_{m}^{(I)}\ket{x_{m}^{(I)}}$. This is motivated by the fact that the
distances involved in the experiment are much greater than the spread of the
position coherent state,
$\braket{x_{m}^{(I)}}\gg\sqrt{\braket{(\sigma_{x_{m}}^{(I)})^{2}}}$. 888This
occurs in quantum optics when the cavity field is off high intensity,
$\hat{x}\ket{\alpha}=\frac{a+a^{\dagger}}{\sqrt{2}}\ket{\alpha}\approx\frac{\alpha+\alpha^{*}}{\sqrt{2}}\ket{\alpha}=\text{Re}(\alpha)\ket{\alpha}=x\ket{\alpha}$
We will further assume that each path is not deformed by the gravitational
interaction, which is reasonable given the small acceleration experienced by
the masses [7]. Thus the wavefunction of the combined $A$ and $B$ system at a
time $\tau$ is given by:

$\displaystyle\ket{\Psi_{AB}(t_{s}+\tau)}$
$\displaystyle=\sum_{m,n=-j}^{j}c_{m}^{(A)}c_{n}^{(B)}e^{i\pi\frac{g^{2}}{\omega_{M}^{2}}(m^{2}+n^{2})}$
$\displaystyle\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;e^{-i\phi_{m,n}}\Ket{m;x_{m}^{(A)}}\Ket{n;x_{n}^{(B)}}\,,$

where the phases are given by:

$\phi_{m,n}^{\text{L}}=\frac{GM_{A}M_{B}\tau}{\hbar}\left|\Delta s+\Delta
x(2j+m-n)\right|^{-1}\,,$
$\phi_{m,n}^{\text{P}}=\frac{GM_{A}M_{B}\tau}{\hbar}\left|\Delta s^{2}+\Delta
x^{2}(m-n)^{2}\right|^{-1}\,,$

for the linear (L) and the parallel (P) set-ups. We took $\Delta x=\Delta y$
in order to give a comparison between the two set-ups. After the recombining
process, the final state of the system (at a final time $t_{f}=2t_{s}+\tau$)
is given by:

$\displaystyle\ket{\Psi_{AB}(t_{f})}$
$\displaystyle=\biggl{[}\sum_{m,n=-j}^{j}c_{m}^{(A)}c_{n}^{(B)}e^{i2\pi\frac{g^{2}}{\omega_{M}^{2}}(m^{2}-n^{2})}$
$\displaystyle\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;e^{-i\phi_{m,n}}\Ket{m}\otimes\ket{n}\biggr{]}\otimes\ket{x_{0}^{A};x_{0}^{B}}\,.$

The spin degrees of freedom are independent of the position ones, as the two
masses are recombined in the center of their respective trap. Noticing the
phases $\phi_{m,n}$, the remaining spin state is not separable and is
entangled999With the considered approximation, the proposed protocol can be
expressed as resulting from an evolution defined by a spin Hamiltonian. For
instance, for the linear set-up:
$\hat{H}_{G}^{\text{L}}=\frac{GM_{A}M_{B}}{\left|\Delta s+\Delta
x(2j+\hat{J}_{z}^{(A)}-\hat{J}_{z}^{(B)})\right|}.$ . Such entanglement will
not arise in the case of a classical c-numbered gravitational Hamiltonian
[10].

##  VI  von Neumann entropy

As a first measure of the entanglement, we will use the von Neumann entropy:

$S(\rho_{A})\coloneqq-\text{Tr}\left[\rho_{A}\ln(\rho_{A})\right],$

where $\rho_{A}$ is the reduced density matrix of sub-system $A$.

We will consider an experimental set-up as proposed in Ref. [7], with $\tau=2$
s, $M_{A}\approx M_{B}\sim 10^{-14}$ kg, and $\Delta x\approx 250\,\mu$m, with
the addition of Casimir screening, which allows $\Delta s=50\,\mu$m. Different
spins state families will be compered, and their parameters will be optimized
seeking to find the Maximal Gravity-induced Entanglement Entropy (MGEE) that
can be generated for a given spin family. For conciseness, some of the plots
of the numerical analysis will be reported in Appendix C.

Let us first study the case of CSS as easily implementable in a laboratory
setting. Each embedded CSS – defined according to Eq.(13) – has two parameters
$\phi_{I}$ and $\theta_{I}$, with $I\in\\{A,B\\}$. Hence, finding the MGEE is
an optimization problem with 4 parameters. The von Neumann entropy has been
found to be independent of the parameters $\phi_{I}$, and we will not further
consider them. The gravitational entanglement as a function of $\theta_{I}$ is
plotted in Fig. 5 for both parallel and linear geometries.

(a) J = 2

(b) J = 10

Figure 5: Gravitational entanglement as a function of the rotation angles
$\theta_{I}$ of two CSS with (a) $J=2$ and (b) $J=10$ in two SGS
interferometers the (i) linear and (ii) parallel set-ups. The symmetry axes
are in green.

As expected, a significant entropy is produced when the mass probability
distributions of the two interferometers are at their closest. In the linear
set-up, this implies that most of the mass of the left interferometer is
placed on the right ($\theta_{A}$ small) while in the right interferometer the
mass is on the left ($\theta_{B}$ large). In the parallel set-up, the mass
probability distributions have to match on both sides, i.e.
$\theta_{A}\approx\theta_{B}$. Higher spins require higher precision on the
rotation angles $\theta_{I}$: for instance, in the parallel set-up, the
requirement $\theta_{A}\approx\theta_{B}$ is more restrictive as $j$ increases
as the entanglement region (yellow band in Fig. 5) shrinks in width.

Hence, it is possible to maximize the entropy over $\theta_{I}$ and find the
optimal rotation angles $\theta_{I}^{(o)}$ of the input CSSs for a given $j$
that gives the MGEE. The found optimal values respect the symmetries of the
interferometers geometries:  $\theta_{A}^{(o)}=\pi-\theta_{B}^{(o)}$ for the
linear set-up, and $\theta_{A}^{(o)}=\theta_{B}^{(o)}$ in the parallel case.
Thus, among the four parameters ($\theta_{A,B},\phi_{A,B}$), only one is
independent in the entanglement.

The gravitational entanglement as a function of $\theta_{A}$ is plotted in
Fig. 6 for different spins values and fixed time ($\tau=2$ s). In the linear
set–up, $\theta_{A}^{(o)}$ increases with the spin value (and the mass is
progressively concentrate on one side of the intefermoeter), while the MGEE
remains approximately constant. In the parallel case and for small spins
values ($j\leq 5$), MGEE is at $\theta_{A}=\pi/2$. As the spin increases
($j>5$), a local minimum is formed at $\theta_{A}=\pi/2$ and two maxima emerge
symmetric with respect to $\theta_{A}=\pi/2$.

Figure 6: Entanglement entropy as a function of the only independent parameter
$\theta_{A}$, for different spins $j$. The parallel set-up is plotted with
dashed lines and the liner set-up with solid lines.

The time evolution of the MGEE can be computed for the optimal CSSs
$\ket{\theta_{A}^{(o)}(\tau)}$, where $\theta_{A}^{(o)}$ becomes a time
dependent parameter. This is plotted in Fig. 7.

Figure 7:  Time evolution of the MGEE generated by optimal CSS,
$\ket{\theta_{A}^{(o)}}$. The parallel set-up is plotted with dashed lines and
the liner set-up with solid lines.

Note that the entanglement is significantly enhanced in the parallel set-up.
However, in the linear set-up, the time evolution of the optimal case does not
differ significantly among different spins. This is also clear from Fig. 11 in
Appendix C, where the MGEE as a function of the spin value $j$ is plotted.

In the linear case, only a few $m$-components of the quantum state (the one
closest to $-j$ for the left interferometer and to $+j$ for the right
interferometer) generate entanglement. Given that there is not a significant
improvement in increasing the Hilbert space in the linear set-up, following,
we will only consider the parallel interferometric geometry.

In the parallel set-up, the MGEE monotonically increases over $j$ (see Fig.
11), tending to an asymptotic value. We will approximate the asymptotic value
to be the MGEE for $j=30$ (as it is the largest spin computed). At short
interaction time, the increment in the entanglement is less prominent for
higher $j$, i.e. the MGEE approaches its asymptotic value at smaller $j$ for
shorter interaction times $\tau$. Hence, to exploit the improvements of large
spins, we require a large enough interaction time, for instance, as previously
proposed, $\tau=2s$.

Let us now consider a superposition of CSS, restricting the discussion to
states of the form:

$\ket{\psi}=\mathcal{N}\left(\Ket{\frac{\pi}{2}+\delta\theta,\delta\phi}+\Ket{\frac{\pi}{2}-\delta\theta,-\delta\phi}\right)\,,$
(15)

with $\delta\theta,\delta\phi\in[0,\pi]$.101010The case of
$\ket{\psi}=\mathcal{N}\left(\ket{\theta_{I}^{(o)}+\delta\theta,\delta\phi}+\ket{\theta_{I}^{(o)}-\delta\theta,-\delta\phi}\right)$
has also been studied, but no enhancement in the entanglement entropy were
found for the cases of $\theta_{I}^{(o)}\neq\pi/2$, i.e. for $j>5$.

The entanglement entropy as a function of $\delta\theta$ and $\delta\phi$ is
plotted in Fig. 12a in Appendix C. The $\delta\phi$ variable leads to
interference-like patterns, and this superposition parameter will be discarded
as it does not significantly increase the entanglement entropy. There is an
interesting dependency on the $\delta\theta$ variable, and the entanglement
entropy as a function of $\delta\theta$ is plotted in Fig. 8.

With a superposition of CSS, the MESS is higher with respect to the
corresponding MGEE given by input CSS, for any given spin (with $j>1/2$, as
the case of $j=1/2$ is trivial). Interestingly, for $j\geq 2$, the MGEE for
superpositions of CSS is larger than the MGEE asymptotic value for CSS.
However, the MGEE is approximately constant for any $j\geq 2$. Hence, by
creating a superposition of CSS, the entropy can be maximized for $j=2$, and
further increment of the spin value does not lead to significant improvements.
The optimal superposition parameter, which maximizes the entropy, tends to
$\pi/2$ as $j$ increases.111111 A notable behaviour is the minima at
$\delta\phi=\pi/2$, as it represents a NOON state case. Contrary to quantum
sensing where such states are the most sensible, since they maximize the
Fisher information, for gravitational entanglement detection, they do not
represent a good choice of states, since they minimize the entanglement.

Figure 8:  Gravitational entanglement at time $\tau=2s$ between two
superpositions of CSS (Eq.(15)) as a function of superposition parameter
$\delta\theta$, describing the deviation from an initial CSS with
$\theta=\pi/2$.

Finally, it is possible to study SSS of the general form given in Eq. (14).
The entanglement entropy can be computed as a function of the squeezing
parameter ($\chi$), for both one-axis twisted and two-axis twisted states and
of the input CSS angle $\theta$ (the spin state before the squeezing). This is
plotted in Fig. 12b and 13 in Appendix C. It is found that using SSS does not
significantly increase the MGEE. The behaviour becomes chaotic for higher
squeezing parameters, where the SSS have higher entanglement entropy with
respect to the corresponding CSSs. Such high squeezing is not reasonable in an
experimental setting and so SSS does not represents a good input spin state.
Hence, we will discard such states from our further analysis.

To summarise the main results, the MGEE for different families of spin state
and spin values is given in Table 1. It seems that in terms of generating
gravitationally induced entanglement for these experimental parameters, a
superposition of coherent spin states with $j=2$ is most favorable.

Maximal Gravity-induced Entanglement Entropy
---
Spin |  $j=1/2$ |  $j=2$ |  $j=5$ |  $j=10$
Coherent |  0.59 |  1.19 |  1.32 |  1.32
Superposition |  0.59 |  1.40 |  1.38 |  1.39
1-axis Squeezed |  0.59 |  1.37 |  1.38 |  1.39
2-axis Squeezed |  0.59 |  1.22 |  1.36 |  1.33
Table 1:  Maximal Gravity-induced Entanglement Entropy generated by optimizing
families of spin states in two GSG interferometer (parallel set-up)
interacting via gravity with $\tau=2$ s, $M_{A}\approx M_{B}\sim 10^{-14}$ kg,
$\Delta x\approx 250\,\mu$m, and $\Delta s=50\,\mu$m.

##  VII Negativity

The von Neumann entropy requires a full knowledge of the quantum state of the
system and does not represent a suitable experimental observable. In
particular, when describing the test masses as an open quantum system
interacting with the environment (i.e. including decoherence), the system will
entangle with the environment, and the entanglement entropy generated by
gravity will be corrupted. Hence, the von Neumann entropy is not a good
entanglement measure for globally mixed states. In order to solve this issue,
the negativity of the quantum state is considered [64, 65]. Similar to the PPT
criterion121212It has to be noted that in our analysis, the PPT criterion
(which only consider the smallest negative eigenvalue) does not represent a
good measure of entanglement, even if in principle it can detect entanglement.
As $j$ increases, the smallest eigenvalue decreases (so the PPT criterion),
while the number of negative eigenvalues increases. With the inclusion of
decoherence such a quantity is then not robust, as it is too small., the
negativity is a sufficient (but not necessary) measure of entanglement. The
negativity is defined as the sum of the negative eigenvalues of the partial
transposed (denoted with the superscript $PT$) density matrix of the system,
i.e. $\braket{\mathcal{W}}=\sum_{\lambda^{PT}<0}\lambda^{PT}$. Note that –
according to our definition – such a quantity is negative (while in other
definitions, its absolute value is considered). It is possible to construct a
witness operator that measures this quantity as

$\hat{\mathcal{W}}=\sum_{\lambda<0}\ket{\lambda}\bra{\lambda}^{PT},$

where $\ket{\lambda}$ is the eigenvector corresponding to the $\lambda$
eigenvalue, and $PT$ represents the partial transposition of one of the
Hilbert spaces.

It is then possible to expand the $\hat{\mathcal{W}}$-operator in terms of a
complete basis of observables which can be implemented in a laboratory
setting. In the case of a spin $j$, these observables are the $2j+1$
dimensional irreducible representation of $SU(2)$, i.e. the generalized Gell-
Mann matrices. The measurements expansion and statistics in terms of Gell-Mann
matrices have already been studied for a qudit, and we refer to the Appendix
of Ref [62] for further explanation.

As in the previous section, the negativity of the presented protocol for
different spin states, spin values, and interaction time has been computed, to
find the optimal spin state which minimize the negativity (recall that the
negativity is a negative quantity). The behaviour of the negativity resembles
the one of the entanglement entropy. For this reason, the plots are presented
in Appendix D, while a summary of the results is given in Table 2.

Minimal Gravity-induced Negativity
---
Spin |  $j=1/2$ |  $j=2$ |  $j=5$ |  $j=10$
Coherent |  -0.44 |  -1.40 |  -2.02 |  -2.43
Superposition |  -0.44 |  -1.76 |  -2.61 |  -3.02
Table 2:  Minimal Negativity generated by optimizing families of spin states
interacting via gravity.

Some qualitative differences between the von Neumann entropy and the
negativity can be found from the plot in Appendix D. First of all, the
dependency on the spin value of the minimal negativity is more prominent. The
negativity monotonically decreases as function of $j$ toward an asymptotic
value, which saturates at higher spin values. Hence, when considering this
observable, the employment of larger spins become more preferable. Secondly,
there is no local maximum in negativity arising for the rotation
$\theta_{A}=\pi/2$ when the spin increases in value, and the optimal parameter
$\theta_{A}^{(o)}=\pi/2$ for any $j$ computed(in contrast with the
entanglement entropy, where for $j\geq 5$, a local minima appeared). Thirdly,
when considering superposition of CSS, the optimal $\delta\theta$ does not
tend to $\pi/2$, but decreases in value.

##  VIII Decoherence

There are many sources of decoherence. The main sources of decoherence were
already mentioned in [7]. Then, the analysis of [59, 44, 66, 62, 63] have
reinforced the constraints along with external jitters and gravity gradient
noise [67], gravitational decoherence [68, 69], dipole-induced decoherence
[70], decoherence/dephasing due to the plate [60], and gravity-induced
decoherence due to the apparatus [71]. These computations are based on earlier
analysis of decoherence in a quantum system [72, 73]. There are also sources
of decoherence due to internal phonon vibration [74, 75], and due to the
effect of rotation [42, 43].

The aim of this section is to give a first comparative analysis of decoherence
between different spins. Overall, sources of decoherence can be divided into
two groups: external decoherence (due to scattering of external particles) and
internal decoherence (due to internal structure of the mass). We will consider
the former but not the latter, as we assume the internal structure of the mass
to be independent of the spin. However, as we will see, larger spins increase
the total superposition and the Hilbert space, hence changing the behaviour of
the system under some sources of external decoherence.

Following Ref. [76], the decoherence due to scattering of external particles
with de Broglie wavelength $\lambda_{B}$ onto a mass in a superposition
$\Delta x$ can be further divided into two subgroups, namely, short-wavelength
and long-wavelength limits. In the case of the former ($\lambda_{B}\ll\Delta
x$), a single scattering event would resolve all the which-path information of
the superstition.131313For instance, at T = 3 K, the de Broglie wavelength of
$\text{O}_{2}$ air molecules is $\lambda_{\text{dB}}\approx 2\cdot 10^{-10}$m
which is much smaller compared to $\Delta x$ considered in the experiment.
While in the case of the latter ($\lambda_{B}\gg\Delta x$), multiple
scattering events are needed to resolve the information.141414For black-body
photons and $T\sim 1$ K, $\lambda_{\text{BB}}\sim 10^{-3}\gg\Delta x$.

Both long and short wavelength decoherence limits are compared for different
spins in regard to how the negativity is affected. To give a meaningful
comparison, the decoherence rates are kept as a variable (as they are
dependent on the specifics of the experiment - material of the mass,
temperature and pressure). The specifics on how scattering decoherence is
included in the model is given in Appendix E. Fig. 9 plots the Negativity as a
function of the rotation angle $\theta_{A}$ of the input CSS for different
decoherence rates.

In the case of short-wavelength, the optimal angle is still
$\theta_{A}^{(o)}=\pi/2$, for any decoherence rate $\gamma_{short}$. However,
in the case of long-wavelength limit, the optimal angle changes according to
the quantity $\Gamma_{long}\Delta x^{2}$.

(a) Short-wavelength

(b) Long-wavelength

Figure 9: Negativity as a function of the rotation angle and decoherence rate
for (a) short-wavelength and (b) long-wavelength limits. The optimal angle as
a function of the decoherence rate is plotted in green.

As was done previously, the negativity can be minimized over input spin states
at different decoherence rates. The results on the minimal negativity
generated by gravity with decoherence are plotted in Fig. 10.

(a) Short-wavelength

(b) Long-wavelength

Figure 10:  Minimal negativity for different spins values as a function of the
decoherence rate for (a) short-wavelength and (b) long-wavelength limits.

In the short-wavelength limit, larger spins are more resilient under
decoherence but, conversely, are more sensitive to long-wavelength decoherence
sources. This is expected as larger spins imply a larger total superposition
distance, which is linear in $j$. On one hand, the dependency on the
superposition distance in the long-wavelength limit leads to a worse
performances of large spins, which have a larger superposition. On the other
hand, in the short-wavelength limit, a larger superposition does not affect
the decoherence, as a scattering event would completely resolve the path taken
even for the smallest superposition ($j=1/2$). In an experimental scenario,
black-body noise (main source of long-wavelength decoherence) can be
significantly reduced by cooling the system and the environment with
temperature $T<1$ K, as the decoherence rate scales as $T^{9}$ [72, 77].
However, for molecules scattering (short-wavelength), a high-vacuum is
required, which may represent a bigger experimental challenge. With this in
mind we conclude that large spins are more resilient under the most
challenging source of external decoherence, namely, air molecules scattering.

##  IX Conclusion

An arbitrary spin embedded in a mesoscopic mass undergoes splitting in $2j+1$
paths when a magnetic gradient is applied, representing a generalization of
the Stern-Gerlach interferometer to an arbitrary spin. The total spatial
superposition size increases linearly in $j$ while the splitting time remains
constant. Therefore, the presented a dynamics exemplifyies a controlled
protocol to entangle any spin state with a spatial degree of freedom. This
creates a joint state which includes a spatial superposition of coherent
states, weighted accordingly to the spin distribution in the quantization axis
representation. Spin coherent states, superposition of these, and squeezed
states have been used in this study.

The entanglement generated by two adjacent GSG interferometers interacting via
gravity as quantified by the von Neumann entropy, witnesses the quantum nature
of the gravitational field. The aforementioned spin state families were varied
in order to find the optimal spin state which maximizes the gravitational
entanglement, and, similarly, which minimize the negativity both with and
without decoherence.

Such a procedure is first applied to the case of coherent spin states, which
are the most easily implementable spin states in experiments. Two
interferometer geometries are considered, namely, the linear and parallel set-
ups. With the former, no improvements are found, as only a small fraction of
the state is involved in the generation of the entanglement. In the latter,
the maximum entanglement entropy monotonically increases as a function of the
spin and reaches an asymptotic value for a given interaction time. A similar
behaviour is found for the minimal negativity, which monotonically decreases
with the size of the Hilbert space. For instance, by considering a spin
$j=10$, the entanglement entropy is increased of a factor of $2.2$ with
respect of the spin-$1/2$ case (and, equivalently, the negativity – an
entanglement monotone – decreased of a factor of $5.5$).

Superposition of coherent spin states can be used to further increase the
optimal entropy (negativity), which is increased of a factor of $2.5$ ($6.9$)
for $j=10$. Squeezed spin states, with both one-axis and two-axis twisted
Hamiltonian, do not significantly increase the generated entanglement entropy.
The experimental challenges of creating these last two families of spin states
may overcome the small increment in the entanglement entropy, compared to the
CSS.

The presented dynamics are of particular interest when considering NV-center
nano-diamonds. It is a known experimental challenge to create a single-center
nano-diamond, as multiple vacancies may be present. Until now, the dynamics of
a multiple NV-centers mass in a magnetic gradient was unknown. With the
assumption that the coupling with the magnetic gradient is the same among
every center and neglecting the spin-spin coupling, the mass will evolve
according to the proposed GSG interferometer.

The challenges of creating a mesoscopic superposition in an experimental
setting are numerous and further studies should be carried out. The
diamagnetic response from a mass in a magnetic gradient behaves as a trap
(which allowed the presented interferometer) but also significantly increases
the splitting time. As such a phenomenon is present both in this study, as
well as in the original protocol, we have, in a sense, compared them here in
an equal footing. Further works on low decoherence environments are encouraged
in order to create such a large and massive superposition. An alternative
solution to this issue may be achieved from a material science prospective
with the design of a low diamagnetic susceptibility mesoscopic mass.

To conclude, it is of great importance to experimentally study quantum aspects
of the gravitational field. The experimental challenge of reaching mesoscopic
masses in a quantum regime demands theoretical developments to find protocols
that are more sensible to quantum phenomena of gravity. The question that we
must pose is: Does the employment of large spins – i.e. a $2j+1$ spatial
superposition of a mass weighted according to a spin distribution – give a
significant improvement to make quantum natured gravity-induced entanglement
detection achievable in an experiment? While improvements in the entanglement
and a witness were found, we conclude that the entanglement cannot be
increased of the orders of magnitude required to conduct the experiment with
the current quantum technologies. However, once significant developments
regarding noise isolation and splitting dynamics are achieved, the present
protocol may lead to the necessary increment which may distinguish detectable
and undetectable entanglement signals.

## Acknowledgements

L.B. would like to aknowledge Engineering and Physical Sciences Research
Council (EPSRC) grants (EP/R513143/1 and EP/W524335/1), supported, in part, by
Perimeter Institute for Theoretical Physics (Research at Perimeter Institute
is supported by the Government of Canada through the Department of Innovation,
Science and Economics Development and by the Province of Ontario through the
Ministry of Colleges and University). S.B. would like to acknowledge EPSRC
grants (EP/N031105/1, EP/S000267/1, and EP/X009467/1) and grant ST/W006227/1.
We thank the Center for Information Technology of the University of Groningen
for their support and for providing access to the Hábrók high performance
computing cluster.

## References

  * Dyson [2013] F. Dyson, Is a graviton detectable?, International Journal of Modern Physics A 28 (2013), 10.1142/S0217751X1330041X.
  * Rothman and Boughn [2006] T. Rothman and S. Boughn, Can Gravitons be Detected?, Foundations of Physics 36 (2006), 10.1007/s10701-006-9081-9.
  * Keil  _et al._ [2021] M. Keil, S. Machluf, Y. Margalit, Z. Zhou, and et al, in  _Molecular Beams in Physics and Chemistry: From Otto Stern’s Pioneering Exploits to Present-Day Feats_ (Springer International Publishing, 2021) Chap. Stern-Gerlach Interferometry with the Atom Chip, pp. 263–301, 10.1007/978-3-030-63963-1-14.
  * Bennett  _et al._ [1996] C. H. Bennett, D. P. DiVincenzo, J. A. Smolin, and W. K. Wootters, Mixed state entanglement and quantum error correction, Phys. Rev. A 54 (1996), 10.1103/PhysRevA.54.3824.
  * Horodecki  _et al._ [2009] R. Horodecki, P. Horodecki, M. Horodecki, and K. Horodecki, Quantum entanglement, Rev. Mod. Phys. 81 (2009), 10.1103/RevModPhys.81.865.
  * Bose [2016] S. Bose, Matter wave ramsey interferometry and the quantum nature of gravity, Available at https://www.youtube.com/watch?v=0Fv-0k13s_k (2016), fundamental Problems of Quantum Physics, ICTS, Bangalore.
  * Bose  _et al._ [2017] S. Bose, A. Mazumdar, G. W. Morley, H. Ulbricht, M. Toroš, M. Paternostro, A. A. Geraci, P. F. Barker, M. S. Kim, and G. Milburn, Spin Entanglement Witness for Quantum Gravity, Phys. Rev. Lett. 119 (2017), 10.1103/PhysRevLett.119.240401.
  * Marletto and Vedral [2017] C. Marletto and V. Vedral, Gravitationally-induced entanglement between two massive particles is sufficient evidence of quantum effects in gravity, Phys. Rev. Lett. 119 (2017), 10.1103/PhysRevLett.119.240402.
  * Marshman  _et al._ [2020] R. J. Marshman, A. Mazumdar, and S. Bose, Locality & Entanglement in Table-Top Testing of the Quantum Nature of Linearized Gravity, Phys. Rev. A 101 (2020), 10.1103/PhysRevA.101.052110.
  * Bose  _et al._ [2022] S. Bose, A. Mazumdar, M. Schut, and M. Toroš, Mechanism for the quantum natured gravitons to entangle masses, Phys. Rev. D 105 (2022), 10.1103/PhysRevD.105.106028.
  * Hensen  _et al._ [2015] B. Hensen _et al._ , Loophole-free Bell inequality violation using electron spins separated by 1.3 kilometres, Nature 526 (2015), 10.1038/nature15759.
  * Christodoulou, Marios and Rovelli, Carlo [2019] Christodoulou, Marios and Rovelli, Carlo, On the possibility of laboratory evidence for quantum superposition of geometries, Phys. Rev. B 792 (2019), 10.1016/j.physletb.2019.03.015.
  * Belenchia  _et al._ [2018] A. Belenchia _et al._ , Quantum Superposition of Massive Objects and the Quantization of Gravity, Phys. Rev. D 98 (2018), 10.1103/PhysRevD.98.126009.
  * Danielson  _et al._ [2022] D. L. Danielson, G. Satishchandran, and R. M. Wald, Gravitationally mediated entanglement: Newtonian field versus gravitons, Phys. Rev. D 105 (2022), 10.1103/PhysRevD.105.086001.
  * Carney, Daniel [2022] Carney, Daniel, Newton, entanglement, and the graviton, Phys. Rev. D 105 (2022), 10.1103/PhysRevD.105.024029.
  * Christodoulou, Marios and Di Biagio, Andrea and Aspelmeyer, Markus and Brukner, Časlav and Rovelli, Carlo and Howl, Richard [2022] Christodoulou, Marios and Di Biagio, Andrea and Aspelmeyer, Markus and Brukner, Časlav and Rovelli, Carlo and Howl, Richard, Locally mediated entanglement through gravity from first principles, Phys. Rev. Lett. 130 (2022), 10.1103/PhysRevLett.130.100202.
  * 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 (2022), q-2022-08-17-779.
  * Christodoulou  _et al._ [2023] M. Christodoulou, A. D. Biagio, R. Howl, and C. Rovelli, Gravity entanglement, quantum reference systems, degrees of freedom, Classical and Quantum Gravity 6 (2023), 10.1088/1361-6382/acb0aa.
  * Chen  _et al._ [2023] L.-Q. Chen, F. Giacomini, and C. Rovelli, Quantum States of Fields for Quantum Split Sources, Quantum 7 (2023), q-2023-03-20-958.
  * Vinckers  _et al._ [2023] U. K. B. Vinckers, A. de la Cruz-Dombriz, and A. Mazumdar, Quantum entanglement of masses with nonlocal gravitational interaction, Phys. Rev. D 107 (2023), 10.1103/PhysRevD.107.124036.
  * Chakraborty  _et al._ [2023] S. Chakraborty, A. Mazumdar, and R. Pradhan, Distinguishing Jordan and Einstein frames in gravity through entanglement (2023), arXiv:2310.06899 [gr-qc] .
  * Biswas  _et al._ [2023] D. Biswas, S. Bose, A. Mazumdar, and M. Toroš, Gravitational Optomechanics: Photon-Matter Entanglement via Graviton Exchange, Phys. Rev. D 108 (2023), 10.1103/PhysRevD.108.064023.
  * the possibility of experimental detection of the discreteness of time [2020] O. the possibility of experimental detection of the discreteness of time, Christodoulou, marios and rovelli, carlo, Front. Phys. 8 (2020), 10.3389/fphy.2020.00207.
  * experiment to test the discreteness of time [2022] A. experiment to test the discreteness of time, Christodoulou, marios and di biagio, andrea and martin-dussaud, pierre, Quantum 6 (2022), 10.22331/q-2022-10-06-826.
  * Bose  _et al._ [2023] S. Bose, A. Mazumdar, M. Schut, and M. Toroš, Entanglement witness for the weak equivalence principle, Entropy 6 (2023), 10.3390/e25030448.
  * Barker  _et al._ [2022] P. F. Barker, S. Bose, R. J. Marshman, and A. Mazumdar, Entanglement based tomography to probe new macroscopic forces, Phys. Rev. D 106 (2022), 10.1103/PhysRevD.106.L041901.
  * Elahi and Mazumdar [2023] S. G. Elahi and A. Mazumdar, Probing massless and massive gravitons via entanglement in a warped extra dimension, Phys. Rev. D 108 (2023), 10.1103/PhysRevD.108.035018.
  * Al  _et al._ [2018] A. Al, W. Cong, and R. B. Mann, Optomechanical quantum cavendish experiment, Phys. Rev. A 98 (2018), 10.1103/PhysRevA.98.043811.
  * Miao  _et al._ [2020] H. Miao, D. Martynov, H. Yang, and A. Datta, Quantum correlations of light mediated by gravity, Phys. Rev. A 101 (2020), 10.1103/PhysRevA.101.063804.
  * Matsumura and Yamamoto [2020] A. Matsumura and K. Yamamoto, Gravity-induced entanglement in optomechanical systems, Phys. Rev. D 101 (2020), 10.1103/PhysRevD.102.106021.
  * Haine [2021] S. A. Haine, Searching for signatures of quantum gravity in quantum gases, New Journal of Physics 23 (2021), 10.1088/1367-2630/abd97d.
  * Howl  _et al._ [2021] R. Howl, V. Vedral, D. Naik, M. Christodoulou, C. Rovelli, and A. Iyer, Non-Gaussianity as a Signature of a Quantum Theory of Gravity, PRX Quantum 2 (2021), 10.1103/PRXQuantum.2.010325.
  * 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), arXiv, arXiv:2307.08133 [gr-qc] .
  * Degen  _et al._ [2017] C. L. Degen, F. Reinhard, and P. Cappellaro, Quantum sensing, Reviews of Modern Physics 89 (2017), 10.1103/RevModPhys.89.035002.
  * Ma  _et al._ [2011] J. Ma, X. Wang, C. P. Sun, and F. Nori, Quantum spin squeezing, Physics Reports 509 (2011), 10.1016/j.physrep.2011.08.003.
  * Maleki  _et al._ [2021] Y. Maleki, M. O. Scully, and A. M. Zheltikov, Quantum metrology with superposition spin coherent states: Insights from Fisher information, Phys. Rev. A 104 (2021), 10.1103/PhysRevA.104.053712.
  * Rahman [2019] A. T. M. A. Rahman, Large spatial Schrodinger cat using a levitated ferrimagnetic nanoparticle, New Journal of Physics 21 (2019), 10.1103/PhysRevA.104.053712.
  * Dicke [1954] R. H. Dicke, Coherence in Spontaneous Radiation Processes, Phys. Rev. Journals Archive 93 (1954), 10.1103/PhysRev.93.99.
  * Schliemann [2015] J. Schliemann, Coherent Quantum Dynamics: What Fluctuations Can Tell, Phys. Rev. A 92 (2015), 10.1103/PhysRevA.92.022108.
  * Wang  _et al._ [2022] Z.-Q. Wang, Y.-P. Wang, J. Yao, R.-C. Shen, W.-J. Wu, J. Qian, J. Li, S.-Y. Zhu, and J. Q. You, Giant spin ensembles in waveguide magnonics, Nature Communications 13 (2022), 10.1038/s41467-022-35174-9.
  * Delić  _et al._ [2020] U. Delić, M. Reisenbauer, K. Dare, D. Grass, V. Vuletić, N. Kiesel, and M. Aspelmeyer, Cooling of a levitated nanoparticle to the motional quantum ground state, Science 367 (2020), 10.1126/science.aba3993.
  * Japha and Folman [2023] Y. Japha and R. Folman, Quantum Uncertainty Limit for Stern-Gerlach Interferometry with Massive Objects, Phys. Rev. Lett. 130 (2023), 10.1103/PhysRevLett.130.113602.
  * Ma  _et al._ [2021] Y. Ma, M. S. Kim, and B. A. Stickler, Torque-free manipulation of nanoparticle rotations via embedded spins, Phys. Rev. B 104 (2021), 10.1103/PhysRevB.104.134310.
  * Pedernales  _et al._ [2020] J. S. Pedernales, G. W. Morley, and M. B. Plenio, Motional Dynamical Decoupling for Interferometry with Macroscopic Particles, Phys. Rev. Lett. 125 (2020), 10.1103/PhysRevLett.125.023602.
  * Marshman  _et al._ [2022] R. J. Marshman, A. Mazumdar, R. Folman, and S. Bose, Constructing nano-object quantum superpositions with a Stern-Gerlach interferometer, Phys. Rev. Res. 4 (2022), 10.1103/PhysRevResearch.4.023087.
  * Zhou  _et al._ [2022a] R. Zhou, R. J. Marshman, S. Bose, and A. Mazumdar, Catapulting towards massive and large spatial quantum superposition, Phys. Rev. Res. 4, 043157 (2022a), arXiv:2206.04088 [quant-ph] .
  * Zhou  _et al._ [2023] R. Zhou, R. J. Marshman, S. Bose, and A. Mazumdar, Mass-independent scheme for enhancing spatial quantum superpositions, Phys. Rev. A 107, 032212 (2023), arXiv:2210.05689 [quant-ph] .
  * Zhou  _et al._ [2022b] R. Zhou, R. J. Marshman, S. Bose, and A. Mazumdar, Gravito-diamagnetic forces for mass independent large spatial superpositions (2022b), arXiv:2211.08435 [quant-ph] .
  * Englert  _et al._ [1988] B.-G. Englert, J. Schwinger, and M. O. Scully, Is spin coherence like Humpty-Dumpty? I. Simplified treatment, Foundations of Physics 18 (1988), 10.1007/BF01909939.
  * Amit  _et al._ [2019] O. Amit _et al._ , ${T}^{3}$ stern-gerlach matter-wave interferometer, Phys. Rev. Lett. 123 (2019), 10.1103/PhysRevLett.123.083601.
  * Margalit  _et al._ [2021] Y. Margalit _et al._ , Realization of a complete Stern-Gerlach interferometer: Towards a test of quantum gravity, Science 7 (2021), 10.1126/sciadv.abg2879.
  * Machluf  _et al._ [2013] S. Machluf, Y. Japha, and R. Folman, Coherent stern–gerlach momentum splitting on an atom chip, Nature communications 4 (2013), 10.1126/sciadv.abg2879.
  * Radcliffe [1971] J. M. Radcliffe, Some properties of coherent spin states, Journal of Physics A: General Physics 4 (1971), 10.1088/0305-4470/4/3/009.
  * world arising out of quantum physics under the restriction of coarse-grained measurements [2007] C. world arising out of quantum physics under the restriction of coarse-grained measurements, Kofler, johannes and brukner, caslav, Phys. Rev. Lett. 99 (2007), 10.1103/PhysRevLett.99.180403.
  * Kitagawa and Ueda [1993] M. Kitagawa and M. Ueda, Squeezed spin states, Phys. Rev. A 47 (1993), 10.1103/PhysRevA.47.5138.
  * Casimir and Polder [1948] H. B. G. Casimir and D. Polder, The Influence of retardation on the London-van der Waals forces, Phys. Rev. 73 (1948), 10.1103/PhysRev.73.360.
  * Casimir [1948] H. B. G. Casimir, On the Attraction Between Two Perfectly Conducting Plates, Indag. Math. 10 (1948).
  * Kim  _et al._ [2005] H.-Y. Kim, J. O. Sofo, D. Velegol, M. W. Cole, and G. Mukhopadhyay, Static polarizabilities of dielectric nanoclusters, Phys. Rev. A 71 (2005), 10.1103/PhysRevA.72.053201.
  * van de Kamp  _et al._ [2020] T. W. van de Kamp, R. J. Marshman, S. Bose, and A. Mazumdar, Quantum Gravity Witness via Entanglement of Masses: Casimir Screening, Phys. Rev. A 102 (2020), 10.1103/PhysRevA.102.062807.
  * Schut  _et al._ [2023] M. Schut, A. Grinin, A. Dana, S. Bose, A. Geraci, and A. Mazumdar, Relaxation of experimental parameters in a quantum-gravity-induced entanglement of masses protocol using electromagnetic screening, Phys. Rev. Res. 5 (2023), 10.1103/PhysRevResearch.5.043170.
  * 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 (2020), 10.1140/epjd/e2020-10077-8.
  * Tilly  _et al._ [2021] J. Tilly, R. J. Marshman, A. Mazumdar, and S. Bose, Qudits for witnessing quantum-gravity-induced entanglement of masses under decoherence, Phys. Rev. A 104 (2021), 10.1103/PhysRevA.104.0524168.
  * Schut  _et al._ [2022] M. Schut, J. Tilly, R. J. Marshman, S. Bose, and A. Mazumdar, Improving resilience of quantum-gravity-induced entanglement of masses to decoherence using three superpositions, Phys. Rev. A 105 (2022), 10.1103/PhysRevA.105.032411.
  * Zyczkowski  _et al._ [1998] K. Zyczkowski, P. Horodecki, A. Sanpera, and M. Lewenstein, Volume of the set of separable states, Phys. Rev. A 58 (1998), 10.1103/physreva.58.883.
  * B. [2005] P. M. B., Logarithmic negativity: A full entanglement monotone that is not convex, Phys. Rev. Lett. 95 (2005), 10.1103/PhysRevLett.95.090503.
  * Rijavec  _et al._ [2021] S. Rijavec, M. Carlesso, A. Bassi, V. Vedral, and C. Marletto, Decoherence effects in non-classicality tests of gravity, New Journal of Physics 23 (2021), 10.1088/1367-2630/abf3eb.
  * Toroš  _et al._ [2021] M. Toroš, T. W. Van De Kamp, R. J. Marshman, M. S. Kim, A. Mazumdar, and S. Bose, Relative acceleration noise mitigation for nanocrystal matter-wave interferometry: Applications to entangling masses via quantum gravity, Phys. Rev. Res. 3 (2021), 10.1103/PhysRevResearch.3.023178.
  * Toroš  _et al._ [2020] M. Toroš, A. Mazumdar, and S. Bose, Loss of coherence of matter-wave interferometer from fluctuating graviton bath (2020), 10.48550/arXiv.2008.08609.
  * Wu  _et al._ [2021] M.-Z. Wu, M. Toroš, S. Bose, and A. Mazumdar, Quantum gravitational sensor for space debris, Phys. Rev. D 107 (2021), 10.1103/PhysRevD.107.104053.
  * Fragolino  _et al._ [2023] P. Fragolino, M. Schut, M. Toroš, S. Bose, and A. Mazumdar, Decoherence of a matter-wave interferometer due to dipole-dipole interactions (2023), 2307.07001.
  * Gunnink  _et al._ [2023] F. Gunnink, A. Mazumdar, M. Schut, and M. Toroš, Gravitational decoherence by the apparatus in the quantum-gravity induced entanglement of masses, Classical and Quantum Gravity 40 (2023), 10.1088/1361-6382/ad0496.
  * Romero-Isart [2011] O. Romero-Isart, Quantum superposition of massive objects and collapse models, Phys. Rev. A 84 (2011), 10.1103/physreva.84.052121.
  * Chang  _et al._ [2009] D. E. Chang, C. A. Regal, S. B. Papp, D. J. Wilson, J. Ye, O. Painter, H. J. Kimble, and P. Zoller, Cavity opto-mechanics using an optically levitated nanosphere, Proceedings of the National Academy of Sciences 107 (2009), 10.1073/pnas.0912969107.
  * Henkel and Folman [2022] C. Henkel and R. Folman, Internal decoherence in nano-object interferometry due to phonons, AVS Quantum Science 4 (2022), 10.1116/5.0080503.
  * Henkel and Folman [2023] C. Henkel and R. Folman, Universal limit on spatial quantum superpositions with massive objects due to phonons (2023), 2305.15230.
  * Schlosshauer [2007] M. Schlosshauer,  _The quantum-to-classical transition and decoherence_ (Springer International Publishing, 2007) 10.1007/978-3-540-35775-9.
  * Schlosshauer [2019] M. Schlosshauer, Quantum Decoherence, Phys. Rept. 831 (2019), 10.1016/j.physrep.2019.10.001.
  * Heremans  _et al._ [1994] J. Heremans, C. H. Olk, and D. T. Morelli, Magnetic susceptibility of carbon structures, Phys. Rev. B 49 (1994), 10.1103/PhysRevB.49.15122.
  * Sakurai [1994] J. J. Sakurai,  _Modern Quantum Mechanics_ (Addison-Wesley, Reading, MA, 1994).
  * Bose  _et al._ [1997] S. Bose, K. Jacobs, and P. L. Knight, Preparation of nonclassical states in cavities with a moving mirror, Phys. Rev. A 56 (1997), 10.1103/PhysRevA.56.4175.

##  Appendix A Diamagnetic Hamiltonian

A uniform magnetic field applied to a diamagnetic material leads to a
diamagnetic response. For a suitable choice of magnetic profile, this can be
treated as a quantum harmonic oscillator trapping the mass, where the
frequency of the magnetic field is related to the frequency of the harmonic
oscillator. Additionally, the spin embedded in the diamagnetic material will
couple to a magnetic field gradient in a Stern-Gerlach-like interaction. For a
general magnetic field $\bf{B}$, the Hamiltonian of the diamagnetic system
reads:

$\hat{H}=\frac{\hat{\bf{p}}^{2}}{2M}-\frac{M}{2}\frac{\chi_{m}}{\mu_{0}}\hat{\bf{B}}^{2}+\tilde{g}\mu_{B}\hat{\bf{J}}\cdot\hat{\bf{B}}\,,$

where $\hat{\bf{p}}$ and $\hat{\bf{J}}$ are the 3-dimensional momentum and
spin operators of the mass M, respectively, $\hat{\bf{B}}$ is the magnetic
field operator, $\chi_{m}$ the mass (density) magnetic susceptibility,
$\mu_{0}$ the vacuum permeability, $\mu_{B}$ is the Bohr magneton and
$\tilde{g}$ is the Landé g-factor.

By choosing the magnetic gradient of the form
$\bm{B}=\left(\partial_{x}B\right)x\bm{e}_{z}$ (following [44]), where
$\partial_{x}B$ is the strength of the magnetic gradient and $\bm{e}_{z}$ is
the unit vector align with the $z$ spin axis, the Hamiltonian on the $x$-axis
(the others axes have free evolution) can be further simplifies to:

$\displaystyle\hat{H}_{x}=\frac{\hat{p}_{x}^{2}}{2M}-\frac{1}{2}M\frac{\chi_{m}\left(\partial_{x}B\right)^{2}}{\mu_{0}}\hat{x}^{2}+\tilde{g}\mu_{B}\left(\partial_{x}B\right)\hat{J}_{z}\hat{x}\,,$

where we restricted the Hamiltonian to the splitting axis, $\hat{x}$.

Treating the trapped systems as a quantum harmonic oscillator, we can
determine the frequency of the oscillation. Consequently we determine the
coupling constant $g$ used in Eq. 3, which was defined as the coefficient of
the term $\hbar\hat{J}_{z}(\hat{a}+\hat{a}^{\dagger})$ in the Hamiltonian,
where the creation and annihilation operators come from the definition
$\hat{x}=\sqrt{\hbar/(2m\omega)}(\hat{a}+\hat{a}^{\dagger})$.

$\omega_{M}=\sqrt{\frac{|\chi_{m}|}{\mu_{0}}}\left(\partial_{x}B\right)\;\;\;\;\implies\;\;\;\;g=\tilde{g}\mu_{B}\sqrt{\frac{\left(\partial_{x}B\right)}{2\hbar
M}\sqrt{\frac{|\chi_{m}|}{\mu_{0}}}}$

and the Hamiltonian becomes the one given in Eq. (2). Finally, the
superposition distance (see Eq. 9) and splitting time can be computed:

$t_{s}\coloneqq\frac{\pi}{\omega_{M}}=\frac{\pi}{\left(\partial_{x}B\right)}\sqrt{\frac{\mu_{0}}{|\chi_{m}|}}\,,\;\;\;\;\;\;\;\;\;\;\Delta
x=\frac{2\tilde{g}\mu_{B}\mu_{0}}{M|\chi_{M}|\left(\partial_{x}B\right)}$

Considering a diamond with $|\chi_{M}|\approx 6\cdot
10^{-9}\,\mbox{m}^{3}/\mbox{kg}$ [78] and $M\approx 10^{-14}$ kg:

$t_{s}\approx
44.7\left|\left(\partial_{x}B\right)\right|^{-1}\mbox{s}\,,\;\;\;\;\;\;\;\;\;\;\Delta
x=7.5\cdot 10^{-7}\left|\left(\partial_{x}B\right)\right|^{-1}\mbox{m}\,.$

In order to reach the wanted superposition of $\Delta x\approx 2.50\cdot
10^{-4}$m, a magnetic gradient $\left(\partial_{x}B\right)\approx 3.0\cdot
10^{-3}\text{T/m}$, which implies a splitting time of $t_{s}\approx 1.4\cdot
10^{4}$s. This is the well-known diamagnetic issue (also present in the
original proposal): a large superposition requires a weak magnetic gradient
which significantly increases the splitting time.

##  Appendix B Splitting Unitary Evolution

Using time in frequency unity, the Hamiltonian in Eq. (2) gives the evolution:

$U_{I}\left(t\right)=e^{-ita^{\dagger}a+itkJ_{z}\left(a^{\dagger}+a\right)}$
(16)

where $r=\frac{\omega_{0}}{\omega_{M}}$ and $k=\frac{g}{\omega_{M}}$. Using
Schwinger’s Oscillator model of angular momentum, the spin algebra can be
expressed in term of the algebra of two uncoupled quantum harmonic oscillators
[79]. Each quantum harmonic oscillator is associated with ladder operators
$b_{+}$ and $b_{+}^{\dagger}$ ($b_{-}$ and $b_{-}^{\dagger}$) creating and
annihilating a particle with $m=\frac{1}{2}$ ($m=-\frac{1}{2}$). In this
formalism, the projection of the angular momentum on the $z$-axis is equal to
half the difference between of quanta in the $+$ state and in the $-$ state,
$J_{z}=\frac{1}{2}\left(b_{+}b_{+}^{\dagger}-b_{-}b_{-}^{\dagger}\right)$.
Applying this to Eq. (16) recovers the unitary evolution studied in [80], with
two commuting quantum harmonic oscillators instead of one. For completeness,
we will derive the simplified unitary evolution. Following [80], let us define
the spin-displacement operator:

$T=e^{-\frac{k}{2}\left(b_{+}^{\dagger}b_{+}-b_{-}^{\dagger}b_{-}\right)\left(a^{\dagger}-a\right)}\,.$

Using the Baker-Cambell-Hausdorf (BCH) expansion, the following identities can
be derived:

$\displaystyle TaT^{\dagger}$
$\displaystyle=a+\frac{k}{2}\left(b_{+}^{\dagger}b_{+}-b_{-}^{\dagger}b_{-}\right)\,,$
$\displaystyle Ta^{\dagger}T^{\dagger}$
$\displaystyle=a^{\dagger}+\frac{k}{2}\left(b_{+}^{\dagger}b_{+}-b_{-}^{\dagger}b_{-}\right)\,,$
$\displaystyle Tb_{\pm}^{\dagger}b_{\pm}T^{\dagger}$
$\displaystyle=b_{\pm}^{\dagger}b_{\pm}\,.$

The last one implies that $TJ_{z}T^{\dagger}=J_{z}$. Recalling that
$Uf\left(\\{X_{i}\\}\right)U^{\dagger}=f\left(\\{UX_{i}U^{\dagger}\\}\right)$,
for any $U$ unitary:

$TU_{I}\left(t\right)T^{\dagger}=e^{itk^{2}J_{z}^{2}}e^{-ita^{\dagger}a}\,.$

Multiplying on the left (right) by $T^{\dagger}$ $(T)$, the unitary evolution
can be rewritten as:

$U_{I}\left(t\right)=e^{itk^{2}J_{z}^{2}}e^{-kJ_{z}^{\dagger}\left(a-a^{\dagger}\right)}e^{-ita^{\dagger}a}e^{-kJ_{z}\left(a^{\dagger}-a\right)}\,.$

Furthermore, noticing that

$e^{-ita^{\dagger}a}\left(J_{z}(a^{\dagger}-a)\right)e^{ita^{\dagger}a}=J_{z}\left(a^{\dagger}e^{-it}-ae^{it}\right)\,,$

the unitary can be further simplified to:

$U_{I}\left(t\right)=e^{itk^{2}J_{z}^{2}}e^{-kJ_{z}\left(a-a^{\dagger}\right)}e^{-kJ_{z}\left(a^{\dagger}e^{-it}-ae^{it}\right)}e^{-ita^{\dagger}a}\,.$

Finally, applying again the BCH expansion:

$U_{I}\left(t\right)=e^{ik^{2}J_{z}^{2}\left(t-\sin(t)\right)}e^{kJ_{z}\left(\eta(t)a^{\dagger}-\eta(t)^{*}a\right)}e^{-ita^{\dagger}a}\,,$

where the function $\eta(t)=1-e^{-it}$ was defined. We note the
$\eta(t)b^{\dagger}-\eta(t)^{*}b$ component of the unitary evolution, is
typical of a time-dependent displacement operators in the position basis,
$\mathcal{D}\left(\eta(t)\right)=e^{\left(\eta(t)b^{\dagger}-\eta(t)^{*}b\right)}$
. Thus, recalling that one of the definitions of position coherent state is
the displacement of the vacuum, i.e.
$\ket{\alpha}\coloneqq\mathcal{D}(\alpha)\ket{0}$, the quantum state in Eq.
(4) follows trivially.

##  Appendix C Additional Entanglement Results

In this Appendix, some of the plots regarding the vonNeumann entropy generated
by gravity are presented. In Fig. 11, the maximized entanglement entropy for
input CSSs is plotted as function of the total anglular momentum $j$. In Fig.
12a, the entanglement entropy is plotted as function of $\delta\theta$ and
$\delta\phi$ for a superposition of CSSs defined according to Eq. (15). In
Fig. 12b and 13, the entanglement entropy for SSS is ploted.

Figure 11:  Maximum entanglement entropy generated by optimal CSS
$\ket{\theta_{A}^{(o)}}$ for (a) the liner and (b) the parallel set-ups, as a
function of the spin $j$, at $\tau=2s$

(a) Superposition of CSS

(b) Squeezed Spin State

Figure 12:  Gravitational entanglement of (a) superpositions of CSS (Eq. (15),
$J=5$) as a function of $\delta\theta$ and $\delta\phi$ with $J=5$ and (B)
squeezed spin state ($J=10$) as a function of the squeezing parameters
$\chi_{i}$ and the rotation angle $\theta$ for both (i) one-axis twisting, and
(ii) two axis twisting. Figure 13:  Gravitational Entanglement as a function
of the squeezing parameter $\chi_{i}$ of CSS input state with $\theta=\pi/2$,
for both (i) one-axis twisting and (ii) two axis twisting.

##  Appendix D Negativity Numerical Results

In this Appendix, the results regarding the negativity generated by the
gravitational interaction for different spins lengths and states are reported.
For the related discussion, see Sc. VII.

(a) J = 5

(b) J = 10

Figure 14: Negativity for Spin Coherent States with two angles parameters

(a) Negativity as a function of $\theta_{A}$

(b) Minimal Negativity as a function of J

(c) Minimal Negativity as a function of $\tau$

Figure 15: Negativity for Spin Coherent States with one angles parameters

(a) Two-parameters Negativity: J = 5

(b) Two-parameters Negativity: J = 10

(c) One-parameter Negativity:

Figure 16: Negativity for Superposition of CSS, as function of the two
parameters $\delta\phi$ and $\delta\theta$ for (a) $J=5$, and (b) $J=10$, (c)
and as function of $\delta\theta$ for different spins values.

##  Appendix E Including Decoherence in the Model and Sources of Decoherence

The general evolution due to decoherence caused by random particle scattering
on a system in a superposition is described by the Master equation: [76, 77],

$\frac{\partial\rho(\mathbf{x},\mathbf{x}^{\prime},t)}{\partial
t}=-F(\mathbf{x}-\mathbf{x}^{\prime})\rho(\mathbf{x},\mathbf{x}^{\prime},t)$

where $F$ plays the role of decoherence factor and, generally, it is given by:
[76, 77]

$F(\mathbf{x}-\mathbf{x}^{\prime})=\int
dq\rho(q)v(q)\int\frac{d\hat{\mathbf{n}}d\hat{\mathbf{n}}^{\prime}}{4\pi}\left(1-e^{iq(\hat{\mathbf{n}}-\hat{\mathbf{n}}^{\prime})(\mathbf{x}-\mathbf{x}^{\prime})/\hbar}\right)\left|f(q\hat{\mathbf{n}},q\hat{\mathbf{n}}^{\prime})\right|^{2}$
(17)

where $q$ is the normalized momentum of the incoming particle, $\rho(q)$ the
momentum-space density, $v(q)$ the speed of particles with momentum $q$,
$f(q\hat{\mathbf{n}},q\hat{\mathbf{n}}^{\prime})$ the scattering
amplitude151515Defined according to:
$\braket{q}{\hat{T}}{q^{\prime}}=\frac{i}{2\pi\hbar
q}\delta\left(q-q^{\prime}\right)f(\mathbf{q},\mathbf{q}^{\prime})$. With $T$
the scattering part of the $S$-matrix.. It is possible to simplify the above
equation taking the short-wavelength limit, where the wavelength of the
environmental particles $\lambda$ is much smaller than the superposition size,
$\lambda\ll\Delta x$ (using the expression of the wavelength in terms of its
momentum this is equivalent to $\frac{q|\mathbf{x}-\mathbf{x}^{\prime}|}{h}\gg
1$), or the long-wavelength limit, where $\lambda\gg\Delta x$ (or equivalently
$\frac{q|\mathbf{x}-\mathbf{x}^{\prime}|}{h}\ll 1$). Within these limits, the
integrals over $\mathbf{n},\mathbf{n}^{\prime}$ in Eq. 17 can be simplified
and evaluated. The evolution of the off-diagonal terms of the density matrix
can be described in the short-and long-wavelength limit, respectively, as:
[76]

$\frac{\partial\rho(\mathbf{x},\mathbf{x}^{\prime},t)}{\partial
t}=-\gamma_{\text{short}}\rho(\mathbf{x},\mathbf{x}^{\prime},t)\,,\;\;\;\;\;\;\;\frac{\partial\rho(\mathbf{x},\mathbf{x}^{\prime},t)}{\partial
t}=-\Gamma_{\text{long}}(\mathbf{x}-\mathbf{x}^{\prime})^{2}\rho(\mathbf{x},\mathbf{x}^{\prime},t)\,.$

where $\gamma_{short}$ and $\Gamma_{long}$ are decoherence rates (the former
in units of Hz, while the latter of Hz/m2), and are computed via scattering
theory according to the process that they describe. In order to compare how
these processes effect the negativity for different spins values, the
decoherence rates will be kept as a variable.

Each subsystem of the bipartite system describing the masses is assumed to
evolve under decoherence uncorrelated to the other subsystem, i.e. the
scattering particles do not scatter with both masses. Furthermore, it is
assumed that the decoherence rate is the same for both subsystems, i.e. both
masses are surrounded by the same environment and have the equivalent
scattering properties. Under these assumptions and after a time $\tau$, the
density matrix of the composite system under decoherence in the short-
wavelenght limit:

$\rho_{m,m^{\prime},n,n^{\prime}}^{(\text{short})}=\rho_{m,m^{\prime},n,n^{\prime}}e^{-(2-\delta_{mm^{\prime}}-\delta_{nn^{\prime}})\gamma_{\text{short}}\tau}.$

For the short-wavelength, and were $m$ and $m^{\prime}$ ($n$ and $n^{\prime}$)
are the spin quantum number of subsystem A (B). For the long-wavelength:

$\rho_{m,m^{\prime},n,n^{\prime}}^{(\text{long})}=\rho_{m,m^{\prime},n,n^{\prime}}e^{-\Gamma_{\text{long}}(\Delta
x)^{2}\left[\left(m-m^{\prime}\right)^{2}+\left(n-n^{\prime}\right)^{2}\right]\tau}$

The decoherence terms cause an exponential decay of the off-diagonal terms of
the density matrix. After a certain time all off-diagonal terms will go to
zero, and the density matrix becomes a classical mixed state.

It is possible to note that in the spin $1/2$ case, the two processes are
described by the density matrix:

$\rho_{1/2}^{(\text{short/long})}=\begin{pmatrix}\rho_{00,00}&\rho_{00,10}e^{-\gamma\tau}&\rho_{10,00}e^{-\gamma\tau}&\rho_{10,10}e^{-2\gamma\tau}\\\
\rho_{00,01}e^{-\gamma\tau}&\rho_{00,11}&\rho_{10,01}e^{-2\gamma
t}&\rho_{10,11}e^{-\gamma\tau}\\\
\rho_{01,00}e^{-\gamma\tau}&\rho_{01,10}e^{-2\gamma\tau}&\rho_{11,00}&\rho_{11,10}e^{-\gamma
t}\\\
\rho_{01,01}e^{-2\gamma\tau}&\rho_{01,11}e^{-\gamma\tau}&\rho_{11,01}e^{-\gamma\tau}&\rho_{11,11}\end{pmatrix}$

where it was possible to define a total decoherence rate
$\gamma_{\text{tot}}=\gamma_{\text{short}}+\Gamma_{\text{long}}\Delta x^{2}$.
This will not be possible in the case of spins with $J\geq 1/2$.