# Satisfying the Einstein-Podolsky-Rosen criterion with massive particles

J. Peise1, I. Kruse1, K. Lange1, B. Lücke1, L. Pezzè2, J. Arlt3, W. Ertmer1,
K. Hammerer4, L. Santos4, A. Smerzi2, C. Klempt1∗ 1Institut für Quantenoptik,
Leibniz Universität Hannover, Welfengarten 1, D-30167 Hannover, Germany
2 QSTAR, INO-CNR and LENS, Largo Enrico Fermi 2, I-50125, Firenze, Italy
3 Institut for Fysik og Astronomi, Aarhus Universitet, Ny Munkegade 120,
DK-8000 Århus C, Denmark
4 Institut für Theoretische Physik, Leibniz Universität Hannover, Appelstraße
2, D-30167 Hannover, Germany
∗ Correspondence should be addressed to C.K<EMAIL_ADDRESS>hannover.de)

In 1935, Einstein, Podolsky and Rosen (EPR) questioned the completeness of
quantum mechanics by devising a quantum state of two massive particles with
maximally correlated space and momentum coordinates. The EPR criterion
qualifies such continuous-variable entangled states, where a measurement of
one subsystem seemingly allows for a prediction of the second subsystem beyond
the Heisenberg uncertainty relation. Up to now, continuous-variable EPR
correlations have only been created with photons, while the demonstration of
such strongly correlated states with massive particles is still outstanding.
Here, we report on the creation of an EPR-correlated two-mode squeezed state
in an ultracold atomic ensemble. The state shows an EPR entanglement parameter
of $0.18(3)$, which is $2.4$ standard deviations below the threshold $1/4$ of
the EPR criterion. We also present a full tomographic reconstruction of the
underlying many-particle quantum state. The state presents a resource for
tests of quantum nonlocality and a wide variety of applications in the field
of continuous-variable quantum information and metrology.

In their original publication 1, Einstein, Podolsky and Rosen describe two
particles A and B with correlated position
$x_{\mathrm{B}}=x_{\mathrm{A}}+x_{0}$ and anti-correlated momentum
$p_{\mathrm{B}}=-p_{\mathrm{A}}$ (see Fig. 1a). When coordinates
$x_{\mathrm{A}}$ and $p_{\mathrm{A}}$ are measured in independent realizations
of the same state, the correlations allow for an exact prediction of
$x_{\mathrm{B}}$ and $p_{\mathrm{B}}$. EPR assumed that such exact predictions
necessitate an ”element of reality” which predetermines the outcome of the
measurement. Quantum mechanics however prohibits the exact knowledge of two
noncommuting variables like $x_{\mathrm{B}}$ and $p_{\mathrm{B}}$, since their
measurement uncertainties are subject to the Heisenberg relation $\Delta
x_{\mathrm{B}}\Delta p_{\mathrm{B}}\geq\frac{1}{2}$. EPR thus concluded that
quantum mechanics is incomplete - under their assumptions which are today
known as ”local realism”. Later, the notion of EPR correlations was
generalized to a more realistic scenario, yielding a criterion 2; 3 for the
uncertainties $\Delta x_{\mathrm{B}}^{\mathrm{inf}}$, $\Delta
p_{\mathrm{B}}^{\mathrm{inf}}$ of the inferred predictions for
$x_{\mathrm{B}}$ and $p_{\mathrm{B}}$. The EPR criterion is met if these
uncertainties violate the Heisenberg inequality for the inferred uncertainties
$\Delta x_{\mathrm{B}}^{\mathrm{inf}}\,\Delta
p_{\mathrm{B}}^{\mathrm{inf}}\geq\frac{1}{2}$. The EPR criterion also
certifies steering, a concept termed by Schrödinger 4; 5 in response to EPR
which has attracted a lot of interest in the past years 6. An experimental
realization of states satisfying the EPR criterion is not only desirable in
the context of the fundamental questions raised by EPR, but also provides a
valuable resource for many quantum information tasks, including dense coding,
quantum teleportation 7 and quantum metrology 8. Some quantum information
tasks specifically require the strong type of entanglement that is tested by
the EPR criterion, as for example one-sided device independent entanglement
verification 9.

Figure 1:  Einstein-Podolsky-Rosen correlations. a EPR’s original work
describes two particles A and B with maximally correlated position and
momentum coordinates $x_{\mathrm{A/B}}$ and $p_{\mathrm{A/B}}$. b Spin
dynamics in a Bose-Einstein condensate can be used to create EPR correlations
between $N_{\mathrm{A/B}}$ atoms in two different Zeeman levels A and B. The
correlations appear in amplitude $x_{\mathrm{A/B}}$ and phase
$p_{\mathrm{A/B}}$ quadratures which are defined as a function of the creation
and annihilation operators $a^{\dagger}_{\mathrm{A/B}}$ and $a_{\mathrm{A/B}}$
in the two modes.  Figure 2:  Quadrature distributions. a Recorded probability
distributions of the quadratures (first row). Distributions of the quadratures
according to the reconstructed state (second row, see below for details on
reconstruction). Ideal distributions of the quadratures of a two-mode squeezed
state with the reconstructed squeezing parameter $\xi_{\mathrm{fit}}=0.63$
(third row). b Two-mode variances $V_{\mathrm{X}(\theta)}^{+}$ and
$V_{\mathrm{X}(\theta)}^{-}$ as a function of the local oscillator phase
$\theta$. c EPR parameter $V_{\mathrm{X}(\theta)}^{+}\times
V_{\mathrm{X}(\theta+\frac{\pi}{2})}^{-}$ as a function of the local
oscillator phase $\theta$. Data points below the dashed line meet the EPR
criterion and therefore violate the Heisenberg inequality for inferred
uncertainties. At the same time, the EPR criterion also certifies
inseparability 12. d The weaker inseparability parameter
$V_{\mathrm{X}(\theta)}^{+}+V_{\mathrm{X}(\theta+\frac{\pi}{2})}^{-}$ as a
function of the local oscillator phase $\theta$. The dotted line indicates
inseparability of the underlying quantum state. The dashed line is a
sufficient condition for the EPR criterion. The error bars indicate the
statistical uncertainty (one standard deviation).

Up to now, the creation of continuous-variable entangled states satisfying the
EPR criterion was only achieved in optical systems. In a seminal publication
10, the EPR criterion was met by a two-mode squeezed vacuum state generated by
optical parametric down-conversion. In this experiment, and in more recent
investigations 11; 12, continuous variables are represented by amplitude
$x_{\mathrm{A/B}}$ and phase $p_{\mathrm{A/B}}$ quadratures, satisfying the
commutation relation $[x_{\mathrm{A/B}},p_{\mathrm{A/B}}]=i$. These
quadratures can be measured accurately by optical homodyning. The correlations
are captured by the four two-mode variances
$V_{x}^{\pm}=\mathrm{Var}(x_{\mathrm{A}}\pm x_{\mathrm{B}})$ and
$V_{p}^{\pm}=\mathrm{Var}(p_{\mathrm{A}}\pm p_{\mathrm{B}})$. These variances
were proven to fulfill a symmetric form of Reid’s inequality 3
$V_{x}^{-}\times V_{p}^{+}<\frac{1}{4}$, which is a sufficient EPR criterion
since $V_{x}^{-}=(\Delta x_{\mathrm{B}}^{\mathrm{inf}})^{2}$ and
$V_{p}^{+}=(\Delta p_{\mathrm{B}}^{\mathrm{inf}})^{2}$. In recent years,
continuous-variable entangled optical states have been applied for proof-of-
principle quantum computation and communication tasks 7. Despite these
advances with optical systems, an experimental realization of EPR correlations
with massive particles is desirable, because of the similarity to the original
EPR proposal and since massive particles may be more tightly bound to the
concept of local realism 2; 3.

Entangled states of massive particles have been generated in neutral atomic
ensembles, promising fruitful applications in precision metrology due to the
large achievable number of entangled atoms 13; 14; 15; 16. They have been
created by atom-light interaction at room-temperature 17; 14, in cold samples
18; 19; 20; 21; 22, and by collisional interactions in Bose-Einstein
condensates 13; 23; 24; 16; 25. For Gaussian states of two collective atomic
modes, the inseparability criterion 26; 27 $V_{x}^{-}+V_{p}^{+}<2$ has been
used to demonstrate entanglement 17; 14; 28, but the strong correlations
necessary to meet the more demanding EPR criterion $V_{x}^{-}\times
V_{p}^{+}<\frac{1}{4}$ have not been achieved so far.

Here we report on the creation of an entangled state from a spinor Bose-
Einstein condensate (BEC) which meets the EPR criterion. We exploit spin-
changing collisions to generate a two-mode squeezed vacuum state in close
analogy to optical parametric down-conversion. The phase and amplitude
quadratures are accessed by atomic homodyning. Their correlations yield an EPR
entanglement parameter of $0.18(3)$, which is $2.4$ standard deviations below
the threshold $1/4$ of the EPR criterion. Finally, we deduce the density
matrix of the underlying many-particle state from a Maximum Likelihood
reconstruction.

##  I Results

Two-mode squeezed vacuum. In our experiments, a BEC with $2\times 10^{4}$ 87Rb
atoms in the Zeeman level $\left(F,m_{F}\right)=\left(1,0\right)$ generates
atom pairs in the levels $\left(1,\pm 1\right)$ due to spin-changing
collisions (see Fig. 1b), ideally yielding the two-mode squeezed state

$\left|\xi\right\rangle=\sum_{n=0}^{\infty}\frac{(-i\tanh{\xi})^{n}}{\cosh{\xi}}\left|n,n\right\rangle,$
(1)

where $\xi=\Omega t$ is the squeezing parameter, which depends on the spin
dynamics rate $\Omega=2\pi\times 5.1$ Hz and the spin dynamics duration $t=26$
ms. The notation $\left|n,m\right\rangle$ represents a two-mode Fock state in
the two Zeeman levels $\left(1,\pm 1\right)$. The generated two-mode squeezed
state can be characterized by the quadratures
$x_{\mathrm{A/B}}=\frac{1}{\sqrt{2}}(a^{\dagger}_{\mathrm{A/B}}+a_{\mathrm{A/B}})$
and
$p_{\mathrm{A/B}}=\frac{i}{\sqrt{2}}(a^{\dagger}_{\mathrm{A/B}}-a_{\mathrm{A/B}})$
for the two levels $\left(1,\pm 1\right)$. These exhibit EPR correlations,
since the variances $V_{x}^{-}=V_{p}^{+}=e^{-2\xi}$ are squeezed, while the
conjugate variances $V_{x}^{+}=V_{p}^{-}=e^{2\xi}$ are anti-squeezed. The
state fulfills Reid’s EPR criterion for $\xi>\frac{1}{2}\ln(2)\approx 0.35$
which corresponds to a spin dynamics duration of more than $11$ ms. In the
limit of large squeezing, our setup presents an exact realization of the
perfect correlations with massive particles envisioned by EPR.

Quadratures and the EPR criterion. The quadratures in the two modes are
simultaneously detected in our experiments by unbalanced homodyne detection
(see methods). Atomic homodyne detection was first demonstrated in Ref. 28,
and later applied to discriminate between vacuum and few-atom states in a
quantum Zeno scenario 29. A small radio-frequency pulse couples $15\%$ of the
BEC in the level $\left(1,0\right)$ (the local oscillator) symmetrically to
the two modes $\left(1,\pm 1\right)$. The local oscillator phase is
represented by the BEC phase relative to the phase sum of the two ensembles in
$\left(1,\pm 1\right)$. It can be varied in our experiments by shifting the
energy of the level $\left(1,-1\right)$ for an adjustable time. From the
measured number of atoms in both levels, we obtain a linear combination of the
quadratures according to
$X_{\mathrm{A/B}}(\theta)=x_{\mathrm{A/B}}\cos(\theta-\frac{\pi}{4})+p_{\mathrm{A/B}}\sin(\theta-\frac{\pi}{4})$.
Figure 2a shows two-dimensional histograms of these measurements for
$\theta=\frac{3}{4}\pi$ and $\theta=\frac{5}{4}\pi$, corresponding to the
$x$\- and $p$-quadratures. The histograms demonstrate the strong correlation
and anticorrelation of these two quadratures, as expected for the EPR case.
The variances along the two diagonals, represented by
$V_{\mathrm{X}(\theta)}^{\pm}=\mathrm{Var}(X_{\mathrm{A}}(\theta)\pm
X_{\mathrm{B}}(\theta))$, are shown in Fig. 2b and reveal the expected two-
mode squeezing behavior. From these measurements, we quantify the EPR
entanglement by Reid’s criterion, yielding $V_{x}^{+}\times
V_{p}^{-}=0.18(3)$, which is $2.4$ standard deviations below the limit of
$\frac{1}{4}$. In addition, the data also fulfills the inseparability
criterion as $V_{x}^{+}+V_{p}^{-}=0.85(8)$, which is $15$ standard deviations
below the classical limit of $2$ (see Fig. 2d), and meets the criterion for a
symmetric (”two-way”) steering between the systems 6. We estimate that the
product value could be reduced to $V_{x}^{+}\times V_{p}^{-}=0.13(3)$ if the
radio frequency intensity noise was eliminated by stabilization or
postcorrection. The experimental creation of entangled massive particles which
satisfy the continuous-variable EPR criterion presents the main result of this
publication.

Figure 3:  Squeezing dynamics. a/b Time evolution of the four two-mode
quadrature variances $V_{x}^{\pm}$, $V_{p}^{\pm}$ during spin dynamics. The
data follows the ideal squeezing/anti-squeezing behavior according to the
independently measured spin dynamics rate (dashed line). c Time evolution of
the EPR parameter $V_{x}^{-}\times V_{p}^{+}$ during spin dynamics. The
squeezed variances and the EPR parameter are well reproduced by a simple noise
model (solid line). The error bars indicate the statistical uncertainty (one
standard deviation). Figure 4:  Result of the state reconstruction. a Elements
$\left\langle
N_{\mathrm{A}},N_{\mathrm{B}}\right|\rho\left|N_{\mathrm{A}},N_{\mathrm{B}}\right\rangle$
of the reconstructed density matrix. The vacuum state $\left|0,0\right\rangle$
and the two-particle twin Fock state $\left|1,1\right\rangle$ show the largest
values. The large contributions of the twin Fock states reveal the creation of
atomic pairs during spin dynamics. For example, the state
$\left|2,2\right\rangle$ has by far the strongest weight among all states with
a total of four particles (highlighted in dark blue). b Reconstructed
distribution of the difference of the number of atoms in the two modes A and B
(solid bars). The distribution is strongly peaked in comparison to a
Poissonian distribution with the same mean number of particles (open bars). c
Reconstructed distribution of the total number of atoms in the two modes A and
B.

Squeezing dynamics. Figure 3 shows the squeezing dynamics due to the spin-
changing collisions. For these measurements, we fix the local oscillator phase
to the values $\theta\approx 3\pi/4$ and $\theta\approx 5\pi/4$ to record only
the $x$\- and $p$-variances. As a function of the evolution time, the
variances $V_{x}^{-},V_{p}^{+}$ are squeezed below the vacuum reference of
$1$, while the variances $V_{x}^{+},V_{p}^{-}$ exhibit an antisqueezing
behavior (Fig. 3a and b). From these data, we extract the EPR parameter
$V_{x}^{-}\times V_{p}^{+}$, as a function of evolution time (see Fig. 3c).
The EPR parameter is quickly pushed below $1$ and follows the prediction for
an ideal squeezed state. It eventually reaches a minimum at the optimal
squeezing time of $26$ ms, as used for the data in Fig. 2. The data is well
reproduced by a simple noise model, which includes a radio-frequency intensity
noise of $0.4\%$ and a local oscillator phase noise of $0.044\pi$ (see
methods).

Full state reconstruction. The total of $2864$ homodyne measurements obtained
for different local oscillator phases at the optimal evolution time allow for
a full reconstruction of the underlying many-particle state. Previously,
tomography of an atomic state was demonstrated either by reconstruction of the
Wigner function 30 or the Husimi $Q$-distribution 25; 21. However, both
methods yield a characterization of the state’s projection on the fully
symmetric subspace only. The well-developed methods in quantum optics 31
allowed for a full reconstruction of an optical two-mode squeezed state by
homodyne tomography 32; 11. Despite the beautiful tomography data, the optical
state reconstruction assumed either Gaussian states or averaged over all phase
relations, such that the coherence properties could not be resolved.

In contrast, we obtain an unbiased, positive semidefinite density matrix by
Maximum Likelihood reconstruction 31; 33 of the experimental data, free of any
a-priori hypothesis. This represents the second major result of this
publication. The recorded data for each local oscillator phase are binned in
two-dimensional histograms (see Fig. 2a) presenting the marginal distributions
for the $x_{\mathrm{A/B}}$ and $p_{\mathrm{A/B}}$ variables. The reconstructed
state is the one which optimally reproduces the measured histograms by a
superposition of harmonic oscillator wave functions 31. The coefficients of
this superposition are estimates of the density matrix elements of the
underlying quantum state (see methods).

Figure 4 shows the result of the reconstruction. The diagonal matrix elements
(Fig. 4a) witness the predominant creation of atom pairs. The two-particle
twin Fock state $\left|1,1\right\rangle$ shows the strongest contribution
besides the vacuum state. Likewise, the twin Fock states
$\left|2,2\right\rangle$ to $\left|4,4\right\rangle$ have the strongest
contribution for a given total number of particles. The strong nonclassicality
of the reconstructed state becomes also apparent in the distributions of the
difference and the sum of the particles (Fig. 4b and c). The distribution of
the number difference is strongly peaked at zero and is much narrower than a
Poissonian distribution with the same number of particles. The distribution of
the total number of atoms shows an indication of the characteristic even/odd
oscillations, which is caused by the pair production in the underlying spin
dynamics.

##  II Discussion

For an evaluation of the created state, we have extracted a logarithmic
negativity of $1.43\pm 0.06$ from the reconstructed density matrix. This value
is above the threshold of zero for separable states and signals
nonseparability of the reconstructed state. The quantum Fisher information 34
$F_{Q}$ for the state projected on fixed-N subspaces reveals that
$F_{Q}/\bar{n}=1.5\pm 0.1$, where $\bar{n}$ is the average number of
particles. Since $F_{Q}/\bar{n}>1$ the state is a resource for quantum
enhanced metrology 34. Furthermore, we fit an ideal two-mode squeezed state
with variable squeezing parameter $\xi$ to the reconstructed two-mode density
matrix with maximal fidelity. With a fidelity of $78.4$%, the experimentally
created state matches a two-mode squeezed state with a squeezing parameter of
$\xi_{\mathrm{fit}}=0.63$. The fidelity increases to $90$% if we include local
oscillator phase noise and statistical noise. The unwanted contributions in
the density matrix, including the off-diagonal terms in Fig. 4a, can be well
explained by four effects. Firstly, the purity of the reconstructed state is
limited by the finite number of homodyne measurements. Secondly, small drifts
in the microwave intensity of the dressing field (on the order of $0.1$ %),
which shifts the level $(1,-1)$, result in a small drift of the local
oscillator phase. Thirdly, a small drift of the radio-frequency coupling
strength during homodyning virtually increases the variance in the
$(x_{\mathrm{A}}+x_{\mathrm{B}})$ and the $(p_{\mathrm{A}}+p_{\mathrm{B}})$
directions. Finally, we did not correct for the detection noise of our
absorption imaging.

Our experimental realization of the EPR criterion demonstrates a strong form
of entanglement intrinsically connected to the notion of local realism. In the
future, the presented atomic two-mode squeezed state allows to demonstrate the
continuous-variable EPR paradox with massive particles. Since the two modes A
and B are Zeeman levels with an opposite magnetic moment, the modes can be
easily separated with an inhomogenous magnetic field to ensure a spatial
separation. The nonlocal EPR measurement could then be realized by homodyning
with two spatially separated local oscillators. These can be provided by
splitting the remaining BEC into the levels $\left(2,\pm 1\right)$ which have
the same magnetic moment as the two EPR modes. Furthermore, this setup can be
complemented by a precise atom number detection to demonstrate a violation of
a Clauser-Horne-Shimony-Holt-type inequality. Such a measurement presents a
test of local realism with continuous-variable entangled states. In this
context, neutral atoms provide the exciting possibility to investigate the
influence of increasingly large particle numbers and possible effects of
gravity.

##  III Methods

Experimental sequence. We start the experiments with an almost pure Bose-
Einstein condensate of $20,000$ 87Rb atoms in an optical dipole potential with
trap frequencies of $2\pi\times(280,220,180)$ Hz. At a homogeneous magnetic
field of $2.6$ G with fluctuations of about $70$ $\mu$G, the condensate is
transferred from the level $\left(F,m_{F}\right)=\left(2,2\right)$ to the
level $\left(1,0\right)$ by a series of three resonant microwave pulses.
During this preparation, two laser pulses resonant to the $F=2$ manifold rid
the system of atoms in unwanted spin states. Directly before spin dynamics is
initiated, the output states $\left(1,\pm 1\right)$ are emptied with a pair of
microwave $\pi$-pulses from $\left(1,+1\right)$ to $\left(2,+2\right)$ and
from $\left(1,-1\right)$ to $\left(2,-2\right)$ followed by another light
pulse. This cleaning procedure ensures that no thermal or other residual
excitations are present in the two output modes, which may destroy the EPR
signal 35.

Figure 5: Schematic overview of the experimental sequence. The state is
produced via a spin dynamics process at time $t_{1}$. The remaining condensate
in the $m_{F}=0$ state acts as the local oscillator and has a phase imprinted
upon it by an energy shift of the $m_{F}=-1$ state relative to the $m_{F}=1$
state. This phase accumulates during time $t_{2}$ with $500\,\mathrm{Hz}$. At
time $t_{3}$ a radio-frequency pulse couples the $m_{F}=0$ condensate to the
$m_{F}=\pm 1$ states. This coupling acts as a three port beam splitter. An
inhomogeneous magnetic field is applied at $t_{4}$ before the absorption
detection is performed at $t_{5}$. A microwave dressing shifts the energy of
the $m_{F}=-1$ state during the spin dynamics time $t_{1}$ and the radio-
frequency coupling $t_{3}$.

Figure 5 shows a schematic overview of the following experimental sequence. A
microwave frequency which is red-detuned to the transition between the levels
$\left(1,-1\right)$ and $\left(2,-2\right)$ by about $208$ kHz is
adiabatically ramped on within $675\,\mu$s. The microwave shifts the level
$\left(1,-1\right)$ by about $500$ Hz, depending on the chosen detuning, to
compensate for the quadratic Zeeman effect of $q=491$ Hz, such that multiple
spin dynamics resonances can be addressed 36; 16. Each resonance condition
belongs to a specific spatial mode of the states $\left(1,\pm 1\right)$ to
which the atoms are transferred. If the energy of the level
$\left(1,-1\right)$ is reduced, more internal energy is released, and higher
excited spatial modes are populated (for details, see Ref. 36). Here, we
choose the first resonance, where spin dynamics leads to a population of the
levels $\left(1,\pm 1\right)$ in the ground state of the effective potential.
This ensures an optimal spatial overlap between the atoms in the three
contributing levels. This resonance condition is reached, when the input state
(two atoms in the BEC in the level $\left(1,0\right)$ at the energy of the
chemical potential) is exactly degenerate with the output state (two atoms in
the levels $\left(1,\pm 1\right)$ including dressing, trap energy and mean-
field shift). Due to this degeneracy, the phase relation between the initial
condensate and the output state stays fixed during the spin dynamics evolution
time. For this configuration, we have checked that spin dynamics is the only
relevant process which produces atoms in the state $\left(1,\pm 1\right)$ (see
Ref. 29, Fig. 3). Subsequently, the microwave dressing field is ramped down
within $675\,\mu$s, stays off for a variable duration between $25$ and
$1150\,\mu$s and is quickly switched on again. The variable hold time allows
for an adjustment of the local oscillator phase relative to the output state.

For the atomic homodyning, a radio-frequency (rf) pulse with a frequency of
$1.834$ MHz and a duration of $\tau=30\,\mu$s couples the level
$\left(1,0\right)$ with the levels $\left(1,\pm 1\right)$. The microwave
dressing field is chosen such that both rf transitions are resonant but the
resonance condition for spin dynamics is not fulfilled. Afterwards, the dipole
trap is switched off to allow for a ballistic expansion. After an initial
expansion of $1.5$ ms to reduce the density, a strong magnetic field gradient
is applied to spatially separate the atoms in the three Zeeman levels.
Finally, the number of atoms in the three clouds is detected by absorption
imaging on a CCD camera with a large quantum efficiency. The statistical
uncertainty of a number measurement is dominated by the shot noise of the
photoelectrons on the camera pixels and amounts to $16$ atoms. We estimate the
uncertainty of the total number calibration to be less than $1$%.

Three-mode unbalanced homodyning. The rf coupling is described by the three-
mode unitary operation $e^{-i{H}\tau/\hbar}$, where

${H}=\frac{\hbar\Omega_{+1}}{2\sqrt{2}}\big{(}{a}_{\rm
A}^{\dagger}{a}_{0}+{a}_{\rm
A}{a}_{0}^{\dagger}\big{)}+\frac{\hbar\Omega_{-1}}{2\sqrt{2}}\big{(}{a}_{\rm
B}^{\dagger}{a}_{0}+{a}_{\rm B}{a}_{0}^{\dagger}\big{)},$

and $\Omega_{\pm 1}$ are Rabi frequencies for the $(1,0)\leftrightarrow(1,\pm
1)$ transition (in general $\Omega_{+1}\neq\Omega_{-1}$). To calculate the
mode transformation, we use $[{H},{a}_{\rm
A}]=-\hbar\Omega_{+1}{a}_{0}/2\sqrt{2}$, $[{H},{a}_{\rm
B}]=-\hbar\Omega_{-1}{a}_{0}/2\sqrt{2}$ and
$[{H},{a}_{0}]=-\hbar(\Omega_{+1}{a}_{\rm A}+\Omega_{-1}{a}_{\rm
B})/2\sqrt{2}$. We have

$\begin{pmatrix}{a}_{\rm A}\\\ {a}_{\rm B}\\\ {a}_{0}\end{pmatrix}_{\rm
out}=\begin{pmatrix}\tfrac{\tilde{\Omega}_{+1}^{2}c+\tilde{\Omega}_{-1}^{2}}{2}&\tfrac{\tilde{\Omega}_{+1}\tilde{\Omega}_{-1}(c-1)}{2}&\tfrac{\tilde{\Omega}_{+1}s}{i\sqrt{2}}\\\
\tfrac{\tilde{\Omega}_{+1}\tilde{\Omega}_{-1}(c-1)}{2}&\tfrac{\tilde{\Omega}_{-1}^{2}c+\tilde{\Omega}_{+1}^{2}}{2}&\tfrac{\tilde{\Omega}_{-1}s}{i\sqrt{2}}\\\
\tfrac{\tilde{\Omega}_{+1}s}{i\sqrt{2}}&\tfrac{\tilde{\Omega}_{-1}s}{i\sqrt{2}}&c\end{pmatrix}\begin{pmatrix}{a}_{\rm
A}\\\ {a}_{\rm B}\\\ {a}_{0}\end{pmatrix},$ (2)

where $c=\cos(\Omega\tau/2)$, $s=\sin(\Omega\tau/2)$, and

$\tilde{\Omega}_{+1}=\frac{\Omega_{+1}}{\Omega},\qquad\tilde{\Omega}_{-1}=\frac{\Omega_{-1}}{\Omega},\qquad{\rm
with}\quad\Omega=\sqrt{\frac{\Omega_{+1}^{2}+\Omega_{-1}^{2}}{2}},$ (3)

are rescaled Rabi frequencies. Below, we illustrate how the measurement of the
number of particles in the $m_{F}=\pm 1$ mode after the rf coupling, $N_{\rm
A}=({a}_{\rm A}^{\dagger}{a}_{\rm A})_{\rm out}$ and $N_{\rm B}=({a}_{\rm
B}^{\dagger}{a}_{\rm B})_{\rm out}$, gives access to the number conserving
quadratures

$x_{\rm A,B}=\frac{\hat{a}_{\rm
A,B}^{\dagger}{a}_{0}+{a}_{0}^{{\dagger}}\hat{a}_{\rm
A,B}}{\sqrt{2n_{0}}},\qquad p_{\rm A,B}=\frac{{a}_{0}^{{\dagger}}\hat{a}_{\rm
A,B}-\hat{a}_{\rm A,B}^{{\dagger}}{a}_{0}}{i\sqrt{2n_{{0}}}},$ (4)

where $n_{0}=\langle{a}_{0}^{{\dagger}}{a}_{0}\rangle$ is the average number
of particles in the condensate before homodyne (similarly, $n_{\rm
A,B}=\langle\hat{a}_{\rm A,B}^{{\dagger}}\hat{a}_{\rm A,B}\rangle$). In our
experiment, $(n_{\rm A}+n_{\rm B})/n_{0}\sim 10^{-4}$. We thus neglect
fluctuations of the number of particles in the $m_{F}=0$ mode, replacing
${a}_{0}^{\dagger}{a}_{0}$ with its mean value $n_{0}\approx n_{\rm
A}+n_{0}+n_{\rm B}=N_{\rm tot}$.

Number difference. The quadrature difference can be experimentally obtained by
measuring the difference of the number of particles in the $\pm 1$ modes. From
Eq. (2) we can directly calculate $N_{\rm A}-N_{\rm B}$. To the leading order
in $n_{0}$, we have

$\displaystyle\frac{N_{\rm A}-N_{\rm B}}{\sqrt{s^{2}N_{\rm tot}}}$
$\displaystyle\approx$ $\displaystyle\frac{s\sqrt{N_{\rm
tot}}(\tilde{\Omega}_{+1}^{2}-\tilde{\Omega}_{-1}^{2})}{2}$
$\displaystyle+\frac{\tilde{\Omega}_{+1}\big{[}c(\tilde{\Omega}_{+1}^{2}-\tilde{\Omega}_{-1}^{2})+2\tilde{\Omega}_{-1}^{2}\big{]}}{2}p_{\mathrm{A}}$
$\displaystyle+\frac{\tilde{\Omega}_{-1}\big{[}c(\tilde{\Omega}_{+1}^{2}-\tilde{\Omega}_{-1}^{2})-2\tilde{\Omega}_{-1}^{2}\big{]}}{2}p_{\mathrm{B}}.$

Since $\Omega_{+1}$ and $\Omega_{-1}$ only differ by $1.7\,\%$ in our
experiments, and $c(\tilde{\Omega}_{+1}^{2}-\tilde{\Omega}_{-1}^{2})\ll
2\tilde{\Omega}_{\pm 1}^{2}$, we can simplify this equation and obtain

$p_{\rm A}-p_{\rm B}=\frac{N_{\rm A}-N_{\rm
B}-s^{2}\big{(}\tilde{\Omega}_{+1}^{2}-\tilde{\Omega}_{-1}^{2}\big{)}N_{\textrm{tot}}/2}{\sqrt{s^{2}N_{\textrm{tot}}}}.$
(6)

Number sum. The quadrature sum is obtained by adding the number of particles
in the $\pm 1$ modes after homodyning:

$\displaystyle\frac{N_{\rm A}+N_{\rm B}}{\sqrt{s^{2}c^{2}N_{\rm tot}}}$
$\displaystyle\approx$ $\displaystyle\frac{s\sqrt{N_{\rm
tot}}}{c}+\tilde{\Omega}_{+1}p_{\rm A}+\tilde{\Omega}_{-1}p_{\rm B}$ (7)

Taking $\tilde{\Omega}_{-1}\approx\tilde{\Omega}_{+1}\approx 1$, we have

$p_{\rm A}+p_{\rm B}=\frac{N_{\rm A}+N_{\rm
B}-s^{2}N_{\textrm{tot}}}{\sqrt{s^{2}c^{2}N_{\textrm{tot}}}}.$ (8)

Finally, the mean transfer of particles from $m_{F}=0$ to $m_{F}=\pm 1$ and
the mean number difference is used to calculate

$\left\langle\frac{N_{\rm A}+N_{\rm B}}{N_{\mathrm{tot}}}\right\rangle\approx
s^{2},\quad{\rm and}\quad\left\langle\frac{N_{\rm A}-N_{\rm
B}}{N_{\mathrm{tot}}}\right\rangle\approx\frac{1}{2}s^{2}\big{(}\tilde{\Omega}_{+1}^{2}-\tilde{\Omega}_{-1}^{2}\big{)}.$
(9)

Observing a transfer of 15% of the atoms from $m_{F}=0$ to $m_{F}=\pm 1$ we
deduce $c^{2}\approx 0.85$.

To summarize, Eqs. (6) and (8) are used to experimentally obtain $p_{\rm A}\pm
p_{\rm B}$ from the measurement of the number of particles in the output
modes. The quadratures $x_{\rm A}\pm x_{\rm B}$ are obtained with the same
method, by applying a relative $\pi/2$ phase between the pump and side modes
before homodyne detection.

Entanglement criteria for continuous variables. Criteria for identifying
continuous-variable entanglement between the systems A and B, with no
assumption on the quantum state of the local oscillator, have been discussed
in Ref. 37.

Separability. For mode-separable states, ${\rho}_{\rm
sep}=\sum_{k}p_{k}{\rho}^{(k)}_{\rm A}\otimes{\rho}^{(k)}_{\rm B}$ ($p_{k}>0$,
$\sum_{k}p_{k}=1$), we have 37; 38

$V_{x}^{\pm}+V_{p}^{\mp}\geq 2-\frac{{n}_{\rm A}+n_{\rm B}}{{n}_{0}},$ (10)

where $V_{x}^{\pm}={\rm Var}(x_{\mathrm{A}}\pm x_{\mathrm{B}})$ and
$V_{p}^{\pm}={\rm Var}(p_{\mathrm{A}}\pm p_{\mathrm{B}})$ are the variances of
quadrature sum and difference, respectively. A violation of Eq. (10) signals
non-separability, i.e. ${\rho}\neq{\rho}_{\rm sep}$. Equation (10) generalizes
the criterion of Refs. 26 and 27 that was derived for standard quadrature
operators (i.e. when the $m_{F}=0$ mode is treated parametrically, the
operator ${a}_{0}$ being replaced by $\sqrt{n_{0}}$ ).

EPR criterion. Reid’s EPR criterion corresponds to a violation of the
Heisenberg uncertainty relation on system B, when measurements are performed
on system A. This requires the two-mode state to be non-separable and to have
strong correlations between the sum and difference of position and momentum
quadratures, $x_{\rm A}\pm x_{\rm B}$ and $p_{\rm A}\mp p_{\rm B}$. We point
out that not all non-separable states fulfill Reid’s criterion. The position-
momentum quadratures for the B mode satisfy the commutation relation $[x_{\rm
B},p_{\rm B}]=-i({a}_{\rm B}^{{\dagger}}{a}_{\rm
B}-{a}_{0}^{{\dagger}}{a}_{0})/{n}_{0}$. The corresponding Heisenberg
uncertainty relation is $\Delta^{2}x_{\rm B}\Delta^{2}p_{\rm
B}\geq\tfrac{1}{4}(1-\tfrac{{n}_{\rm B}}{{n}_{0}})^{2}$. Let us introduce the
quantities $x_{\rm ext}(x_{\rm A})$ and $p_{\rm ext}(p_{\rm A})$, which are
the estimate of $x_{\rm B}$ and $p_{\rm B}$ on system B, respectively, given
the results of quadrature measurements on the system A. We then indicate as
$\Delta^{2}x_{\rm B}^{\rm inf}$ the squared deviation of the estimate from the
actual value, averaged over all possible results $x_{\rm A}$,

$\Delta^{2}x_{\rm B}^{\rm inf}=\int dx_{\rm B}\int dx_{\rm A}P(x_{\rm
A},x_{\rm B})[x_{\rm B}-x_{\rm ext}(x_{\rm A})]^{2},$ (11)

and similarly for $\Delta^{2}p_{\rm B}^{\rm inf}$, where $P(x_{\rm A},x_{\rm
B})$ is the joint probability. Reid’s criterion thus reads 37
$\Delta^{2}x_{\rm B}^{\rm inf}\Delta^{2}p_{\rm B}^{\rm
inf}<\tfrac{1}{4}(1-\tfrac{{n}_{\rm B}}{{n}_{0}})^{2}$. Taking $x_{\rm
ext}(x_{\rm A})=x_{\rm A}-(\bar{x}_{\rm A}-\bar{x}_{\rm B})$ and $p_{\rm
ext}(p_{\rm A})=-p_{\rm A}+(\bar{p}_{\rm B}+\bar{p}_{\rm A})$, where the bar
indicates statistical average, Reid’s criterion translates into a condition
for the product of quadrature variances:

$V_{x}^{-}\times V_{p}^{+}<\frac{1}{4}\bigg{(}1-\frac{{n}_{\rm
B}}{{n}_{0}}\bigg{)}^{2}.$ (12)

In our case, ${n}_{\rm A}/{n}_{0}\sim{n}_{\rm B}/{n}_{0}\approx 10^{-4}$.
Therefore corrections in Eqs. (10) and (12) due to finite number of particles
in the $m_{F}=0$ are negligible. We are thus in a continuous-variable limit.

We point out that the above EPR criterion – consistent with the analysis of
the experimental data presented in Figs. 2b and 3 – uses quadrature variances
with symmetric contributions from A and B. In this case the EPR threshold is
1/4. The above inequalities and entanglement criteria can be generalized (and
optimized) for asymmetric contributions, see Refs. 3 and 12.

Quantum state tomography. Here we discuss the protocol used for quantum state
tomography and, very briefly, its theoretical basis. A more detailed
discussion can be found in Refs. 31 and 33. We point out that our state
reconstruction is performed without assumption neither on the state nor on the
experimental quadrature distribution, in particular we do not assume our state
to be a Gaussian state.

We have collected a total of $N=2864$ measurements of the quadratures $x_{\rm
A}$ and $x_{\rm B}$ at different values of the local oscillator phase $\theta$
relative to the side modes. The measurement results are binned in 2D
histograms (see Fig. 2a, where the typical bin width is $dx=0.25$) such that
we take $x_{\rm A,B}$ to have a discrete spectrum. To simplify the notation,
let us indicate as $x$ the square bin $[x_{\rm A},x_{\rm A}+dx]$, $[x_{\rm
B},x_{\rm B}+dx]$. Given a quantum state ${\rho}$ (unknown here), the
probability to observe a certain sequence of results ($n_{x,\theta}$
measurement in the bin $x$, when the phase value is set to $\theta$, with
$\sum_{x,\theta}n_{x,\theta}=N$) is

$\mathcal{L}(\\{n_{x,\theta}\\}|{\rho})=\frac{N!}{\prod_{x,\theta}n_{x,\theta}!}\prod_{x,\theta}P_{{\rho}}(x,\theta)^{n_{x,\theta}},$
(13)

indicated as likelihood function. In Eq. (13),
$P_{{\rho}}(x,\theta)=P_{{\rho}}(x|\theta)P(\theta)$ is the joint probability,
$P_{{\rho}}(x|\theta)=\langle x|U_{\theta}^{\dagger}{\rho}U_{\theta}|x\rangle$
is the conditional probability, with $|x\rangle=|x_{\rm A},x_{\rm B}\rangle$,
$U_{\theta}=e^{-i\theta({n}_{A}+{n}_{B})}$, and $P(\theta)$ is the fraction of
measurements done when the phase is equal to $\theta$. The maximum likelihood
(ML)

${\rho}_{\rm
ML}=\arg\big{[}\max_{{\rho}}\mathcal{L}(\\{n_{x,\theta}\\}|{\rho})\big{]}$
(14)

is the state that maximizes $\mathcal{L}(\\{n_{x,\theta}\\}|{\rho})$ on the
manifold of density matrices. To find the ML we use the chain of inequalities
31; 33

$\displaystyle\mathcal{L}(\\{n_{x,\theta}\\}|{\rho})^{\frac{1}{N}}$
$\displaystyle\leq$ $\displaystyle{\rm
Tr}\big{[}{\rho}{R}(\bm{f},\bm{a})\big{]}\prod_{x,\theta}(a_{x,\theta})^{f_{x,\theta}}\leq\lambda(\bm{f},\bm{a})\prod_{x,\theta}(a_{x,\theta})^{f_{x,\theta}},$
(15)

where $a_{x,\theta}$ are arbitrary positive numbers
($\bm{a}=\\{a_{x,\theta}\\}$ is the corresponding vector),
$f_{x,\theta}=n_{x,\theta}/N$ are relative frequencies
($\bm{f}=\\{f_{x,\theta}\\}$ is the corresponding vector), and

${R}=\frac{1}{N}\sum_{x,\theta}\frac{n_{x,\theta}}{P(x|\theta)}U_{\theta}|x\rangle\langle
x|U_{\theta}^{\dagger},$ (16)

is a non-negative operator with largest eigenvalue $\lambda(\bm{f},\bm{a})$.
The second inequality is saturated by taking ${\rho}={\rho}_{\rm ML}$ with
support on the subspace corresponding to the maximum eigenvalue of ${R}$:
${R}(\bm{f},\bm{a}){\rho}_{\rm ML}=\lambda(\bm{f},\bm{a}){\rho}_{\rm ML}$. The
first inequality is a Jensen’s inequality between the geometric and the
arithmetic average (which follows from the concavity of the logarithm). It is
saturated if and only if $a_{x,\theta}=P(x,\theta)$, which also implies ${\rm
Tr}[{R}(\bm{f},\bm{a}){\rho}_{\rm ML}]=1$ and thus $\lambda(\bm{f},\bm{a})=1$.
In conclusion, the search for the ML can be recast in the operator equation
${R}{\rho}_{\rm ML}={\rho}_{\rm ML}$ or, equivalently (since $R$ and
$\rho_{\rm ML}$ are Hermitian operators),

${R}{\rho}_{\rm ML}{R}={\rho}_{\rm ML}.$ (17)

Numerically, this equation is solved iteratively: we start the protocol from a
unit matrix ${\rho}^{(0)}=\mathcal{N}[\mathbbm{1}]$ and apply repetitive
iterations according to Eq. (17),
${\rho}^{(k+1)}=\mathcal{N}[{R}({\rho}^{(k)}){\rho}^{(k)}{R}({\rho}^{(k)})]$
being the $k$th step of the algorithm, where $\mathcal{N}$ denotes
normalization to unit trace. The convergence (which is not guaranteed in
general) is checked. The method guarantees that ${\rho}_{\rm ML}$ is a non-
negative definite operator. In practical implementations, it is most
convenient to work in the atom-number representation and write
${\rho}=\sum_{n_{\rm A},n_{\rm B},m_{\rm A},m_{\rm B}=0}^{n_{\rm
cut}}\rho_{n_{\rm A},n_{\rm B},m_{\rm A},m_{\rm B}}|n_{\rm A},n_{\rm
B}\rangle\langle m_{\rm A},m_{\rm B}|$, where $n_{\rm cut}$ is a cutoff (in
our case $n_{\rm cut}\approx 10$). We use $\langle n|x_{\rm
A,B}\rangle=\tfrac{e^{-x_{\rm A,B}/2}}{\pi^{1/4}\sqrt{2^{n}n!}}H_{n}(x_{\rm
A,B})$, where $H_{n}$ is the Hermite polynomial of order $n$.

Simulation of ideal state reconstruction. To check the consistency of the used
tomography method, we have simulated the reconstruction of an ideal two-mode
squeezed vacuum state $|\xi\rangle$, Eq. (1). The simulation follows three
steps: i) we generate distributions for the quadratures $x_{\rm A,B}$ at
different values of $\theta$, according to the probability $P(x_{\rm
A,B}|\theta)=|\langle\xi|{U}_{\theta}|x_{\rm A,B}\rangle|^{2}$; ii) we
generate $p$ random quadrature data for each $\theta$ value (for a total of
$N=p\times n_{\theta}$, where $n_{\theta}$ is the number of $\theta$ values
considered). This simulates, via Monte Carlo sampling, the acquisition of
experimental data. iii) We perform a ML reconstruction, using the same
numerical code and method used for the analysis of the experimental data. In
Fig. 6 we plot the quantum fidelity between the reconstructed state,
${\rho}_{\rm ML}$, and the two-mode squeezed vacuum state,
$\mathcal{F}=\sqrt{\langle\xi|{\rho}_{\rm ML}|\xi\rangle}$. When the number of
measurements $p$ per $\theta$ value is increased, the fidelity converges to an
asymptotic value $\mathcal{F}_{\infty}\lesssim 1$. The asymptotic fidelity
$\mathcal{F}_{\infty}$ tends to $1$ when decreasing the bin size $dx$.

Figure 6: Simulated quantum state tomography of an ideal squeezed state. Here
we show the quantum fidelity between the reconstructed state ${\rho}_{\rm ML}$
and the ideal two-mode squeezed vacuum state, as a function of the number of
measurements $p$ per $\theta$ value (here we consider $n_{\theta}=29$ phase
values). The dark blue circles are obtained for a binning of $dx=0.25$, the
light blue diamonds for $dx=0.1$. The dashed lines are a guide to the eye, the
error bars (one standard deviation) are obtained with a bootstrap method.

Furthermore, to characterize the entanglement of the reconstructed state, we
have evaluated the logarithmic negativity and the quantum Fisher information
(QFI). The logarithmic negativity is defined as
$E[{\rho}]=\log_{2}\sqrt{{\rho}_{\rm ppt}{\rho}_{\rm ppt}^{\top}}$, where
${\rho}_{\rm ppt}$ is the partial transpose of ${\rho}_{\rm ML}$. Mode-
entanglement is obtained for39 $E[{\rho}]>0$. The QFI for the state projected
over subspaces of a fixed number of particles $n$,
$\rho=\sum_{n}Q_{n}\rho_{n}$, is given by 40

$F_{Q}[{\rho}]=2\sum_{n}Q_{n}\sum_{k_{n},k_{n}^{\prime}}\frac{(p_{k_{n}}-p_{k^{\prime}_{n}})^{2}}{p_{k_{n}}+p_{k^{\prime}_{n}}}|\langle
k_{n}|{J}_{\bm{r}}|k_{n}^{\prime}\rangle|^{2}$ (18)

where ${\rho}_{n}=\sum_{k_{n}}p_{k_{n}}|k_{n}\rangle\langle k_{n}|$ is in
diagonal form and ${J}_{\bm{r}}$ is the collective pseudo-spin operator
(pointing along an arbitrary direction $\bm{r}$ in the three-dimensional
space). The QFI is then maximized over $\bm{r}$, for further details see Refs.
8. Particle entanglement, useful for sub-shot-noise metrology, is obtained
for40 $F_{Q}[{\rho}]\geq\bar{n}$, where $\bar{n}$ corresponds to the average
number of particles in the two-mode state. Similarly to the results of
simulations shown in Fig. 6 we obtain that, in the limit $p\to\infty$ and
$dx\to 0$, the logarithmic negativity and QFI converge to
$E[|\xi\rangle]=2\xi/\log 2$ and
$F_{Q}[|\xi\rangle]=4\sinh^{2}\xi\cosh^{2}\xi$, respectively, which are
analytical values calculated for the two-mode squeezed vacuum state.

Figure 7: Simulated quantum state tomography of a noisy squeezed state.
Tomography reconstruction of a two-mode squeezed state affected by phase noise
and a systematic shift in the variance of the quadrature sum, as in the
experiment. Here we have used the experimental parameters: $p\approx 100$
measurements, $n_{\theta}=29$ phase values and a bin size of $dx=0.25$. The
reconstructed state agrees very well with Fig. 4 of our manuscript, both
qualitatively (showing both the presence of off-diagonal terms and asymmetry)
and quantitatively (with a quantum a fidelity of about 90%).

Noise model and simulation of noisy state reconstruction. The main sources of
noise in our apparatus are phase fluctuations and noise of the rf coupling
strength. The phase noise is assumed to have a Gaussian distribution
$P_{\sigma}(\theta)=e^{-\theta^{2}/2\sigma^{2}}/\sqrt{2\pi\sigma^{2}}$ and we
estimate a width $\sigma\approx 0.36$. Correlated fluctuations of
$\Omega_{+1}$ and $\Omega_{-1}$ affect (to first order) only measurements of
the quadrature sum. We have evaluated that this effect systematically
increases the variance by $0.12$. Both these effects are included in the solid
line of Fig. 3.

We have simulated the state reconstruction in presence of these noise effects.
We model the state in presence of phase noise as

${\rho}_{\rm pn}(\sigma)=\int_{-\pi}^{\pi}d\theta
P_{\sigma}(\theta)|\xi,\theta\rangle\langle\xi,\theta|$ (19)

where
$|\xi,\theta\rangle=\sum_{n=0}^{+\infty}e^{-in\theta}\frac{(\tanh\xi)^{n}}{\cosh\xi}|n,n\rangle$.
The systematic shift of the quadrature sum is included in the calculation of
the quadrature distributions used to generate random data. Results are shown
in Fig. 7. We see that statistical noise (i.e. the limited sample size) and
phase noise are responsible for the appearance of off-diagonal terms, very
similar to the ones observed in Fig. 4. Note that phase noise alone is not
responsible for the appearance of off-diagonal terms in the density matrix.
This can be seen by rewriting Eq. (19) as ${\rho}_{\rm
pn}(\sigma)=\sum_{n,m=0}^{+\infty}\tilde{P}_{\sigma}(n-m)\tfrac{(\tanh\xi)^{n+m}}{\cosh^{2}\xi}|n,n\rangle\langle
m,m|$, where $\tilde{P}_{\sigma}(n-m)=\int_{-\pi}^{\pi}d\theta
P_{\sigma}(\theta)e^{i(n-m)\theta}$.

Figure 7 shows a slight asymmetry of the reconstructed state due to the
systematic shift of the variance sum: this effect is also observed in Fig. 4.
The quantitative agreement between the simulated density matrix $\rho_{\rm
sim}$ and the experimental density matrix $\rho_{\rm exp}$ is excellent, with
a quantum fidelity ${\rm Tr}[\sqrt{\sqrt{\rho_{\rm sim}}\rho_{\rm
exp}\sqrt{\rho_{\rm sim}}}]\approx 0.9$.

## References

  * Einstein  _et al._ (1935) A. Einstein, B. Podolsky, and N. Rosen, “Can quantum-mechanical description of physical reality be considered complete?” Phys. Rev. 47, 777–780 (1935).
  * Furry (1936) W. H. Furry, “Note on the quantum-mechanical theory of measurement,” Phys. Rev. 49, 393–399 (1936).
  * Reid (1989) M. D. Reid, “Demonstration of the Einstein-Podolsky-Rosen paradox using nondegenerate parametric amplification,” Phys. Rev. A 40, 913–923 (1989).
  * Schrödinger (1935) E. Schrödinger, “Discussion of probability relations between separated systems,” Proc. Cambridge Philos. Soc. 31, 555–563 (1935).
  * Schrödinger (1936) E. Schrödinger, “Probability relations between separated systems,” Proc. Cambridge Philos. Soc. 32, 446–452 (1936).
  * Wiseman  _et al._ (2007) H. M. Wiseman, S. J. Jones, and A. C. Doherty, “Steering, entanglement, nonlocality, and the Einstein-Podolsky-Rosen paradox,” Phys. Rev. Lett. 98, 140402 (2007).
  * Braunstein and van Loock (2005) S. L. Braunstein and P. van Loock, “Quantum information with continuous variables,” Rev. Mod. Phys. 77, 513–577 (2005).
  * Pezzè and Smerzi (2014) L. Pezzè and A. Smerzi, “Quantum theory of phase estimation,” in  _Atom Interferometry_, Proceedings of the International School of Physics ”Enrico Fermi”, Course 188, Varenna, edited by G.M. Tino and M.A. Kasevich (IOS Press, Amsterdam, 2014) pp. 691–741.
  * Opanchuk  _et al._ (2014) B. Opanchuk, L. Arnaud, and M. D. Reid, “Detecting faked continuous-variable entanglement using one-sided device-independent entanglement witnesses,” Phys. Rev. A 89, 062101 (2014).
  * Ou  _et al._ (1992) Z. Y. Ou, S. F. Pereira, H. J. Kimble, and K. C. Peng, “Realization of the Einstein-Podolsky-Rosen paradox for continuous variables,” Phys. Rev. Lett. 68, 3663–3666 (1992).
  * D’Auria  _et al._ (2009) V. D’Auria, S. Fornaro, A. Porzio, S. Solimeno, S. Olivares, and M. G. A. Paris, “Full characterization of Gaussian bipartite entangled states by a single homodyne detector,” Phys. Rev. Lett. 102, 020502 (2009).
  * Reid  _et al._ (2009) M. D. Reid, P. D. Drummond, W. P. Bowen, E. G. Cavalcanti, P. K. Lam, H. A. Bachor, U. L. Andersen, and G. Leuchs, “Colloquium : The Einstein-Podolsky-Rosen paradox: From concepts to applications,” Rev. Mod. Phys. 81, 1727–1751 (2009).
  * Estève  _et al._ (2008) J. Estève, C. Gross, A. Weller, S. Giovanazzi, and M. K. Oberthaler, “Squeezing and entanglement in a Bose-Einstein condensate,” Nature 455, 1216–1219 (2008).
  * Wasilewski  _et al._ (2010) W. Wasilewski, K. Jensen, H. Krauter, J. J. Renema, M. V. Balabas, and E. S. Polzik, “Quantum noise limited and entanglement-assisted magnetometry,” Phys. Rev. Lett. 104, 133601 (2010).
  * Berrada  _et al._ (2013) T. Berrada, S. van Frank, R. B cker, T. Schumm, J.-F. Schaff, and J. Schmiedmayer, “Integrated Mach-Zehnder interferometer for Bose-Einstein condensates,” Nature Commun. 4, 2077 (2013).
  * Lücke  _et al._ (2014) B. Lücke, J. Peise, G. Vitagliano, J. Arlt, L. Santos, G. Tóth, and C. Klempt, “Detecting multiparticle entanglement of Dicke states,” Phys. Rev. Lett. 112, 155304 (2014).
  * Julsgaard  _et al._ (2001) B. Julsgaard, A. Kozhekin, and E. S. Polzik, “Experimental long-lived entanglement of two macroscopic objects,” Nature 413, 400–403 (2001).
  * Appel  _et al._ (2009) J. Appel, P. J. Windpassinger, D. Oblak, U. B. Hoff, N. Kærgaard, and E. S. Polzik, “Mesoscopic atomic entanglement for precision measurements beyond the standard quantum limit,” Proc. Natl. Acad. Sci. U. S. A. 106, 10960–10965 (2009).
  * Chen  _et al._ (2011) Z. Chen, J. Bohnet, S. Sankar, J. Dai, and J. Thompson, “Conditional spin squeezing of a large ensemble via the vacuum Rabi splitting,” Phys. Rev. Lett. 106, 133601 (2011).
  * Schleier-Smith  _et al._ (2010) M. H. Schleier-Smith, I. D. Leroux, and V. Vuletić, “States of an ensemble of two-level atoms with reduced quantum uncertainty,” Phys. Rev. Lett. 104, 073604 (2010).
  * Haas  _et al._ (2014) F. Haas, J. Volz, R. Gehr, J. Reichel, and J. Estève, “Entangled states of more than 40 atoms in an optical fiber cavity,” Science 344, 180–183 (2014).
  * Behbood  _et al._ (2014) N. Behbood, F. M. Ciurana, G. Colangelo, M. Napolitano, G. Tóth, R. J. Sewell, and M. W. Mitchell, “Generation of macroscopic singlet states in a cold atomic ensemble,” Phys. Rev. Lett. 113, 093601 (2014).
  * Riedel  _et al._ (2010) M. Riedel, P. Böhi, Y. Li, T. Hänsch, A. Sinatra, and P. Treutlein, “Atom-chip-based generation of entanglement for quantum metrology,” Nature 464, 1170–1173 (2010).
  * Hamley  _et al._ (2012) C. D. Hamley, C. S. Gerving, T. M. Hoang, E. M. Bookjans, and M. S. Chapman, “Spin-nematic squeezed vacuum in a quantum gas,” Nature Phys. 8, 305–308 (2012).
  * Strobel  _et al._ (2014) H. Strobel, W. Muessel, D. Linnemann, T. Zibold, D. B. Hume, L. Pezzè, A. Smerzi, and M. K. Oberthaler, “Fisher information and entanglement of non-Gaussian spin states,” Science 345, 424–427 (2014).
  * Duan  _et al._ (2000) L.-M. Duan, G. Giedke, J. I. Cirac, and P. Zoller, “Inseparability criterion for continuous variable systems,” Phys. Rev. Lett. 84, 2722–2725 (2000).
  * Simon (2000) R. Simon, “Peres-Horodecki separability criterion for continuous variable systems,” Phys. Rev. Lett. 84, 2726–2729 (2000).
  * Gross  _et al._ (2011) C. Gross, H. Strobel, E. Nicklas, T. Zibold, N. Bar-Gill, G. Kurizki, and M. K. Oberthaler, “Atomic homodyne detection of continuous-variable entangled twin-atom states,” Nature 480, 219–223 (2011).
  * Peise  _et al._ (2015) J. Peise, B. Lücke, L. Pezzè, F. Deuretzbacher, W. Ertmer, J. Arlt, A. Smerzi, L. Santos, and C. Klempt, “Interaction-free measurements by quantum Zeno stabilisation of ultracold atoms,” Nat. Commun. 6, 6811 (2015).
  * Schmied and Treutlein (2011) R. Schmied and P. Treutlein, “Tomographic reconstruction of the Wigner function on the Bloch sphere,” New J. Phys. 13, 065019– (2011).
  * Lvovsky and Raymer (2009) A. I. Lvovsky and M. G. Raymer, “Continuous-variable optical quantum-state tomography,” Rev. Mod. Phys. 81, 299–332 (2009).
  * Vasilyev  _et al._ (2000) M. Vasilyev, S.-K. Choi, P. Kumar, and G. M. D’Ariano, “Tomographic measurement of joint photon statistics of the twin-beam quantum state,” Phys. Rev. Lett. 84, 2354–2357 (2000).
  * Hradil  _et al._ (2004) Z. Hradil, J. Řeháček, J. Fiurášek, and M. Ježek, “3 maximum-likelihood methods in quantum mechanics,” in  _Quantum State Estimation_, Lecture Notes in Physics, Vol. 649, edited by Matteo Paris and Jaroslav Řeháček (Springer Berlin Heidelberg, 2004) pp. 59–112.
  * Pezzè and Smerzi (2009) L. Pezzè and A. Smerzi, “Entanglement, nonlinear dynamics, and the Heisenberg limit,” Phys. Rev. Lett. 102, 100401 (2009).
  * Lewis-Swan and Kheruntsyan (2013) R. J. Lewis-Swan and K. V. Kheruntsyan, “Sensitivity to thermal noise of atomic Einstein-Podolsky-Rosen entanglement,” Phys. Rev. A 87, 063635 (2013).
  * Klempt  _et al._ (2010) C. Klempt, O. Topic, G. Gebreyesus, M. Scherer, T. Henninger, P. Hyllus, W. Ertmer, L. Santos, and J. J. Arlt, “Parametric amplification of vacuum fluctuations in a spinor condensate,” Phys. Rev. Lett. 104, 195303 (2010).
  * Ferris  _et al._ (2008) Andrew J. Ferris, Murray K. Olsen, Eric G. Cavalcanti, and Matthew J. Davis, “Detection of continuous variable entanglement without coherent local oscillators,” Phys. Rev. A 78, 060104 (2008).
  * Raymer  _et al._ (2003) M. G. Raymer, A. C. Funk, B. C. Sanders, and H. de Guise, “Separability criterion for separate quantum systems,” Phys. Rev. A 67, 052104 (2003).
  * Gühne and Tóth (2009) Otfried Gühne and Géza Tóth, “Entanglement detection,” Physics Reports 474, 1–75 (2009).
  * Pezzè  _et al._ (2015) Luca Pezzè, Philipp Hyllus, and Augusto Smerzi, “Phase-sensitivity bounds for two-mode interferometers,” Phys. Rev. A 91, 032103 (2015).

##  IV Acknowledgments

We thank G. Tóth for inspiring discussions. We thank W. Schleich for a review
of our manuscript. We acknowledge support from the Centre for Quantum
Engineering and Space-Time Research (QUEST) and from the Deutsche
Forschungsgemeinschaft (Research Training Group 1729). We acknowledge support
from the European Metrology Research Programme (EMRP). The EMRP is jointly
funded by the EMRP participating countries within EURAMET and the European
Union. L.P. and A.S. acknowledge financial support by the EU-STREP project
QIBEC. J.A. acknowledges support by the Lundbeck Foundation.