Krzysztof Pawlowski

Matteo Fadel

Philipp Treutlein

Yvan Castin

Alice Sinatra

K Pawłowski

Mesoscopic quantum superpositions in bimodal Bose-Einstein condensates: decoherence and strategies to counteract it

Keywords: numbers: 03.75.Gg, 03.75.Mn, 42.50.Dv

come

I. INTRODUCTION

While mesoscopic superpositions of coherent states of light with up to one hundred photons [1][2][3], and Greenberger-Horne-Zeilinger states with up to 14 trapped ions [4,5] have been generated and observed in experiments, Schrödinger cat states with atomic gases are still out of reach [6]. Bose-Einstein condensates of ultracold atoms, confined in conservative potentials made by light or magnetic fields, are excellent candidates to take up the challenge as they are to a good approximation isolated systems. Characterized by a macroscopic population of a single particle state, condensates offer in principle the unprecedented possibility of generating large orbito-spinorial Schrödiger cats, that is superpositions of two coherent spin states with opposite phases, each in a single and well controlled quantum state concerning the orbital degrees of freedom. Nevertheless, decoherence originating from particle losses and total atom number fluctuations, as well as from the intrinsic multimode nature of the atomic field and nonzero initial temperature, is usually not negligible. The aim of this work is to present strategies to counteract decoherence, within the possibilities and constraints of specific experiments on bimodal condensates. Other proposals starting with monomode condensates, see e.g. [7], are not discussed here.

In analogy to the well-known optical proposal of Yurke and Stoler of 1986 in reference [8], in bimodal condensates the entanglement stems from the interactions between atoms that introduce a nonlinearity of the Kerr type for the atomic field. A state where all the atoms are in a superposition of the two modes with a well defined relative phase, a so-called "phase state" or "coherent spin state", dynamically evolves into a Schrödinger cat state, superposition of two phase states with opposite relative phases [9,[START_REF] Castin | Bose-Einstein Condensates in Atomic Gases[END_REF]. At twice the cat-state time the system returns into a single phase state giving rise to a revival peak in the contrast of the interference pattern between the two modes [START_REF] Castin | [END_REF].

The influence of particle losses on the revival peak amplitude has been studied analytically in reference [12]. In section II of the present paper we show that there is a simple quantitative relation between the amplitude of the revival peak and the cat-state fidelity. As an application, for N = 300 rubidium 87 atoms in two separated spatial modes, with three-body losses but no fluctuations of the total number of particles and at zero temperature, we calculate that a Schrödinger cat is obtained in time t cat = 128 ms with a fidelity F ≃ 0.8.

In section III we concentrate on the use of two internal states of a hyperfine transition in sodium or rubidium atoms. For Rb we consider the states |F = 1, m F = -1 and |F = 2, m F = 1 that have been used to generate spin squeezing in state-dependent potentials on a chip [13,14]. A particular interest in these states resides in the fact that they form the clock transition in atomic clock experiments with trapped atoms on a chip [15]. We present a strategy to obtain mesoscopic superpositions using these systems, despite severe intrinsic and experimental constraints, including particle losses and Poissonian fluctuations of the total particle number. We first perform a numerical study where we optimize the Fisher information of the obtained state by exploring systematically the parameter space for experimentally accessible configurations. In contrast to reference [16], where the Husimi function of the macroscopic superposition was considered (this distribution does not exhibit fringes), we look at the interference fringes of the Wigner function to quantify the survival of quantum correlations in the presence of decoherence (see endnote [17]). We also calculate the Fisher information of the state after averaging over many stochastic realizations in order to quantify its usefulness for metrology. The numerical study is followed by an analytical part that gives a limpid interpretation of the results, see section IV.

Finally in section V we give up the two-mode approximation for a truly multimode description of the bosonic system and we estimate what are the constraints on the temperature of the Bose-condensed gas used in the preparation of the initial phase state, in order to obtain the desired mesoscopic superposition with a good fidelity and a significant revival in the phase contrast. We conclude in section VI.

II. CAT-STATE FIDELITY VERSUS CONTRAST REVIVAL

In this section we show that there is a simple relation between the fidelity of the state obtained at t cat and the amplitude of the contrast revival peak at t rev = 2t cat . For simplicity, we consider in this section two spatially separated components, with the same scattering length and loss rates, as one would have by using two symmetric Zeeman sub-levels as internal states, or by using two spatially separated BECs in the same internal state.

We neglect fluctuations of the total particle number assuming that an initial state with a fixed number of particles can be prepared, for example by melting a Mott insulator phase in an optical lattice, or by non-destructive detection of the atoms with an optical cavity.

Since the bosonic field populates two orthogonal modes with corresponding annihilation operators â and b, one can attribute an effective spin 1/2 to the bosons and introduce the usual single-spin Bloch representation and the usual dimensionless collective spin operators [6]:

Ŝx = â †b + b † â 2 ; Ŝy = â †b -b † â 2i ; Ŝz = â † â -b †b 2 . (1)
We consider an initial phase state with N particles, on the equator of the Bloch sphere, with a relative phase ϕ = 0 between the two modes:

|ψ(0) = 1 √ N ! â † + b † √ 2 N |0 ≡ π 2 ; 0 N . (2)
Here the phase state with N atoms is defined as

|θ; ϕ N ≡ 1 √ N ! cos θ 2 e i ϕ 2 â † + sin θ 2 e -i ϕ 2 b † N |0 .
(3) The relative phase ϕ ∈ [-π, π] has a meaning modulo 2π and the polar angle θ ∈ [0, π]. The initial state (2) evolves under the influence of a nonlinear spin Hamiltonian resulting from the elastic s-wave interactions inside each mode [START_REF] Castin | [END_REF]12],

H = χ Ŝ2 z = χ 2 N 2 a + N 2 b - N 2 2 , (4)
and in the presence of particle losses (one-, two-and three-body) within each spatial component. The whole evolution is governed by the master equation for the density operator [12,18]:

d dt ρ = 1 i Ĥ, ρ + L 1 [ρ] + L 2 [ρ] + L 3 [ρ] , (5)
where the Liouvilian operators are L m = L

and similarly for the mode b. Note that there are no collisions between modes a and b because they are spatially separated. The rates γ (m) are related to the loss rate constants K m and to the (in practice Gross-Pitaevskii) normalised condensate wavefunction φ(r) in one of the modes by mγ (m) = K m d 3 r|φ(r)| 2m [18].

The loss rate constants are such that, in the spatially homogeneous zero-temperature Bose gas with N particles and mean density ρ, the m-body losses lead to a decay

d dt N = -K m ρ m-1 N .
In the absence of losses, at the time t cat = π 2χ , the system is in a Schrödinger cat state given by

|ψ(t cat ) = e -i π 2 Ŝ2 z π 2 ; 0 N = e i π 8 (N 2 -2) e i π 2 N Ŝz π 2 ; 0 N + ie i π 2 N π 2 ; π N √ 2 . (7)
This results from the identity exp(-iπn

2 /2) = exp(iπ/4)[exp(iπn) -i]/ √ 2,
for n integer, and from the expansion of the initial state over Fock states.

|ψ(0) = 1 2 N/2 N Na=0 N ! N a !N b ! 1/2 |N a , N b (8)
with N b = N -N a . Equation ( 7) agrees with equation (19) in reference [19] up to a global phase factor but it disagrees with reference [9]. By using the relation exp(iα Ŝz ) π 2 ; ϕ N = π 2 ; ϕ + α N , it can be rewritten as

|ψ(t cat ) N even = e -iπ/4 π 2 ; 0 N + i π 2 ; π N √ 2 (9) |ψ(t cat ) N odd = e -iπ/8 π 2 ; π 2 N -π 2 ; 3π 2 N √ 2 .( 10
)
In presence of losses, we introduce the fidelity F (t) of the state ρ at time t, that is its overlap with the "target" state that one would obtain in the lossless case:

F (t) ≡ Tr{ ρ(t)|ψ 0 (t) ψ 0 (t)|} (11)
The normalized first order correlation function between the two modes gives the contrast of the interference pattern if the two modes are made to interfere:

g (1) (t) = 2 N Ŝx (t) (12)
Its maximum value is one, realized at t = 0 when the system is in a phase state. In the lossless case g (1) (t) = ±1 is recovered at multiples of the revival time t rev = 2t cat . We show here that for weak losses, one has to a very good approximation

F (t cat ) = |g (1) (t rev )| 1/2 (13)
A. Proof in the constant loss rate approximation

The Monte Carlo wave function method [20][21][22] provides us with a stochastic formulation of the master equation (5). In this point of view the density matrix is seen as a statistical mixture of pure states | ψ(t) , each of which evolves under the influence of a non-hermitian effective Hamiltonian H eff and of random quantum jumps. In terms of the jump operators Ĵm ǫ that annihilate m particles in component ǫ = a, b:

Ĵm a = γ (m) âm , Ĵm b = γ (m)b m (14)
the effective Hamiltonian takes the form

H eff = H - i 2 ǫ=a,b 3
m=1 ( Ĵ † ǫ ) m Ĵm ǫ . ( 15
)
For a Monte Carlo wave function | ψ(t) normalized to unity, quantum jumps occur with a (total) rate

ǫ=a,b 3
m=1 ψ(t)|( Ĵ † ǫ ) m Ĵm ǫ | ψ(t)
. The so-called "constant loss rate approximation", introduced in reference [12], consists of the replacement

( Ĵ † ǫ ) m Ĵm ǫ → γ (m)
N m ǫ in the effective Hamiltonian, where Nǫ = N/2 is the mean initial number of particles in each component. Under this approximation, which can be used when the mean fraction of lost particles is small, the probability that n quantum jumps have occurred at time t is given by a Poisson law with parameter n = λt where

λ = 2 3 m=1 γ (m) N 2 m . (16)
In this approximation the effective Hamiltonian indeed reduces to H eff = H -i 2 λ so that the probability that no jump occurs during a time delay

τ is ||e -iH eff τ / | ψ || 2 = e -λτ .
In expression [START_REF] Castin | [END_REF] of the fidelity, only the realizations where no particles were lost at the cat-state time contribute. In this subspace the density matrix evolves only under the influence of the effective Hamiltonian and remains in a pure state P N ρ(t)P N = e -λt |ψ 0 (t) ψ 0 (t)|, where P N projects onto the subspace with N atoms. The fidelity at the cat-state time is then

F (t) = e -λtcat . (17)
On the other hand, we have shown in reference [12] that the peak in the contrast at the revival time in the presence of losses is g (1) (t rev ) = (-1) N e -λtrev = (-1) N e -2λtcat , which concludes the proof. A more detailed analysis, beyond the constant loss rate approximation, is performed analytically and numerically in Appendix A for one-and three-body losses (see endnote [START_REF]The contribution of two-body losses, whose origin is spin changing collisions, can be made negligible by choosing extremum Zeeman sub-levels[END_REF]). We show there that, in the interesting regime in which the number of lost atoms at the revival time is smaller than one, mλt rev < 1 (the fidelity and the revival would be killed by the losses otherwise), the relative correction to (13) ||g (1)

(t rev )| -F (t rev /2) 2 | F (t rev /2) 2 ≈ 2 (mπ) 2 mλt rev × mλt rev N (18
) can be interpreted, up to a factor 2/(mπ) 2 , as the product between the number of lost atoms and the fraction of lost atoms at the revival time. It is hence ≪ 1.

B. Numerical example

We show in Fig. 1 an example where we solve numerically the master equation in the presence of three-body losses (see endnote [START_REF]Extremely long lifetimes of the order of hours might be obtained in cryogenic environment[END_REF]) and compare the evolutions of the g (1) function and of the fidelity, confirming that the relation (13) approximately holds also beyond the constant loss rate approximation. The height of the first revival peak is 0.63 and the fidelity of the cat state is 0.79. The conclusion of this section is twofold. First, losses should be limited to less than one particle on average at the cat-state time to preserve a high fidelity. Second, we have shown that there is a simple quantitative relation (13) between the amplitude of the revival peak in the contrast and the cat-state fidelity. The physical reason is that each loss event introduces a random shift of the relative phase between the modes (see section IV), corresponding to a rotation of the state around the z axis, by an angle of order χt where t is the time at which the loss occurred. As χt is of the order of π at the cat-state time or the revival time, one particle lost on average is sufficient to kill both the cat state and the phase revival.

III. REALISTIC ANALYSIS FOR RUBIDIUM OR SODIUM ATOMS ON A HYPERFINE TRANSITION

This section gives a description of the two-mode dynamics as close as possible to the experimental state of the art, including losses and particle number fluctuations. The two condensed modes correspond to two different atomic internal sub-levels, already used and coupled in cold atom experiments by a hyperfine transition. As N fluctuates, we take a different perspective on the catstate formation: the goal is no longer to prepare with highest fidelity the pure cat state ( 9) or [START_REF] Castin | Bose-Einstein Condensates in Atomic Gases[END_REF], it is rather to produce a mixed cat state with maximal usefulness for precision measurements, that is maximal Fisher information. The "catiness" of the mixed state is revealed by fringes in the Wigner distribution function. We now concentrate on the two internal states of ru-bidium 87 |F = 1, m F = -1 and |F = 2, m F = 1 that have been used to generate spin squeezing in state dependent potentials on a chip [13,14]. The experimental constraints that we consider are (i) large two-body losses in |F = 2, m F = 1 due to spin changing collisions, (ii) limited background lifetime (we take 1/K 1 = 5 s) in both states due to imperfect vacuum, (iii) fluctuations of the total number of atoms N . Concerning this last effect, we remark that, even in the absence of losses, the orientation of the cat state depends on N modulo 4. This is apparent from equation (7) where a N -dependent rotation around the z-axis acts on a state (the state between parentheses) with N -independent coefficients in the Fock basis. If N fluctuates with a standard deviation ≫ 1 as in regular experiments, the interference fringes at the cat-state time are then completely washed out when averaging over N . In these conditions one might think that it would be difficult, if not impossible, to create a cat state under the experimental constraint mentioned above. We will show that this is not the case. However, in order to counteract decoherence we will have to consider a more general situation than the one described in Section II. We will (i) de-symmetrize the initial mixture by performing a large pulse instead of a π/2-pulse, (ii) de-symmetrize the two trapping potentials, and (iii) allow for an overlap between the two spatial modes. This is schematized in Fig. 2.

After rubidium we consider the two internal states of sodium 23 |F = 1, m F = 0 and |F = 2, m F = -2 in more general, cigar-shaped or pancake-shaped statedependent potentials. In this case spin changing collisions between atoms in F = 2 are suppressed and a-b losses are negligible [28]. The losses can then be significantly lowered provided a very good vacuum is achieved, which allows us to push up further the atom number in the quantum superposition.

B. Numerical calculations

We first performed a numerical study to determine the optimal experimental conditions within the given constraints.

The system state is supposed to be initially in a statistical mixture of phase states

ρ(0) = ∞ N =0 p(N ) |θ; ϕ N N θ; ϕ| (19)
where the phase state |θ; ϕ N with N atoms is given in (3), and p(N ) is the distribution of the total number of atoms, assumed to be Poissonian of average N . The master equation obeyed by ρ(t) is still of the form of Eq.( 5), but with non-symmetric m-body loss rates γ

γ (m) ǫ = K (m) ǫ m d 3 r |φ ǫ (r)| 2m , (20)
γ ab = K ab 2 d 3 r |φ a (r)| 2 |φ b (r)| 2 , (21)
where

K (m) ǫ
and K ab are loss rate constants, and γ (m) ǫ and γ ab are calculated using the stationary normalized condensate wave functions φ ǫ (r) for N a = Na and N b = Nb . As now the modes can spatially overlap, we also include two-body processes, with rate γ ab , where one atom in a and one atom in b are lost at the same time [START_REF] Egorov | [END_REF] (see endnote [29]).

The unitary part of the evolution in the master equation is calculated with the zero-temperature mean-field model Hamiltonian

H GP = ∞ Na,N b =0 E GP (N a , N b ) |N a , N b N a , N b | (22)
where E GP is the Gross-Pitaevskii energy

E GP (N a , N b ) = ǫ=a, b N ǫ φ * ǫ h ǫ φ ǫ + g ǫǫ 2 N ǫ |φ ǫ | 4 + g ab N a N b |φ a | 2 |φ b | 2 . ( 23
)
The single particle Hamiltonians h a and h b include the kinetic energy and the trapping potential. The stationary condensate wave functions φ ǫ and the Gross-Pitaevskii energy E GP have been computed numerically for different pairs (N a , N b ) (in practice a few thousands) to construct the Hamiltonian (22).

In order to find the optimal conditions, we were scanning the experimental parameters space, each time performing the evolution starting from the initial condition (19), optimizing entanglement witnesses that are sensitive to the presence of a Schrödinger cat. To avoid extreme parameters that would make the experimental realization more difficult, we have restricted the search to trapping frequencies ratios smaller than 20. Details of our procedure are given in Appendix B, and two examples of results are shown in the next subsection.

C. Fisher information and Wigner function of the cat state

For optimized conditions issued by our search algorithm (see Appendix B), in Fig. 3 and Fig. 4 we show the resulting time evolution of the Fisher information, and the Wigner distribution at the cat-state time, obtained respectively for rubidium 87 and sodium 23, for the hyperfine transitions mentioned above.

The corresponding cuts through the atomic density distribution along the z-axis of the trap for the two states are shown in Fig. 5. The Fisher information F ( Ŝz ) that we plot quantifies the sensitivity of the state to a small rotation around axis n lying in the equator of the Bloch sphere, when a measurement of the observable Ŝz is performed: Hz, ∆z = 1.620 a ⊥ (distance between the trap centers). Scattering lengths aaa = 100.4a0, a bb = 95.44a0, a ab = 98.13a0 [START_REF] Egorov | [END_REF]. One-body, two-body, and three-body loss rate constants

F ( Ŝz ) = Max n F ( Ŝz , n) , (24)
F ( Ŝz , n) = lim α→0 ∞ k=0 1 p(k|α, n) d p(k|α, n) dα 2 , (25)
K (1) a = K (1) b = 0.2 s -1 , K (2) b = 8.1 × 10 -20 m 3 /s, K ab = 1.51 × 10 -20 m 3 /s, K (3)
a = 6 × 10 -42 m 6 /s [START_REF] Egorov | [END_REF]. In the general two-mode model of section IV A, these parameters lead to χ = 12.895 s -1 , χ = 12.888 s -1 , γ

(2) b = 6.436 × 10 -3 s -1 , γ ab = 1.032 × 10 -3 s -1 , γ (3)
a = 5.15 × 10 -6 s -1 .
where p(k|α, n) is the probability of finding k atoms in the minority component a in the rotated state ρα = e -iα Ŝ•n ρe iα Ŝ•n . We have chosen Ŝz as the observable with respect to which we define the Fisher information, because one can show that in the ideal lossless case, even for a large pulse as in Fig. 2, F ( Ŝz ) reaches the quantum Fisher information obtained by maximizing F ( Ô, m) = 0.01 s -1 . In the general two-mode model of section IV A, these parameters lead to χ = 8.763 s -1 , χ = 8.729 s -1 . We expect no relevant a-b or b-b two-body losses here [28], and we have checked that for the considered parameters the contribution of three-body losses is negligible.

both with respect to the measured observable Ô and to the rotation axis m of the state. In an experiment one has to consider in addition the finite resolution of the atom number counting in the modes a and b. If the detection system does not quite reach single atom resolution, one can still detect the cat state and determine the Fisher information if one chooses a spin observable in the xy plane, oriented along the direction of the fringes in Figs. 3 and4 respectively. The optimized results in Fig. 3 and Fig. 4 include Poissonian fluctuations of the total particle number, finite lifetime and particle losses for both states, that is one-body losses and, for 87 Rb, three-body and two-body losses including inter-component a-b losses.

In the lossless case, the maximal Fisher information achievable when starting with an initial phase state (3) with a total atom number N and a pulse angle θ is

F ( Ŝz ) = N 2 sin 2 θ + N cos 2 θ . (26)
In Fig. 3 and Fig. 4 we show as a function of time this maximal Fisher information with N = N a (t) + N b (t) and θ = arccos Na(t)-N b (t) Na(t)+N b (t) corresponding to the timedependent mean atom numbers in our system in the presence of losses. It is remarkable that, although 30 particles are lost on average in the majority component in Fig. 3 for 87 Rb atoms (see Appendix B), high contrast fringes are obtained in the Wigner function at the cat-state time, and the corresponding Fisher information is reduced by a factor less than half with respect to its maximal possible value with the same number of atoms. From equation [START_REF] Willems | [END_REF] it is apparent that a non-symmetric pulse θ = π 2 reduces the maximal Fisher information in the lossless case. The situation is however different in the presence of losses, where the best pulse angle, as well as the best trap parameters can only be derived from an optimization procedure that is specific to the selected transition and the atomic species plus the experimental constraints (see Appendix B).

IV. PHYSICAL INTERPRETATION OF THE REALISTIC-ANALYSIS OPTIMUM

We provide here a simple physical interpretation of the mixed cat state with maximal Fisher information numerically determined in section III for experimentally realistic conditions including losses and particle number fluctuations.

A. Analytical model in the general case

As in the numerical simulations of section III B, in the general non-symmetric case, we use a master equation of the form (5), with non-symmetric m-body loss rates (20) and ( 21), and with an initial condition (19) representing a statistical mixture of phase states with Poissonian fluctuations of the total atom number with average N . The difference here, in order to perform an analytical study and develop some intuition, is that the unitary evolution of the density matrix is not calculated with the fully nonlinear Hamiltonian (22) but with a non-symmetric S 2 z Hamiltonian [30]

Ĥ = χ N Ŝz + χ Ŝ2 z , (27)
where we omitted terms that give a constant phase drift or a global phase shift. The microscopic expressions of χ and χ are [30]

χ = 1 2 (∂ Na µ a + ∂ N b µ b -∂ Na µ b -∂ N b µ a ) Na , Nb , (28)
χ = 1 2 (∂ Na µ a -∂ N b µ b ) Na , Nb , (29)
where µ a and µ b are the chemical potentials of the a and b condensate respectively. We calculate χ and χ by solving the stationary two-component Gross-Pitaevskii equation in the considered trap geometry for different atom numbers around the averages Na = N cos 2 θ/2 and Nb = N sin 2 θ/2, while the loss parameters ( 20) and ( 21) are calculated as mentioned earlier for N a = Na and N b = Nb . As can be noted from Fig. 3 and Fig. 4 (see the captions), for all the optimal configurations found by our algorithm, one has |(χχ)/(χ + χ)| ≪ 1. We explain here the reason.

B. Compensation of random phase shifts

By evolving an initial phase state |θ; 0 N with the Hamiltonian (27) in the absence of losses, a cat state appears at the time t cat = π 2χ . We define the "unrotated cat" with N particles as

|unrot cat N ≡ 1 √ 2 |θ; 0 N + i e i π 2 N |θ; π N . ( 30
)
By using the result for the evolution with a pure χ Ŝ2

z

Hamiltonian [19] and by including the effect of the additional N -dependent drift term χ N Ŝz in the Hamiltonian, we obtain at the cat-state time:

|ψ(t cat ) N = e i N π 2 Ŝz(1-χ χ ) |unrot cat N . (31)
This shows that for χ = χ the dependence on N of the cat orientation, impossible to avoid when χ = 0 as in equation (7), is now eliminated.

Let us now consider the effect of one-body losses (with a rate γ) in one of the two components. Starting from a phase state, the trajectory with one atom lost at time t 1 in mode a or b, can be expressed in terms of a Hamiltonian evolution starting from a state with initially N -1 atoms plus a random t 1 -dependent shift of the relative phase:

| ψ1 a,b (t) = N a,b e -i ( χ±χ) Ŝz t1 t1-dependent drift e -iHt/ |ψ(t = 0) N -1 (32)
where N a,b includes a global phase and a normalization factor, and the plus or minus sign in the t 1 -dependent shift refers to a loss in component a or b, respectively. This shows that the random shift due to losses that comes from the quantum jump and from the N -dependent drift velocity in (27), can be set to zero in one of the two components by adjusting χ to ∓χ, the two effects compensating each other. In particular, for χ = χ both the random shift due to losses in b and the deterministic Ndependent rotation of the cat state that is present even without losses in (31) are suppressed [31]. This conclusion, based here on the analysis of a conditional state with a single lost atoms (32), holds also in the case of two-and three-body losses [32].

We now have to distribute the roles of a-mode and bmode to the two hyperfine states of 87

Rb |F = 1, m F = -1 and |F = 2, m F = 1 .
The key point is that the dominant loss process is two-body losses in |F = 2, m F = 1 . These are the ones that should be compensated. In addition there will be unavoidable one-body losses in the majority component. All the significant losses should be concentrated in a single component where they can be compensated. This explains the (at first sight) counterintuitive choice of taking |b = |F = 2, m F = 1 as the majority component, done in Fig. 2.

C. Coherent state description

A particularly simple interpretation of our results is obtained in the coherent state description that we will adopt in this subsection. To this aim we note that a Poissonian mixture of phase states for two modes is identical to a statistical mixture of Glauber coherent states with random total phase and a fixed relative phase

ρ = N N N e - N N ! |θ; ϕ N N θ; ϕ| = 2π 0 dΘ 2π |α, β α, β| (33)
where |α, β is a two-mode coherent state α = Na e iϕa and β = Nb e iϕ b , with ϕ = ϕ aϕ b the relative phase between the coherent states, Θ = ϕa+ϕ b 2 the total phase and N = Na + Nb the mean total atom number. To show the equality (33) one expands the phase states and the coherent states over Fock states |N a , N b . The integral over Θ suppresses coherences between Fock states with different total numbers of particles.

The next step is to remark that the Hamiltonian ( 27) can be elegantly written as the sum of two independent Hamiltonians plus a term that depends on N only,

Ĥ = χ a 2 N 2 a + χ b 2 N 2 b -χ N 2 4 , (34)
χ a = χ + χ ; χ b = χ -χ , (35)
despite the fact that the two modes overlap and interact with each other. The term that depends on N only is irrelevant because there are no coherences between states of different N . For each state |α, β appearing in the statistical mixture (33), the evolution of a and b modes under the influence of the Hamiltonian (34) and of losses other than a-b losses, is decoupled.

Evolution of the coherent states in the presence of losses

In the remainder of this section we consider the evolution of the two-mode coherent state |α, β under the influence of the Hamiltonian (34) and one-body losses. Although strictly speaking these states are not physical and the integral in (33) randomizing the total phase should be taken into account, the analysis gives some insight into the compensation condition, and it allows to introduce a fidelity that is not trivially zero in a case in which the total number of particles is not fixed.

Perfect compensation case -Let us consider the effect of one-body losses first in the case χ = χ that is χ b = 0. In section IV B we refer to this condition as "compensation" because in the Monte Carlo wave function approach, the random phase shifts coming from the losses and the Ndependent drift of the relative phase compensate. After the transformation (34) we can call it as well "no effective interactions in b". In this case, even in the presence of losses, the state of mode b remains a pure state: it is an exponentially decreasing coherent state

|ψ b (t) = | β where β = βe -γ b 2 t . (36)
This can be seen in the Monte Carlo wave function method, where, after renormalisation, we obtain | ψb (t) = | β for any quantum trajectory evolving under the influence of the non-hermitian Hamiltonian

H eff = -i 2 γ b b † b and k jumps with jump operator C = √ γ b b.
Since the mode b is effectively non-interacting (χ b = 0), it constitutes a perfect phase reference even in the presence of losses. Only its amplitude decreases in time. A similar conclusion was already reached in references [31,32].

In the absence of losses in a, γ a = 0, with χ a = 0, the mode a evolves as described in reference [8], going through a Schrödinger cat at time t cat = π χa = π 2χ and a revival at time t rev = 2π χa = π χ . In particular, for

|ψ a (0) = |α , we have |ψ 0 a (t cat = π χ a ) = 1 √ 2 e -i π 4 |α + e i π 4 | -α (37)
and

|ψ 0 a (t rev = 2π χ a ) = | -α . (38)
The exponent on ψ 0 a recalls that this is the ideal, lossless case in mode a.

What happens in the presence of one-body losses of rate γ a in mode a ? Something close to a cat state can only be obtained if these losses are very weak (less than one atom lost on average at the cat-state time). This means |α| 2 γ a t cat < 1 and hence γ a t cat ≪ 1. Within the Monte Carlo wave function approach, we introduce the non-normalized state vector | ψa (t) , corresponding to a trajectory for mode a where no atoms were lost in that mode at time t. Noting that the effective non-hermitian Hamiltonian can be written in a form equivalent to (34), as the sum of commuting parts, and introducing Ĥa ≡ χ a N 2 a /2, we have

| ψa (t) = e -γa Nat/2 e -i Ĥat |α = Ãe -i Ĥat |α (39)
with

α = αe -γa 2 t and à = e -|α| 2 2 e | α| 2 2 . ( 40
)
In the coherent state description, and before taking the integral over Θ, we define the fidelity of the state resulting from the evolution with losses as

F ≡ | ψa (t)|ψ 0 a (t) | 2 . ( 41
)
From the previous equations, at the cat-state time (neglecting for |α| ≫ 1 the vanishing overlap -α|α ) one then has

F (t cat ) = | Ã α|α | 2 = |e -|α| 2 (1-e -γa 2 t cat ) | 2 ≃ e -|α| 2 γatcat
(42) and similarly at the revival time

F (t rev ) ≃ e -|α| 2 γatrev ≃(F (t cat )) 2 . ( 43
)
Let us now look at the relative amplitude of the revival peak of the normalized g (1) function. For |α| 2 ≫ 1 we obtain:

g (1) (t rev ) = â (t rev ) â (0) ≃ -e -|α| 2 γatrev , (44)
showing that the amplitude of the revival peak directly gives informations on the cat-state fidelity

|g (1) (t rev )| = (F (t cat )) 2 . ( 45
)
This is again the relation (13), this time for coherent states and in the more general asymmetric case.

Note -The fact that one can restrict to the zero-loss subspace to define the cat-state fidelity is less clear when the target state is a coherent superposition of Glauber coherent states: the action of the jump operator â describing the loss of a particle does not render the coherent state |α orthogonal to itself (contrarily to the case of states with well defined particle numbers). The zero-loss subspace restriction performed in equations ( 42)-( 44) is however exact for the defined quantities F(tcat), F(trev) and â (trev) in the limit γatcat → 0 at fixed |α| 2 γatcat. This can be checked from the exact expressions, obtained using equations ( 6) and (7) of reference [18] to calculate the density operator ρa of mode a in presence of one-body losses:

F(t) = ψ 0 a (t)|ρa(t)|ψ 0 a (t) = k∈N γ k a k! [0,t] k dt1 . . . dt k e -γa k j=1 t j × exp -2|α| 2 1 -e -γat/2 cos(χa k j=1 tj ) (46)
â (trev) = Tr[âρa(trev)] = -e -γatrev /2 e -|α| 2 (1-e -γa trev ) × exp γa γa + iχa 1 -e -γatrev (47)
For the fidelity, one is helped by the fact that, in this large |α| 2 limit, the randomness of the particle-loss time, combined with the evolution with the quartic Hamiltonian Ĥa, effectively results at times of order 1/χa into a large random phase shift of the coherent state amplitude α.

What is actually measured in an experiment is | ab † |, where the expectation value is taken in (33) and the integral over Θ must be performed. This experimental contrast then reads

R exp ≡ | âb † |(t rev ) ≃ |α|e -|α| 2 γatrev |β|e -γ b
2 trev . (48) If the fraction of atoms lost in b at t = t rev is small, then e -γ b 2 trev ≃ 1 and one essentially recovers (44).

Imperfect compensation of the lossy mode -If χ b ≪ χ a but χ b = 0, there are some residual effective interactions in the mode b. As a consequence our phase reference starts to undergo a phase collapse. This modifies the contrast as follows: (49) plus small corrections due to losses in b. The compensation constraint becomes stringent for large atom numbers as one must have χ b /χ a ≪ 2/ √ πN b . If no compensation is done at all, that is χ b ≃ χ a , there will be no revival at all in R exp . Indeed as the mode b is lossy with |β| 2 γ b t rev > 1 it has a phase collapse with no revival.

R exp ≃ |α|e -|α| 2 γatrev |β|e -χ 2 b |β| 2 t 2 rev 2

V. MULTIMODE ANALYSIS OF THE CAT-STATE FORMATION: NONZERO TEMPERATURE EFFECTS

In this paper, up to now, we have analysed the quantum dynamics of the bosonic field in a two-mode model.

Reality is however multimodal, and there is always a nonzero thermal component in the initial state of the system, which can endanger the cat-state production even in the absence of losses. In this section we discuss nonzero temperature effects, both on the cat-state fidelity and on the contrast revival, in the Bogoliubov approximation.

A. Proposed experimental procedure

In the multimode case, one must revisit the definition of the initial state (2) and explain how to prepare it. In order to avoid any excitation induced by the π/2 pulse (see equation (28) in reference [33]), we assume as in reference [34] that the gas is initially non-interacting, g aa (0 -) = g ab (0 -) = g bb (0 -) = 0, and prepared at thermal equilibrium at the lowest accessible temperature T with all the N bosons in internal state |a . At time zero, to obtain a phase state, one applies an instantaneous π/2 pulse between the |a and |b states, which transforms the atomic field operators in the Heisenberg picture as follows:

ψa (r, 0 + ) = 1 √ 2 ψa (r, 0 -) -ψb (r, 0 -) (50)
ψb (r, 0

+ ) = 1 √ 2 ψa (r, 0 -) + ψb (r, 0 -) (51)
To obtain the nonlinear spin dynamics required to get a cat state, one adiabatically increases the interaction strength g aa (t) = g bb (t) = g(t) up to the final value g f in a time t ramp , while keeping g ab = 0, and one lets the system evolve until the much longer cat-state production time t cat or contrast revival time t rev .

Note -Experimentally, to suppress interactions, one can start with a condensate at low enough atomic density, perform the π/2 pulse and spatially separate the components a and b. The total number of Bogoliubov excitations created by the pulse in each component in the homogeneous case, for gaa(0 -) = g bb (0 -) = 4π 2 a(0 -)/m and g ab = 0, is N exc σ (0 + ) ≃ 0.395N ρa(0 -) 3 (from equation (38) and Appendix C of reference [35]) and should be ≪ 1. The condition g ab = 0 is ensured by spatial separation of the a and b components right after the pulse using state dependent potentials [14,36]. Note that the interaction dynamics are much slower than the π/2 pulse and the subsequent spatial separation. Once the components are split, the effective interaction strength is increased by adiabatically reducing the volume of the trapping potentials of the two components. Alternatively, an atomic species with a Feshbach resonance in state |a could be used, which allows tuning the interaction strength gaa = 0 [37]. The π/2 pulse could then be performed in real space (rather than on the spin degrees of freedom) with all atoms in |a by adiabatically ramping up a barrier in the trapping potential to split the atomic cloud. Subsequently, the Feshbach resonance is used to tune the interactions in both wells of the resulting double well potential to a nonzero value. Finally, if one prefers to avoid barrier splitting and interaction suppression by decompression, a possibility is to use spin-1 bosonic particles, with |a and |b the internal states of maximal spin ± along the quantization axis Oz as for example |F = 1, mF = ±1 . The spinor symmetry then imposes equal coupling constants gaa = g bb in the two states. Unfortunately, the internal scattering lengths of |F = 1, mF = ±1 are expected to have a magnetic Feshbach resonance at opposite values ±B0 of the magnetic field along the quantization axis Oz. Generically one thus cannot achieve gaa = g bb = 0 for a given value of such a magnetic field Bez. A first solution is to make B rapidly oscillate in time between opposite values such that on average gaa = g bb = 0. A second solution is to rapidly and coherently transfer back the b atoms into the internal state a after the π/2 pulse and the spatial separation of the two spin components, e.g. with a spatially resolved laser-induced Raman transition. -For simplicity, and taking into account recent experiments on degenerate gases in flat bottom potentials [38,39], we assume in this section that each spin component is trapped in a cubic box of volume V = L 3 with periodic boundary conditions. One can then take advantage of the fact that the Bogoliubov mode functions are plane waves with known amplitudes, which makes explicit calculations straightforward. As an immediate illustration, we give an adiabaticity condition for the interaction switching in Appendix C, for the Hann ramp

g(t) = g f 2 1 -cos πt t ramp for 0 < t < t ramp . (52)
B. Analysis at zero temperature

In the ideal limit of T = 0, the system is initially prepared in its ground state, with the N bosons in internal state |a with a vanishing wavevector k = 0. Just after the π/2 pulse, due to (50,51), the system is in the state

|ψ(0 + ) = 1 (N !) 1/2 2 N/2 [ĉ † a,0 (0 -) + ĉ † b,0 (0 -)] N |0 (53)
where the bosonic operator ĉσ,k annihilates a particle in internal state |σ with wavevector k and |0 is the vacuum. The binomial expansion gives

|ψ(0 + ) = 1 2 N/2 N Na=0 N ! N a !N b ! 1/2 × |N a : a, k = 0; N b : b, k = 0 . (54)
In this form, each Fock state is the ground state of the system at the considered fixed values of N a and N b = N -N a . Under adiabatic switching of the interaction strength, it is transformed into the instantaneous ground state |ψ 0 (N a , N b ; t) of the interacting system (taken with a real wavefunction), with instantaneous energy E 0 (N a , N b ; t). The global state of the system is then at time t:

|ψ adiab (t) = 1 2 N/2 N Na=0 N ! N a !N b ! 1/2 × e -i t 0 dτ E0(Na,N b ;τ )/ |ψ 0 (N a , N b ; t) . (55)
This defines the equivalent of the phase state and its evolution in the multimode theory. In the large N limit, one recovers a Ŝ2 z spin dynamics as in Eq. ( 4) by expanding E 0 (N a , N b ; τ ) around ( Na , Nb ) = (N/2, N/2) up to second order in N a -Na = -(N b -Nb ) = (N a -N b )/2. At t > t ramp , this gives (see the note in the next paragraph)

|ψ adiab (t) ≃ e -i t 0 dτ E0( Na, Nb ;τ )/ 2 N/2 e -iχ Ŝ2 z (t-t0) × N Na=0 N ! N a !N b ! 1/2 |ψ 0 (N a , N b ; t) (56)
with the collective spin operator Ŝz = ( Na -Nb )/2 and the spin nonlinearity coefficient

χ = 1 ∂ 2 E 0 ∂N 2 σ ( Na , Nb ; g = g f ) ( 57
)
where σ is any of the a, b. The phenomenology of catstate formation and contrast revival of the two-mode model is straightforwardly recovered, up to a retardation time t 0 due to the adiabatic ramping of the interaction,

t 0 dτ 1 ∂ 2 E 0 ∂N 2 σ ( Na , Nb ; τ ) t>tramp = χ(t -t 0 ) (58)
The pure state (56), and the resulting cat state at the appropriate time, exhibits entanglement between the external orbital degrees of freedom and the internal spin degrees of freedom. This entanglement can be eliminated by adiabatically ramping down the interaction strength to zero, to transform back each |ψ 0 (N a , N b ; t) into the Fock state |N a : a, k = 0; N b : b, k = 0 with spin-state independent orbital modes.

Note -If one expands E0(Na, N b ; τ ) in equation (55) up to fourth order in Na -Na = -(N b -Nb ) in the spirit of figures 3 and 4 (red curve vs green curve), one finds a state |ψquart(t) that differs from the state |ψ quad (t) resulting from the second order expansion (as given by Eq. ( 56)), because the Bogoliubov ground-state energy of the uniform gas is not purely quadratic in N (contrarily to the Gross-Pitaevskii approximation). At the first time tcat where |ψ quad (t) is the target cat state, we find an overlap of the form ψquart(t chat )|ψ quad (t chat ) = N π/2; ϕ = 0| exp(iα Ŝ4 z )|π/2; ϕ = 0 N with α ≃ π 24

∂ 4 Nσ E 0 ( Na, Nb ;g=g f ) ∂ 2
Nσ E 0 ( Na, Nb ;g=g f ) . In the phase state, Ŝ4 z = N (3N -2)/16. The small parameter controlling the expansion is thus ε = N π/2; ϕ = 0|α Ŝ4 z |π/2; ϕ = 0 N = -3(2π) 1/2 (ρa 3 f ) 1/2 /48 in the thermodynamic limit, where ρ = N/L 3 is the total density. For the parameters of Fig. 7 we find the very small value ǫ ≃ -0.001. This legitimates the quadratic expansion of E0 for the uniform gas. In reality, cubic box potentials correspond to hard walls rather than to periodic boundary conditions. For the parameters of Fig. 7, the healing length in a given spin component ξσ = /(mµσ) 1/2 is significantly smaller than the box size L, so to calculate the Gross-Piatevskii chemical potential µσ, we use the approximate condensate wavefunction φσ(r) = [µσ/(Nσg f )] 1/2 α=x,y,z tanh(rα/ξ) tanh[(L-rα)/ξ], knowing that the hyperbolic tangent form is exact for a single wall. The normalisation of φσ to unity leads to the equation of state in the box Nσ = µσ(L-2ξσ) 3 /g f . We obtain ξσ/L ≃ 0.17 and ǫ ≃ 0.01, which again validates the quadratisation of E0. -

C. Fidelity at nonzero temperature

In practice, the system is prepared at a nonzero temperature T . The fidelity of the cat-state preparation, less than one, can be obtained by the following general reasoning. Let us call Û the unitary evolution operator during t cat mapping the initial zero-temperature system state |ψ 0 (0 -) onto the cat state |cat (with fidelity one):

Û |ψ 0 (0 -) = |cat . ( 59
)
If the system is prepared in an initial state |ψ(0 -) orthogonal to |ψ 0 (0 -) , for example in an excited eigenstate, then the state produced at time t cat by the same preparation procedure will be orthogonal to the target state |cat , which corresponds to a zero fidelity. If the system is prepared in the density operator ρ, the cat state is obtained with a fidelity

F = ψ 0 (0 -)|ρ|ψ 0 (0 -) = P 0 ( 60
)
where P 0 is the probability that the system is initially in its ground state. In practice, ρ corresponds to the canonical ensemble at temperature T for an ideal gas in internal state |a . If T is small enough as compared to the critical temperature T c , one can consider that the condensate is never empty and one can relax the condition that the number of noncondensed particles is less than or equal to the total particle number [40,41]. In a given single-particle mode of wavevector k, the number of excitations n k then follows the usual exponential law

P k (n k = n) = [1 -exp(-βE k )]e -βE k n , n ∈ N (61) with β = 1/(k B T ), E k = 2 k 2 /(2m) and k ∈ 2π L Z 3 * .
The system is in its ground state if all modes are in their ground state. This leads to the cat-state fidelity

F = k =0 (1 -e -βE k ) . ( 62
)
This result is a universal function of k B T /∆, where ∆ = E 2π/L = 2 (2π/L) 2 /(2m) is the minimal excitation energy, that is the energy gap. It is plotted as a black solid line in Fig. 6. This shows that one must have initially a small number of excitations in the system in order to have a fidelity close to one, hence the stringent requirement on temperature:

k B T < ∆ 4 = 2 (2π) 2 8mL 2 (63)
We take for the box size L = 1 µm to have about the same chemical potential as in the harmonic trap, for the parameters of Fig. 1. The energy ∆/4 then corresponds to a temperature of 28 nK. This temperature is close to the range ≈ 40 -30 nK already accessed by direct in situ evaporative cooling [42,43]. It might also be reached by using as a coolant the subnanokelvin gases prepared in very weak traps [44]. , which are universal functions of kBT /∆, Fspin and F Bog spin depend on the particle number N and on the interaction strength, adiabatically ramped up and down between 0 and the a-a and b-b scattering length a f , see subsection V E. Here the parameters are the ones of figure 7: N = 300, 4πa f /L = 0.0667 and tramp = 20t adiab ramp . The temperature is expressed in units of ∆/kB, where ∆ = 2 (2π) 2 2mL 2 is the minimal excitation energy.

0

D. Contrast at nonzero temperature

We now calculate the g (1) (t) function of the condensate in the multimode case (see endnote [45]) within the Bo-goliubov approximation, in its U (1)-symmetry preserving version [46,47]. The number-conserving noncondensed fields Λσ (r, t) in spin state σ can be generally expanded over Bogoliubov modes, which are here plane waves. Due to g ab = 0, the two spin components decouple. Due to the adiabatic interaction ramp, the mode amplitudes are, up to a global phase factor, given by the instantaneous Bogoliubov steady state expressions. After the π/2 pulse, we thus get

Λσ (r, t) Λ † σ (r, t) adiab = k =0 bσ,k (0 + )e -i t 0 dτ ǫ k (τ )/ U k (t) V k (t) e ik•r V 1/2 + b † σ,k (0 + )e i t 0 dτ ǫ k (τ )/ V k (t) U k (t) e -ik•r V 1/2 (64)
with the instantaneous real amplitudes and energies

U k (t) + V k (t) = 1 U k (t) -V k (t) = E k E k + 2µ σ (t) 1/4 (65) ǫ k (t) = [E k (E k + 2µ σ (t))] 1/2 . ( 66
)
Here the instantaneous chemical potential µ σ (t) in internal state σ is given in the mean-field approximation µ σ (t) = g(t) Nσ /V where Nσ = N/2 is the mean number of particles in that spin state. The quasi-particle annihilation and creation operators bσ,k and b † σ,k obey the usual bosonic commutation relations at equal times. The operators for the numbers of quasi-particles nσ,k = b † σ,k bσ,k are constants of motion in the Bogoliubov approximation, which neglects the quasi-particle interactions, and coincide here with the particle number operators at time 0 + since the gas is still non interacting immediately after the pulse:

nσ,k (0 + ) = (ĉ † σ,k ĉσ,k )(0 + ) ∀k = 0 . ( 67
)
The first-order coherence function of the condensate in the multimode case is defined similarly to equation (12) as

g (1) (t) = ĉ † a,0 (t)ĉ b,0 (t) ĉ † a,0 (0 + )ĉ b,0 (0 + ) . ( 68
)
We use the usual modulus-phase representation ĉσ,0 = exp(i θσ )[ĉ † σ,0 ĉσ,0 ] 1/2 where θσ is the condensate phase operator in spin state σ, canonically conjugated to the operator number of particles in the condensate mode of that spin state, [ θσ , ĉ † σ,0 ĉσ,0 ] = -i, and we perform the usual approximation replacing the weakly fluctuating moduli by constants, so that g (1) (t) ≃ e -i[ θa(t)-θb (t)] ,

(69) which expresses the fact that the loss of contrast is due to the condensate phase spreading dynamics. At the Bogoliubov order, and neglecting rapidly oscillating terms of negligible contribution at long times, the phase evolution after the π/2 pulse is given by [48] -

d θσ (t) dt = µ 0 ( Nσ ) + k =0 dǫ k dN σ ( Nσ )n σ,k (0 + ) (70)
with µ 0 (N σ ) the zero-temperature Bogoliubov chemical potential of a single-component gas with N σ particles [49]:

µ 0 (N σ ) = g V N σ - 1 2 + g V k =0 V k (N σ )[U k (N σ ) + V k (N σ )] + N σ g 2V E k + g 2 N σ V   d 3 k (2π) 3 1 2E k - 1 V k =0 1 2E k   . ( 71
)
The last contribution is a finite size effect; the difference between the integral and the sum between square brackets was evaluated in reference [50] to be mC(0)/[(2π

) 2 L]
with C(0) ≃ 8.91364. Linearising the dependence of µ 0 and dǫ k dNσ in Nσ around Nσ and integrating over time, we obtain (see endnote [51])

( θa -θb )(t) = ( θa -θb )(0 + ) -A(t)( Na -Nb )(0 + ) - k =0 γ k (t)(n a,k -nb,k )(0 + ) . (72)
The time-dependent, dimensionless coefficients are given by

A(t) = t 0 dτ dµ 0 dN σ ( Nσ , τ ) + k =0 d 2 ǫ k dN 2 σ ( Nσ , τ ) nσ,k (0 + ) , (73)
γ k (t) = t 0 dτ dǫ k dN σ ( Nσ , τ ) . ( 74
)
They are affine functions of t for t > t ramp . In particular, one has

A(t) t>tramp = χ T (t -t 0,T ) (75)
where the spin nonlinearity coefficient and the retardation time, contrarily to the coefficients γ k , are now temperature dependent (see endnote [52]):

χ T = dµ 0 dN σ ( Nσ , g = g f )+ k =0 d 2 ǫ k dN 2 σ ( Nσ , g = g f ) nσ,k (0 + ) .
(76) To calculate the expectation value in (69), we first exponentiate the relation (72), separating the various contributions in three mutually commuting groups: from left to right, a first group containing the operators θa (0 + ) and Na (0 + ), a second one containing the quasi-particle number operators nσ,k (0 + ) and a third one containing θb (0 + ) and Nb (0 + ). The Baker-Campbell-Hausdorff formula for two operators X and Ŷ applied to the first group and to the third group reduces to exp( X + Ŷ ) = exp( X) exp( Ŷ ) exp(-1 2 [ X, Ŷ ]), hence to e -i θa(0 + )+iA(t) Na(0 + ) = e -i θa(0 + ) e iA(t) Na(0 + ) e iA(t)/2 (77)

e -iA(t) Nb (0 + )+i θb (0 + ) = e -iA(t) Nb (0 + ) e i θb (0 + ) e -iA(t)/2 (78) since the commutator of θσ (0 + ) and Nσ (0 + ) is proportional to the identity. It remains to express from equations ( 50) and ( 51) the various post-pulse operators in terms of the pre-pulse operators:

( Na -Nb )(0 + ) = -(ĉ † a,0 ĉb,0 + ĉ † b,0 ĉa,0 )(0 -) , + k =0 (n a,k -nb,k )(0 + ) (79) (n a,k -nb,k )(0 + ) = -(ĉ † a,k ĉb,k + ĉ † b,k ĉa,k )(0 -) , (80) θσ (0 + ) = θa (0 -) + O(N -1/2 ) , (81)
using Eq. ( 67) for the first two identities and reference [35] for the third one, and to take the expectation value first in the vacuum state of the ĉb,k (0 -) (see the note in the next paragraph) and then in the thermal state of the ĉa,k =0 (0 -), eliminating the mode (a, k = 0) using conservation of particle number, ĉ † a,0 ĉa,0 = Nk =0 ĉ † a,k ĉa,k . We also need nσ,k (0 We finally obtain the Bogoliubov prediction for the first order coherence function for N bosons prepared at temperature T , with an interaction ramped up after the π/2 pulse from its initial value 0 to its final value g f :

+ ) = 1 2 [exp(βE k ) -1] -1 for k = 0,
g (1) (t) ≃ cos N -1 [A(t)] k =0 1 -e -βE k 1 -cos[γ k (t)+A(t)] cos[A(t)] e -βE k (82)
keeping in mind that, after the ramp, that is at times t > t ramp , A(t) = χ T (tt 0,T ) as in (75). In figure 7a, we plot this prediction as a function of time for various temperatures. As it is apparent from expression (82), g (1) (t)

results for large N from a narrow function cos N -1 [A(t)] selecting thin temporal windows in a slowly varying envelope function. At low temperature, when only a few noncondensed modes are populated, the envelope function oscillates in time, which results in a nonmonotonic behavior of the height of the successive revival peaks as one can see in the figure. In particular for k B T /∆ = 0.25 the third revival is almost perfect. Indeed, as | cos N -1 [A(t)]| is very small for large N , except for the A(t) integer multiple of π, one can to a good approximation replace A(t) by such an integer multiple in the product over k in equation ( 82), which results in the envelope function

G(t) = k =0 1-e -βE k
1-cos[γ k (t)]e -βE k . At k B T ≪ ∆, only the ground noncondensed mode degenerate multiplicity k = 2π/L is significantly populated, leading to an almost periodic function G(t) oscillating between ≃ 1 and ≃ F 2 with an angular frequency d dt γ 2π L (t), where F is the catstate fidelity (62). In figure 7b, we show the value of g (1) (t) at the first revival time of the cosine prefactor, χ T (tt 0,T ) = π, as a function of temperature. Both plots show that thermal excitations essentially destroy the revival, except at temperatures below the first excited mode energy ∆.

E. Spin fidelity at nonzero temperature

The fidelity considered in subsection V C is an orbitospinorial fidelity: it measures the overlap of the actual physical state of the system with the target state (56), which is an orbito-spinorial cat state. In pratical applications, however, one mainly measures pure spin observables, that do not act on the orbital part of the manybody state. This is the case for the collective spin operator Ŝz used as a reference observable in the Fisher information of section III. It is then more appropriate to consider a spin fidelity F spin . The main question is whether or not this spin fidelity is significantly less sensitive to nonzero temperature effects than the full fidelity F . This question is answered in this subsection.

For simplicity we assume in this subsection that N is an integer multiple of 4 so that the cat spin state is, according to equation (9),

|spin cat N = |+ N + i|-N √ 2 (83)
where |ǫ N , with ǫ = ±1, is the collective spin state with all the N spins in the same state (|a + ǫ|b )/ √ 2. In a Bose gas with an orbito-spinorial density operator ρ(t) at time t, the cat spin state is realized with a spin fidelity

F spin (t) = Tr [ρ(t)|spin cat N N spin cat|] . (84)
We assume that at the initial time t = 0 -, the Bose gas in a single realisation occupies in the internal state |a a N -boson Fock state |ψ(0 -) . This state samples the ideal gas thermal equilibrium density operator and is

FIG. 7: First order coherence function g (1) (t) of the condensate, as given by equation (68) in the multimode Bogoliubov approximation (82). N = 300 lossless 87 Rb atoms are initially prepared in internal state |a = |F = 1, mF = -1 in a box [0, L] 3 with periodic boundary conditions, L = 1 µm, at temperature T in the absence of interactions. At t = 0 + they are subjected to a π/2-pulse towards the internal state |b = |F = 1, mF = 1 and to an adiabatic Hann ramping (52) of the interaction strength with a duration tramp = 20 t adiab ramp = 0.21 ms, where t adiab ramp is defined in (C10), up to the final value g f = 4π 2 a f /m, with the a-a and b-b scattering lengths a f = 100.4 a0. In (a): g (1) as a function of time for temperatures kBT = 0 (black solid line), kBT = ∆/4 (red solid line), kBT = ∆/2 (blue solid line), from bottom to top at the first and third revival times and from top to bottom at the second revival time. In (b): its absolute value as a function of temperature at the first revival time (at which |g (1) (t)| has a maximum). We recall that ∆ = 2 (2π) 2 2mL 2 is the energy of the first single-particle excited state in the box.

characterised by the occupation numbers (n k ) k∈ 2π L Z 3 of the single-particle modes of wavevectors k in the quantisation volume [0, L] 3 :

|ψ(0 -) = k (ĉ † a,k ) n k (n k !) 1/2 |0 (85)
We use here the Schrödinger picture. After the instan-taneous a ↔ b π/2 pulse, the state of the Bose gas is, according to equations (50,51),

|ψ(0 + ) = k ĉ † a,k +ĉ † b,k √ 2 n k (n k !) 1/2 |0 = (n a,k ) k k P n k (n a,k ) 1/2 |a : (n a,k ) k , b : (n b,k ) k (86)
In the second form, obtained from the first one by using the binomial theorem, the ket is the Fock state with mode occupation numbers n a,k in internal state a and n b,k in internal state b, P n (n a ) = n!2 -n na!n b ! is the classical binomial probability that n incoming particles are split into n a particles in the output channel a and n b = nn a in the output channel b, one sets n b,k = n kn a,k and the sum runs over all the occupation numbers (n a,k ) k such that 0 ≤ n a,k ≤ n k .

The system then evolves as follows during the time t. One switches on adiabatically the interaction strength g aa = g bb from 0 to g f during t ramp with the Hann halframp (66) (we recall that g ab = 0 at all times). The interaction strength remains constant during the time t -2t ramp . It is then switched off adiabatically over the time interval [tt ramp , t] with the time-reversed Hann half-ramp. In this process, the ideal gas Fock state |a : (n a,k ) k , b : (n b,k ) k is adiabatically turned into a Fock state of Bogoliubov quasi-particles, with the same occupation numbers for k = 0 (see endnote [53]), and is turned back into itself when the interactions are switched off, up to a global phase shift given by the time integral of the instantaneous eigenenergy E divided by , with

E((n a,k ) k , (n b,k ) k , t) = σ=a,b   E 0 (N σ , t) + k =0 ǫ k (N σ , t)n σ,k   (87)
The Bogoliubov eigenenergies ǫ k as functions of the total number of particles N σ = k n σ,k in internal state σ are given by equation (66). The ground-state energy E 0 (N σ ) of N σ interacting bosons in internal state σ is obtained by integrating µ 0 (N σ ) over N σ in equation (71) (knowing that E 0 (N σ = 0) = 0). So at time t the Bose gas is in the state

|ψ(t) = (n a,k ) k k P n k (n a,k ) 1/2 × e -i t 0 dτ E((n a,k ) k ,(n b,k ) k ,τ )/ |a : (n a,k ) k , b : (n b,k ) k (88)
As shown in Appendix D, the spin fidelity corresponding to the single-realisation density operator ρ(t) =

|ψ(t) ψ(t)| is then F single spin (t) = (n a,k ) k 1 -i(-1) Sz √ 2 k P n k (n a,k ) × e -i t 0 dτ E((n a,k ) k ,(n b,k ) k ,τ )/ 2 (89)
where S z = (N a -N b )/2 (see endnote [54]). It remains to average this result over the thermal canonical distribution of the (n k ) k in the initial ideal Bose gas to obtain the sought spin fidelity,

F spin (t) = (n k ) k =0 F single spin (t) k =0 (1-e -βE k )e -βE k n k (90)
where

E k = 2 k 2
2m and the number of condensate particles n 0 is adjusted in each realisation to have a fixed total number N of particles, n 0 = Nk =0 n k . This average can in practice be taken with a Monte Carlo simulation, and the local maximum of F spin (t) close to the expected cat-state formation time t cat can be found numerically. The resulting spin fidelity for the physical parameters of Fig. 7 is plotted as symbols with error bars in Fig. 6, as a function of temperature. As expected, it is larger than the orbito-spinorial fidelity F , plotted as a black (lower) solid line in that figure. Unfortunately, over the temperature range where it is larger than 1/2, the spin fidelity is only slightly larger than the orbito-spinorial fidelity F . This means that the stringent temperature requirement (63) also applies to the spin cat-state formation.

We now go through a sequence of approximations to get a more inspiring analytical result and some physical explanation of the sensitivity of F spin to temperature. First, as we did in section V D, we take advantage of the fact that, in the large N limit, N σ has weak relative fluctuations around its mean value Nσ = N/2, (N σ -Nσ )/N ≈ N -1/2 . Expanding the energy E((n a,k ) k , (n b,k ) k , t) in equation (87) to second order in N σ -Nσ and replacing the coefficients of the quadratic terms by their thermal averages we obtain

t 0 dτ E((n a,k ) k , (n b,k ) k , t)/ ≃ t 0 dτ [2E 0 ( Nσ , τ ) + k =0 ǫ k ( Nσ , τ )n k ]/ +   k =0 γ k (t)(n a,k -n b,k )   S z + A(t)S 2 z (91)
The first contribution in equation (91) does not depend on S z nor on the (n a,k ) k =0 and it will not contribute at all to the spin fidelity. In the other contributions, the time-dependent coefficients γ k (t) and A(t) are given by equations (73,74). Second, in the spirit of the particlenumber-conserving Bogoliubov methods [46,47], we take as independent variables in each internal state σ the total number of particles N σ and the occupation numbers (n σ,k ) k =0 of the Bogoliubov modes. This is here an approximation as the π/2 pulse introduces a small correlation between the difference of the total particle numbers N a -N b and the difference of the noncondensed particle numbers N a,nc -N b,nc , of the order of f 1/2 nc , where f nc is the initial noncondensed fraction (see endnote [55]). In practice, in expression (89), we perform to leading order in f nc the substitution n0 na,0=0

P n0 (n a,0 ) → N Na=0 P N (N a ) (92)
Finally, using the identity

1 -i(-1) Sz √ 2 = e -i π 4 e i π 2 S 2 z ( 93
)
valid for N integer multiple of 4, we obtain the Bogoliubov approximation for the single realisation spin fidelity

F single,Bog spin (t) = N Na=0 P N (N a )e iS 2 z [ π 2 -A(t)] × (n a,k ) k =0   k =0 P n k (n a,k )e -iSzγ k (t)(n a,k -n b,k )   2 (94)
The sum over each n a,k from 0 to n k can be evaluated analytically with the binomial theorem. The modulus square of the result can be averaged analytically over the thermal distribution of the (n k ) k =0 using α n k = [1 + nk (1α)] -1 here with |α| ≤ 1, which gives

F Bog spin (t) = N Na,N ′ a =0 P N (N a )P N (N ′ a ) cos[(S 2 z -S ′ z 2 )(A(t) -π 2 )] k =0 {1 + nk [1 -cos(γ k (t)S z ) cos(γ k (t)S ′ z )]} (95)
where nk = (e βE k -1) -1 is the initial mean occupation number of the mode k in the internal state a. The approximate result (95) is readily evaluated numerically as a function of time for the physical parameters of Fig. 7.

It is found that the sought spin fidelity peak is located extremely close to the expected cat-state formation time t cat , such that A(t cat ) = π 2 , and its value, plotted as a red (upper) solid line in Fig. 6, is in very good agreement with the Monte Carlo results (red circles) resulting from the full expression (89).

Note -The approximation made in (91), consisting of the replacement of coefficients of the quadratic terms by their thermal averages, is not essential. Without it, one obtains

F Bog spin (t) = N Na,N ′ a =0 PN (Na)PN (N ′ a )e -i(S 2 z -S ′2 z )[A 0 (t)-π 2 ] k =0 [1 + nk (1 -D k )]
where A0(t) is the zero temperature value of A(t),

D k = e -iα k (S 2 z -S ′2 z ) cos(γ k (t)Sz) cos(γ k (t)S ′ z ) and α k (t) = 1 2 t 0 dτ ∂ 2 Nσ ǫ k ( Nσ, τ ) .
We have verified that this result is very close to the less refined approximation (95) for the parameters of Fig. 6. -A physical insight in the temperature sensitivity of the spin fidelity is obtained by rewriting the single realisation spin fidelity at the cat-state time t cat from equation ( 94) as

F single,Bog spin (t cat ) = cos N ∆θ th 2 partition 2 (96)
where the average is taken over the partition noise in the noncondensed modes accompanying the π/2 pulse, that is with the binomial probability distribution P n k (n a,k ) for each n a,k , and

∆θ th = k =0 γ k (t cat )(n a,k -n b,k ) (97
that is in a coherent superposition of rotated cat states (what appears here is the operator Ŝz ). As the cat-state time scales as N , the coefficients γ k (t cat ) scale as N 0 and are of order unity. This shows that the presence of a single thermal excitation in the initial state of the system, by activating the partition noise, will give quantum fluctuations of ∆θ th of order unity, sufficient to compromise the cat-state fidelity (see endnote [57]). This high sensitivity to thermal excitations was anticipated in reference [58]. Equation ( 96) is not only physically appealing, it also provides a lower bound to the peak spin fidelity F Bog spin (t cat ) that is very close to the actual value for large N . Indeed, when N ≫ 1 in a fixed volume and at a fixed temperature, cos N (∆θ th /2) is a very narrow function of ∆θ th with a width smaller than the discreteness of the distribution of ∆θ th . As a consequence, only the realisations with ∆θ th = 0 contribute significantly. As the coefficients γ k for different wave numbers k are in general incommensurable, this imposes that for all allowed wave numbers k = 0 in the quantization volume:

k / ||k||=k (n a,k -n b,k ) = 0 . ( 99
)
For a given realisation of the initial thermal occupation numbers n k , this occurs for a given k with the binomial probability P N k (N k /2) where N k is the total number of initial thermal excitation in the degenerate manifold of wave number k, N k = k / ||k||=k n k . Note that this probability is zero for odd N k . As N is here even, cos N is nonnegative and we obtain after thermal average the inequality F Bog spin (t cat ) ≥ F Bog, minor spin with the lower bound

F Bog, minor spin ≡ k =0 [P N k (N k /2)] 2 therm = F k =0 3 F 2 1 2 , d k 2 , 1 + d k 2 ; 1, 1; e -2βE k ( 100
)
where

F

VI. CONCLUSION

We have studied the interaction-induced formation of mesoscopic quantum superpositions in bimodal Bose-Einstein condensates including limiting effects such as particle losses and fluctuations of the total number of particles. We have explained how these effects can be compensated, giving two examples for sodium and rubidium Bose-Einstein condensates. To quantify the survival of quantum correlations in the presence of decoherence, and their usefulness for metrology, we have calculated the Fisher information and the Wigner function of the obtained state, and we have also shown that, in the presence of losses, there is a simple quantitative relation between the cat-state fidelity and the amplitude of the revival peak in phase contrast. Finally, giving up the twomode description in a last multimode section, we have described a possible procedure to prepare the initial state, and we have studied the influence of a nonzero initial temperature on the amplitude of the phase revival and on the cat-state fidelity. Two fidelities are introduced: the full orbito-spinorial fidelity F and a purely spinorial fidelity F spin , defined once the orbital degrees of freedom are traced out. We find that macroscopic superpositions can be obtained at nonzero temperature with a high fidelity, with no substantial gain of F spin with respect to F , provided the temperature is lower than about one quarter of the energy of the first single-particle excited state.

FIG. 10: The algorithm searches for good configurations in the parameter space by maximizing the Fisher information and verifies the presence of a small-amplitude cat state signaled by high contrast fringes in the Wigner function. : Budget of losses at time tcat . Parameters are the same as in Fig. 3.

with the ideal gas initial conditions U k (0

+ ) = 1, V k (0 + ) = 0.
Here ρ σ = N σ /V is the density in component σ.

In the quasi-adiabatic regime it is convenient to project (U k (t), V k (t)) onto the instantaneous stationary Bogoliubov mode real amplitudes (U k (t), V k (t)) of (65) of energy ǫ k (t) and on the corresponding mode (V k (t), U k (t)) of energy -ǫ k (t):

U k (t) V k (t) = A k (t) U k (t) V k (t) + B k (t) V k (t) U k (t) (C3) with A k (t) = U k (t)U k (t) -V k (t)V k (t) (C4) B k (t) = -V k (t)U k (t) + U k (t)V k (t) (C5) leading to the differential system i d dt A k (t) B k (t) = ǫ k (t) -i Ω k (t) -i Ω k (t) -ǫ k (t) A k (t) B k (t) (C6) with the initial conditions A k (0 + ) = 1, B k (0 + ) = 0. The symplectic symmetry imposes |U k (t)| 2 -|V k (t)| 2 = |A k (t)| 2 -|B k (t)| 2 = 1. The Rabi angular frequency Ω k (t) = U k (t) d dt V k (t) -V k (t) d dt U k (t) = d dt [U k (t) + V k (t)] U k (t) + V k (t) = - 1 2 ρ σ d dt g(t) E k + 2ρ σ g(t) (C7)
constitutes the nonadiabatic coupling. The number of quasi-particle excitations in the stationary Bogoliubov mode k present at the end of the ramp is given by

n exc k = |B k (t > t ramp )| 2 . (C8)
The evolution enters the adiabatic regime when the Rabi coupling is much weaker than the Bohr frequency:

|Ω k (t)| ≪ 2ǫ k (t) . (C9)
This is most stringent at the minimal wavenumber k

= 2π/L. For 2ρ σ g f ≫ ∆ = 2 (2π/L) 2 2m
, this is then most stringent at times ≪ t ramp , where the Hann expression (52) can be quadratised (see endnote [62]). One finally gets from the adiabaticity condition (C9):

t ramp ≫ t adiab ramp = L 3 2 (ρ σ g f ) 1/2 m 3/2 24π 2 √ 3 . (C10)
The corresponding time scale is much shorter than the cat-state formation time t cat ≃ π/(2χ) since the gas is weakly interacting:

χt adiab ramp ≃ (ρ σ a 3 f ) 1/2 3(3π) 1/2 ≪ 1 (C11)
with χ ≃ g f /( L 3 ) and g f = 4π 2 a f /m. In the quasi-adiabatic regime, one can treat the Rabi coupling to first order in time dependent perturbation theory, replacing in the equation for B k (t) the amplitude A k (t) by its zeroth-order, adiabatic expression Na on the x-axis is in the initial mean number of atoms in the minority component. For a given Na the different points correspond to successive configurations explored by the algorithm in its convergence process. Restricting to configurations with trap aspect ratios smaller than 20 we select out the darker points (in blue). Scattering lengths and loss rates are as in Fig. 3 for 87 Rb.

n exc k ≃ tramp 0 dt Ω k (t)e -2i t 0 dτ ǫ k (τ )/ 2 . (C12)
As the Hann ramp (52) leads to vanishing derivatives d dt g at t = 0 and t = t ramp , the number of excitations drops rapidly with k:

n exc k ∼ k→+∞ ρ σ g f 8E 3 k π t ramp 2 2 cos 2 [(E k + ρ σ g f /2)t ramp ] . (C13)
A numerical calculation of (C12) for the parameters of Fig. 7 (N σ = N/2 = 150 and 4πa f /L = 0.0667) confirms the condition (C10): for t ramp = t adiab ramp , the total number of created excitations in each spin component is ≃ 0.5; it drops to ≃ 0.01 for t ramp = 20 t adiab ramp .

Appendix D: Details on the calculation of the spin fidelity

In this appendix we derive expression (89) for the spin fidelity of the state |ψ(t) in equation ( 88) with respect to the cat state (83). To this end, it suffices to calculate the matrix element of a purely spinorial physical observable Ôspin between Fock states with occupation numbers

(n σ,k ) k and (n ′ σ,k ) k : X = a : (n ′ a,k ) k , b : (n ′ b,k ) k | Ôspin |a : (n a,k ) k , b : (n b,k ) k (
D1) This is conveniently evaluated in the first quantization formalism, where the Fock state reads

|a : (n a,k ) k , b : (n b,k ) k = N ! s j=1 n a,kj !n b,kj ! 1/2 × Ŝ|a, k 1 ⊗n a,k 1 |b, k 1 ⊗n b,k 1 . . . |a, k s ⊗n a,ks |b, k s ⊗n b,ks (D2)
where we labeled the populated wave vectors as k 1 , . . . , k s and used the notation |u ⊗n = |u ⊗ . . . ⊗|u (n factors). We have introduced the symmetrisation operator

Ŝ = 1 N ! σ∈SN Pσ (D3)
where the sum runs over all permutations σ of N elements and Pσ is the usual permutation operator representing σ in the Hilbert space. As Ôspin commutes with Ŝ, due to the indistinguishability of the particles, and as Ŝ2 = Ŝ, it is enough to symmetrise the ket only, which gives

X = σ∈SN n ′ a,k 1 ⊗ a, k 1 | n ′ b,k 1 ⊗ b, k 1 | . . . n ′ a,ks ⊗ a, k s | n ′ b,ks ⊗ b, k s | Ôspin Pσ |a, k 1 ⊗n a,k 1 |b, k 1 ⊗n b,k 1 . . . |a, k s ⊗n a,ks |b, k s ⊗n b,ks s j=1 n a,kj !n b,kj !n ′ a,kj !n ′ b,kj ! 1/2 (D4)
In the matrix element of the numerator, one can move the orbital part k j | of the bras through Ôspin to contract them with the orbital part of the kets. As the quantum state (88) results from the a-b partition of an initial Fock state in internal state a with occupation numbers

(n k ) k , one has n a,k + n b,k = n ′ a,k + n ′ b,k = n k , ∀ k.
As a consequence, the only permutations σ that can give a nonzero contribution are those who leave stable (or setwise invariant) the subsets corresponding to a given k, that is {1, . . . , n k1 }, {1 + n k1 , . . . , n k2 + n k1 }, . . ., {1 + Nn ks , . . . , N }. This gives the purely spinorial expression:

X = n ′ a,k 1 ⊗ a| n ′ b,k 1 ⊗ b| . . . n ′ a,ks ⊗ a| n ′ b,ks ⊗ b| Ôspin σ1 Pσ1 |a ⊗n a,k 1 |b ⊗n b,k 1 . . . σs Pσs |a ⊗n a,ks |b ⊗n b,ks s j=1 n a,kj !n b,kj !n ′ a,kj !n ′ b,kj ! 1/2 (D5)
where in the sums the permutation σ j runs over S n k j the permutation group of n kj elements. It remains to take for Ôspin the orthogonal projector on the spin cat state (83), Ôspin = |spin cat N N spin cat|, to obtain (see endnote [63]):

X = [1+i(-1) N ′ b ] √ 2 [1-i(-1) N b ] √ 2 k n k ! 2 N k n a,k !n b,k !n ′ a,k !n ′ b,k ! 1/2 (D6)
where N σ = k n σ,k . Since this is factorisable in a function of the (n σ,k ) k times a function of the (n ′ σ,k ) k , it finally leads to the desired expression (89) of the spin fidelity of the single realisation (88), knowing that (-1) N b = (-1) Sz for N/2 even integer.

In the remaining part of this appendix, we give a justification to the writing (98) of the spin state vector, which led to an enlightening interpretation of expression (96) for the single realization peak fidelity in the Bogoliubov approximation. To this aim we rewrite equation (D5), where the orbital degrees of freedom have been traced out, in the form × Ŝpartial |a ⊗n a,k 1 |b ⊗n b,k 1 . . . |a ⊗n a,ks |b ⊗n b,ks (D8) and |χ ′ defined in the same way with (n σ,k ) k replaced by (n ′ σ,k ) k . The projector Ŝpartial performs a partial symmetrization restricted to the aforementioned permutations, forming a subgroup G of S N , that leave setwise invariant the subsets corresponding to a given k, that is {1, . . . , n k1 }, {1 + n k1 , . . . , n k2 + n k1 }, . . ., {1 + Nn ks , . . . , N }: where P n (n a ) = 2 -n n! na!n b ! (with n b = nn a ) is the binomial probability distribution. Note that |ψ spin (t) is not bosonic as it is only partially symmetrized. However, if we are interested in the spin dynamics in the phase space bosonic sector, which is enough to study the spin cat-state formation, we can perform the full symmetriza-tion and consider Ŝ|ψ spin (t) which amounts to replacing Ŝpartial with Ŝ. In the spirit of the Bogoliubov approximation, we further perform the substitution (92) and quadratize the energy around Na and Nb as in equation (91) to finally obtain Ŝ|ψ Bogol spin (t) = e -iC(t) N Na=0

X = χ ′ | Ôspin |χ (
Ŝpartial = 1 k n k ! σ∈G
[P N (N a )] 1/2 e -iA(t)S z e -i k =0 γ k (t)(n a,k -n b,k )Sz |N a : a, N b : b partition (D12)

where |N a : a, N b : b is a spin Fock state, C(t) is the partition-independent integral contribution on the righthand side of equation (91), and the brackets indicate the average over the partition noise in the noncondensed modes. The value of Ŝ|ψ Bogol spin (t) at the cat-state time (such that A(t) = π/2) reproduces equation (98), and its scalar product with the target state (30) reproduces the form (96). The writing (D12) of the state shows that the effect of finite temperature is captured by a two-mode model, see equation (8), supplemented by a dephasing environment. This kind of model was already used in the context of spin squeezing [19,64] in particular to predict the optimum spin squeezing at finite temperature [64]. Contrarily to references [19,64], the stochastic element enters here in two different ways: the average over the partition noise in the noncondensed modes results in a coherent superposition of kets, while the average over the initial thermal excitations in component a is a classical average at the level of the density matrix which results in a statistical mixture.

To be complete, we give the expression of the mean value of the total spin, which is along the x axis for the considered initial state of the system. To this end, we give another writing of equation (D11): )n b,k1 ] 1/2 |n a,k1 +1 : a, n b,k1 -1 : b , and so forth for the n k2 , n k3 , . . . subsequent terms. By performing the energy quadratization (91) but not the Bogoliubov substitution (92), we obtain the single realization result as a sum of the contributions of the various single particle modes:

ψ(t)| Ŝ+ |ψ(t) = q C q (t) . (D14)
The condensate contribution is

C 0 (t) = n 0 2 cos N -1 [A(t)] k =0
cos[γ k (t) + A(t)] cos A(t)

n k (D15)
while the noncondensed q = 0 mode contribution is

C q (t) = n q 2 cos N -1 [A(t) + γ q (t)] × k =0
cos[γ k (t) + γ q (t) + A(t)] cos[γ q (t) + A(t)]

n k -δ k,q (D16)
where δ k,q is a Kronecker delta. If in equation (D15) one approximates n 0 by its mean value and one performs the thermal average over the n k , one recovers exactly the result (82) for the condensate first order coherence function. Interestingly the contributions of the noncondensed modes to the mean spin have different revival times than the condensate. As a consequence they do not contribute to the major peaks in Ŝx (t), they contribute to side peaks of very small relative amplitudes O(f nc ) (f nc is the initial noncondensed fraction). mean value, which introduces an error O(N -1 tfnc) where fnc is the noncondensed fraction; at the revival time this introduces a small O(fnc) error.

[52] The nonzero temperature correction to χ was missed in reference [35].

[53] The number of particles in the condensate modes is no longer well defined, but what matters in the number conserving Bogoliubov theory is the total number of particles Nσ in each spin state, which is well defined.

[54] The result (89) holds under the assumption that the single-particle orbital states are the same in the internal states a and b. In practice, this means that a and b experience the same trapping potential. If the a and b traps were spatially translated to ensure an effective g ab = 0 coupling constant, they must be translated back to the same location at the cat-state time. where Var stands for the variance and Cov for the covariance.

[56] This state has a norm less than one because it corresponds to the restriction of the spin density operator to the bosonic sector. This spin density operator, obtained as a trace of the full density operator over the orbital variables, can indeed populate various irreducible representations of the permutation group SN , not simply the bosonic one, when the orbital variables do not occupy a single mode of the field. We give more details in Appendix D.

[57] This reasoning fails if γ k = 2π for the considered single excitation. This is not possible at the cat-state time such that A(tcat) = π/2, as in the weakly interacting regime and for tcat ≫ tramp, γ k /A ≃ [E k /(E k + 2g f ρσ)] 1/2 < 1. This would be possible for the third cat-state time A(t

cat ) = 5π/2 for well chosen parameters such that γ k /A ≃ 4/5. For the geometry and interaction parameters of Fig. 6, this is approximately realized in the first excited manifold k = 2π/L for a total particle number N = 168 (if one restricts to N = 0 modulo 4). As the first excited manifold then does not contribute to the dephasing, the spin fidelity will be limited by the second excited manifold at energy 2∆. This clever tuning allows to increase the required temperature by a factor of two.

  âm , ρ â † m + âm ρ, â † m

FIG. 1 :

 1 FIG.1: Fidelity and absolute value of contrast versus time for a split Bose-Einstein condensate of N = 300 |F = 1, mF = -1 87 Rb atoms in two identical and spatially separated harmonic potentials in the presence of three-body losses. Scattering length a = 100.4a0, trapping frequency ω/2π = 500 Hz, with a three-body loss constant rate K3 = 6 × 10 -42 m 6 /s[START_REF] Egorov | [END_REF]. This gives χ = 12 s -1 and γ (3) = 2.6 × 10 -7 s -1 .

FIG. 2 :

 2 FIG. 2: Top: Trapping configuration for rubidium 87 atoms: two cigar-shaped harmonic traps displaced by ∆z along the "long" trap axis. Bottom: Representation of the initial state on the Bloch sphere: |b = |F = 2, mF = 1 is the majority component and |a = |F = 1, mF = -1 is the minority component spin state. The initial state is close to the south pole (θ close to π).

  m = 1, 2, 3 and ǫ = a, b:

40 FIG. 3 :

 403 FIG. 3: Optimal cat state that we predict for realistic experimental conditions with the two rubidium 87 states |a = |F = 1, mF = -1 and |b = |F = 2, mF = 1 . Top: Fisher information (24) as a function of time, calculated with the Gross-Pitaevskii Hamiltonian (22) (green dashdotted line), and with the general two-mode model of section IV A (red solid line); for comparison, the blue dotted curve equation (26) gives the maximal Fisher information that one could obtain for the time-dependent mean atom numbers. Bottom: Wigner function at tcat = 112 ms calculated with the Gross-Pitaevskii Hamiltonian (22). The Wigner function, in the south hemisphere of the Bloch sphere, is projected onto the x -y plane. Parameters: The total atom number of average N = 150 has Poissonian fluctuations, Na = 5.71, Nb = 144.29, ω ⊥ = 2π × 1000 Hz, ωza = 2π × 850 Hz, ω zb = 2π × 50Hz, ∆z = 1.620 a ⊥ (distance between the trap centers). Scattering lengths aaa = 100.4a0, a bb = 95.44a0, a ab = 98.13a0[START_REF] Egorov | [END_REF]. One-body, two-body, and three-body loss rate constants K

40 FIG. 4 :

 404 FIG. 4: Optimal cat state that we predict for realistic experimental conditions with the two sodium 23 states |a = |F = 1, mF = 0 and |b = |F = 2, mF = -2 . Top: Fisher information (24) as a function of time, calculated with the Gross-Pitaevskii Hamiltonian (22) (green dash-dotted line), and with the general two-mode model of section IV A (red solid line); for comparison, the blue dotted curve equation (26) gives the maximal Fisher information that one could obtain for the time-dependent mean atom numbers. Bottom: Wigner function at tcat = 178 ms calculated with the Gross-Pitaevskii Hamiltonian (22). The Wigner function, in the south hemisphere of the Bloch sphere, is projected onto the x -y plane. Parameters: The total atom number of average N = 150 has Poissonian fluctuations, Na = 22, N b = 128, ω ⊥a = 2π × 1415 Hz, ω ⊥b = 2π × 115 Hz, ωza = 2π × 612 Hz, ω zb = 2π × 772 Hz, ∆z = 0.62 × 10 -6 m (distance between the trap centers). Scattering lengths aaa = 52.91a0, a bb = 64.25a0, a ab = 64.25a0 [28]. One-body loss rate constants K (1) a = K (1) b

F=1FIG. 5 :

 5 FIG. 5: Density cuts along the z-axis of the majority (blue solid line) and the minority (red dashed line) component for parameters of Fig. 3 for rubidium 87 (left), and of Fig. 4 for sodium 23 (right).

2 ; 2 ;

 22 and the fact that exp[-i θa (0 -)][cos A(t)] N exp[i θa (0 -)] = [cos A(t)] N -1 . Note -To perform the average on the vacuum in mode b, one uses the operatorial relation 0 : b| exp[iγ(â † b + âb † )]|0 : b = (cos γ) â † â where â, b are two bosonic annihilation operators with standard commutation relations, γ is a real number and the expectation value is taken in the vacuum state of b. As â † b + âb † conserves the total boson number, it suffices to prove the relation in a Fock state |na : a , that is to evaluate Sz = na/2| exp(2iγ Ŝx)|Sz = na/2 according to equation (1). Up to a rotation of angle π/2 around Oy, this is also Sx = na/2| exp(-2iγ Ŝz)|Sx = na/2 = na π 0| exp(-2iγ Ŝz)| π 2 ; 0 na . The sought relation then results from the known property exp(-2iγ Ŝz)| π ϕ na = | π 2 ; ϕ -2γ na of the phase states. -

FIG. 11 :

 11 FIG.11: Initial (red peak on the right) and the final (blue peak on the left) probability distribution of the total number of atoms. Around 20% of the atoms are lost, but practically only in the majority component b. Table:Budget of losses at time tcat . Parameters are the same as in Fig.3.

FIG. 12 :

 12 FIG.12:Example of output of the optimization program described in Appendix B for the Fisher information at the catstate time. Na on the x-axis is in the initial mean number of atoms in the minority component. For a given Na the different points correspond to successive configurations explored by the algorithm in its convergence process. Restricting to configurations with trap aspect ratios smaller than 20 we select out the darker points (in blue). Scattering lengths and loss rates are as in Fig.3for 87 Rb. Fig.3corresponds to one of the blue points with maximal Fisher information around 1500.

Fig. 3

 3 FIG.12:Example of output of the optimization program described in Appendix B for the Fisher information at the catstate time. Na on the x-axis is in the initial mean number of atoms in the minority component. For a given Na the different points correspond to successive configurations explored by the algorithm in its convergence process. Restricting to configurations with trap aspect ratios smaller than 20 we select out the darker points (in blue). Scattering lengths and loss rates are as in Fig.3for 87 Rb. Fig.3corresponds to one of the blue points with maximal Fisher information around 1500.

  n a,kj !n b,kj !

  k n k ! is the cardinality of G = {σ 1 • . . . • σ s }, one has indeed Ŝ2 partial = Ŝpartial . Also Ŝpartial commutes with Ôspin . Using expression (88) for the state wave vector |ψ(t) in a single realization and (D7), we obtainψ(t)| Ôspin |ψ(t) = ψ spin (t)| Ôspin |ψ spin (t) (D10)with the vector |ψ spin (t) defined as follows:|ψ spin (t) = 2 N/2 (n a,k ) k k P n k (n a,k ) e -i t 0 dτ E((n a,k ) k ,(n b,k ) k ,τ )/Ŝpartial |a ⊗n a,k 1 |b ⊗n b,k 1 . . . |a ⊗n a,ks |b ⊗n b,ks , (D11)

2 e

 2 |ψ spin (t) = (n a,k ) k k P n k (n a,k ) 1/-i t 0 dτ E((n a,k ) k ,(n b,k ) k ,τ )/|n a,k1 : a, n b,k1 : b ⊗ . . . ⊗ |n a,ks : a, n b,ks : b (D13) where |n a,k : a, n b,k : b are spin Fock states. We calculate the action of Ŝ+ = Ŝx + i Ŝy = N i=1 |a i i b| on equation (D13). The first n k1 terms of Ŝ+ act on the first Fock state |n a,k1 : a, n b,k1 : b and give [(1+n a,k1

  [55] At time t = 0 + , for a single realisation, one hasCov (Na-N b ,Na,nc-N b,nc ) [Var (Na-N b )Var (Na,nc-N b,nc )] 1/2 = k =0 n k N 1/2

  [58] E. Mimoun, PhD thesis, Université Pierre et Marie Curie, Paris, France (2010), https://tel.archives-ouvertes.fr/tel-00527457 [59] I. S. Gradshteyn and I. M. Ryzhik, in Tables of Inte-

  ) is a random thermal shift of the a-b condensate relative phase, already present in operatorial form in equation (72). The form (96) is obtained by summing over N a in equation (94), taking into account the fact that A(t cat ) = π/2. This is exactly the single realisation spin fidelity that one would obtain if the collective spin was in the quantum state (see endnote[56]) |ψ single spin (t cat ) = e i∆θ th Ŝz |spin cat N partition

  is the orbito-spinorial fidelity (62), d k is the degeneracy of the manifold k, and 3 F 2 is a generalised hypergeometric function, see §9.100 in reference [59]. The lower bound (100) is represented as a dashed red line in Fig. 6. Remarkably this universal function of k B T /∆ is almost indistinguishable from F Bog spin (t cat ) at the scale of the figure (see endnote [60]).

Acknowledgments

We acknowledge useful discussions with Dominique Spehner and Anna Minguzzi. Dominique Spehner made

analytical calculations about how to choose parameters to satisfy the compensation condition in the Gaussian regime, that will be exploited in a future work. During an internship at ENS in 2000, Uffe Poulsen obtained with a different technique for a one-dimensional harmonically trapped Bose gas in the Hartree-Fock limit similar results for the multimode contrast (82) (unpublished work). K. P. acknowledges support from Polish National Science Center Grant No UMO-2014/13/D/ST2/01883. M.F. and P.T. acknowledge financial support from the Swiss National Science Foundation through NCCR QSIT.

Appendix A: Fidelity and revival beyond the constant loss rate approximation

In section II we derived a simple relation ( 13) between the cat-state fidelity and the revival amplitude using the constant loss rate approximation. In this appendix we calculate the first correction to this approximation. As in section II, we restrict to the case in which the two components are symmetric and spatially separated. All analytical results are derived in the frame of the stochastic wavefunction approach, while the numerical results come from the exact diagonalization method described in Appendix A of reference [32], applied to the master equation.

The initial state is the phase state placed on the equator of a pseudo-Bloch sphere with exactly N atoms, i.e. |ψ(t = 0) = θ = π 2 ; φ = 0 N . Due to particle losses the state evolves into a mixed state,

where ρn is the unnormalized density matrix corresponding to the restriction of ρ to the subspace with exactly n atoms. The trace of the state ρn is the probability that the total number of atoms is equal to n:

In the stochastic wavefunction approach there is at time t only one stochastic wavefunction with the initial number of atoms, the one that has not experienced any quantum jump:

It means that within this subspace the (unnormalized) state remains pure, i.e. ρN = |ψ N (t) ψ N (t)|.

In the lossless case, the total number of atoms is fixed to N . Hence the time-dependent fidelity between the state in the lossless case, denoted with ψ (0) (t) , and the density matrix (A1) depends only on the state restricted FIG. 8: Relative deviations between g (1) (trev) and F 2 (tcat) in the presence of one-body losses. The approximate formula (A9) is represented as a full line, while the dashed lines are exact solutions for N = 100 and N = 300. The values of the abscissa 2N (γ (1) /χ) 2 corresponding to the trapping angular frequency ω = 2π×500 Hz and the scattering length a = 100.4 Bohr radii and one-body loss rate equal to K1 = 0.01Hz are marked as dotted vertical lines for N = 100 (left line) and N = 300 (right line). Note that as in section II, we restrict here to the case in which the two components are symmetric and spatially separated.

to the subspace with the N atoms:

We relate the fidelity to the normalized first order correlation function:

Tr Ŝx ρn = N n=0 g (1) n (t), (A5) where g

N Tr Ŝx ρn is the contribution to g (1) (t) of the subspace with n atoms.

In what follows we use the notations γ ≡ γ (1) t rev for one-body losses and γ ≡ γ (3) t rev for three-body losses, where t rev = π/χ is the first revival time.

One-body losses

We now look at corrections to the constant loss rate approximation in the presence of one-body losses. In this case there are two jump operators: Ĵa,1 = γ (1) â and Ĵb,1 = γ (1) b. The fidelity evaluated from Eq. (A4) is equal to F (t) = e -N γ (1) t (exactly as in the constant loss rate approximation).

In the case of one-body losses the full g (1) function at the time t rev can be calculated exactly:

We quantify the discrepancy between g (1) (t rev ) and (F (t cat )) 2 with the relative deviation |g (1) -F 2 |/F 2 , plotted for N = 300 (green dashed line) and N = 100 (blue dotted line) in Fig. 8.

The contribution to g (1) from the subspace with the initial number of atoms reads

Thus, if one restricts to the subspace with N atoms, the fidelity-contrast relation ( 13) becomes exact. The small discrepancy is due to the contributions g n (t) from the other subspaces n < N . The leading one is

(A8) By including this correction we obtain the approximate formula ||g (1)

N -1 (t rev )

In Fig. 8 we compare the approximate expression (A9) and the exact value of the relative correction calculated from (A6). We note that g

N -1 ≤ 8 (πe) 2 , the equality holding for γ = 2 N .

Three-body losses

Let us now consider the case of three-body losses. As the two components do not overlap, there are only two associated jump operators: Ĵa,3 = γ (3) â3 and Ĵb,3 = γ (3) b3 . From Eq. (A4) we obtain the fidelity

where

In the case of three-body losses we cannot compute analytically the first order correlation functions. Using the stochastic wavefunction approach we can however calculate the contributions to g (1) of subspaces with N and FIG. 9: Relative deviations between g (1) (trev) and F 2 (tcat) in the presence of three-body losses. The approximate formula (A21) is represented as a full line, while the symbols linked by dashed lines are exact solutions for N = 100 and N = 300. The values of the abscissa N 5 (γ (3) /χ) 2 /8 corresponding to the trap and loss parameters of Fig. 1, are marked as dotted vertical lines for N = 100 (left line) and N = 300 (right line). Note that as in section II, we restrict here to the case in which the two components are symmetric and spatially separated.

N -3 atoms:

N -3 (t rev ) = (-1)

where

In the limit of large atom numbers, the binomial distribution can be approximated with a Gaussian distribution and the sums over k with integrals:

where

. Using this continuous limit we approximate Eqs. (A10)-(A11) with

N (t rev ) = -

.

(A19)

The contribution of the subspace with N -3 atoms to g (1) (t rev ), in the limit γN 2 → 0, N → ∞ with γN 3 fixed, is g

N -3 (t rev ) ≈

Corrections to the relation (13), stating that F (t cat ) = |g (1) (t rev )| 1/2 , then come from two sources: from the difference

N | and from the difference g (1) g

(1)

N . In the limit of a small lost fraction γN 2 ≪ 1 we obtain:

The leading corrections come from (ii), as confirmed by Fig. 9, which compares the approximate analytical result

with an exact numerical calculation. We note that, both for three-body and one-body losses, up to a numerical factor 2/(mπ) 2 , the relative corrections (A9) and (A21) can be interpreted as the product between the number of lost atoms and the fraction of lost atoms. In the interesting regime in which the number of lost atoms at the revival time is smaller than one (the fidelity and the revival would be killed by the losses otherwise), the corrections are then smaller than the dominant contribution of the losses coming from the N particles subspace, by a factor 1/N .

Appendix B: Search algorithm

The numerical algorithm we use to find optimal parameters within experimental constraints, to create the cat in the case of a hyperfine transition in rubidium or sodium, is described in Fig. 10. We fix the average total number of atoms and the pulse preparing the initial phase state, and the code varies some parameters to maximize the Fisher Information (24) at the cat-state time. For example, for rubidium, the variational parameters are the radial trap frequencies ω ⊥ assumed to be equal for the two species, the longitudinal frequencies ω az and ω bz and the distance between trap centers along z, denoted ∆z. Once a favor-able configuration is found by the algorithm, we proceed to a verification step by calculating both the Fisher Information and the Wigner function beyond the χ Ŝ2 z approximation, meaning that instead of the Hamiltonian (27) in the master equation, we use (22). The Wigner function is defined as W (θ, ϕ) = ∞ N =0 p(N )W N (θ, ϕ) where W N (θ, ϕ) is normalized to unity [61]. We show a result in figure 3 for a particular configuration. For this particular configuration, corresponding to the rubidium 87 case, in Fig. 11 we show the probability distributions for the total number of atoms and we summarize the loss budget. Note that although 30 particles are lost on average in the majority component, high contrast fringes are obtained in the Wigner function. Finally in Fig. 12 we show an example of output of the optimization program, where the Fisher information at the cat-state time is maximized for different initial values of the average number of atoms in the minority component, for scattering lengths and loss rates of 87 Rb as in Fig. 3. As the initial atom number in the minority component is increased, first the optimal Fisher information increases, as one expects it from [START_REF] Willems | [END_REF] in the absence of losses, then it decreases due to the non-compensated, one-body losses in the minority component. By restricting to non extreme configurations, with ratios between the trap frequencies smaller than 20, we select out the points in blue. Fig. 3 corresponds to one of the blue points with maximal Fisher information around 1500.

Appendix C: Adiabaticity of the interaction ramp

In the multimode analysis of the cat-state production scheme at zero temperature in the box [0, L] 3 , after the π/2 pulse, one ramps the interaction strength g(t) in each spin state σ from 0 to a final positive value g f according to the Hann semi-window (52). We determine here, in the Bogoliubov approximation, the number of quasi-particles created by the ramp. Requiring that this number is ≪ 1 ensures adiabaticity of the process and its compatibility with cat-state production.

For a general time dependence of the coupling amplitude, the expansion (64) of the noncondensed field in spin state σ takes the form [46] Λσ (r, t)

where the complex Bogoliubov modal amplitudes obey the equations of motion

grals, Series, and Products, 5th ed., edited by A. Jeffrey (Academic Press, San Diego 1994).

[60] The last expression in equation (100) results from the average of [PN k (N k /2)] 2 over the N k with their probability distribution π k (N k ) = (1-q k ) d k q N k k (N k +d k -1)! N k !(d k -1)! where q k = exp(-βE k ). This results from the fact that N k is the sum of d k independent random variables with the same law n → (1 -q k )q n k , the binomial coefficient giving the number of ways to arrange N k identical objects in d k boxes.

[61] J. P. Dowling, G. S. Agarwal, W. P. Schleich, Phys. Rev. A 49, 4101 (1994).

[62] In the opposite limit 2ρσg f ≪ ∆, the adiabaticity con-dition is most stringent at time t = tramp/2, where the derivative of g(t) is maximal. One then gets the condition tramp ≫ t adiab ramp = π 8 ρσ g f ∆ 2 . For the parameters of Fig. 7, one is in the intermediate regime 2ρσg f ≃ ∆ but the adiabaticity times t adiab ramp issued from the two limiting cases agree within 30%.