Atomic four-wave mixing via condensate collisions

# Atomic four-wave mixing via condensate collisions

A. Perrin, C. M. Savage, D. Boiron, V. Krachmalnicoff, C. I. Westbrook and K.V. Kheruntsyan 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
###### Abstract

We perform a theoretical analysis of atomic four-wave mixing via a collision of two Bose-Einstein condensates of metastable helium atoms, and compare the results to a recent experiment. We calculate atom-atom pair correlations within the scattering halo produced spontaneously during the collision. We also examine the expected relative number squeezing of atoms on the sphere. The analysis includes first-principles quantum simulations using the positive -representation method. We develop a unified description of the experimental and simulation results.

###### pacs:
03.75.Kk, 34.50.-s, 03.75.Nt
: New J. Phys.

## 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, 2, 3, 4, 5, 6, 7, 8, 9, 10, 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, 13, 14]).

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 [15] 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 [15]. Figure 1 is a schematic momentum space diagram of these experiments. Two condensates, whose atoms have approximately equal but opposite momenta, and , interact by four-wave mixing, while they spatially overlap, to produce correlated atomic pairs with approximately equal but opposite momenta, and , satisfying momentum conservation, . Figure 1 corresponds to the experimental data shown in figure 2 of [15], 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 -representation method [16, 17, 18]. 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 [19] as the simulation time increases, eventually leading to diverging results.

An empirically estimated upper bound for the positive- simulation time (with controllable sampling error) of condensates with -wave scattering interactions is given approximately by [20]

 tsim≲2.5m(ΔV)1/3/[4πℏaρ2/30], (1)

where is the atom mass, is the -wave scattering length, is the condensate peak density, and is the volume of the elementary cell of the computational lattice, with lattice spacings of , , and . 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 s), which are about times shorter than the experimental interaction time of [15]. Despite this, our positive- 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 [21, 22, 23, 24]) to establish the range of their validity.

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 -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 Sec. 2 of this paper we will summarize the experimental results we wish to analyze. In Sec. 3 we discuss order of magnitude estimates. In Sec. 4 we describe simulations using the positive -representation method, and in Sec. 5 we discuss the results of our simulations. Sec. 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 to atoms confined in a magnetic trap whose frequencies are: Hz and Hz. A sudden Raman outcoupling drives the trapped He* from the Zeeman sublevel into the magnetic field insensitive state .  [15]. The Raman transition also splits the initial () condensate into two roughly equally populated condensates with opposite velocities along the direction. The magnitude of each velocity is equal to the recoil velocity cm/s, defined by the momentum of the photons used to create the colliding condensates , m. The relative velocity of the two condensates is about times higher than the speed of sound of the initial condensate, ensuring that 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 () of the total initial atom number into the halo. The system is shown in three-dimensions in an accompanying video of the experimental results after a  ms time of flight 111A 3 dimensional, animated rendition of the atomic positions 320 ms after release from the trap. The vertical positions are derived from the arrival times as described in [15]. 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 4 wave mixing (see [15]). The 4 condensates are excluded from the analysis given in the text.. For the purposes of this paper, the experiment consists in the acquisition of the three dimensional 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 and collinearly. The back-to-back correlation results from binary, elastic collisions between atoms, whereas the collinear correlation is a two particle interference effect involving members of two different pairs: a Hanbury Brown-Twiss correlation [25]. 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,

 G(2)(k1,k2)=⟨:^n(k1)^n(k2):⟩. (2)

Here, is the momentum density operator, are 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

 ⟨:^n(k1)^n(k2):⟩=⟨^a†(k1)^a†(k2)^a(k2)^a(k1)⟩. (3)

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

 g(2)CL(Δki)=∫Dd3k G(2)(k,k+eiΔki)∫Dd3k ⟨^n(k)⟩⟨^n(k+eiΔki)⟩, (4)

where is the unit vector in the direction. The normalization of ensures that for uncorrelated densities . 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 is lost.

The average back-to-back (BB) correlation function between two diametrically opposite points, one of which is additionally displaced by in the -direction, is defined similarly to :

 g(2)BB(Δki)=∫Dd3k G(2)(k,−k+eiΔki)∫Dd3k ⟨^n(k)⟩⟨^n(−k+eiΔki)⟩. (5)

The experimental observations can be summarized as follows. The width of both correlation functions along the axial direction of the condensate, the -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 and for collinear and back-to-back cases respectively. The experimental results are summarized in the following table

 σBBy,z/kr σCLy,z/kr σCLy,z/σBBy,z 0.081±0.004 0.069±0.008 0.85±0.15
(6)

One can also use the data to deduce the averaged radial width 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 -axis, but along the accessible directions we find .

## 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 [15], the width of the back-to-back and collinear 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 , an approximate theoretical treatment based on a simple ansatz for the pair wavefunction predicts a Gaussian dependence of the back-to-back (BB) and collinear (CL) correlation functions on the relative wavevectors [25]:

 G(2)(k,−k+niΔki) ∝ exp(−Δk2i2(σBBi)2), (7) G(2)(k,k+niΔki) ∝ exp(−Δk2i2(σCLi)2). (8)

The widths of the back-to-back () and collinear () correlations are related to the momentum-space width of the initial (source) condensate via [25]

 σBBi/σi = √2, (9) σCLi/σi = 2, (10)

and therefore the width of the back-to-back correlation is times smaller than the width of the collinear correlation.

In Sec. 5.1 the initial momentum-space widths are found to be and , assuming atoms. Expressing the experimentally measured widths in units of , we can rewrite (6) as

 σBBy,z/σy,z σCLy,z/σy,z σCLy,z/σBBy,z 1.47±0.07 1.25±0.15 0.85±0.15
(11)

and therefore, (9) is in agreement with the measured width of the radial back-to-back correlation function, whereas (10) overestimates the width of the collinear correlation function by almost . As we show below, first-principles simulations using the positive- 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 collinear correlations than for the back-to-back 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- 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 [26], 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 (in , or direction) of the colliding condensates imposes a minimum width

 δki≳σi. (12)

This limit suggests that the halo could be anisotropic. As noted above however, the experiment in [15] is not highly sensitive to such an anisotropy, and measures the width chiefly in the -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 . Here we discuss two mechanisms that impose a finite radial width on the halo.

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

 δk≃mℏkrΔt, (13)

where we assumed . 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- simulations however, numerical stability problems limit the maximum collision time that can be simulated, as discussed in Sec. 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 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 [27]. Under this approximation, the present system is equivalent to the dissociation of a condensate of molecular dimers studied in detail in [13, 28, 29]. The latter system in turn is analogous to parametric down-conversion in optics [30]. The details of the approximate solutions, common to condensate collisions and molecular dissociation, and the relationship between them are given in C. The resulting width of the halo found from this approach is

 δk≃4πaρ0kr. (14)

We see that in this regime, the width is proportional to the scattering length 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) , 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 can be understood via the quadratic dependence of the energy on momentum : to get the same excitation at a given energy offset , (16ar), one requires smaller absolute momentum offset at larger than at small .

Positive- simulations covering the transition from the spontaneous to stimulated regimes are available for Na condensate collisions as in [14]. 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 s and s the actual numerical results for the width of the spherical halo vary, respectively, between and , whereas Equation (14) predicts .

For He, on the other hand, the small mass and the larger scattering length of He atoms limit the maximum simulation time to 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 . Adding this width in quadrature to the momentum width of the initial condensate, , gives , not far from the experimentally observed radial momentum width of . 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 BECs is given by

 ^H=∫dx{ℏ22m|∇^Ψ|2+ℏU02^Ψ†^Ψ†^Ψ^Ψ}, (15)

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 –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 assumes a UV momentum cutoff . In our numerical simulations the momentum cutoff is imposed explicitly via the finite computational lattice. If the lattice spacings (, , ) in each spatial dimensions are chosen to be much larger than the –wave scattering length , then the respective momentum cutoffs satisfy . In this case the coupling constant is given by the familiar expression [31] without the need for explicit renormalization.

To model the dynamics of quantum fields describing the collision of two BECs, we use the positive –representation approach [16]. In this approach the quantum field operators and are represented by two complex stochastic –number fields and whose dynamics is governed by the following stochastic differential equations [14]:

 ∂Ψ(x,t)∂t = iℏ2m∇2Ψ−iU0~ΨΨΨ+√−iU0Ψ2ζ1(x,t), (16a) ∂~Ψ(x,t)∂t = −iℏ2m∇2~Ψ+iU0Ψ~Ψ~Ψ+√iU0~Ψ2ζ2(x,t). (16b)

Here, and are real independent noise sources with zero mean, , and the following nonzero correlation:

 ⟨ζj(x,t)ζk(x′,t′)⟩=δjkδ(3)(x−x′)δ(t−t′). (16q)

The stochastic fields and are independent of each other [] except in the mean, , where the brackets refer to stochastic averages with respect to the positive –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 and :

 ⟨[^Ψ†(x,t)]m[^Ψ(x′,t)]n⟩=⟨[~Ψ(x,t)]m[Ψ(x′,t)]n⟩. (16r)

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

 Ψ(x,0)=⟨^Ψ(x,0)⟩=√ρ0(x)/2(eikrx+e−ikrx), (16s)

with . Here, 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 sublevel into the magnetic field insensitive state , using two horizontally counter-propagating lasers and a third vertical laser [15]. In this geometry, the Raman transitions split the initial () condensate into two equally populated condensates and simultaneously impart velocities of 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 , having the –wave scattering length of nm ([15] and references therein), while the initial density profile refers to that of the trapped atoms in the state having the scattering length of nm [32]. The same distinction in terms of the scattering length in question applies to the definition of the interaction strength , in which has to be understood as for the trapped condensate or as for the outcoupled cloud.

In our simulations we assume for simplicity that the outcoupling from the trapped state is efficient, in which case the entire population is transferred into the state and therefore we have to only model -wave scattering interactions between the atoms in the state. In the experiment, on the other hand, the transfer efficiency is only about and therefore the collisions between the atoms in the and are not completely negligible and maybe 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- numerical simulations of collisions of two condensates of He atoms ( kg) as in the experiment of [15]. The key parameters in our main numerical example are the collision velocity, cm/s, and the peak density of the initial trapped condensate, m. The trap frequencies are matched exactly with the experimental values, Hz and Hz. The -wave scattering length for the magnetically trapped atoms in the sublevel is nm; the -wave scattering length for the outcoupled atoms in the sublevel is  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 . With these parameters, the average kinetic energy of colliding atoms is K, which is about times larger than the mean-field energy of the initial condensate K.

The duration of simulation in the main example is s. This is considerably smaller than the estimated duration of collision in the experiment, s (see A). The number of scattered atoms in our numerically simulated example at s is , representing % 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 , 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 s with the experimentally measured fraction of % at the end of collision at s, we first note that these time scales 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 % can be extrapolated to about % to correspond to the expected fraction at s. Next, one has to scale this value by a factor to account for the fact that in the experiment only of the initial atom number was transferred to the state of the colliding condensates. Accordingly, our theoretical estimate of % should be proportionally scaled down to conversion, in good agreement with the experimentally estimated fraction of % (see also 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 . The numerical result is fitted with a Gaussian (dashed line), centered at and having the radial width of m, where . The fitted radial width of of the numerical simulation is in reasonable agreement with the simple estimate of Equation (13), which gives for s.

Figure 5 shows the numerical results for the back-to-back and collinear correlations (solid lines with circles), defined in Equations (4) and (5). Due to the symmetry of the and 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 now refers to the actual wave-vector component ). 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:

 g(2)BB(Δki)−1 = 9.2exp{−Δk2i/[2(σBBi)2]}, (16t) g(2)CL(Δki)−1 = exp{−Δk2i/[2(σCLi)2]}, (16u)

where the correlation widths and are shown in the table (16v) below. The Gaussian fits (dashed lines) for the slices of the initial momentum distribution are scaled to the same peak value as and have and .

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

 σBBx/σx σBBy,z/σy,z σCLx/σx σCLy,z/σyz 1.18 1.39 1.27 1.57
(16v)

The ratios between the collinear and back-to-back correlation widths are and . The errors due to stochastic sampling on all quoted values of the correlation widths are smaller than %.

The values for and 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 s and the experimentally measured values after a s interaction time may be due to the evolution of the condensates past s, not attainable within the positive- method. The above numerical results for the correlation widths can be also 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 back-to-back and collinear widths by and , respectively, in the present example.

The amplitude of the correlation functions can also be inferred by simple models. In fact, the collinear 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  [25]. This is in agreement with the positive- simulations. The back-to-back correlation amplitude, on the other hand, can be substantially higher and display super-bunching () [13, 22] 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 [15], the amplitude of the back-to-back correlation can be linked to the inverse population of the atomic modes on the halo. As we show in 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 Sec. 5.1, except that the data are analyzed at s, which is half the previous interaction time. We found in Sec. 5.1 that and the width of the halo are all nearly the same. In Sec. 3, however, we argue that the widths of the correlation functions and the halo are governed by different limits [Equations (10),(9) and (13) or (14), respectively]. The example in this section illustrates this point.

Figure 6 shows two orthogonal slices of the -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 s and the fitted Gaussian gives the radial width of . This is precisely twice the width in Figure 4 and is in agreement with the simple qualitative estimate of Equation (13).

The back-to-back and collinear correlation functions after s collision time are qualitatively very similar to those shown in figure 5, except that the Gaussian fits are

 g(2)BB(Δki)−1 = 35.6exp{−Δk2i/[2(σBBi)2]}, (16w) g(2)CL(Δki)−1 = exp{−Δk2i/[2(σCLi)2]}, (16x)

with the correlation widths given by

 σBBx/σx σBBy,z/σy,z σCLx/σx σCLy,z/σyz 1.16 1.28 1.27 1.48
(16y)

The ratios between the widths are and .

For the correlation functions, the main difference compared to the case for s is that the peak value of the back-to-back 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- 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 , which is approximately half the number at 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, cm/s ( m), 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- simulations for the momentum density distribution at 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 Sec. 5.1. The width of the Gaussian function fitted to the numerically calculated radial momentum distribution is given by . 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 ) in the present example is slightly shifted towards the centre of the halo, , which is a feature predicted in [27] to occur when the ratio of the kinetic energy to the interaction energy per particle is reduced.

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

 g(2)BB(Δki)−1 = 9.0exp{−Δk2i/[2(σBBi)2]}, (16z) g(2)CL(Δki)−1 = exp{−Δk2i/[2(σCLi)2]}, (16aa)

with the correlation widths given by

 σBBx/σx σBBy,z/σy,z σCLx/σx σCLy,z/σy,z 1.16 1.35 1.31 1.51
(16ab)

where and . The ratios between the collinear and back-to-back correlation widths are and .

As we see from these results, the absolute widths of the correlation functions are practically unchanged compared to the main numerical example (16v). 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 Sec. 5.1.

The number of scattered atoms in this example is about , which is approximately times smaller than in Sec. 5.1 and corresponds to conversion. This scaling is in agreement with the rate equation approach [22], 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 Sec. 5.1, except that the scattering lengths and are artificially halved, i.e. nm and 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, m. The total number of atoms is now smaller, . 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 [25].

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 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 s is , corresponding to conversion of the initial total number . The fraction of itself corresponds approximately to a scaling law of , 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 , where and (cf. with and in figure 5, which are smaller). The Gaussian fits to the correlation functions in figure 12 are

 g(2)BB(Δki)−1 = 49exp{Δk2i/[2(σBBi)2]}, (16ac) g(2)CL(Δki)−1 = 0.94exp{Δk2i/[2(σCLi)2]}, (16ad)

where the widths and are given by

 σBBx/σx σBBy,z/σy,z σCLx/σx σCLy,z/σy,z 1.18 1.53 1.42 1.81
(16ae)

We see that the relative widths are practically unchanged, implying that the absolute widths are broadened. The ratios between the collinear and back-to-back correlation widths are slightly increased and are given by and .

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 noninteracting 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 and ,

 Vi−j=⟨[Δ(^Ni−^Nj)]2⟩⟨^Ni⟩+⟨^Nj⟩=1+⟨:[Δ(^Ni−^Nj)]2:⟩⟨^Ni⟩+⟨^Nj⟩, (16af)

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 , which corresponds to the level of fluctuations in the absence of any correlation between and . Variance smaller than one, , 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 (%) squeezing of the relative number fluctuations corresponds to .

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 -wave scattering sphere, we assign the indices in Equation (16af) to one of the four quadrants as illustrated in figure 13. The total atom number operator in each quadrant within the -wave scattering sphere is defined after elimination of the regions in momentum space occupied by the two colliding condensates

 ^Ni(t)=∫Didkxdky∫+∞−∞dkz^n(k,t). (16ag)

Operationally, this is implemented by discarding the data points beyond . In addition, the quadrants are defined on a 2D plane after integrating the momentum distribution along the -direction, which in turn only takes into account the 3D data points satisfying , i.e. lying in the narrow spherical shell with . 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, 33]. Due to strong back-to-back pair correlations, we expect the relative number fluctuations in the diametrically opposite quadrants to be squeezed, , while the relative number variance in the neighboring quadrants, such as and , is expected to be larger than, or equal to, one. The positive- simulations confirm these expectations and are shown in figure 14, where we see strong () relative number squeezing for the diametrically opposite quadrants, .

These results assume a uniform detection efficiency of , whereas if the efficiency is less than % (), then the second term in Equation (16af) should be multiplied by . This implies, that for as an example, the above prediction of % relative number squeezing will be degraded down to a much smaller but still measurable value of % squeezing (). Even with perfect detection efficiency, our simulations do not lead to ideal (%) 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 [33]).

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

 Vi−j=1+⟨^Ni⟩[g(2)ii−g(2)ij], (16ah)

where the second-order correlation function is defined according to

 g(2)ij=⟨:^Ni^Nj:⟩⟨^Ni⟩⟨^Nj⟩. (16ai)

Equation (16ah) helps to relate the relative number squeezing, , to the violation of the classical Cauchy-Schwartz inequality , studied extensively in quantum optics with photons [34, 30]. The analysis presented here (see also [33] 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 Secs. 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 [25] 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 Sec. 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 back to back and collinear 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 collinear 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 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, 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 [35, 36]. This effect would be a kind of atomic analogue of superradiance observed when off-resonant light is shined on a condensate [37, 38]. In addition, our results may be useful beyond the cold atom community: theoretical descriptions of correlation measurements in heavy ion collisions [39] may benefit from some of our insights.

The authors acknowledge stimulating discussions with A. Aspect, K. Mølmer, M. Trippenbach, P. Deuar, and M. Davis, and thank the developers of the XMDS software [40]. CS and KK acklowledge support from the Australian Research Council and the Queensland State Government. Part of this work was done at the Institut Henri Poincare – 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 [22]. Denoting by and the density distributions of the two condensates, the number of scattered atoms at a given time can be written

 Nsc(t)=2∫t0dt′∫d3x 2σ0vrρ1(x,t′)ρ2(x,t′) (16aj)

where is the cross section for a collision of two particles. In this latter formula nm is the scattering length between atoms [15].

The time-dependent density of the two condensates can be calculated from the expansion of a condensate in the Thomas-Fermi regime described in [41]. This approach suggests two different time scales 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 . Taking for the Thomas-Fermi radius of the initial condensate, one can show that is on the order of  ms. At the same time, the condensates expand during their separation on a time scale s. This latter effect appears to be predominant in the evaluation of Equation (16aj) and can be taken as a definition of the collision duration . The numerical evaluation of Equation (16aj) gives and the estimated total number of scattered atoms corresponds to the experimentally observed of the initial total number of atoms in the trapped condensate.

## Appendix B Occupation number of the scattering modes and amplitude of the back-to-back correlation

In order to estimate the occupation number of the scattering modes one needs to compare the number of scattered atoms to the number of scattering modes . To achieve this one has to first consider the volume of a scattering mode , given by the first-order coherence volume (also dubbed as “phase grain”in [12, 14]). 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 collinear correlation function as one expects in a Hanbury Brown-Twiss experiment. For simplicity we match the scattering mode volume to the coherence volume of the source condensate in momentum space,

 Vm≃βσx(σyz)2, (16ak)

where is a geometrical factor which depends on the geometry of the modes. Approximating the source condensate in momentum space by a Gaussian , one has .

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

 Nm=VVm, (16al)

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

 V = ∫d3k exp[−(k−kr)2/(2δk2)] (16am) ≃ 4π√2πk2rδk, (16an)

for . If we apply this estimate to the results of the main numerical example (see Sec. 5.1), we find . As , this implies an occupation number per mode of . Such an estimate confirms that the system is indeed in the spontaneous regime and that bosonic stimulation effects are negligible.

The simple model of [15] for the back-to-back correlation predicts that its height is given by

 g(2)BB(0)=1+Nm/Nsc (16ao)

Using the above estimate of and the actual value of found from the numerical simulations, we obtain that the height of the back-to-back correlation peak should be approximately given by . This compares favorably with the actual numerical result of . Similarly, we obtain the back-to-back correlation peak of: in the example with the shorter collision time (compare with the numerical result of ); in the example with the smaller collision velocity (compare with ); and in the example with the smaller scattering length (compare with ).

## 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 -wave scattering sphere as we show here.

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

 nk(t)=¯¯¯g2¯¯¯g2−Δ2ksinh2(√¯¯¯g2−Δ2kt). (16ap)

Here, the constant is given by

 ¯¯¯g=2U0ρ0=8πℏa00ρ0m, (16aq)

where corresponds to the coupling constant of [27], and we note that the results of [27] contain typographical errors and have to be corrected as follows [42]: given the Hamiltonian of (1), with , the coupling in (2), (7), (9), and (10), as well as in the definition of