Atomic four-wave mixing via condensate collisions

A Perrin¹, C M Savage², D Boiron¹, V Krachmalnicoff¹, C I Westbrook¹ and K V Kheruntsyan^3,4

¹ Laboratoire Charles Fabry de l'Institut d'Optique, CNRS, Univ Paris-Sud, Campus Polytechnique, RD128, 91127 Palaiseau Cedex, France
² ARC Centre of Excellence for Quantum-Atom Optics, Department of Physics, Australian National University, Canberra ACT 0200, Australia
³ ARC Centre of Excellence for Quantum-Atom Optics, School of Physical Sciences, University of Queensland, Brisbane, QLD 4072, Australia

⁴ Author to whom any correspondence should be addressed.

E-mail: christoph.westbrook@institutoptique.fr and karen.kheruntsyan@uq.edu.au

Received 11 December 2007
Published 30 April 2008

Contents

• 1. Introduction
• 2. Summary of experimental results
  • 2.1. Overview of the experiment
  • 2.2. Main results
• 3. Qualitative analysis
  • 3.1. Width of the pair correlation functions
  • 3.2. Width of the scattered halo
• 4. Model
• 5. Results and discussion
  • 5.1. Main numerical example
  • 5.2. Shorter collision time
  • 5.3. Smaller collision velocity
  • 5.4. Smaller scattering length
  • 5.5. Relative number squeezing and violation of Cauchy-Schwartz inequality
• 6. Summary
• Acknowledgments
• Appendix A. Duration of the collision
• Appendix B. Occupation number of the scattering modes and amplitude of the BB correlation
• Appendix C. Width of the s-wave scattering sphere in the undepleted `pump' approximation
• Appendix D. Positive-P simulation parameters
• References

1. Introduction

Recent years have seen the introduction of powerful new tools for studying degenerate quantum gases. For example, on the experimental side correlation measurements offer a new experimental probe of many-body effects [1]–[11]. On the theoretical side, the challenges posed by the new experimental techniques are being met by quantum dynamical simulations of large numbers of interacting particles in realistic parameter regimes. These are becoming possible due to the advances in computational power and improvements in numerical algorithms (for recent examples, see [12]–[15]).

In this paper, we study metastable helium (⁴He^*), which is currently unique in quantum atom optics in that it permits a comparison of experimentally measured [16] and theoretically calculated quantum correlations. This is one of the first examples in which experimental measurements can be considered in the context of first-principles calculations. Our goal in this paper is to confront a theoretical analysis with the results of recent experiments on atomic four-wave mixing via a collision of two Bose–Einstein condensates (BECs) of metastable ⁴He* atoms [16]. Figure 1 is a schematic momentum space diagram of these experiments. Two condensates, whose atoms have approximately equal but opposite momenta, k₁ and k₂≊–k₁, interact by four-wave mixing, while they spatially overlap, to produce correlated atomic pairs with approximately equal but opposite momenta, k₃ and k₄, satisfying momentum conservation, k₁ + k₂ = k₃ + k₄. Figure 1 corresponds to the experimental data shown in figure 2 of [16], since after time-of-flight expansion, atomic momentum is mapped into atomic position.

We perform first-principles quantum simulations of the collision dynamics using the positive P-representation method [17]–[20]. The advantage of this method is that given the Hamiltonian of the interacting many-body system, no additional approximations are imposed to simulate the quantum dynamics governed by the Hamiltonian. The drawback on the other hand, is that it usually suffers from large sampling errors and the boundary term problem [21] as the simulation time t_sim increases, eventually leading to diverging results.

An empirically estimated upper bound for the positive-P simulation time (beyond which the stochastic trajectories start to make large excursions in phase space, leading to boundary term problem and uncontrollable sampling errors [21]) for the evolution of condensates with s-wave scattering interactions is given approximately by [22]

where m is the atom mass, a is the s-wave scattering length, ρ₀ is the condensate peak density, and ΔV = ΔxΔyΔz is the volume of the elementary cell of the computational lattice, with lattice spacings of Δx, Δy and Δz. Applying this formula to metastable helium, we see that this is a particularly challenging case among commonly condensed atoms due to its small atomic mass and relatively large scattering length. Our simulations are restricted to short interaction times (typically   25 μs), which are about six times shorter than the experimental interaction time of [16]. Despite this, our positive-P simulations provide useful insights into the experimental observations; in addition, they can serve as benchmarks for approximate theoretical methods (such as the Hartree–Fock–Bogoliubov method [23]–[26] or the truncated Wigner method [12, 13]) to establish the range of their validity.

We note that the simulations of BEC collisions of [12, 13] using the truncated Wigner method are in a different regime than the ones carried out here using the positive-P method. More specifically, the authors of [12, 13] simulate condensates at much higher densities, for which the approximations of the method [13, 27] are satisfied. The advantage of the Wigner method over the positive-P is that it does not suffer from boundary term problem and can be used to simulate condensate collisions for longer durations.

In the present paper, we calculate atom pair correlations within the scattering halo produced spontaneously during the collision (see figure 1). The scattering halo is a spherical shell in momentum space. In the limit of small occupation of the scattered modes, the s-wave nature of the collisions ensures an approximately uniform atom population over the halo. We consider the strength and the width of the correlation signal, as well as the momentum width of the halo. We also analyze relative atom number squeezing and the violation of the classical Cauchy–Schwartz inequality.

In section 2 of this paper, we will summarize the experimental results we wish to analyze. In section 3, we discuss order of magnitude estimates. In section 4, we describe simulations using the positive P-representation method, and in section 5, we discuss the results of our simulations. Section 6 summarizes our work.

2. Summary of experimental results

2.1. Overview of the experiment

The starting point of the experiment is a ⁴He* condensate of 10⁴–10⁵ atoms confined in a magnetic trap whose frequencies are: ω_x/2π = 47 Hz and ω_y/2π = ω_z/2π = 1150 Hz. A sudden Raman outcoupling drives the trapped ⁴He* from the m_x = 1 Zeeman sublevel into the magnetic field insensitive state m_x = 0 [16]. The Raman transition also splits the initial (m_x = 1) condensate into two roughly equally populated condensates with opposite velocities along the x-direction. The magnitude of each velocity is equal to the recoil velocity v_r = 9.2 cm s^–1, defined by the momentum of the photons used to create the colliding condensates , k_r = 5.8×10⁶ m^–1. The relative velocity 2v_r of the two condensates is about eight times higher than the speed of sound of the initial condensate, ensuring that the relevant elementary excitations of the condensates correspond to free particles.

During the separation of the condensates, elastic collisions occurring between atoms with opposite velocities scatter a small fraction (5%) of the total initial atom number into the halo. The system is shown in three dimensions (3D) in an accompanying video of the experimental results after a 320 ms time of flight^Note5 . For the purposes of this paper, the experiment consists in the acquisition of the 3D positions of the particles scattered into the collision halo after the time of flight. This information permits the reconstruction of the 3D momentum vectors of the individual particles after they have ceased interacting with each other.

2.2. Main results

Knowledge of the momentum vectors in turn permits the construction of two-particle correlation functions in momentum space. The correlation function shows features for particles traveling both back-to-back (BB) and collinearly. The BB correlation results from binary, elastic collisions between atoms, whereas the collinear (CL) correlation is a two particle interference effect involving members of two different pairs: a Hanbury Brown–Twiss correlation [28]. Both correlation functions are anisotropic because of the anisotropy of the initial colliding condensates.

To quantify these correlations, we first introduce the unnormalized normally-ordered second-order correlation function between the densities at two points in momentum space,

Here, is the momentum density operator, and are the Fourier transforms of the field creation and annihilation operators and , and the colons :: stand for normal ordering of the operators according to which all creation operators stand to the left of the annihilation operators, so that

Because of a low data rate, the correlation measurements must be averaged over the entire collision sphere to get statistically significant results. The average CL second-order correlation as a function of the relative displacement Δk_i in the k_i-direction (i = x,y, z) is defined as

where e_i is the unit vector in the k_i-direction. The normalization of g_CL⁽²⁾(Δk_i) ensures that for uncorrelated densities g_CL⁽²⁾(Δk_i) = 1. The integration domain in (4) selects a certain region of interest in momentum space and can be defined, for example, to contain only a narrow spherical shell and to eliminate the initial colliding condensates. Due to the averaging, the dependence of the correlation functions on the direction k is lost.

The average BB correlation function g_BB⁽²⁾(Δk_i) between two diametrically opposite points, one of which is additionally displaced by Δk_i in the k_i-direction, is defined similarly to g_CL⁽²⁾(Δk_i):

The experimental observations can be summarized as follows. The width of both correlation functions along the axial direction of the condensate, the x-axis, is limited by the resolution of the detector and hence contains little information about the collision. In the radial direction (with respect to the symmetry of the colliding condensates), one observes a peak which can be fitted to a Gaussian function with rms widths σ_{y, z}^CL and σ_{y, z}^BB for the CL and BB cases, respectively. The experimental results are summarized in the following table

One can also use the data to deduce the averaged radial width δk of the scattering halo. Figure 2 shows a cross-section of the halo, averaged over all accessible scattering angles. The presence of the unscattered condensates prevents observation of the shell along the x-axis, but along the accessible directions we find δk≊0.067k_r.

3. Qualitative analysis

In this section, we discuss some simple, analytical estimates of the measured quantities. In later sections, we will do more precise, numerical calculations which will verify the results of this section.

3.1. Width of the pair correlation functions

As discussed in [16], the width of the BB and CL correlation functions should be on the order of the momentum width of the initial condensate, which in turn is proportional to the inverse width of its spatial profile. For a Gaussian density profile of the initial condensate in position space , and therefore a Gaussian density distribution in momentum space, , with σ_i = 1/2w_i, an approximate theoretical treatment based on a simple ansatz for the pair wavefunction predicts a Gaussian dependence of the BB and CL correlation functions on the relative wavevectors Δk_i [28]:

The widths of the BB (σ_i^BB) and CL (σ_i^CL) correlations are related to the momentum-space width σ_i of the initial (source) condensate via [28]

and therefore the width of the BB correlation is times smaller than the width of the CL correlation. Similar predictions of correlation widths have been made and discussed in [13, 24].

In section 5.1, the initial momentum-space widths are found to be σ_x = 0.0025k_r and σ_{y, z} = 0.055k_r, assuming N = 9.84×10⁴ atoms. Expressing the experimentally measured widths in units of σ_{y, z}, we can rewrite (6) as

and therefore, equation (9) is in agreement with the measured width of the radial BB correlation function, whereas (10) overestimates the width of the CL correlation function by almost 60%. As we show below, first-principles simulations using the positive-P method and incorporating atom–atom interactions result in widths which are closer to the experimental values.

The discrepancy between the two theoretical approaches (which apparently is larger for the CL correlations than for the BB ones) comes mostly from the fact that the above calculation is made for a Gaussian shape of the initial BEC density profile, whereas in practice and in the positive-P simulations the spatial density of a harmonically trapped condensate is closer to an inverted parabola (as in the Thomas–Fermi limit) rather than to a Gaussian. An alternative theoretical model [29], based on the undepleted source condensate approximation and a numerical solution to the linear operator equations of motion for scattered atoms, also confirms that for short times the momentum-space correlation widths are narrower if the source condensate has a parabolic spatial density profile, compared to the case of a Gaussian density profile.

3.2. Width of the scattered halo

A second, experimentally accessible quantity in a BEC collision is the width in momentum space of the halo on which the scattered atoms are found. Clearly the momentum spread σ_i (in i = x-, y- or z-direction) of the colliding condensates imposes a minimum width

This limit suggests that the halo could be anisotropic. As noted above, however, the experiment in [16] is not highly sensitive to such an anisotropy, and measures the width chiefly in the y- and z-directions.

Other physical considerations also affect this width, and suggest that the halo should rather be isotropic, in which case we can drop the index from δk. Here, we discuss two mechanisms that impose a finite radial width on the halo.

If the pairs are produced during a finite time interval Δt, the total energy of the pair is necessarily broadened by /Δt. This is true even if the relative momentum is well defined. For a mean k-vector k_r, the finite interaction time between the colliding BECs results in a broadening of

where we assumed δk/k_r  1. In the experiment, the collision time is sufficiently long that the above effect does not impose a limitation on the width of the sphere. In the positive-P simulations, however, numerical stability problems limit the maximum collision time that can be simulated, as discussed in section 5, and this time does indeed impose a width on the halo. For short collision times, where the scattering is in the spontaneous regime, our numerical results for the width δk are in good agreement with the simple estimate of equation (13).

For long collision times, it can happen that so many atoms are scattered that Bose enhancement and stimulated effects become important. In this case, the width of the scattering shell can be estimated by a slightly more involved approximate approach based on analytic solutions for the uniform system within the undepleted `pump' (source condensate) approximation [30]. Under this approximation, the present system is equivalent to the dissociation of a condensate of molecular dimers studied in detail in [14, 31, 32]. The latter system in turn is analogous to parametric down-conversion in optics [33]. The details of the approximate solutions, common to condensate collisions and molecular dissociation, and the relationship between them are given in appendix C. The resulting width of the halo found from this approach is

We see that in this regime, the width is proportional to the scattering length a and the peak density ρ₀, but it no longer depends on the collision duration.

The physical interpretation of equation (14) is that with the stronger effective coupling (or nonlinearity) aρ₀, one can excite and amplify spectral components that are further detuned from the exact resonance condition (or further `phase mismatched'). The inverse dependence on collision momentum k_r can be understood via the quadratic dependence of the energy on momentum k: to get the same excitation at a given energy offset , (C.3), one requires smaller absolute momentum offset δk at larger k_r than at small k_r.

Positive-P simulations covering the transition from the spontaneous to stimulated regimes are available for ²³Na condensate collisions as in [15]. The numerical results in this case are in agreement with the simple analytic estimate of equation (14). More specifically, we find that for collision durations between 300 and 640 μs the actual numerical results for the width of the spherical halo vary, respectively, between δk/k_r≊0.13 and δk/k_r≊0.087, whereas equation (14) predicts δk/k_r≊0.096.

For ⁴He*, on the other hand, the small mass and the larger scattering length of ⁴He* atoms limit the maximum simulation time to t_sim  25 μs. This is far from the stimulated regime, and therefore we do not have a direct comparison of the numerical results with equation (14). The experiment is also not in the stimulated regime. We are nevertheless tempted by the numerical ²³Na result to extrapolate equation (14) to ⁴He* BEC collisions in the long time limit and we obtain δk/k_r≊0.05. Adding this width in quadrature to the momentum width of the initial condensate, σ_{y, z}≊0.055k_r, gives , not far from the experimentally observed radial momentum width of δk≊0.067k_r. We thus suggest that the mechanism leading to equation (14) may play a role in the experiment.

4. Model

The effective field theory Hamiltonian governing the dynamics of the collision of BECsis given by

where is the field operator with the usual commutation relation is the atomic mass, the first term is the kinetic energy, and the second term describes the s-wave scattering interactions between the atoms. The trapping potential for preparing the initial condensate before the collision is omitted since we are only modeling the dynamics of the outcoupled condensates in free space. The use of the effective delta function interaction potential U(x–y) = U₀δ(x–y) assumes a UV momentum cutoff k^max. In our numerical simulations, the momentum cutoff is imposed explicitly via the finite computational lattice. If the lattice spacings (Δx, Δy and Δz) in each spatial dimension are chosen to be much larger than the s-wave scattering length a, then the respective momentum cutoffs satisfy k_{x, y, z}^max  1/a. In this case the coupling constant U₀ is given by the familiar expression  [34] without the need for explicit renormalization.

To model the dynamics of quantum fields describing the collision of two BECs, we use the positive P-representation approach [17]. In this approach, the quantum field operators and are represented by two complex stochastic c-number fields Ψ(x,t) and whose dynamics is governed by the following stochastic differential equations [15]:

Here, ζ₁(x, t) and ζ₂(x, t) are real independent noise sources with zero mean, ζ_j(x,t) = 0, and the following non-zero correlation:

The stochastic fields Ψ(x,t) and are independent of each other ] except in the mean, , where the brackets  ... refer to stochastic averages with respect to the positive P-distribution function. In numerical realizations, this is represented by an ensemble average over a large number of stochastic realizations (trajectories). Observables described by quantum mechanical ensemble averages over normally ordered operator products have an exact correspondence with stochastic averages over the fields Ψ(x, t) and :

The initial condition for our simulations is a coherent state of a trapped condensate, modulated with a standing wave that imparts initial momenta ±k_r (where and v_r is the collision velocity) in the x-direction,

with . Here, ρ₀(x) is the density profile given by the ground state solution to the Gross–Pitaevskii equation in imaginary time. The above initial condition models a sudden Raman outcoupling of a BEC of trapped ⁴He* atoms in the m_x = 1 sublevel into the magnetic field insensitive state m_x = 0, using two horizontally counter-propagating lasers and a third vertical laser [16]. In this geometry, the Raman transitions split the initial (m_x = 1) condensate into two equally populated condensates and simultaneously impart velocities of ±v_r onto the two halves. As a result, the two outcoupled condensates undergo a collision and expand in free space. Accordingly, in our dynamical simulations, the field represents the atoms in the untrapped state m_x = 0, having the s-wave scattering length of a₀₀ = 5.3 nm ([16] and references therein), while the initial density profile ρ₀(x) refers to that of the trapped atoms in the m_x = 1 state having the scattering length of a₁₁ = 7.51 nm [35]. The same distinction in terms of the scattering length in question applies to the definition of the interaction strength , in which a has to be understood as a₁₁ for the trapped condensate or as a₀₀ for the outcoupled cloud.

In our simulations, we assume for simplicity that the outcoupling from the trapped m_x = 1 state is 100% efficient, in which case the entire population is transferred into the m_x = 0 state and therefore we have only to model s-wave scattering interactions between the atoms in the m_x = 0 state. In the experiment, on the other hand, the transfer efficiency is only about 60% and therefore the collisions between the atoms in the m_x = 0 and m_x = 1 are not completely negligible and may be responsible for some of the deviations between the present theoretical results and the experimental observations.

5. Results and discussion

5.1. Main numerical example

Here, we present the results of positive-P numerical simulations of collisions of two condensates of ⁴He* atoms (m≊6.65×10^–27 kg) as in the experiment of [16]. The key parameters in our main numerical example are the collision velocity, v_r = 9.2 cm s^–1, and the peak density of the initial trapped condensate, ρ₀ = 2.5×10¹⁹ m^–3. The trap frequencies are matched exactly with the experimental values, ω_x/2π = 47 Hz and ω_y/2π = ω_z/2π = 1150 Hz. The s-wave scattering length for the magnetically trapped atoms in the m_x = 1 sublevel is a₁₁ = 7.5 nm; the s-wave scattering length for the outcoupled atoms in the m_x = 0 sublevel is a₀₀ = 5.3 nm. Other simulation parameters are given in appendix D.

The initial state of the trapped condensate is found via the solution of the Gross–Pitaevskii equation in imaginary time. Given the above trap frequencies and the peak density as a target, we find that the total number of atoms in the main example is N = 9.84×10⁴. With these parameters, the average kinetic energy of colliding atoms is E_kin/k_B = mv_r²/2k_B≊2.0×10^–6 K, which is about 7.4 times larger than the mean-field energy of the initial condensate  K.

The duration of simulation in the main example is t_f = 25 μs. This is considerably smaller than the estimated duration of collision in the experiment, 140 μs (see appendix A). The number of scattered atoms in our numerically simulated example at t_f = 25 μs is ~1750, representing ~1.8% of the total number of atoms in the initial BEC. Operationally, the fraction of scattered atoms is determined as the total number of atoms contained within the scattering halo (see figure 3 showing two orthogonal slices through the momentum density distribution) after eliminating the regions of momentum space occupied by the two colliding condensates. We implement the elimination by simply discarding the data points corresponding to |k_x| > 0.99k_r, which fully contain the colliding condensates. This cuts off a small fraction of the scattered atoms as well, but the procedure is simple to implement operationally and is unambiguous.

In order to compare our calculated fraction of scattered atoms at t_f = 25 μs with the experimentally measured fraction of 5% at the end of collision at ~140 μs, we first note that these timescales are relatively short and correspond to the regime of spontaneous scattering. The number of scattered atoms increases approximately linearly with time, therefore our calculated fraction of 1.8% can be extrapolated to about 10% to correspond to the expected fraction at ~140μs. Next, one has to scale this value by a factor 0.6² to account for the fact that in the experiment only 60% of the initial number of atoms was transferred to the m_x = 0 state of the colliding condensates. Accordingly, our theoretical estimate of 10% should be proportionally scaled down to 4% conversion, in good agreement with the experimentally estimated fraction of 5% (see also appendix A).

In figure 4, we plot the radial momentum distribution of scattered atoms (solid line), obtained after angle averaging of the full 3D distribution within the region |k_x|≤0.8k_r. The numerical result is fitted with a Gaussian exp[–(k–k₀)²/(2δk²)] (dashed line), centered at k₀ = 0.98k_r and having the radial width of δk = 0.10k_r≊5.8×10⁵ m^–1, where k = |k|. The fitted radial width of δk = 0.10k_r of the numerical simulation is in reasonable agreement with the simple estimate of equation (13), which gives δk≊0.075k_r for Δt = 25 μs.

Figure 5 shows the numerical results for the BB and CL correlations (solid lines with circles), defined in equations (4) and (5). Due to the symmetry of the y- and z-directions, the results in these directions are practically the same. In order to verify the hypothesis that the shape and therefore the width of the pair correlation functions is governed by the width of the momentum distribution of the source condensate, we also plot the actual initial momentum distributions of the source condensate in the two orthogonal directions (with the understanding that the horizontal axis Δk_i now refers to the actual wavevector component k_i). The actual data points for the correlation functions and for the momentum distribution of the source are shown by the circles and squares, respectively, and are fitted with Gaussian curves for simplicity and to guide the eye. The Gaussian fits for the correlation functions (solid lines) give:

where the correlation widths σ_i^BB and σ_i^CL are shown in the table (22) below. The Gaussian fits (dashed lines) for the slices of the initial momentum distribution n₀(k_i)exp {–k_i²/[2(σ_i)²]} are scaled to the same peak value as g_BB/CL⁽²⁾(0)–1 and have σ_x = 0.0025k_r and σ_{y, z} = 0.055k_r.

By comparing the solid and the dashed lines, we see that the shape of the correlation functions indeed closely follow the shape of the momentum distribution of the source. More specifically, we find that the following results provide the best fit to our numerical data:

The ratios between the CL and BB correlation widths are σ_x^CL/σ_x^BB≊1.08 and σ_{y, z}^CL/σ_{y, z}^BB≊1.13. The errors due to stochastic sampling on all quoted values of the correlation widths are smaller than 3%.

The values for σ_{y, z}^CL/σ_{y, z} and σ_{y, z}^BB/σ_{y, z} can be compared with the respective experimentally measured values of table (11) and we see reasonably good agreement, even though the numerical data are for a much shorter collision time. The remaining discrepancy between the numerical data at t_f = 25 μs and the experimentally measured values after a ~140 μs interaction time may be due to the evolution of the condensates past 25 μs, not attainable within the positive-P method. The above numerical results for the correlation widths can also be compared with the simple analytic estimate based on the Gaussian ansatz treatment of equations (9) and (10). We find that the approximate analytic results overestimate the BB and CL widths by ~20 and 40%, respectively, in the present example.

The amplitude of the correlation functions can also be inferred by simple models. In fact, the CL correlation function is a manifestation of the Hanbury Brown and Twiss effect since it involves pairs from two independent spontaneous scattering events and we expect an amplitude of g⁽²⁾_CL(0) = 2  [28]. This is in agreement with the positive-P simulations. The BB correlation amplitude, on the other hand, can be substantially higher and display super-bunching (g_BB⁽²⁾(0)  1) [14, 24] since the origin of this correlation is a simultaneous creation of a pair of particles in a single scattering event.

In a simple qualitative model [16], the amplitude of the BB correlation can be linked to the inverse population of the atomic modes on the halo. As we show in appendix B, this model follows the trends observed in our first-principles numerical simulations.

5.2. Shorter collision time

Here, we present the results of numerical simulation for the same parameters as in our main numerical example from section 5.1, except that the data are analyzed at t_f = 12.5 μs, which is half the previous interaction time. We found in section 5.1 that σ_yz^BB, σ_yz^CL and the width of the halo δk are all nearly the same. In section 3, however, we argue that the widths of the correlation functions and the halo are governed by different limits (equations (9), (10) and (13) or (14), respectively). The example in this section illustrates this point.

Figure 6 shows two orthogonal slices of the s-wave scattering sphere in momentum space (cf figure 3), whereas figure 7 is the corresponding radial distribution after angle averaging. The most obvious feature of the distribution is that it is broader than at t_f = 25 μs and the fitted Gaussian gives the radial width of δk = 0.20k_r. This is precisely twice the width in figure 4 and is in agreement with the simple qualitative estimate of equation (13).

The BB and CL correlation functions after t_f = 12.5 μs collision time are qualitatively very similar to those shown in figure 5, except that the Gaussian fits are

with the correlation widths given by

The ratios between the widths are σ_x^CL/σ_x^BB≊1.09 and σ_{y, z}^CL/σ_{y, z}^BB≊1.16.

For the correlation functions, the main difference compared to the case for 25 μs is that the peak value of the BB correlation is now larger, reflecting the lower atomic density on the scattering halo. The correlation widths, on the other hand, are practically unchanged, at least within the numerical sampling errors of the positive-P simulations; the errors are at the level of the third significant digit in the quoted values, which we suppress. The number of scattered atoms in this example is about 850, which is approximately half the number at 25 μs, confirming the approximately linear dependence on time in the spontaneous scattering regime.

5.3. Smaller collision velocity

In this example, we present the results of simulations in which the collision velocity is smaller by a factor than before, v_r ' = 6.5 cm s^–1 (k_r ' = 4.1×10⁶ m^–1), while all other parameters are unchanged. In practice, this can be achieved by changing the propagation directions of the Raman lasers that outcouple the atoms from the trapped state. As in the previous example, the halo width illustrates equation (13).

The results of positive-P simulations for the momentum density distribution at t_f = 25 μs are shown in figures 8 and 9. The most obvious feature of the distribution is again the fact that it is now broader than in our main example of section 5.1. The width of the Gaussian function fitted to the numerically calculated radial momentum distribution is given by δk = 0.21k_r '. This is again in excellent agreement with the simple analytic estimate of equation (13), which predicts the broadening to be inversely proportional to the collision velocity. We also note that the peak momentum (relative to k_r ') in the present example is slightly shifted towards the center of the halo, k₀ = 0.92k_r ', which is a feature predicted in [30] to occur when the ratio of the kinetic energy to the interaction energy per particle is reduced.

Figure 8. Same as in figure 3, except for times smaller collision velocity, vr ' = 6.46 cm s-1 (kr ' = 4.09×106 m- 1). The axes for the momentum components ki (i = x, y and z) are in units of smaller recoil momentum than in figure 3, and therefore the absolute radius of the s-wave scattering sphere is smaller in the present example.

Figure 9. Same as in figure 4 except for times smaller collision velocity vr ' (kr ' = 4.1×106 m- 1). The width and the peak of the fitted Gaussian are δk = 0.21kr ' = 8.6×105 m- 1 and k0 = 0.92kr '.

The BB and CL correlation functions in this example are again qualitatively very similar to those shown in figure 5, except that the Gaussian fits are

with the correlation widths given by

where σ_x/k_r '≊0.0035 and σ_x/k_r '≊0.078. The ratios between the CL and BB correlation widths are σ_x^CL/σ_x^BB≊1.13 and σ_{y, z}^CL/σ_{y, z}^BB≊1.12.

As we see from these results, the absolute widths of the correlation functions are practically unchanged compared to the main numerical example (22). This provides further evidence that, at least for short collision times, the correlation widths are governed by the momentum width of the source condensate, which is unchanged in the present example compared to the case of section 5.1.

The number of scattered atoms in this example is about 1270, which is approximately times smaller than in section 5.1 and corresponds to ~1.3% conversion. This scaling is in agreement with the rate equation approach [24], according to which the number of scattered atoms is proportional to the square root of the collision energy and hence to the collision momentum, which is times smaller here.

5.4. Smaller scattering length

Finally, we present the results of numerical simulations for the same parameters as in our main numerical example from section 5.1 , except that the scattering lengths a₁₁ and a₀₀ are artificially halved, i.e. a₀₀ = 2.65 nm and a₁₁ = 3.75 nm. The trap frequencies are unchanged and we modify the chemical potential to arrive at the same peak density of the initial BEC in the trap, ρ₀≊2.5×10¹⁹ m^–3. The total number of atoms is now smaller, N≊3.5×10⁴. One effect of changing the scattering length is that it changes the size and shape of the trapped cloud, and therefore also its momentum distribution. The shape is slightly closer to a Gaussian and therefore also to the treatment in [28].

Due to the smaller scattering length, the density distribution in position space of the initial trapped condensate is now narrower and conversely the momentum distribution of the colliding condensates is broader. On the other hand, the width of the halo (see figures 10 and 11 at t_f = 25 μs) of scattered atoms is practically unchanged compared to the example of figure 4, as it is governed by the energy–time uncertainty consideration (13), for the spontaneous scattering regime. The only quantitative difference is the lower peak density on the scattering sphere, which is due to the weaker strength of atom–atom interactions resulting in a slower scattering rate. The number of scattered atoms at 25 μs is ~180, corresponding to 0.51% conversion of the initial total number N≊3.5×10⁴. The fraction 0.51% itself corresponds approximately to a scaling law of ~a^3/2, which is the same as the scaling of the total initial number of trapped atoms in the Thomas–Fermi limit for a fixed peak density.

Since the widths of the correlation functions are governed by the width of the momentum distribution of the initial colliding condensates, we expect corresponding broadening of the correlation functions as well (see figure 12). To quantify this effect, we fit the momentum distribution of the initial BEC by a Gaussian exp {–k_i²/[2(σ_i)²]}, where σ_x = 0.0036k_r and σ_{y, z} = 0.068k_r (cf with σ_x = 0.0025k_r and σ_{y, z} = 0.055k_r in figure 5, which are smaller). The Gaussian fits to the correlation functions in figure 12 are

where the widths σ_i^BB and σ_i^CL are given by

We see that the relative widths are practically unchanged, implying that the absolute widths are broadened. The ratios between the CL and BB correlation widths are slightly increased and are given by σ_x^CL/σ_x^BB≊1.20 and σ_{y, z}^CL/σ_{y, z}^BB≊1.18.

These numerical results make the present example—with the diminished role of atom–atom interactions—somewhat closer to the simple analytic predictions of equations (9) and (10) based on a Gaussian ansatz for non-interacting condensates.

5.5. Relative number squeezing and violation of Cauchy–Schwartz inequality

Another useful measure of atom–atom correlations is the normalized variance of the relative number fluctuations between atom numbers and in a pair of counting volume elements denoted via i and j,

where is the fluctuation. This definition uses the conventional normalization with respect to the shot-noise level characteristic of Poissonian statistics, such as for a coherent state, . In this case, the variance V_i–j = 1, which corresponds to the level of fluctuations in the absence of any correlation between and . Variance smaller than one, V_i–j < 1, implies reduction (or squeezing) of fluctuations below the shot-noise level and is due to quantum correlation between the particle number fluctuations in and . Perfect (100%) squeezing of the relative number fluctuations corresponds to V_i–j = 0.

In the context of the present model for the BEC collision experiment and possible correlation measurements between atom number fluctuations on diametrically opposite sides of the s-wave scattering sphere, we assign the indices i,j = A,B, C and D in equation (32) to one of the four quadrants as illustrated in figure 13. The total atom number operator in each quadrant within the s-wave scattering sphere is defined after elimination of the regions in momentum space occupied by the two colliding condensates

Operationally, this is implemented by discarding the data points beyond |k_x| > 0.8k_r. In addition, the quadrants are defined on a 2D plane after integrating the momentum distribution along the z-direction, which in turn only takes into account the 3D data points satisfying |1–k²/k_r²| < 0.28, i.e. lying in the narrow spherical shell k_r±δk with δk≊0.14k_r. The elimination of the inner and outer regions of the halo is done to minimize the sampling error in our simulations, since these regions have vanishingly small population and produce large noise in the stochastic simulations.

The choice of the quadrants as above is a particular implementation of the procedure of binning, known to result in a stronger correlation signal and larger relative number squeezing [11, 36]. Due to strong BB pair correlations, we expect the relative number fluctuations in the diametrically opposite quadrants to be squeezed, V_A–C, V_B–D < 1, while the relative number variance in the neighboring quadrants, such as V_A–B and V_C–D, is expected to be larger than or equal to one. The positive-P simulations confirm these expectations and are shown in figure 14, where we see strong (~80%) relative number squeezing for the diametrically opposite quadrants, V_A–C,B–D≊1–0.8 = 0.2.

These results assume a uniform detection efficiency of η = 1, whereas if the efficiency is less than 100% (η < 1), then the second term in equation (32) should be multiplied by η. This implies, that for η = 0.1 as an example, the above prediction of ~80% relative number squeezing will be degraded down to a much smaller but still measurable value of ~8% squeezing (V_A–C,B–D≊1–0.08 = 0.92). Even with perfect detection efficiency, our simulations do not lead to ideal (100%) squeezing. This can be understood in terms of a small fraction of collisions that take place with a center-of-mass momentum offset that is (nearly) parallel to one of the borders between the quadrants. As a result, the respective scattered pairs fail to appear in diametrically opposite quadrants during the (finite) propagation time (see also [36]).

For the symmetric case with and , the variance V_i–j can be rewritten as

where the second-order correlation function g_ij⁽²⁾ is defined according to

Equation (34) helps to relate the relative number squeezing, V_i–j < 1, to the violation of the classical Cauchy–Schwartz inequality g₁₂⁽²⁾ > g₁₁⁽²⁾, studied extensively in quantum optics with photons [33, 37]. The analysis presented here (see also [36] on molecular dissociation) shows that the Cauchy–Schwartz inequality, and its violation, is a promising area of study in quantum atom optics as well.

6. Summary

An important conclusion that we can draw from the numerical simulations is that the predicted widths of the correlation functions are remarkably robust against the parameter variations we were able to explore (in section 5.1 through 4). This gives us confidence in our physical interpretation of the width as being chiefly due to the initial momentum width of the condensate. The discrepancy with the analytical calculation of [28] seems to be primarily due to the different cloud shapes used. The width of the halo varies with the parameters we tested in a predictable way and also confirms the discussion in section 3.

As for comparison with the experiment, the numerically calculated widths of the scattering halo and the correlation functions coincide with the experimental ones to within better than 20% in most cases. The main discrepancy with the experiment is in the ratio of the BB and CL correlation widths. From the experimental point of view, these ratios are more significant than the individual widths since some sources of uncertainty, such as the number of atoms and the size of the condensates, cancel. The discrepancy may mean that the CL correlations are not sufficient to characterize the size and momentum distribution in the source at this level of accuracy. The discrepancies may of course also be due to the numerous experimental imperfections, especially the fact that the Raman outcoupling was only 60% efficient, and therefore an appreciable trapped m_x = 1 condensate was left behind. This defect may be remedied in future experiments. On the other hand, the current simulations neglect the unavoidable interaction of the scattered atoms with unscattered, m_x = 0 condensates as they leave the interaction region. This interaction could alter the trajectories of the scattered atoms in a minor, but complicated way. Future numerical work must examine this possibility further.

Still, the overall message of this work is that a first principles quantum field theory approach can quantitatively account for experimental observations of atomic four-wave mixing experiments. This work represents the first time that this sort of numerical simulation has been carefully confronted with an experiment. An interesting extension would be to examine the regime of stimulated scattering. It has been predicted that a highly anisotropic BEC could lead to an anisotropic population of the scattering halo [38, 39]. This effect would be a kind of atomic analogue of superradiance observed when off-resonant light is shone on a condensate [40, 41]. In addition, our results may be useful beyond the cold atom community: theoretical descriptions of correlation measurements in heavy ion collisions [42] may benefit from some of our insights.

Acknowledgments

We acknowledge stimulating discussions with A Aspect, K Mølmer, M Trippenbach, P Deuar and M Davis, and thank the developers of the XMDS software [43]. CS and KK acknowledge support from the Australian Research Council and the Queensland State Government. Part of this work was done at the Institut Henri Poincaré—Centre Emile Borel; KK thanks the institute for hospitality and support. The atom optics group is supported by the SCALA program of the EU and by the IFRAF institute.

Appendix A. Duration of the collision

In order to estimate the collision duration one can consider a simple classical model of the collision [24]. Denoting by ρ₁(x, t) and ρ₂(x, t) the density distributions of the two condensates, the number of scattered atoms N_sc(t) at a given time can be written

where σ₀ = 8πa₀₀² is the cross-section for a collision of two particles. In this latter formula a₀₀≊5.3 nm is the scattering length between m_x = 0 atoms [16].

The time-dependent density of the two condensates can be calculated from the expansion of a condensate in the Thomas–Fermi regime described in [44]. This approach suggests two different timescales for the collision duration. First, the separation of the two condensates occurs in a time defined by the ratio of the longitudinal size of the condensates and their relative velocity t_sep = R_x/v_r. Taking for R_x the Thomas–Fermi radius of the initial condensate, one can show that t_sep is on the order of 1 ms. At the same time, the condensates expand during their separation on a timescale t_exp = 1/ω_y = 1/ω_z≊140 μs. This latter effect appears to be predominant in the evaluation of equation (A.1) and t_exp can be taken as a definition of the collision duration Δt. The numerical evaluation of equation (A.1) gives N_sc(Δt)≊0.66N_sc(∞) and the estimated total number of scattered atoms corresponds to the experimentally observed 5% of the initial total number of atoms in the trapped condensate.

Appendix B. Occupation number of the scattering modes and amplitude of the BB correlation

In order to estimate the occupation number of the scattering modes one needs to compare the number of scattered atoms N_sc to the number of scattering modes N_m. To achieve this one has to first consider the volume of a scattering mode V_m, given by the first-order coherence volume (also dubbed `phase grain' in [12, 15]). Such a volume corresponds in fact to the coherence volume of the source condensate, and in practice it can also be deduced from the measurement of the width of the CL correlation function g_CL⁽²⁾(Δk_i) as one expects in a Hanbury Brown–Twiss experiment. For simplicity, we match the scattering mode volume V_m to the coherence volume of the source condensate in momentum space,

where β is a geometrical factor which depends on the geometry of the modes. Approximating the source condensate in momentum space by a Gaussian exp [–x²/(2σ_x²)–(y² + z²)/(2σ_{y, z}²)], one has β = (2π)^3/2.

The number of scattering modes N_m can in turn be estimated from knowledge of the total volume of the scattering shell V,

where the volume V is determined from the value of the width of the scattering shell δk:

for δk  k_r. If we apply this estimate to the results of the main numerical example (see section 5.1), we find N_m≊26 400. As N_sc = 1750, this implies an occupation number per mode of N_sc/N_m≊0.066. Such an estimate confirms that the system is indeed in the spontaneous regime and that bosonic stimulation effects are negligible.

The simple model of [16] for the BB correlation predicts that its height is given by

Using the above estimate of N_m and the actual value of N_sc found from the numerical simulations, we obtain that the height of the BB correlation peak should be approximately given by ~16. This compares favorably with the actual numerical result of 10.2. Similarly, we obtain the BB correlation peak of: ~62 in the example with the shorter collision time (compare with the numerical result of 36.6); ~18 in the example with the smaller collision velocity (compare with 10); and ~70 in the example with the smaller scattering length (compare with 50).

Appendix C. Width of the s-wave scattering sphere in the undepleted `pump' approximation

To estimate the width of the halo of scattered atoms beyond the spontaneous regime we use the analytic solutions for a uniform system in the so-called undepleted `pump' approximation in which the number of atoms in the colliding condensates are assumed constant. This approximation is applicable to short collision times. Nevertheless, it formally describes the regime of stimulated scattering and can be used to estimate the width of the s-wave scattering sphere as we show here.

The problem of BEC collisions in the undepleted `pump' approximation was studied in [30]; the solutions for the momentum distribution of the s-wave scattered atoms are formally equivalent to those obtained for dissociation of a BEC of molecular dimers in the undepleted molecular condensate approximation [14, 31]. For a uniform system with periodic boundary conditions, one has the following analytic solution for momentum mode occupation numbers:

Here, the constant is given by

where corresponds to the coupling constant of [30], and we note that the results of [30] contain typographical errors and have to be corrected as follows [45]: given the Hamiltonian of (1), with , the coupling g in (2), (7), (9) and (10), as well as in the definition of Δ(p) after (9), should be replaced by 2g. In the problem of molecular dissociation, the constant corresponds to  [14], where χ is the atom–molecule coupling and ρ₀ is the molecular BEC density.

The parameter Δ_k in equation (C.1) corresponds to the energy offset from the resonance condition

where is the collision momentum; in molecular dissociation, corresponds to the effective dissociation energy , using the notation of [14].

From equation (C.1), we see that modes with experience Bose enhancement and grow exponentially with time, whereas the modes with oscillate at the spontaneous noise level. The absolute momenta of the exponentially growing modes lie near the resonant momentum , and therefore we can use the condition to define the approximate width of the s-wave scattering sphere. First, we write k = k_r + Δk and assume for simplicity that k_r is large enough so that Δk  k_r. Then the condition can be approximated by

This can be solved for Δk and used to define the width δk = Δk/2 of the s-wave scattering sphere as

The reason for defining it as half of Δk is to make δk closer in definition to the half-width at half maximum and to the rms width around k_r.

The above simple analytic estimate (C.5) gives δk/k_r≊0.05 for the present ⁴He^* parameters. For comparison, the actual width of the analytic result (C.1) varies between δk/k_r≊0.12 and δk/k_r≊0.027 for durations between and , corresponding, respectively, to t≊20 μs and t≊140 μs in the present ⁴He^* example.

Appendix D. Positive-P simulation parameters

The positive-P simulations in our main numerical example of section 5 are performed on a computational lattice with 1400×50×70 points in the (x, y, z)-directions, respectively. The length of the quantization box along each dimension is L_x = 252 μm, L_y = 20.52 μm and L_z = 30.76 μm. The computational lattice in momentum space is reciprocal to the position space lattice and has the lattice spacing of Δk_i = 2π/L_i, giving Δk_x = 2.49×10⁴ m^–1, Δk_y = 3.06×10⁵ m^–1 and Δk_z = 2.04×10⁵ m^–1. The momentum cutoffs are k_x^(max) = 1.75×10⁷ m^–1, k_y^(max) = 7.66×10⁶ m^–1 and k_z^(max) = 7.15×10⁶ m^–1.

The momentum cutoff in the collision direction, k_x^(max), is more than 3 times larger than the collision momentum k_r, and hence it captures all relevant scattering processes of interest, including the energy non-conserving scatterings (k_r) + (k_r)→(3k_r) + (–k_r) and (–k_r) + (–k_r)→(–3k_r) + (k_r) [15]. In all our figures, the regions of momentum space covering k_x≊ ± 3k_r are not shown for the clarity of presentation of the main halo. These scattering processes, which produce a weak but not negligible signal at k_x≊± 3k_r, i.e. outside the main halo, are enhanced by Bose stimulation due to the large population of the colliding condensate components at k_x≊ k_r, respectively. In the remaining y- and z-directions, such processes are absent and therefore the number of lattice points and the momentum cutoffs can be smaller.

Since the momentum distribution of the initial condensate is narrowest in the k_x-direction, one may question whether the resolution of Δk_x = 2.49×10⁴ m^–1 with 1400 lattice points is sufficient. We check this by repeating the simulations with 4200×40×40 lattice points and quantization lengths of L_x = 753 μm and L_y = L_z = 15.4 μm, which give smaller lattice spacing Δk_x = 8.24×10³ m^–1, together with Δk_y = Δk_z = 4.08×10⁵ m^–1, k_x^(max) = 1.75×10⁷ m^–1 and k_y^(max) = k_z^(max) = 8.16×10⁶ m^–1. Our results on the new lattice reproduce the previous ones, within the sampling errors of the stochastic simulations. We typically average over 2800 stochastic trajectories, and take 128 time steps in the simulations over 25 μs collision time. A typical simulation of this size takes about 100 h on 7 CPUs running in parallel at 3.6 GHz clock speed.

References

[1]

Jeltes T et al 2007 Nature 445 402 
CrossRefPubMed

[2]

Hellweg D, Cacciapuoti L, Kottke M, Schulte T, Sengstock K, Ertmer W and Arlt J J 2003 Phys. Rev. Lett. 91 010406 
CrossRefPubMed

[3]

Esteve J, Trebbia J B, Schumm T, Aspect A, Westbrook C I and Bouchoule I 2006 Phys. Rev. Lett. 96 130403 
CrossRefPubMed

[4]

Rom T, Best T, van Oosten D, Schneider U, Fölling S, Paredes B and Bloch I 2006 Nature 444 733 
CrossRefPubMed

[5]

Öttl A, Ritter S, Köhl M and Esslinger T 2005 Phys. Rev. Lett. 95 090404 
CrossRefPubMed

[6]

Ritter S, Öttl A, Donner T, Bourdel T, Köhl M and Esslinger T 2007 Phys. Rev. Lett. 98 090402 
CrossRefPubMed

[7]

Hofferberth S, Lesanovsky I, Schumm T, Schmiedmayer J, Imambekov A, Gritsev V and Demler E 2007 Preprint 0710.1575
Preprint

[8]

Schellekens M, Hoppeler R, Perrin A, Viana Gomes J, Boiron D, Westbrook C I and Aspect A 2005 Science 310 648 
CrossRefPubMed

[9]

Yasuda M and Shimizu F 1996 Phys. Rev. Lett. 77 3090 
CrossRefPubMed

[10]

Fölling S, Gerbier F, Widera A, Mandel O, Gerike T and Bloch I 2005 Nature 434 481 
CrossRefPubMed

[11]

Greiner M, Regal C A, Stewart J T and Jin D S 2005 Phys. Rev. Lett. 94 110401 
CrossRefPubMed

[12]

Norrie A A, Ballagh R J and Gardiner C W 2005 Phys. Rev. Lett. 94 040401 
CrossRefPubMed

[13]

Norrie A A, Ballagh R J and Gardiner C W 2006 Phys. Rev. A 73 043617 
CrossRef

[14]

Savage C M, Schwenn P E and Kheruntsyan K V 2006 Phys. Rev. A 74 033620 
CrossRef

[15]

Deuar P and Drummond P D 2007 Phys. Rev. Lett. 98 120402 
CrossRefPubMed

[16]

Perrin A, Chang H, Krachmalnicoff V, Schellekens M, Boiron D, Aspect A and Westbrook C I 2007 Phys. Rev. Lett. 99 150405 
CrossRefPubMed

[17]

Drummond P D and Gardiner C W 1980 J. Phys. A: Math. Gen. 13 2353 
IOPscience

[18]

Gardiner C W and Zoller P 2000 Quantum Noise  (Berlin: Springer) 

[19]

Steel M J, Olsen M K, Plimak L I, Drummond P D, Tan S M, Collett M J, Walls D F and Graham R 1998 Phys. Rev. A 58 4824 
CrossRef

[20]

Drummond P D and Corney J F 1999 Phys. Rev. A 60 R2661 
CrossRef

[21]

Gilchrist A, Gardiner C W and Drummond P D 1997 Phys. Rev. A 55 3014 
CrossRef

[22]

Deuar P and Drummond P D 2006 J. Phys. A: Math. Gen. 39 1163 
IOPscience

[23]

Ziń P, Chwedeńczuk J, Veitia A, Rzażcedil;ewski K and Trippenbach M 2005 Phys. Rev. Lett. 94 200401 
CrossRefPubMed

[24]

Ziń P, Chwedeńczuk J and Trippenbach M 2006 Phys. Rev. A 73 033602 
CrossRef

[25]

Morgan S A, Rusch M, Hutchinson D A W and Burnett K 2003 Phys. Rev. Lett. 91 250403 
CrossRefPubMed

[26]

Yurovsky V A 2002 Phys. Rev. A 65 033605 
CrossRef

[27]

Sinatra A, Lobo C and Castin I 2002 J. Phys. B: At. Mol. Opt. Phys. 35 3599 
IOPscience

[28]

Molmer K, Perrin A, Krachmalnicoff V, Leung V, Boiron D, Aspect A and Westbrook C I 2008 Phys. Rev. A 77 033601 
CrossRef

[29]

Ögren H M and Kheruntsyan K V (unpublished)

[30]

Bach R, Trippenbach M and Rzażcedil;ewski K 2002 Phys. Rev. A 65 063605 
CrossRef

[31]

Kheruntsyan K V and Drummond P D 2002 Phys. Rev. A 66 031602 
CrossRef

[32]

Davis M J, Thwaite S J, Olsen M K and Kheruntsyan K V 2008 Phys. Rev. A 77 023617 
CrossRef

[33]

Walls D and Milburn G 1991 Quantum Optics  (New York: Springer) 

[34]

Abrikosov A A, Gorkov L P and Dzyaloshinski I E 1963 Methods of Quantum Field Theory in Statistical Physics  (New York: Springer) 

[35]

Moal S, Portier M, Kim J, Dugué J, Rapol U D, Leduc M and Cohen-Tannoudji C 2006 Phys. Rev. Lett. 96 023203 
CrossRefPubMed

[36]

Savage C M and Kheruntsyan K V 2007 Phys. Rev. Lett. 99 220404 
CrossRefPubMed

[37]

Zou X, Wang L and Mandel L 1991 Opt. Commun. 84 351 
CrossRef

[38]

Vardi A and Moore M G 2002 Phys. Rev. Lett. 89 090403 
CrossRefPubMed

[39]

Pu H and Meystre P 2000 Phys. Rev. Lett. 85 3987 
CrossRefPubMed

[40]

Inouye S, Chikkatur A, Stamper-Kurn D M, Stenger J, Pritchard D and Ketterle W 1999 Science 285 571 
CrossRefPubMed

[41]

Gross M and Haroche S 1982 Phys. Rep. 93 301 
CrossRef

[42]

Wong C Y and Zhang W N 2007 Phys. Rev. C 76 034905 
CrossRef

[43]

Collecut G R and Drummond P D 2001 Comput. Phys. Commun. 142 219 
CrossRef

[44]

Castin Y and Dum R 1996 Phys. Rev. Lett. 77 5315 
CrossRefPubMed

[45]

Trippenbach M, private communication

Notes

Note5  A 3D, animated rendition of the atomic positions 320 ms after release from the trap, available from stacks.iop.org/NJP/10/045021/mmedia. The vertical positions are derived from the arrival times as described in [16]. Each point corresponds to the detection of one atom and the animation shows the sum of 50 separate runs. The ellipsoids at the sides are the colliding condensates. The ellipsoids at the top and bottom result from imperfect Raman polarizations and stimulated atomic four-wave mixing (see [16]). The four condensates are excluded from the analysis given in the text.