Introduction

In their original publication1, Einstein, Podolsky and Rosen describe two particles A and B with correlated position x B=x A+x 0 and anti-correlated momentum p B=−p A (Fig. 1a). When coordinates x A and p A are measured in independent realizations of the same state, the correlations allow for an exact prediction of x B and p 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 B and p B, since their measurement uncertainties are subject to the Heisenberg relation . EPR thus concluded that quantum mechanics is incomplete—under their assumptions that are today known as 'local realism'. Later, the notion of EPR correlations was generalized to a more realistic scenario, yielding a criterion2,3 for the uncertainties , of the inferred predictions for x B and p B. The EPR criterion is met if these uncertainties violate the Heisenberg inequality for the inferred uncertainties . The EPR criterion also certifies steering, a concept termed by Schrödinger4,5 in response to EPR, which has attracted a lot of interest in the past years6. 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 teleportation7 and quantum metrology8. 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 verification9.

Figure 1: Einstein–Podolsky–Rosen correlations.
figure 1

(a) EPR's original work describes two particles A and B with maximally correlated position and momentum coordinates x A/B and p A/B. (b) Spin dynamics in a Bose–Einstein condensate can be used to create EPR correlations between N A/B atoms in two different Zeeman levels A and B. The correlations appear in amplitude x A/B and phase p A/B quadratures that are defined as a function of the creation and annihilation operators and a A/B in the two modes.

Full size image

Up to now, the creation of continuous-variable entangled states satisfying the EPR criterion was only achieved in optical systems. In a seminal publication10, 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 investigations11,12, continuous variables are represented by amplitude x A/B and phase p A/B quadratures, satisfying the commutation relation [x A/B, p A/B]=i. These quadratures can be measured accurately by optical homodyning. The correlations are captured by the four two-mode variances and . These variances were proven to fulfil a symmetric form of Reid's inequality3 , which is a sufficient EPR criterion since and . In recent years, continuous-variable entangled optical states have been applied for proof-of-principle quantum computation and communication tasks7. 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 realism2,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 atoms13,14,15,16. They have been created by atom–light interaction at room temperature14,17, in cold samples18,19,20,21,22 and by collisional interactions in Bose–Einstein condensates13,16,23,24,25. For Gaussian states of two collective atomic modes, the inseparability criterion26,27 has been used to demonstrate entanglement14,17,28, but the strong correlations necessary to meet the more demanding EPR criterion 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 s.d. 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.

Results

Two-mode squeezed vacuum

In our experiments, a BEC with 2 × 104 87Rb atoms in the Zeeman level (F, m F )=(1, 0) generates atom pairs in the levels (1, ±1) due to spin-changing collisions (Fig. 1b), ideally yielding the two-mode squeezed state

where ξ=Ωt is the squeezing parameter, which depends on the spin dynamics rate Ω=2π × 5.1 Hz and the spin dynamics duration t=26 ms. The notation represents a two-mode Fock state in the two Zeeman levels (1, ±1). The generated two-mode squeezed state can be characterized by the quadratures and for the two levels (1, ±1). These exhibit EPR correlations, since the variances are squeezed, while the conjugate variances are anti-squeezed. The state fulfills Reid's EPR criterion for which corresponds to a spin dynamics duration of >11 ms. In the limit of large squeezing, our set-up 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 scenario29. A small radiofrequency pulse couples 15% of the BEC in the level (1, 0) (the local oscillator) symmetrically to the two modes (1, ±1). The local oscillator phase is represented by the BEC phase relative to the phase sum of the two ensembles in (1, ±1). It can be varied in our experiments by shifting the energy of the level (1, −1) for an adjustable time. From the measured number of atoms in both levels, we obtain a linear combination of the quadratures according to . Figure 2a shows two-dimensional histograms of these measurements for and , 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 , are shown in Fig. 2b and reveal the expected two-mode squeezing behaviour. From these measurements, we quantify the EPR entanglement by Reid's criterion, yielding , which is 2.4 s.d. below the limit of . In addition, the data also fulfil the inseparability criterion as , which is 15 s.d. below the classical limit of 2 (Fig. 2d), and meets the criterion for a symmetric ('two-way') steering between the systems6. We estimate that the product value could be reduced to if the radiofrequency intensity noise was eliminated by stabilization or postcorrection. The experimental creation of entangled massive particles that satisfy the continuous-variable EPR criterion presents the main result of this publication.

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 θ≈3π/4 and θ≈5π/4 to record only the x- and p-variances. As a function of the evolution time, the variances are squeezed below the vacuum reference of 1, while the variances exhibit an anti-squeezing behaviour (Fig. 3a,b). From these data, we extract the EPR parameter , as a function of evolution time (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 were well reproduced by a simple noise model, which includes a radiofrequency intensity noise of 0.4% and a local oscillator phase noise of 0.044π (see Methods).

Figure 3: Squeezing dynamics.
figure 3

(a,b) Time evolution of the four two-mode quadrature variances , during spin dynamics. The data follow the ideal squeezing/anti-squeezing behaviour according to the independently measured spin dynamics rate (dashed line). (c) Time evolution of the EPR parameter 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 s.d.).

Full size image

Full state reconstruction

The total of 2,864 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 function30 or the Husimi Q-distribution21,25. However, both methods yield a characterization of the state's projection on the fully symmetric subspace only. The well-developed methods in quantum optics31 allowed for a full reconstruction of an optical two-mode squeezed state by homodyne tomography11,32. 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 reconstruction31,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 (Fig. 2a) presenting the marginal distributions for the x A/B and p A/B variables. The reconstructed state is the one that optimally reproduces the measured histograms by a superposition of harmonic oscillator wave functions31. 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 shows the strongest contribution besides the vacuum state. Likewise, the twin Fock states to 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,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.

Discussion

For an evaluation of the created state, we have extracted a logarithmic negativity of 1.43±0.06 from the reconstructed density matrix. This value is above the threshold of zero for separable states and signals non-separability of the reconstructed state. The quantum Fisher information34 F Q for the state projected on fixed-N subspaces reveals that , where is the average number of particles. Since the state is a resource for quantum enhanced metrology34. Furthermore, we fit an ideal two-mode squeezed state with variable squeezing parameter ξ 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 ξ 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. First, the purity of the reconstructed state is limited by the finite number of homodyne measurements. Second, 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. Third, a small drift of the radiofrequency coupling strength during homodyning virtually increases the variance in the (x A+x B) and the (p A+p 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 (2, ±1) which have the same magnetic moment as the two EPR modes. Furthermore, this set-up 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.

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π × (280, 220, 180) Hz. At a homogeneous magnetic field of 2.6 G with fluctuations of about 70 μG, the condensate is transferred from the level (F, m F )=(2, 2) to the level (1, 0) 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 (1, ±1) are emptied with a pair of microwave π-pulses from (1, +1) to (2, +2) and from (1, −1) to (2, −2) 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 signal35.

Figure 5 shows a schematic overview of the following experimental sequence. A microwave frequency which is red-detuned to the transition between the levels (1, −1) and (2, −2) by about 208 kHz is adiabatically ramped on within 675 μs. The microwave shifts the level (1, −1) 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 addressed16,36. Each resonance condition belongs to a specific spatial mode of the states (1, ±1) to which the atoms are transferred. If the energy of the level (1, −1) is reduced, then 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 (1, ±1) 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 (1, 0) at the energy of the chemical potential) is exactly degenerate with the output state (two atoms in the levels (1, ±1) 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 (1, ±1) (see ref. 29, Fig. 3). Subsequently, the microwave dressing field is ramped down within 675 μs, stays off for a variable duration between 25 and 1,150 μ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.

Figure 5: Schematic overview of the experimental sequence.
figure 5

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 on 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 Hz. At time t 3 a radiofrequency pulse couples the m F =0 condensate to the m F =±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.

Full size image

For the atomic homodyning, a radiofrequency pulse with a frequency of 1.834 MHz and a duration of τ=30 μs couples the level (1, 0) with the levels (1, ±1). The microwave dressing field is chosen such that both radiofrequency 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 charge-coupled device 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 <1%.

Three-mode unbalanced homodyning

The radiofrequency coupling is described by the three-mode unitary operation , where

and Ω±1 are Rabi frequencies for the (1, 0)↔(1, ±1) transition (in general Ω+1≠Ω−1). To calculate the mode transformation, we use , and . We have

where c=cos(Ωτ/2), s=sin(Ωτ/2), and

are rescaled Rabi frequencies. Below, we illustrate how the measurement of the number of particles in the m F =±1 mode after the radiofrequency coupling, and , gives access to the number conserving quadratures

where is the average number of particles in the condensate before homodyne (similarly, ). In our experiment, . We thus neglect fluctuations of the number of particles in the m F =0 mode, replacing with its mean value n 0≈n A+n 0+n B=N tot.

Number difference

The quadrature difference can be experimentally obtained by measuring the difference of the number of particles in the ±1 modes. From equation (2) we can directly calculate N A−N B. To the leading order in n 0, we have

Since Ω+1 and Ω−1 only differ by 1.7% in our experiments, and , we can simplify this equation and obtain

Number sum

The quadrature sum is obtained by adding the number of particles in the ±1 modes after homodyning:

Taking , we have

Finally, the mean transfer of particles from m F =0 to m F =±1 and the mean number difference is used to calculate

Observing a transfer of 15% of the atoms from m F =0 to m F =±1 we deduce c 2≈0.85.

To summarize, equations (6) and (8) are used to experimentally obtain p A±p B from the measurement of the number of particles in the output modes. The quadratures x A±x B are obtained with the same method, by applying a relative π/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, we have37,38

where and are the variances of quadrature sum and difference, respectively. A violation of equation (10) signals non-separability, that is, ρ≠ρ sep. Equation (10) generalizes the criterion of refs 26, 27 that was derived for standard quadrature operators (that is, when the m F =0 mode is treated parametrically, the operator a 0 being replaced by ).

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 A±x B and . We point out that not all non-separable states fulfil Reid's criterion. The position–momentum quadratures for the B mode satisfy the commutation relation . The corresponding Heisenberg uncertainty relation is . Let us introduce the quantities x ext(x A) and p ext(p A), which are the estimate of x B and p B on system B, respectively, given the results of quadrature measurements on the system A. We then indicate as the squared deviation of the estimate from the actual value, averaged over all possible results x A,

and similarly for , where P(x A, x B) is the joint probability. Reid's criterion thus reads37 . Taking and , where the bar indicates statistical average, Reid's criterion translates into a condition for the product of quadrature variances:

In our case, . Therefore, corrections in equations (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, 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, 33. We point out that our state reconstruction is performed without any assumtpions 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=2,864 measurements of the quadratures x A and x B at different values of the local oscillator phase θ 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 A,B to have a discrete spectrum. To simplify the notation, let us indicate as x the square bin [x A, x A+dx], [x B, x B+dx]. Given a quantum state ρ (unknown here), the probability to observe a certain sequence of results (n x, θ measurement in the bin x, when the phase value is set to θ, with ) is

indicated as likelihood function. In equation (13), is the joint probability, is the conditional probability, with , , and P(θ) is the fraction of measurements done when the phase is equal to θ. The maximum likelihood (ML)

is the state that maximizes on the manifold of density matrices. To find the ML we use the chain of inequalities31,33

where a x,θ are arbitrary positive numbers ( a ={a x, θ } is the corresponding vector), f x, θ =n x, θ /N are relative frequencies ( f ={f x, θ } is the corresponding vector), and

is a non-negative operator with largest eigenvalue λ( f , a ). The second inequality is saturated by taking ρ=ρ ML with support on the subspace corresponding to the maximum eigenvalue of R: R( f , a )ρ ML=λ( f , a )ρ 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,θ =P(x, θ), which also implies Tr[R( f , a )ρ ML]=1 and thus λ( f , a )=1. In conclusion, the search for the ML can be recast in the operator equation Rρ ML=ρ ML or, equivalently (since R and ρ ML are Hermitian operators),

Numerically, this equation is solved iteratively: we start the protocol from a unit matrix and apply repetitive iterations according to equation (17), being the kth step of the algorithm, where denotes normalization to unit trace. The convergence (which is not guaranteed in general) is checked. The method guarantees that ρ ML is a non-negative definite operator. In practical implementations, it is most convenient to work in the atom-number representation and write , where n cut is a cutoff (in our case n cut≈10). We use , 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 , equation (1). The simulation follows three steps: (i) we generate distributions for the quadratures x A, B at different values of the phase shift, according to the probability ; (ii) we generate p random quadrature data for each θ value (for a total of N=p × n θ , where n θ is the number of θ 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, ρ ML, and the two-mode squeezed vacuum state, . When the number of measurements p per θ value is increased, the fidelity converges to an asymptotic value . The asymptotic fidelity tends to 1 when decreasing the bin size dx.

Figure 6: Simulated quantum state tomography of an ideal-squeezed state.
figure 6

Here we show the quantum fidelity between the reconstructed state ρ ML and the ideal two-mode squeezed vacuum state, as a function of the number of measurements p per θ value (here we consider n θ =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 s.d.) are obtained with a bootstrap method.

Full size image

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 , where ρ ppt is the partial transpose of ρ ML. Mode-entanglement is obtained for39 E[ρ]>0. The QFI for the state projected over subspaces of a fixed number of particles n, , is given by40

where is in diagonal form and J r is the collective pseudo-spin operator (pointing along an arbitrary direction r in the three-dimensional space). The QFI is then maximized over r , for further details see ref. 8. Particle entanglement, useful for sub-shot-noise metrology, is obtained for40 , where 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→∞ and dx→0, the logarithmic negativity and QFI converge to and , respectively, which are analytical values calculated for the two-mode squeezed vacuum state.

Noise model and simulation of noisy state reconstruction

The main sources of noise in our apparatus are phase fluctuations and noise of the radiofrequency coupling strength. The phase noise is assumed to have a Gaussian distribution and we estimate a width σ≈0.36. Correlated fluctuations of Ω+1 and Ω−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

where . 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 (that is, 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 equation (19) as , where .

Figure 7: Simulated quantum state tomography of a noisy squeezed state.
figure 7

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≈100 measurements, n θ =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%).

Full size image

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 ρ sim and the experimental density matrix ρ exp is excellent, with a quantum fidelity .

Additional information

How to cite this article: Peise, J. et al. Satisfying the Einstein–Podolsky–Rosen criterion with massive particles. Nat. Commun. 6:8984 doi: 10.1038/ncomms9984 (2015).