Multimode analysis of non-classical correlations in double well Bose-Einstein condensates

The observation of non-classical correlations arising in interacting two to size weakly coupled Bose-Einstein condensates was recently reported by Esteve et al. [Nature 455, 1216 (2008)]. In order to observe fluctuations below the standard quantum limit, they utilized adiabatic passage to reduce the thermal noise to below that of thermal equilibrium at the minimum realizable temperature. We present a theoretical analysis that takes into account the spatial degrees of freedom of the system, allowing us to calculate the expected correlations at finite temperature in the system, and to verify the hypothesis of adiabatic passage by comparing the dynamics to the idealized model.


Introduction
Interacting quantum systems are incredibly complex, allowing a range of states that span an enormous Hilbert space. Some of these quantum states, such as squeezed and entangled states, can be utilized for quantum information [1,2] and precision measurement [3] purposes, potentially providing practical improvements over classical techniques. For instance, many common measurements based on quantum mechanical principles are usually limited to a much lower precision than allowed by the fundamental Heisenberg limit. The so-called standard quantum limit can be obtained by using the quantum analogue of classical states, such as the coherent state of light [4], to probe quantities ranging from positions and velocities, to magnetic fields and atom numbers. One can improve the precision by using entangled or squeezed states as inputs to interferometric measurement procedures [3,5].
Experimentally realizing such quantum states presents many technical challenges. Macroscopic superpositions (or Schrödinger cat states) can in principle be used to perform Heisenberg-limited measurements [5], but in practice are extremely difficult to engineer. There has been much success in creating squeezed states of light [6] and the spin of atomic samples [7] that can be used to make measurements more accurate than the standard quantum limit. Interactions between samples of photons and/or atoms can create these squeezed states, but excess technical and thermal noise place practical limitations on the amount of squeezing generated. In order for quantum information and precision measurement techniques based on squeezed states to be practical, these excess fluctuations must be reduced.
Several methods of creating relative number squeezing and entanglement in ultra-cold gases have been proposed, and some realized, such as four-wave mixing [8,9,10] and the Kerr effect [11], molecular dissociation [12,13], Josephson junction analogies [14,15] and coupling with squeezed light [16]. Recently, Estève et al. reported number difference squeezing and entanglement between multiple ultra-cold atomic clouds [15]. The experiment involved loading repulsive rubidium-87 atoms into two or more potential wells at low temperatures. The repulsion between the atoms reduces atomic bunching, and at zero temperature causes the fluctuations in the population differences between the wells to be smaller than the 'classical' binomial distribution. This experiment was the first to perform accurate atom counting within each well, while simultaneously having access to phase correlations between the wells through expansion and interference measurements. Combining these procedures, they observed relative number squeezing that could be used in principle to perform measurements at a precision 3.8 dB better than the standard quantum limit.
One of the cornerstone features of this experiment is a technique employed to reduce the level of thermal fluctuations of the quantum state. Producing multiple condensates by evaporative cooling into multiple wells produced results close to the standard quantum limit. Estève et al. found improved results by first cooling the atoms in a single trap, before modifying the potential to split the condensate between several wells. They postulated that this improvement is due to an adiabatic following of the quasiparticles representing the (quantized) fluctuations. Briefly, the energy of the excitations causing fluctuations is higher in the single condensate than after the splitting; therefore, the number of quasiparticle excitations at a given temperature is lower prior to splitting. As the potential is modified, the system is driven out of thermal equilibrium while the number of quasiparticles remains approximately fixed, and the level of fluctuations becomes lower than what would result from a thermal state at the minimum realizable temperature. In reference [15], the authors test this hypothesis by comparing with a simple, two-mode model, with reasonable agreement between the results.
Here we present a critique of the two-mode model and the hypothesis of adiabatic passage. Going beyond the two mode model reveals regimes where two modes are sufficient, or when the approximation fails. Higher-order spatial modes contribute to the physics when the wells are not well separated, if the atom number is not large, or at higher energy or temperature scales. We employ the perturbative Bogoliubov method [17] taking into account the full multi-mode structure of the system and find the two-mode description is quite accurate for well-separated clouds. We compare the assumption of adiabatic passage to a model of full rethermalization, where the system remains in thermal equilibrium as the clouds are separated while the entropy of the isolated quantum system is constant. As the entropy of the total system depends on contributions from each spatial mode, the full multi-mode statistics need to be calculated for this comparison. We find that only the adiabatic model produces the dramatic improvement in squeezing observed in the experiment. Finally, we perform truncated Wigner [18] simulations of the dynamics and find they agree well with the adiabatic model for a sufficiently large system. Thus we can conclude that thermalization between the quasiparticle modes in this system is slow, and adiabatic passage is a good model for experiments performed on these timescales.

System model
The experiment begins by evaporatively cooling a cloud of 87 Rb atoms in an optical double well potential. The double-well potential is created by superimposing an optical lattice onto Potential (a.u.) m Density (a.u.) Figure 1. A schematic of the double well potential (dashed red) and ground-state meanfield density (solid black) when the chemical potential µ is close to the barrier height. The ground state can be interpreted as a superposition of the atoms being in the left and right wells, sketched by the respective dotted-line "wave-functions". a harmonic dipole trap [19,15]. The potential in the region of the condensate is given by where the atomic mass of 87 Rb is m ≈ 1.44 × 10 −25 kg. The trap is cylindrically symmetric and elongated in the direction of the lattice, ω z < ω ⊥ . The trap has radial trapping frequency ω ⊥ ≈ 2π × 425 Hz and, for the two-well experiment, a longitudinal trapping frequency of ω z ≈ 2π × 60 Hz. The optical lattice creates a barrier between two (or more) low energy regions, or wells, such as that depicted in figure 1 and is generated by interfering two phaselocked lasers at an angle producing a lattice spacing of a ≈ 5.7 µm, where the lattice wavevector is k L = π/a. The peak-to-peak depth V L begins at 2π × 430 Hz before being ramped up. The atoms are described by the Hamiltonian, The interaction constant is U 0 = 4π 2 a s /m, where the s-wave scattering length is a s ≈ 5.39 nm.
We are interested in the low temperature limit, at or near the ground state of the system. To the lowest order of approximation, the solution is given by the ground-state mean-field wave-function ψ 0 , which we find by integrating the imaginary-time Gross-Pitaevskii equation, until the state reaches equilibrium. The chemical potential µ is adjusted to give the required atom number N 0 = |ψ 0 | 2 d 3 r = 1600. If µ is close to the barrier height, the solution is that of two weakly-coupled condensates, sketched in figure 1. The peak-to-peak depth of the optical lattice begins at 2π ×430 Hz. As the gas is cooled, a Bose-Einstein condensate forms in the trap with the potential barrier between the wells significantly smaller than the chemical potential µ -thus at this stage it can be considered a single condensate. The evaporative cooling method is unable to reduce the temperature T significantly below the level of the chemical potential [15], so k B T ∼ µ.
The potential barrier is slowly ramped up by increasing the intensity of the lasers. In the experiment, this was achieved adiabatically (i.e. slow enough that additional heating was not observed) at a rate where the peak-to-peak lattice potential V L increases by approximately 2π × 2 Hz per ms. The BEC is coherently split in two, whereupon the atoms configure themselves in a low energy state possessing relative number squeezing.
The experiment was then able to directly detect the atoms in situ, and count the atoms in each well, N i to a greater accuracy than the shot-noise limit (i.e. ∆N i < √ N i ). Defining the relative number squeezing, or normalized variance, as their results indicate squeezing of up to -6.6 dB, or a reduction in the number difference variance by a factor of 4.5 compared to the binomial distribution expected from ideal, uncorrelated condensates. To quantatively determine the number variance we must go beyond the Gross-Pitaevskii description, which does not include either quantum or thermal fluctuations in the cloud.

Bogoliubov analysis
The Bogoliubov approach is a perturbative expansion valid in the weakly-interacting or large atom number limit [17,20]. In either case, the ground state of the system is close to a coherent state given by the Gross-Pitaevskii equation. The Bogoliubov approach treats fluctuations about the mean-field perturbatively, by writingψ( and assumingδψ is 'small' compared to ψ 0 . The linearized solutions forδψ arising from equation (2) obey One can diagonalize the above linear equation (which couplesδψ withδψ † ) into their respective eigenstates. The annihilation operators for the eigenstates have the form The eigenstates with eigenvalue ǫ i are defined by the functions u i (r) and v i (r). Combining the above yields the Bogoliubov-de Gennes (BdG) equations where for brevity we have defined the Gross-Pitaevskii operator We implement the procedure presented in reference [20] to reliably solve equation (7).
To efficiently solve the BdG equations, it is important to minimize the size of our three-dimensional spatial basis. Here we implement the harmonic oscillator basis [21] as it effectively represents the trapped BEC in the lowest energy states while allowing efficient calculations of the matrix elements between the basis states φ i (r), for example can be found accurately using the Gaussian quadrature sum rules [22,23]. We further optimize the calculations by taking advantage of x, y and z reflection symmetries, allowing us to reduce the computation cost by roughly a factor of 8 2 . The total system is truncated to the 2413 states having energy less than 90 ω z .
We display the condensate density and lowest energy quasiparticle mode solutions in figure 2 for lattice depths of 430 Hz and 1650 Hz. This lowest energy mode can be   (9) and (10). The lowest energy state contributes most strongly to the number difference variance.
interpreted as transferring particles from one well to the other and corresponds to Josephsonstyle oscillations in population. For deep lattices the lowest energy mode has a significant overlap with the antisymmetrized ground state wave-function; that is u(r) and v(r) are approximately proportional to ψ 0 (r) × P (z), where P (z) is −1 for z < 0 and +1 otherwise (i.e. P (z) = Θ(z) − Θ(−z), where Θ is the Heaviside step function). In figure 3 we see that the energy of this lowest mode decreases as the lattice depth is increased, while the energy of the other modes increase. For V L 2π × 1500 Hz, the energy of the second excited mode is several times that of the lowest excitations, and at low temperatures one would expect the majority of quasiparticles in the lowest mode.

Calculating correlations
Based on the solutions to the Bogoliubov-de Gennes equations, one can calculate the variance of the atom number difference between the two wells as a function of the total atom number, temperatures and the potential. We begin by defining the regions of interest around each well, Ω 1 and Ω 2 . For the double well we choose Ω 1 to be the region with z > 0 and Ω 2 to be the region with z < 0. The quantum operator describing the total number of atoms in the ith region isN i = Ωiψ † (r)ψ(r) d 3 r. Expanding the field operator about the mean field using the solutions to the Bogoliubov-de Gennes equations,ψ(r) = ψ 0 (r) + ibi u i (r) −b † i v i (r), and inserting into the above giveŝ The system is symmetric under the transformation z → −z, and so we expect the wells to be evenly balanced, N 1 −N 2 = 0. We take an ansatz for the many-body density matrix by assuming the Bogoliubov modes are in thermal (chaotic) states, with population The variance of the number difference is therefore where we used the symmetry of the system to simplify the expression, and defined the following integrals for brevity: remembering P (z) is −1 for z < 0 and +1 otherwise. We note that some the terms in equation (9) go beyond the accuracy of the first-order Bogoliubov method used to derive them, and many authors therefore choose to disregard these terms. Here these terms are retained, without loss of accuracy, to facilitate a better comparison with the truncated Wigner simulations in section 4, where the full density matrix resulting from thermally-occupied Bogoliubov modes is implemented as initial conditions. For large N 0 and small fluctuations, the dominant contributions to the number variance come from the first sum in equation (9), and in particular the lowest energy (Josephson) mode has the largest contribution. In figure 3 (b), we see that for deep lattices |A 1 | 2 is much larger than those corresponding to other modes, and simultaneously the quasiparticle population in the lowest energy mode is expected to be greatest. Together these might justify the two mode  approximation for lattice intensities V L 2π × 1500 Hz, whereas higher energy excitations are important for weaker lattices.
This would indeed be the case for systems of larger particle number. However, in this system the additional contributions to equation (9) that a two mode model would fail to capture are significant. Even for well separated wells (V L 2π × 1500 Hz) these extra terms contribute approximately one third of the total variance at zero temperature. In this regime, we expect a small loss of precision in our calculations due to the fact that those terms proportional to B ij , C ij and D ik go beyond the accuracy of the first-order Bogoliubov analysis used, and may be poorly approximated.
In figure 4 we compare the results of the full multi-mode analysis presented here, and a simplified model including just the lowest Bogoliubov excitation. This simplified model is somewhat akin to a two-mode approximation, and we see that the two models show reasonable agreement for large barrier heights. In this regime, the two atomic clouds are well-separated and the two-mode approximation is expected to be accurate.
In figure 5 we plot the number squeezing ξ [equation (4)] as a function of barrier height and temperature. Clearly, best squeezing is obtained at lower temperatures and with larger barrier heights.

Adiabatic Passage
We now compare our results to the double-well experiment performed at Universität Heidelberg and examine the hypothesis of adiabatic passage presented in reference [15]. The experiment found improved number squeezing when slowly raising the barrier after evaporation compared to performing evaporative cooling with the barrier at fixed height. If the lattice ramp is slow compared to the energy gaps in the Bogoliubov spectrum, ǫ 1 , ǫ 2 − ǫ 1 , et cetera, one might expect the population in each quasiparticle mode to be fixed during the evolution. If this were true, then the system would no longer be in thermal equilibrium, as the ratio of the excitation energies changes during the ramp (see figure 3). The ramp ends with fewer quasiparticle excitations in the lowest mode than one would expect at realizable temperatures.
However, for sufficiently slow ramps one might expect the system to rethermalize, so that the population in each Bogoliubov changes to conform to a new global temperature as the energy levels move, and we also consider this possibility. The evolution of the potential, and therefore the quasiparticle energies, is slow compared to the gaps in the Bogoliubov spectrum, but it is unclear whether the evolution is adiabatic with respect to the full many-body Hamiltonian. The Bogoliubov description ignores interactions between the quasiparticles that can lead to a redistribution of excitations and a return to thermal equilibrium.
An isolated system undergoing quasi-static evolution will have constant entropy. The entropy of mode i in a thermal state, in terms of the mean occupation, is [24] The entropy of the total system is S = i S i . As the lattice height, and therefore excitation spectrum, varries, we calculate the new temperature required to keep the system in thermal equilibrium with fixed entropy. As many modes make significant contributions to the total entropy, and entropy is exchanged between modes, this analysis is not possible with a twomode model. We compare the two models for thermalization as the lattice potential is increased in figure 6. In all cases we begin in thermal equilibrium with a lattice strength of 430 Hz at a temperature of k B T ≈ µ/6, chosen to roughly correspond to the experimentally observed fluctuations.
As can be seen from figure 6, the adiabatic process predicts a level of fluctuations significantly lower than the fully thermalizing model. In a realistic experiment one might expect the entropy to increase slightly, resulting in even greater fluctuations. As the experiment observes strongly improved number difference squeezing as the lattice is ramped, one can rule out the possibility of rethermalization on the timescales of the experiment. On the other hand, the assumption of adiabatic passage agrees well with the experimental results, and we conclude that the mechanism of effective cooling proposed in [15] is correct.

Dynamical simulation
In this section we perform a direct dynamical simulation using the truncated Wigner [18,25,10,21] method and compare the results to perfect adiabatic passage. As pointed out in the previous section, the number fluctuations of the double-well system realized by [15] are sensitive to second-order interactions between the linearized Bogoliubov modes. For the experimental parameters, the truncated Wigner simulations show dynamics that is inconsistent with the linearized Bogoliubov-de Gennes solutions -in particular, the quasiparticle populations and number difference variance are not constant in time even for a fixed Hamiltonian. The Bogoliubov fluctuations are too large to be considered perturbative, and the linearized approximation breaks down.
Nevertheless, we can probe adiabatic passage dynamically by investigating a slightly different system where the perturbationsδψ are small compared to ψ 0 , and linearization is accurate. To minimize the changes to the system we investigate a system of more atoms with smaller interactions, such that the product N 0 U 0 is unchanged. In principle, such a change could be realized with a Feshbach resonance [17], although this would be challenging. The truncated Wigner method is an approximate phase space method that takes advantage of the Wigner representation of a quantum field. Truncating the higher-order derivatives from the equation of motion for the Wigner function results in a classical Liouville equation that can be efficiently sampled stochastically. Thus we reduced the problem of solving a nonlinear equation for a quantum fieldψ(r) to that of the non-linear motion of an ensemble of classical fields ψ(r). The resulting trajectories obey the (real-time) Gross-Pitaevskii equation, where the initial conditions are sampled from the initial Wigner distribution. Expectation values of symmetrically ordered operators are equal to the ensemble average of their classical counterparts; however care must be taken when using a finite, non-uniform basis [21].
The system we simulate has 1.6 × 10 5 atoms and a scattering length a s = 5.39 × 10 −11 m, and begins in thermal equilibrium at V L = 2π × 430 Hz with a temperature T ≈ µ/6. This system has 100 times more atoms than the experiment, but identical mean-field and Bogoliubov-de Gennes solutions. The initial state is populated by a coherent state condensate surrounded by Bogoliubov modes in thermal (chaotic) states with populations given by the Bose-Einstein statistics [17]. The Wigner distribution for such a state is Gaussian and therefore straightforward to randomly sample. We evolve 4000 independent trajectories with the Gross-Pitaevskii equation, while the lattice is ramped up at a rate 2π × 2 Hz/ms. As we use the non-uniform Hermite-Gauss basis for our simulation, care must be taken to extract the correct quantum expectation values from the averaged data. For arbitrary basis {φ i (r)} such that ψ(r) = c i φ i (r) and matrix P defined by we find the number difference and the number difference squared These relations allow us to efficiently calculate the number difference variance without transforming to the spatial basis, while simplifying the symmetric Wigner corrections for a non-uniform basis [21]. Our dynamical results are plotted in figure 7, along with the predictions for perfect adiabatic following and full thermalization obtained by the same procedure as in the previous section. We see that the number difference variance agrees with adiabatic model, validating the hypothesis of adiabatic following on this timescale.

Systems with different geometry
It is worth discussing changes that occur in multi-well systems of different geometry. First, double well traps can be created in a variety of ways, such as the fully optical approach analysed here [15] or with magnetic traps on atom chips [26]. The latter setup results in two parallel, elongated condensates (with ω x ≪ ω y and ω z ). In such a system, one would expect long wavelength excitations along the transverse dimension x, some of which will have energy less than the Josephson mode. Similarly, one expects multiple Josephson-type modes within a small energy band, coupling the typical double-well excitation with low lying transverse modes (the limit of two infinite condensates can be found in [27]). In this case, a multimode description such as is presented here is essential to analysing these systems. It is unclear, however, if such a geometry is advantageous for creating sub-shot noise number correlations. The increased density of states about the Josephson mode may make adiabatic following more difficult. Furthermore, a clear definition of the relative phase of the condensates, as investigated in reference [15], may be problematic as the phase may vary along the long axis of the system.
In reference [15] the authors also investigated systems of multiple wells. They presented improved results for the inner pair of six coupled wells, which was theoretically investigated using a six-mode Bogoliubov theory. We performed a multi-mode analysis of this system, which involved an increased computational difficulty. In the multi-mode picture, we observed a band of five Josephson-type states of decreasing energy as the barrier height is increased -one mode for each inter-well coupling -while the remaining excitations had increasing energy. Thus, we can conclude that the adiabatic following technique should also be benificial in this arrangement, as was indeed verified experimentally [15]. We were unfortunately unable to repeat the full analysis for that system because of the increased numberical difficulty.

Conclusion
We have performed a multi-mode Bogoliubov analysis of coupled Bose-Einstein condensate systems to find the low temperature physics and calculate the population fluctuations between the two wells. At very low temperatures, repulsive interactions induce relative number squeezing and even entanglement between the wells. We present a critical analysis of the adiabatic passage hypothesis employed in reference [15] as the optical lattice separating the wells is ramped up. Our results indicate that adiabatic evolution of the quasiparticles is possible, causing a decrease in number fluctuations as the lattice height is increased. We conclude that the isolated system can be driven out of thermal equilibrium, to reduce fluctuations to below what is possible at thermal equilibrium. It remains to be seen if this procedure can be improved to reach the quantum regime, where the vacuum fluctuations dominate the physics.
There are several limitations to the Bogoliubov approximation; we are unable to investigate either the low atom number or high temperature limits. The transition between anti-bunching due to repulsion and boson bunching at higher temperatures could be investigated by non-perturbative methods. Quantum Monte Carlo techniques [28] provide the most accurate and reliable results for high temperature bosonic systems. The Positive-P phase space method can be implemented to find thermal statistics of a Bose gas [29]. The Popov approximation [17] is similar to the approach used here, but considers the interactions between the excited modes and back-action on the condensate, and may be an appropriate method for exploring this regime. Related to these, Sinatra et al. [30] have proposed a method for sampling the statistics of a Bose gas described by a Gaussian density matrix.