EPR-entangled Bose-Einstein condensates in state-dependent potentials: a dynamical study

TL;DR

This paper investigates how atomic interactions in bi-modal Bose-Einstein Condensates within state-dependent potentials can generate non-local correlations, enabling EPR steering experiments despite spatial dynamics.

Contribution

It introduces a combined Fock state and time-dependent Hartree-Fock approach to account for spatial and quantum correlations in BECs, demonstrating feasible entanglement generation.

Findings

Successful generation of non-local entanglement in BECs

Potential for EPR steering with spatially separated condensates

Spatial dynamics do not hinder entanglement creation

Abstract

We study generation of non-local correlations by atomic interactions in a pair of bi-modal Bose-Einstein Condensates in state-dependent potentials including spatial dynamics. The wave-functions of the four components are described by combining a Fock state expansion with a time-dependent Hartree-Fock Ansatz, so that both the spatial dynamics and the local and non-local quantum correlations are accounted for. We find that despite the spatial dynamics, our protocole generates enough non-local entanglement to perform an EPR steering experiment with two spatially separated con-densates of a few thousands of atoms.

Equations

$$\hat{\psi }(\text{r})=\sum _{\alpha \in A}{\hat{\psi }}_{\alpha }(\text{r})$$

$$[{\hat{\psi }}_{\alpha }(\text{r}),{\hat{\psi }}_{{\alpha }^{\prime }}^{†}({\text{r}}^{\prime })]=\frac{{\delta }_{\alpha ,{\alpha }^{\prime }}{\delta }_{\text{r},{\text{r}}^{\prime }}}{l^3}$$

$$\begin{gathered} \hat{H}=l^3\sum _{\text{r},\alpha \in A}{\hat{\psi }}_{\alpha }^{†}(\text{r})\left(-\frac{{\hslash }^2}{2m}{\Delta }_{\text{r}}+V_{\alpha }(\text{r},t)\right){\hat{\psi }}_{\alpha }(\text{r}) \\ +l^3\sum _{\begin{array}{c}\text{r}\\ \alpha ,{\alpha }^{\prime }\in A\end{array}}\frac{g_{\alpha {\alpha }^{\prime }}}{2}{\hat{\psi }}_{\alpha }^{†}(\text{r}){\hat{\psi }}_{{\alpha }^{\prime }}^{†}(\text{r}){\hat{\psi }}_{\alpha }(\text{r}){\hat{\psi }}_{{\alpha }^{\prime }}(\text{r}) \end{gathered}$$

$${\varphi }_{\alpha }={\varphi }_{\alpha }(\vec{N},\text{r},t),\vec{N}=(N_{{\alpha }^{\prime }}{)}_{{\alpha }^{\prime }\in A}$$

$$\forall \vec{N},\forall t,l^3\sum _{\text{r}}{\left|{\varphi }_{\alpha }(\vec{N},\text{r},t)\right|}^2=1$$

$${\hat{a}}_{\alpha }(\vec{N},t)\equiv l^3\sum _{\text{r}}{\varphi }_{\alpha }^{\ast }(\vec{N},\text{r},t){\hat{\psi }}_{\alpha }(\text{r})$$

$$|\{N_{\alpha }:{\varphi }_{\alpha }(\vec{N},\text{r},t)\}⟩\equiv \frac{{\prod }_{\alpha }{\left({\hat{a}}_{\alpha }^{†}(\vec{N},t)\right)}^{N_{\alpha }}}{{\left({\prod }_{\alpha }{N}_{\alpha }!\right)}^{1/2}}|\mathrm{vac}⟩$$

$$S\equiv {\int }_{t_{\mathrm{i}}}^{t_{\mathrm{f}}}\mathrm{d}t\left\{\frac{\text{i}\hslash }{2}\left(⟨\{N_{\alpha }:{\varphi }_{\alpha }(\vec{N},\text{r},t)\}|\frac{\mathrm{d}}{\mathrm{d}t}|\{N_{\alpha }:{\varphi }_{\alpha }(\vec{N},\text{r},t)\}⟩-c.c.\right)-E(\vec{N},t)\right\}$$

$$⟨\{N_{\alpha }:{\varphi }_{\alpha }(\vec{N},\text{r},t)\}|\frac{\mathrm{d}}{\mathrm{d}t}|\{N_{\alpha }:{\varphi }_{\alpha }(\vec{N},\text{r},t)\}⟩-c.c.=2l^3\sum _{\text{r},\alpha }{N}_{\alpha }{\varphi }_{\alpha }^{\ast }(\vec{N},\text{r},t)\frac{\mathrm{d}}{\mathrm{d}t}{\varphi }_{\alpha }(\vec{N},\text{r},t)$$

$$\begin{gathered} E(\vec{N},t)\equiv ⟨\{N_{\alpha }:{\varphi }_{\alpha }(\vec{N},\text{r},t)\}|\hat{H}|\{N_{\alpha }:{\varphi }_{\alpha }(\vec{N},\text{r},t)\}⟩ \\ =l^3\sum _{\text{r},\alpha }{N}_{\alpha }{\varphi }_{\alpha }^{\ast }(\vec{N},\text{r},t)h_{\alpha }{\varphi }_{\alpha }(\vec{N},\text{r},t)+l^3\sum _{\text{r},\alpha }\frac{g_{\alpha \alpha }}{2}{N}_{\alpha }(N_{\alpha }-1)|{\varphi }_{\alpha }(\vec{N},\text{r},t){|}^4+l^3\sum _{\text{r},\alpha ,{\alpha }^{\prime }\ne \alpha }\frac{g_{\alpha {\alpha }^{\prime }}}{2}{N}_{\alpha }{N}_{{\alpha }^{\prime }}|{\varphi }_{\alpha }(\vec{N},\text{r},t){|}^2|{\varphi }_{{\alpha }^{\prime }}(\vec{N},\text{r},t){|}^2 \end{gathered}$$

$$h_{\alpha }(\text{r},t)=-\frac{{\hslash }^2}{2m}{\Delta }_{\text{r}}+V_{\alpha }(\text{r},t)$$

$$\text{i}\hslash \frac{\mathrm{d}{\varphi }_{\alpha }(\vec{N},\text{r},t)}{\mathrm{d}t}=\left(h_{\alpha }(\text{r},t)+g_{\alpha \alpha }(N_{\alpha }-1){\left|{\varphi }_{\alpha }(\vec{N},\text{r},t)\right|}^2+\sum _{{\alpha }^{\prime }\ne \alpha }{g}_{\alpha {\alpha }^{\prime }}{N}_{{\alpha }^{\prime }}{\left|{\varphi }_{{\alpha }^{\prime }}(\vec{N},\text{r},t)\right|}^2\right){\varphi }_{\alpha }(\vec{N},\text{r},t)$$

$$|\psi (0)⟩=|\{N_{\alpha }:{\varphi }_{\alpha }(\vec{N},\text{r},0)\}⟩$$

$${\text{e}}^{-\text{i}\hat{H}t/\hslash }|\psi (0)⟩\simeq {\text{e}}^{-\text{i}A(\vec{N},t)}|\{N_{\alpha }:{\varphi }_{\alpha }(\vec{N},\text{r},t)\}⟩$$

$$\begin{gathered} {\text{e}}^{-\text{i}A(\vec{N},t)}= \\ ⟨\{N_{\alpha }:{\varphi }_{\alpha }(\vec{N},\text{r},t)\}|{\text{e}}^{-\text{i}\hat{H}t/\hslash }|\{N_{\alpha }:{\varphi }_{\alpha }(\vec{N},\text{r},0)\}⟩ \end{gathered}$$

$$\hslash \frac{\mathrm{d}A(\vec{N},t)}{\mathrm{d}t}=-l^3\sum _{\text{r},\alpha }\frac{g_{\alpha \alpha }}{2}{N}_{\alpha }(N_{\alpha }-1){\left|{\varphi }_{\alpha }(\vec{N},\text{r},t)\right|}^4-l^3\sum _{\text{r},\alpha ,{\alpha }^{\prime }\ne \alpha }\frac{g_{\alpha {\alpha }^{\prime }}}{2}{N}_{\alpha }{N}_{{\alpha }^{\prime }}{\left|{\varphi }_{\alpha }(\vec{N},\text{r},t)\right|}^2{\left|{\varphi }_{{\alpha }^{\prime }}(\vec{N},\text{r},t)\right|}^2$$

$$\alpha =\sigma ,ϵ\text{mbox}with\sigma =a,b\text{mbox}andϵ=0,1$$

$$|\Psi (0)⟩$$

$$=$$

$$\mathcal{N}=\{\vec{N}\in {\mathbb{N}}^4,N_{a,0}+N_{a,1}=N_a,\text{and}{N}_{b,0}+N_{b,1}=N_b\}$$

$$\begin{gathered} |\Psi (t)⟩=\sum _{\vec{N}\in \mathcal{N}}{\left(\frac{N_a!N_b!}{{\prod }_{\alpha }{N}_{\alpha }!}\right)}^{1/2}{\text{e}}^{-\text{i}A(\vec{N},t)} \\ \times \prod _{\alpha }{C}_{\alpha }^{N_{\alpha }}|\left\{N_{\alpha }:{\varphi }_{\alpha }(\vec{N},\text{r},t)\right\}⟩ \end{gathered}$$

$${\hat{S}}_x^{\sigma }$$

$${\hat{S}}_y^{\sigma }$$

$${\hat{S}}_z^{\sigma }$$

$${\hat{S}}_{x\varphi }^{\sigma }$$

$${\hat{S}}_{y\varphi }^{\sigma }$$

$$\text{mbox}tan{\varphi }_{\sigma }=\frac{⟨{\hat{S}}_y^{\sigma }⟩}{⟨{\hat{S}}_x^{\sigma }⟩}$$

$${\hat{S}}_{\alpha }^a$$

$${\hat{S}}_{\beta }^b$$

$$E_{\mathrm{EPR}}^2\equiv \frac{4\left({\Delta }^2{\hat{S}}_{\alpha }^a{\Delta }^2{\hat{S}}_{\beta }^b-{\text{Covar}}^2({\hat{S}}_{\alpha }^a,{\hat{S}}_{\beta }^b)\right)\left({\Delta }^2{\hat{S}}_{\alpha +\pi /2}^a{\Delta }^2{\hat{S}}_{\beta +\pi /2}^b-{\text{Covar}}^2({\hat{S}}_{\alpha +\pi /2}^a,{\hat{S}}_{\beta +\pi /2}^b)\right)}{\left({\Delta }^2{\hat{S}}_{\alpha }^a{\Delta }^2{\hat{S}}_{\alpha +\pi /2}^a\right){\left|⟨{\hat{S}}_{x\varphi }^b⟩\right|}^2}$$

$${\Delta }^2({\hat{S}}_{\beta }^b-{\hat{S}}_{\beta }^{b,\inf }){\Delta }^2({\hat{S}}_{\beta +\pi /2}^b-{\hat{S}}_{\beta +\pi /2}^{b,\inf })<\frac{1}{4}{\left|⟨{\hat{S}}_{x\varphi }^b⟩\right|}^2$$

$${\hat{S}}_{\beta }^{b,\inf }=⟨{\hat{S}}_{\beta }^b⟩+\frac{\text{mbox}Covar({\hat{S}}_{\alpha }^a,{\hat{S}}_{\beta }^b)}{{\Delta }^2{\hat{S}}_{\alpha }^a}\left({\hat{S}}_{\alpha }^a-⟨{\hat{S}}_{\alpha }^a⟩\right)$$

$${\varphi }_{\alpha }(\vec{N},\text{r},t)$$

Full text

Introduction

The non-local and non-deterministic nature of quantum mechanics is a subject of experimental investigation since a few decades. The most significant result in this direction was the violation of the so-called Bell inequalities with entangled pairs of photon Aspect1982 . With all loopholes now filled in Hanson2015 , this experiment ruled out the possibility that the phenomena described by quantum mechanics could by explained by an underlying local deterministic “hidden variable” theory. More recently, the observation of non-local entanglement between macroscopic massive objects has emerged as a new objective. This result would be an important step in further ruling out local realism, and in pushing back the boundary of observation of quantum effects in the macroscopic world. Experimentally, this is a challenging objective for two reasons: $(i)$ large systems decohere usually faster than individual systems and therefore retain their non-local entanglement only for a short amount of time, $(ii)$ many-body systems have a complex internal dynamics which is likely to affect their non-local entanglement.

Cold atoms experiments provide a promising platform to tackle these limitations and observe non-local entanglement between reasonably large objects Treutlein2016 . They provide experimentalists with clean isolated systems where the impact of decoherence can be fought more efficiently than in other many-body systems ; in particular, the dominant decoherence effect, atomic losses, can be made very weak for a careful choice of the atomic states and can be accurately evaluated theoretically CastinSinatra2008 . At low temperatures, one can prepare a system of bosonic atoms into a few macroscopically occupied modes, a situation which enables a natural generalization of the bipartite system imagined in a Bell-experiment to a macroscopic system. Finally, the van der Walls interactions between the atoms are a powerful source of entanglement creation, able to produce highly entangled many-body states at short times, when the decoherence effects are limited.

The correlations between the particles in a condensate of weakly-interacting bosons have recently been shown to be strong enough Sangouard2016 to allow in principle for a violation of multipartite Bell inequality Acin2014 . However a proper violation would require to perform measurements on each pair separately, which is experimentally challenging. Instead, it was recently proposed Drummond2011 ; Kurizki2011 to violate a weaker form of non-locality usually referred to as “EPR steering”, by reference to the seminal article Rosen1935 which introduced the EPR paradox. In an EPR situation, the measurements done by Alice on her half of an entangled non-local quantum state are shown to apparently “steer” Bob’s second half of the state, eventually leading to a violation of the Heisenberg uncertainty relation which Bob’s observables should satisfy if the quantum state were purely local, the so-called “EPR inequality” Reid1989 ; Doherty2007 . Such an experiment shows that the non-deterministic nature of Quantum Mechanics (in particular the Heisenberg uncertainty relation) is incompatible with locality, thus ruling out all hidden-variable theories locally compatible with Quantum Mechanics.

Already a violation of an EPR inequality Doherty2007 was obtained in a condensate of weakly-interacting bosons using state-changing collisions to create entangled subsystems of a few atoms Klempt2015 . An alternative route relying on light-matter interaction to create entanglement was explored in Ref.Polzik2001 and could be used in an EPR experiment. Here we consider the entanglement scheme proposed by Refs.KPS2013 and Byrnes2013 where two Bose-Einstein condensates $a$ and $b$ of atoms in a superposition of two internal states [math] and $1$ are entangled via atomic interactions after a state-dependent transport (see Fig.1). As it was shown in Refs.KPS2013 and Byrnes2013 , this scheme allows to generate a variety of entangled states by unitary evolution with the Hamiltonian $\chi_{ab}\hat{S}_{z}^{a}\hat{S}_{z}^{b}$, where $\hat{S}_{z}^{a(b)}$ are the $z$ components of the collective spins of the two Bose-Einstein condensates $a$ and $b$. At remarquable times of the order of $\pi/\chi_{ab}$, one observes sharp dips in the entanglement entropy, associated with the formation of non-local macroscopic superpositions of coherent spin states which evoke the two-qubit Bell-state but in a $N$-body system. This Hamiltonian was recently realized with microwave cavities bridged by a superconducting artificial atom and used to generate entangled Schrödinger cat states of the electromagnetic field Schoelkopf2016 . The EPR correlations which are the focus of the present article appear on the other hand at much shorter times $1/N\ll\chi_{ab}t\ll 1/\sqrt{N}$ KPS2013 . While no striking features appear in the entanglement entropy and phase space representation of the state at these times, the non-local correlations are large enough to violate an EPR steering inequality Doherty2007 and the state, which can be seen as a non-local equivalent of a spin-squeezed state, is much more robust against decoherence. Using atomic interactions, spin-squeezed states have been already obtained in clouds whose size can reach a few thousand atoms KetterleSqueezing2007 ; Oberthaler2008 ; Treutlein2010 , and experiments along the same lines to demonstrate non-local entanglement could be envisaged.

However, fast transport of atomic clouds is likely to excite the spatial dynamics of the gas and one should check whether any non-local entanglement remains visible under these conditions. This is what we investigate in this paper by performing realistic numerical simulations of the intertwined spatial and entanglement dynamics. Our approach, adapted from Refs.CastinSinatra2000 ; TreutleinSinatra2009 , describes the state of the system as a coherent superposition of states with a fixed number of particles in a set of distinguishable modes (four in our case) representing the condensates. In each configuration, the wave functions evolve according to a time-dependent variational principle. This formalism is able to describe both the out-of-equilibrium dynamics of the gas, even with important excitations of the spatial eigenmodes, and the dynamics of entanglement between the distinguishable modes. Using this powerful numerical approach, we investigate the effect of spatial dynamics on an EPR entanglement witness and provide experimental parameters for which a significant violation of an EPR inequality could be observed.

I The dynamical model

I.1 Time-dependent description of a system of multimode bosons

Multimodal decomposition

We consider a gas of bosonic atoms distributed in a number of modes of order unity. In practice, we will apply our model to the EPR situation of Ref.KPS2013 using four modes. We imagine a situation in which the atoms in different modes never have the chance to exchange and can be considered as belonging to distinguishable components, so that we can expand the bosonic field operator $\hat{\psi}$ in the form:

[TABLE]

where $\hat{\psi}_{\alpha}$ annihilates a boson in component $\alpha$ at site r, and $A$ is the set of distinguishable components. We work with a space discretized into a cubic lattice of step $l$. This prefigures the numerical implementation of the model and removes any short-distance divergence. The commutation relations of the field operator are given by:

[TABLE]

We assume that the interactions between the atoms take place in the low-energy dilute regime usual in cold atoms, so that they can be represented by contact interactions. We note $g_{\alpha\alpha^{\prime}}$ the coupling constant between one atom of the component $\alpha$ and one of $\alpha^{\prime}$. The Hamiltonian of this system reads:

[TABLE]

with the convention $g_{\alpha^{\prime}\alpha}=g_{\alpha\alpha^{\prime}}$. Note that we allow for a dependence of the external trapping potential $V_{\alpha}$ not only on time but also on the mode $\alpha$. This can be realized experimentally with state-dependent potentials Treutlein2010 .

Variational principle

To describe the many-body dynamics arising from the Hamiltonian (3) we rely on a classical field theory based on the variational principle Ripka1985 . We assume that for each distinguishable bosonic field $\hat{\psi}_{\alpha}$, there is one macroscopically populated single particle mode $\phi_{\alpha}$ to be determined. This wave function depends on position r, time $t$ and number of atoms $\vec{N}=(N_{\alpha^{\prime}})_{\alpha^{\prime}\in A}$ in each mode

[TABLE]

We chose $\phi_{\alpha}$ to be normalized to unity:

[TABLE]

Next, we introduce the operator which annihilates a boson in the wave function $\phi_{\alpha}$:

[TABLE]

This is a bosonic operator since $[\hat{a}_{\alpha}(\vec{N},t),\hat{a}_{\alpha^{\prime}}^{\dagger}(\vec{N},t)]=\delta_{\alpha,\alpha^{\prime}}$. We use it to form the “Fock state” with $N_{\alpha}$ bosons in the wave function $\phi_{\alpha}$:

[TABLE]

where $|\rm vac\rangle$ is the vacuum of bosons, and we use the short-hand notation $\sum_{\alpha}=\sum_{\alpha\in A}$, $\prod_{\alpha}=\prod_{\alpha\in A}$ and $\{X_{\alpha}\}=\{X_{\alpha}\}_{\alpha\in A}$. Note that on state (7), the field operator $\hat{\psi}_{\alpha}(\textbf{r})$ acts as $\phi_{\alpha}(\vec{N},\textbf{r},t)\hat{a}_{\alpha}(\vec{N},t)$. This can be seen by completing $\phi_{\alpha}$ into an orthonormal basis of the Hilbert space of mode $\alpha$ and by expanding the field operator over this basis.

We are now ready to apply the variational principle. The action between the initial time $t_{\rm i}$ and the final time $t_{\rm f}$ reads:

[TABLE]

The first term between curly brackets in (8) can be expressed as

[TABLE]

and the energy $E(\vec{N},t)$ is given by

[TABLE]

where we have defined

[TABLE]

The extremalisation of the action, $\forall\textbf{r},\forall t,\,\delta S/\delta\phi^{*}_{\alpha}(\vec{N},\textbf{r},t)=0$, combined with Eqs.(9) and (10) leads to the equation of evolution of $\phi_{\alpha}$, also known as the “time-dependent Gross-Pitaevskii equation”:

[TABLE]

Accumulated phase of a Fock state

We now want to describe the evolution of the quantum state describing the gas. We imagine that initially this state coincides with a Fock state:

[TABLE]

Without approximations, it is difficult to follow the evolution of this state because the wave functions $\phi_{\alpha}$ are not the exact eigenmodes of the Hamiltonian (3), therefore the state after an evolution time $t$ is no longer a Fock state. We overcome this difficulty by assuming that the coupling to the modes orthogonal to $\phi_{\alpha}$ remains weak all along the evolution. This is justified as long as the dynamics does not create a significant depletion in the condensates. Under this assumption, the state at time $t$ can be approximated by an Hartree-Fock Ansatz:

[TABLE]

where the accumulated phase $A(\vec{N},t)$ is self-consistently defined by

[TABLE]

Deriving this relation with respect to time we obtain its equation of motion

[TABLE]

which completes our description of the dynamics of the gas.

I.2 EPR experiment with two bimodal condensates

Coherent superposition of Fock states

We now place ourselves in the experimental situation of Ref.KPS2013 (represented in Fig.1) where two initially uncorrelated Bose-Einstein condensates of atoms in their internal ground state denoted [math] are prepared in a double well potential in an initial wave function $\phi_{a0}(\textbf{r},t=0)$ for the $a$ (or left) well and $\phi_{b0}(\textbf{r},t=0)$ for the $b$ (or right) well. In this study, we assume that the number of particles is fixed to $N_{a}$ in the $a$ well and to $N_{b}$ in the $b$ well. In each well, an electromagnetic pulse brings the atoms in a superposition of two internal states noted ${0}$ and ${1}$. We assume that the four distinguishable components thus formed

[TABLE]

will remain so during the entire evolution. This is the case in the scheme proposed in Fig.1. Just after the pulses which prepares a coherent superposition of internal states [math] and $1$ in each well, the state is $a$-$b$-factorized and can be written in the compact form:

[TABLE]

where the coefficients $C_{\alpha}$ depend on the strength and on the phase of the pulses (for two $\pi/2$-pulses, $\forall\alpha,\,C_{\alpha}=1/\sqrt{2}$) and the atom number conservation in each well constrains the vector $\vec{N}=(N_{\alpha})_{\alpha\in A}$ to belong to the subset of $\mathbb{N}^{4}$:

[TABLE]

The operators $\hat{a}_{\sigma,0}^{\dagger}(0)$ and $\hat{a}_{\sigma,1}^{\dagger}(0)$ create atoms in the same wave functions $\phi_{\sigma 0}(\textbf{r},t=0)=\phi_{\sigma 1}(\textbf{r},t=0)$, which do not depend on the population $\vec{N}\in\mathcal{N}$, hence the omission of this variable in the list of arguments. The decomposition (19) suggests a two-fold naive procedure to compute the state $|\Psi(t)\rangle$ of the gas at time $t$:

we compute the evolved wave function $\phi_{\alpha}(\vec{N},\textbf{r},t)$ for each mode $\alpha$ and for all vectors $\vec{N}$ such that the prefactor $\left({N_{a}!N_{b}!}/{\prod_{\alpha}N_{\alpha}!}\right)^{1/2}\prod_{\alpha}C_{\alpha}^{N_{\alpha}}$ is significantly non-zero. In the case of $\pi/2$-pulses this prefactor is maximal for $\vec{\bar{N}}=(N_{a}/2,N_{a}/2,N_{b}/2,N_{b}/2)$ and decreases around this maximum with a typical width $N_{a,0}-N_{a}/2=O(N_{a}^{1/2})$ and $N_{b,0}-N_{b}/2=O(N_{b}^{1/2})$. With this procedure, a typical number of order $O(N_{a}^{1/2}N_{b}^{1/2})$ of wave functions must be computed and evolved. 2. 2.

we evolve the phase factor $A(\vec{N},t)$ using Eq.(16), again for all the values of $\vec{N}$ which contribute significantly to the superposition (19).

Doing so, we obtain the state

[TABLE]

Although analytically attractive, this naive procedure seems difficult to implement numerically for large atom numbers for two reasons: $(i)$ it requires us to compute and store a senselessly large number of wave functions as soon as $N_{a},N_{b}$ reach mesoscopic values and $(ii)$ it supposes that we calculate the phase factor $A$, which is of order $O(N)$ (see Eq.(16)), with an accuracy of order $O(1/N)$ as it will be clear in section I.3. In that section, we will explain how to overcome these difficulties with the help of the modulus-phase approximation. Before we do so, let us explain in the following paragraph how a generalized EPR experiment can be performed in our system of multimode bosons.

Collective spins and entanglement witness

To perform a generalized EPR experiment on our many-body system, we need to define generalized conjugated variables in each subsystem $a$ and $b$. To this end, we view each two-level atom as an effective spin and define the collective spins

[TABLE]

which form a spin algebra in each well $\sigma=a,b$. The average spins precess around the $z$ direction at an angular frequency given by the chemical potential difference between states [math] and $1$. We get rid of this effect by unrotating the spin in the $x,y$ plane

[TABLE]

where the rotation angle is defined from the average spin

[TABLE]

The interesting entanglement dynamics takes place in the $y\phi,z$ plane. Because of the non-linearity of the interaction term of the Hamiltonian (3), non-classical correlations build up between the spin components $\hat{S}_{y\phi}^{a}$, $\hat{S}_{z}^{a}$, $\hat{S}_{y\phi}^{b}$, and $\hat{S}_{z}^{b}$. The correlations within one well are responsible for spin-squeezing ; they are not our focus in this article. We are looking for the highest non local correlations between components of the $a$ and $b$ wells. These are obtained for some well chosen quadratures in each well

[TABLE]

Our EPR entanglement witness Drummond2012 ; KPS2013 is constructed from these optimally correlated quadratures and their conjugated ones

[TABLE]

When this quantity is below unity, it shows a violation of the Heisenberg relation which Bob should expect if his quantum state were local, namely:

[TABLE]

where $\Delta^{2}\hat{S}$ is the variance of $\hat{S}$ and $\hat{S}_{\beta}^{b,{\rm inf}}$ is the best guess that Alice can do on the outcome of the measurement of $\hat{S}_{\beta}^{b}$, knowing the outcome of her measurement of $\hat{S}_{\alpha}^{a}$:

[TABLE]

This is the so-called “steering” effect: the violation of the Heisenberg inequality shows that Bob’s measurement acts instantaneously and at a distance on the quantum state of Alice. If one obtains a value of $E_{\rm EPR}$ below unity, one not only proves the entanglement between the $a$ and $b$ subsystems but also rules out the possibility of reconciling Quantum Mechanics (namely the Heisenberg uncertainty principle) with local realism Doherty2007 ; KPS2013 . However, $E_{\rm EPR}<1$ does not prove that no hidden variable theory can account for the results of the experiment, and is therefore strictly weaker than the violation of a Bell inequality.

I.3 Modulus-phase approximation

We now explain how to overcome the limitation mentioned in §I.2.$a$ by performing a linearization of the phase of the wave functions around the central Fock state.

Modulus-phase decomposition of the wave functions

The modulus-phase approximation consists in a linearization of $\phi_{\alpha}(\vec{N})$ around its central and average value for $\vec{N}=\vec{\bar{N}}$, where $\vec{\bar{N}}$ is the configuration that has the greatest weight in the superposition (19). This approximation further assumes that the variation of the modulus of $\phi_{\alpha}$ can be entirely neglected CastinSinatra2000 ; TreutleinSinatra2009 . It is valid when $N_{a},N_{b}$ are large enough such that the typical width of the distribution of Fock states in state (19) (which is of order $O(\sqrt{N})$) is small compared to the central values $\bar{N}_{a,0},\bar{N}_{b,0}$. We write the wave function $\phi_{\alpha}$ for $\vec{N}$ such that $N_{\alpha}-\bar{N}_{\alpha}=O(N^{1/2})$ in the form:

[TABLE]

where

[TABLE]

is the central wave function of the mode $\alpha$, and the vector $\vec{\nabla}$ gathers the derivatives with respect to the number of atoms in each mode:

[TABLE]

Thanks to this approximation, we need only the central wave functions $\bar{\phi}_{\alpha}$ and the derivatives $\vec{\nabla}\theta_{\alpha}$ of the phases in $\vec{N}=\vec{\bar{N}}$ to describe the state $|\Psi(t)\rangle$ of the gas. In practice, we compute $\vec{\nabla}\theta_{\alpha}$ using, besides the central wave function $\bar{\phi}_{\alpha}$, the slightly displaced wave functions $\phi_{\alpha}(\vec{\bar{N}}+\vec{\beta})$ where $\vec{\beta}=(\pm\beta_{\alpha})_{\alpha\in A}$ and $\beta_{\alpha}\ll\bar{N}_{\alpha}$. In the $4$-mode situation considered in (17), this means evolving a total of only nine different Fock states, hence storing only $36$ wave functions.

Simplification of the accumulated phase

The objective of this section is to show how the computation of the accumulated phase $A$ can be avoided. In practice, rather than with the factor $A$ directly, we deal with differences of the kind $A(\vec{N}+\vec{\beta})-A(\vec{N})$, where $\vec{\beta}$ represents a displacement of order unity of the population of each mode. We intend to show that this factor can be reabsorbed by a product of wave function overlaps. The quantity which appears naturally when one computes the averages of the spin components (see below Sec. I.4 and in particular Eqs.(40) and (46)) is the following combination:

[TABLE]

In contrast with $A$, the phase $\Theta$ has a small time derivative and is therefore easy to compute numerically. We give the real part of this time-derivative in the modulus-phase approximation:

[TABLE]

To obtain this result, we use Eq.(16) to express the time-derivative of the first term between parenthesis in Eq.(36), and the Gross-Pitaevskii equation (12) to express the time-derivative of the second term. Note that, consistently with the fact that we neglected the variations of the modulus of $\phi_{\alpha}$, we neglect the imaginary part of $\mathrm{d}\Theta/\mathrm{d}t$ which (like its real part) is of order $O(1/N)$.

I.4 Quantum averages of field operators

We conclude this section by deriving a general formula for the two-position densities of the quantum averages. Let us write an arbitrary products of field operators taken at different sites r et $\textbf{r}^{\prime}$ as:

[TABLE]

where the vectors $\vec{\gamma},\vec{\gamma^{\prime}},\vec{\delta},\vec{\delta^{\prime}}\in\mathbb{N}^{4}$ are of order unity. For the sake of readability, we shall use the short-hand notations

[TABLE]

With the knowledge of $\langle\hat{W}(\textbf{r},\textbf{r}^{\prime})\rangle$ for all $\vec{\gamma},\vec{\gamma^{\prime}},\vec{\delta},\vec{\delta^{\prime}}$ in hand, all the spin averages needed to compute $E_{\rm EPR}$ in Eq.(30) can be reconstructed by integrating over r and $\textbf{r}^{\prime}$ and doing simple linear combinations. We compute the average value of $\hat{W}(\textbf{r},\textbf{r}^{\prime})$ using the expansion (21) of $|\Psi(t)\rangle$ over the Fock states:

[TABLE]

The scalar product appearing on the last line of (40) is non-zero if and only if $\vec{N}^{\prime}=\vec{N}+\vec{\gamma}^{+}-\vec{\delta}^{+}$ ; it can be computed without difficulty:

[TABLE]

We now use the modulus-phase approximation in order to $(i)$ compute the product of wave functions appearing on the second line of (40) and to $(ii)$ eliminate the phase factor $A(\vec{N}+\vec{\gamma}^{+}-\vec{\delta}^{+})-A(\vec{N})$ in favor of the quantity $\Theta(\vec{N},\vec{\gamma}^{+}-\vec{\delta}^{+},t)$ defined in Eq.(36). We obtain the slightly long but very general formula

[TABLE]

where $\vec{\theta}=(\theta_{\alpha})_{\alpha\in A}$. In the numerical simulation, all the relevant average values can be calculated from this formula after a double summation over the spatial indices r and $\textbf{r}^{\prime}$. This result generalizes and gives a compact form to the equations of appendix B of Ref.TreutleinSinatra2009 .

II Numerical simulation

II.1 Experimental and numerical sequence

The experiment we simulated is inspired by the situation available on an atom chip Treutlein2010 , and is depicted on Fig.1. It is realized by the following sequence:

$(i)$ We populate a double well potential with $N_{a}$ atoms of Rubidium $87$ in the hyperfine state $|F=1,m_{\rm F}=-1\rangle$ in the many-body ground state of the $a$ well and similarly in the $b$ well with $N_{b}$ atoms. Different simulations were performed varying the initial trap distance and the atom number in the range $N_{a}=N_{b}=100$ to $N_{a}=N_{b}=5000$. Each well is composed of an isotropic harmonic trap whose trapping frequency is fixed to

[TABLE]

The system is cylindrically symmetric around the axis of the two trap centers. We choose this axis as the $z$ axis of a cylindrical frame and write the spatial dependence of the wave functions as $\phi(\textbf{r})=\phi(r,z)$.

$(ii)$ A $\pi/2$-pulse drives the atoms in both wells into a superposition of the two long-lived hyperfine states $|0\rangle=|F=1,m_{F}=-1\rangle$ and $|1\rangle=|F=2,m_{F}=1\rangle$. The scattering lengths characterizing the interactions between these hyperfine states are known experimentally within an accuracy of about $0.5\%$ Hall2007 ; Sidorov2013 , and depend on the external magnetic field Verhaar2002 . The results presented in this article are for the values

[TABLE]

in units of the Bohr radius $R_{\rm B}\simeq 52.918\times 10^{-12}\,\mbox{m}$ 111We checked that a relative change of the order of 0.5% of the scattering lengths does not significantly affect the results in our scheme by performing test simulations for the sets of values $a_{00}=100.0\,R_{\rm B}\quad a_{11}=95.0\,R_{\rm B}\quad a_{01}=98.0\,R_{\rm B}$ and $a_{00}=100.44\,R_{\rm B}\quad a_{11}=95.47\,R_{\rm B}\quad a_{01}=98.28\,R_{\rm B}$ for $N_{a}=N_{b}=5000$. .

$(iii)$ The trapping potential of atoms in state $|0\rangle$ is displaced by a distance exactly equal to the initial separation between the two wells using a ramp of fixed duration $\omega t_{\rm R}=10$. In order to minimize the creation of spatial excitations, the ramp should be slow in the beginning and in the end, when the atomic clouds are being mixed or demixed, and faster in between. This is achieved by varying the displacement $\delta z$ of the trapping potential of state $1$ using an hyperbolic tangent function of time

[TABLE]

where $\delta z_{\rm max}$ is the initial separation between the two wells and $t_{\rm R}$ is the duration of the ramp.

$(iv)$ After the end of the ramp, component 0 of the $a$ well interacts with component 1 of the $b$ well during an adjustable interaction time $t_{\rm int}$.

$(v)$ Components 0 are ramped back to their initial positions using the ramp (45) backwards.

$(vi)$ Averages and correlations between spin components are computed from Eq.(42).

$(vii)$ The minimal value of $E_{\rm EPR}$ is found by numerically optimizing over the angles $\alpha$ and $\beta$ (see Eqs.(28) and (29)).

The main free parameters of this experiment are the initial spatial separation $\delta z_{\rm max}$ between the two wells and the ramping time $t_{\rm R}$. The separation $\delta z_{\rm max}$ should be large enough to prevent any overlap between the two wells, initially and after ramping back, thus ensuring the spatial separation required for the EPR experiment. The ramping time $t_{\rm R}$ should be long to avoid too violent excitations of the spatial dynamics of the clouds, which would degrade the creation of EPR correlations. There is also an upper bound on $t_{\rm R}$: during the ramp, correlations build up between intrawell spin components (e.g. between $\hat{S}_{y\phi}^{a}$ and $\hat{S}_{z}^{a}$) like in a spin-squeezing experiment TreutleinSinatra2009 ; these correlations are unfavorable for our purpose because they are a source of uncertainty for non local observables such as $\hat{S}_{y\phi}^{a}\hat{S}_{z}^{b}$.

II.2 Discussion of the results

As a first example, in Fig.2 we test the effect of spatial dynamics on EPR correlations by considering a relatively fast transport of the atomic clouds with $\delta z_{\rm max}=10a_{0}$ in unit of the oscillatory length $a_{0}=\sqrt{\hbar/m\omega}$ and the ramping time $\omega t_{\rm R}=10$, and two atom numbers $N_{a}=N_{b}=100$ and $N_{a}=N_{b}=500$. The results of the numerical simulation for $E_{\rm EPR}$ and for the renormalized spin length $\langle\hat{S}_{x\phi}^{\sigma}\rangle$ including spatial dynamics in the modulus-phase approximation are given by symbols respectively in the top and bottom row of Figs.2 and 3. In addition, we show the normalized density overlap

[TABLE]

between wavefunctions $\phi_{\sigma 0}$ and $\phi_{\sigma 1}$, as dash-dotted lines in the bottom row. Each point is obtained by a numerical simulation corresponding to a different interaction time, always including the initial phase in which the potentials are displaced and the final one in which they are ramped back to the initial position before the measurement. We compare these results to an analytical 4-mode model represented in solid-line ; this model assumes that the whole evolution (the $\pi/2$-pulse and the onward and backward ramps) is adiabatic so that the 4 wave functions remain in their instantaneous ground state at all time. With this assumption, the system can be mapped onto a 4-mode model whose non-linearity parameters $\chi_{a}$, $\chi_{b}$ and $\chi_{ab}$, coefficients respectively of $(\hat{S}_{z}^{a})^{2}$, $(\hat{S}_{z}^{b})^{2}$ and $-\hat{S}_{z}^{a}\hat{S}_{z}^{b}$ in the 4-mode Hamiltonian KPS2013 , are time-dependent and given by

[TABLE]

where $\mu_{\alpha}$ is the ground state chemical potential of mode $\alpha$. In practice, for the values (44) of the scattering lengths, we have $\chi_{a}\simeq\chi_{b}\equiv\chi$. In the formulas of Appendix B of Ref.KPS2013 , valid for a purely stationary 4-mode model, we perform the substitutions $\chi t\to\int_{0}^{t}\chi(t^{\prime})\mathrm{d}t^{\prime}$ and $\chi_{ab}t\to\int_{0}^{t}\chi_{ab}(t^{\prime})\mathrm{d}t^{\prime}$, where $t=2t_{\rm R}+t_{\rm int}$ is the total evolution time. The optimized EPR entanglement witness obtained through this approach is shown on the top row of Fig.2 in solid line.

The analytical $4$-mode model captures well the general trend of variation of $E_{\rm EPR}$ and of the renormalized spin length $\langle\hat{S}_{x\phi}^{\sigma}\rangle$. $E_{\rm EPR}$ reaches a minimum for an optimal time KPS2013 , which is the result of a competition between the simultaneous creation of local and non-local correlations and the loss of coherence. Notice that this minimum is lower for 500 atoms ($E_{\rm EPR}\simeq 0.65$) than for 100 ($E_{\rm EPR}\simeq 0.85$) (see also KPS2013 ), indicating that the EPR inequality is violated more strongly for larger atom number.

On top of this trend, a strong spatial dynamics is visible: at low atom number, the dynamics is dominated by regular oscillations of the centers-of-mass of the wave functions in the first excited mode of the harmonic trap. This is mostly an ideal gas behavior, caused by the finite ramping time; the weak interactions manifest themselves in the synchronization of the oscillations of $\phi_{a0}$ and $\phi_{b1}$, caused by the weak coupling between them (see Supplementary Material). This dynamics is responsible for the regular drops in the spin length $\langle\hat{S}_{x\phi}^{\sigma}\rangle$ visible on Fig.2.$(a)$ with a frequency that matches the trap frequency. Remark that the spin length is close to its maximal value only when the density overlap is close to one. However, the converse is not true because perfectly overlapping wavefunctions may still be dephased. Ultimately, the dips of the spin length causes the regular peaks in the EPR entanglement witness. At large atom number $N_{a}=N_{b}=500$ (Fig.2.$(b)$) the spatial dynamics is both damped and rendered more chaotic by the stronger interactions, and the generation of non-local correlations is faster.

We found that the results are improved by reducing the initial distance between the two traps. This increases the ratio between the interaction time and the total time and minimizes the excitation of spatial dynamics. We show an example for $\delta z_{\rm max}=6a_{0}$ and $N_{a}=N_{b}=500$ in Fig.3 where the minimum value of $E_{\rm EPR}\simeq 0.54$ is reached at $t\simeq 0.17$ s. The corresponding plot for $\langle\hat{S}_{x\phi}^{\sigma}\rangle$ is shown in the lower panel of Fig.3. A slight improvement ($E_{\rm EPR}\simeq 0.45$ at $t\simeq 0.16\,$s) is obtained when increasing the number of atoms to $N_{a}=N_{b}=5000$ but the spatial dynamics is more strongly excited.

To complete our investigation, we have evaluated the total loss rate in the interaction configuration (middle row in Fig.1) using stationary wave functions for $N_{a}=N_{b}=500$, scattering lengths and trap frequencies as in Fig.2, and the loss rate constants SidorovDrummond2012 $\kappa_{11}=81\times 10^{-21}\,\mbox{m}^{3}\ \mbox{s}^{-1}$, $\kappa_{01}=15\times 10^{-21}\,\mbox{m}^{3}\ \mbox{s}^{-1}$ for two-body losses and $\kappa_{000}=5.4\times 10^{-42}\,\mbox{m}^{6}\ \mbox{s}^{-1}$ for three-body losses. We find that about $7$ particles are lost on average in $0.2$ s due to two-body losses, while three-body losses are negligible. If we add one-body losses corresponding to a lifetime of $1$ min, we obtain a total lost fraction $N_{\rm lost}/(N_{a}+N_{b})$ of about $10^{-2}$ at $0.2$ s. In analogy with what happens for spin-squeezing, where the losses are negligible as long as the lost fraction is much smaller than the squeezing parameter $\xi^{2}$ (see Eq.(21) in Ref.YunSinatra2008 ), we expect that the effect of losses should not significantly affect the predicted value of $E_{\rm EPR}$ for the considered parameters. We come to a similar conclusion for larger atom numbers $N_{a}=N_{b}=5000$ where we estimate a lost fraction of $2\times 10^{-2}$ at the optimal time for EPR steering $t\simeq 0.16\,$s.

Our results show that an EPR steering experiment can be performed using interactions in Bose-Einstein condensates in state-dependent potential even accounting for the spatial dynamics. This completes our previous study KPS2013 which was limited to stationary modes.

Conclusion

We studied a system of Bose-Einstein Condensates in 4 distinguishable components entangled by collective transport in state-dependent potentials. Far beyond a simple 4-mode approach, we worked out a realistic description of the spatial structure of the system where both the spatial excitations and the entanglement between the distinguishable components are accounted for. From the methodological point of view, we generalized the two-component spatial dynamical model of Refs.CastinSinatra2000 ; TreutleinSinatra2009 to a many-component system, and we worked out a general formula to compute all the relevant correlations. Our compact formulation makes our theory an easily adaptable tool to study the effect of spatial dynamics on entanglement, and provides the experimentalists with a numerical recipe to benchmark their results.

Based on this description, we performed a set of numerical simulations to investigate the effects of spatial dynamics, and in particular of the excitations created by the collective transport of atomic clouds, on the non-local correlations which we expect in our system. We observed that strong excitations of the spatial modes were created, whose behavior became more chaotic as the number of atoms increased. However, despite these excitations, the non-local correlations remain strong enough to perform an EPR steering experiment. If the parameters, and in particular of the collective transport, are properly adjusted, our EPR entanglement witness can reach below the critical value of $1$, and down to $E_{\rm EPR}=0.54$ for $500$ atoms in each well and $0.45$ for $5000$. The lowest values of our entanglement witness were reached at relatively short times, of the order of $0.2$  s. This is longer than the times ($\approx 15$  ms) at which spin squeezing was observed with Bose Einstein condensates Treutlein2010 but remains within reach of cold atomic experiments Dalibard2002 .

Acknowledgments

H. Kurkjian is a FWO [PEGASUS]2 Marie Skłodowska-Curie fellow. This project has received funding from the FWO and the European Union’s Horizon 2020 research and innovation program under the Marie Skłodowska-Curie grant agreement number 665501. K. Pawłowski acknowledges support by the (Polish) National Science Center Grant No. 2014/13/D/ST2/01883.