Atomic four-wave mixing via condensate collisions

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 P-representation method. We develop a unified description of the experimental and simulation results.


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 ( 4 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 4 He * atoms [15]. Figure 1 is a schematic momentum space diagram of these experiments. Two condensates, whose atoms have approximately equal but opposite momenta, k 1 and k 2 ≃ −k 1 , interact by four-wave mixing, while they spatially overlap, to produce correlated atomic pairs with approximately equal but opposite momenta, k 3 and k 4 , satisfying momentum conservation, k 1 + k 2 = k 3 + k 4 . 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 P -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 t sim increases, eventually leading to diverging results.
An empirically estimated upper bound for the positive-P simulation time (with controllable sampling error) of condensates with s-wave scattering interactions is given approximately by [20] t sim 2.5m(∆V ) 1/3 /[4π aρ where m is the atom mass, a is the s-wave scattering length, ρ 0 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 6 times shorter than the experimental interaction time of [15]. 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 [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 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 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 P -representation method, and in Sec. 5 we discuss the results of our simulations. Sec. 6 summarizes our work.

Overview of the experiment
The starting point of the experiment is a 4 He* condensate of 10 4 to 10 5 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 4 He* from the m x = 1 Zeeman sublevel into the magnetic field insensitive state m x = 0. [15]. 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, defined by the momentum of the photons used to create the colliding condensates k r , k r = 5.8 × 10 6 m −1 . The relative velocity 2v r of the two condensates is about 8 times higher than the speed of sound c s = µ/m 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 (5%) 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 320 ms time of flight ‡. 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.

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 normallyordered second-order correlation function between the densities at two points in momentum space, Here,n(k) =â † (k)â(k) is the momentum density operator,â † (k) areâ(k) are the Fourier transforms of the field creation and annihilation operatorsΨ † (x) andΨ(x), 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 collinear (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 ‡ A 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. where e i is the unit vector in the k i direction. The normalization of g (2) CL (∆k i ) ensures that for uncorrelated densities g (2) CL (∆k i ) = 1. The integration domain D 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 back-to-back (BB) correlation function g (2) BB (∆k i ) between two diametrically opposite points, one of which is additionally displaced by ∆k i in the k idirection, is defined similarly to g (2) 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 σ CL y,z and σ BB y,z for collinear and back-to-back cases respectively. The experimental results are summarized in the following table σ BB y,z /k r σ CL y,z /k r σ CL y,z /σ BB y,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 δ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 .

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.

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, n(k) = n(k) ∝ exp[− i=x,y,z k 2 i /(2σ 2 i )], with σ i = 1/w i , 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 ∆k i [25]: The widths of the back-to-back (σ BB i ) and collinear (σ CL i ) correlations are related to the momentum-space width σ i of the initial (source) condensate via [25] and therefore the width of the back-to-back correlation is √ 2 times smaller than the width of the collinear correlation.
In Sec. 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 4 atoms. Expressing the experimentally measured widths in units of σ y,z , we can rewrite (6) as σ BB y,z /σ y,z σ CL y,z /σ y,z σ CL y,z /σ BB y,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 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 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-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 [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.

Width of the scattered halo
A second, experimentally accessible quantity in a BEC collision is the width δk i 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 [15] is not highly sensitive to such an anisotropy, and measures the width chiefly in the y, 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 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 δ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 [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 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 ρ 0 , 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ρ 0 , one can excite and amplify spectral components that are further detuned from the exact resonance condition ∆ k = 0 (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 ∆ k , (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 23 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 300 µs 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 4 He * , on the other hand, the small mass and the larger scattering length of 4 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 23 Na result to extrapolate Equation (14) to 4 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 (0.05k r ) 2 + (0.055k r ) 2 = 0.074k r , 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.

Model
The effective field theory Hamiltonian governing the dynamics of the collision of BECs is given byĤ whereΨ(x, t) is the field operator with the usual commutation relation , m 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 0 δ(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, ∆z) in each spatial dimensions are chosen to be much larger than the s-wave scattering length a, then the respective momentum cutoffs satisfy k max x,y,z ≪ 1/a. In this case the coupling constant U 0 is given by the familiar expression U 0 ≃ 4π a/m [31] 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 [16]. In this approach the quantum field operatorsΨ(x, t) andΨ † (x, t) are represented by two complex stochastic c-number fields Ψ(x, t) andΨ(x, t) whose dynamics is governed by the following stochastic differential equations [14]: Here, ζ 1 (x, t) and ζ 2 (x, t) are real independent noise sources with zero mean, ζ j (x, t) = 0, and the following nonzero correlation: The stochastic fields Ψ(x, t) andΨ(x, t) are independent of each other [Ψ(x, t) = Ψ * (x, t)] except in the mean, Ψ (x, t) = Ψ * (x, t) , 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Ψ(x, t): 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 k r = mv r / and v r is the collision velocity) in the x direction, withΨ(x, 0) = Ψ * (x, 0). Here, ρ 0 (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 4 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 [15]. 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Ψ(x, t) represents the atoms in the untrapped state m x = 0, having the s-wave scattering length of a 00 = 5.3 nm ( [15] and references therein), while the initial density profile ρ 0 (x) refers to that of the trapped atoms in the m x = 1 state having the scattering length of a 11 = 7.51 nm [32]. The same distinction in terms of the scattering length in question applies to the definition of the interaction strength U 0 ≃ 4π a/m, in which a has to be understood as a 11 for the trapped condensate or as a 00 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 to only 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 maybe responsible for some of the deviations between the present theoretical results and the experimental observations.

Main numerical example
Here we present the results of positive-P numerical simulations of collisions of two condensates of 4 He * atoms (m ≃ 6.65 × 10 −27 kg) as in the experiment of [15]. The key parameters in our main numerical example are the collision velocity, v r = 9.2 cm/s, and the peak density of the initial trapped condensate, ρ 0 = 2.5 × 10 19 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 11 = 7.5 nm; the s-wave scattering length for the outcoupled atoms in the m x = 0 sublevel is a 00 = 5.3 nm. Other simulation parameters are given in Appendix 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 4 . With these parameters, the average kinetic energy of colliding atoms is E kin /k B = mv 2 r /2k B ≃ 2.0 × 10 −6 K, which is about 7.4 times larger than the mean-field energy of the initial condensate E M F /k B = 4π 2 a 11 ρ 0 /mk B ≃ 2.7 × 10 −7 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 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 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 2 to account for the fact that in the experiment only 60% of the initial atom number 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).
BB/CL (∆k i ) − 1 as a function of the displacement ∆k i (i = x, y, z) in units of the collision momentum k r , after t f = 25 µs collision time. The circles are the numerical results, angle-averaged over the halo of scattered atoms after elimination of the regions occupied by the two colliding condensates; the solid lines are simple Gaussian fits to guide the eye (see text). For comparison, we also plot the initial momentum distribution n 0 (k i ) of the colliding condensates; the actual data points are shown by the squares and are fitted by a dashed-line Gaussian.
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 0 ) 2 /(2δk 2 )] (dashed line), centered at k 0 = 0.98k r and having the radial width of δk = 0.10k r ≃ 5.8 × 10 5 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 back-to-back and collinear 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 wave-vector 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 σ BB i and σ CL i are shown in the table (22) below. The Gaussian fits (dashed lines) for the slices of the initial momentum distribution n 0 (k i ) ∝ exp{−k 2 i /[2(σ i ) 2 ]} 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 momentum distribution of the source. More specifically, we find that the following results provide the best fit to our numerical data: σ BB x /σ x σ BB y,z /σ y,z σ CL x /σ x σ CL y,z /σ yz 1. 18 1.39 1.27 1.57 (22) The ratios between the collinear and back-to-back correlation widths are σ CL x /σ BB x ≃ 1.08 and σ CL y,z /σ BB y,z ≃ 1.13. The errors due to stochastic sampling on all quoted values of the correlation widths are smaller than 3%.
The values for σ CL y,z /σ y,z and σ BB y,z /σ 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 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 ∼ 20% and 40%, 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 g (2) CL (0) = 2 [25]. This is in agreement with the positive-P simulations. The back-to-back correlation amplitude, on the other hand, can be substantially higher and display super-bunching (g (2) BB (0) ≫ 1) [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 Appendix B, this model follows the trends observed in our first-principles numerical simulations.

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 t f = 12.5 µs, which is half the previous interaction time. We found in Sec. 5.1 that σ BB yz , σ CL yz and the width of the halo δk 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 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 back-to-back and collinear 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 σ BB x /σ x σ BB y,z /σ y,z σ CL x /σ x σ CL y,z /σ yz 1. 16 1.28 1.27 1.48 (25) The ratios between the widths are σ CL x /σ BB x ≃ 1.09 and σ CL y,z /σ BB y,z ≃ 1.16. For the correlation functions, the main difference compared to the case for 25 µ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-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.

Smaller collision velocity
In this example, we present the results of simulations in which the collision velocity is smaller by a factor √ 2 than before, v ′ r = 6.5 cm/s (k ′ r = 4.1 × 10 6 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 Sec. 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 centre of the halo, k 0 = 0.92k ′ r , 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 with the correlation widths given by 16 1.35 1.31 1.51 (28) Figure 8. Same as in figure 3, except for √ 2 times smaller collision velocity, v ′ r = 6.46 cm/s (k ′ r = 4.09 × 10 6 m −1 ). The axis for the momentum components k i (i = x, y, z) are in units of the smaller recoil momentum than in figure 3, and therefore the absolute radius of the s-wave scattering sphere is smaller in the present example. where σ x /k ′ r ≃ 0.0035 and σ x /k ′ r ≃ 0.078. The ratios between the collinear and backto-back correlation widths are σ CL x /σ BB x ≃ 1.13 and σ CL y,z /σ BB y,z ≃ 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 Sec. 5.1.
The number of scattered atoms in this example is about 1270, which is approximately √ 2 times smaller than in Sec. 5.1 and corresponds to ∼ 1.3% 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 √ 2 times smaller here.

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 a 11 and a 00 are artificially halved, i.e. a 00 = 2.65 nm and a 11 = 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, ρ 0 ≃ 2.5 × 10 19 m −3 . The total number of atoms is now smaller, N ≃ 3.5 × 10 4 . 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 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 4 . The fraction of 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 2 i /[2(σ i ) 2 ]}, 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 ∼ √ 2 smaller). The Gaussian fits to the correlation functions in where the widths σ BB i and σ CL i are given by σ BB x /σ x σ BB y,z /σ y,z σ CL x /σ x σ CL y,z /σ y,z 1. 18 1.53 1.42 1.81 (31) 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 σ CL x /σ BB x ≃ 1.20 and σ CL y,z /σ BB y,z ≃ 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 noninteracting condensates.

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 numbersN i andN j in a pair of counting volume elements denoted via i and j, where ∆X =X − X 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, N i + N i . In this case the variance V i−j = 1, which corresponds to the level of fluctuations in the absence of any correlation betweenN i andN j . 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   elimination of the regions in momentum space occupied by the two colliding condensateŝ Operationally, this is implemented by discarding the data points beyond |k x | > 0.8k r . In addition, the quadrants D i 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 2 /k 2 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,33]. Due to strong back-to-back 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, 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 [33]).
For the symmetric case with N i = N j and N 2 i = N 2 j , the variance V i−j can be rewritten as where the second-order correlation function g (2) 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 (2) 12 > g (2) 11 , 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.

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 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 [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. 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 back-to-back 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 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 g (2) 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 /(2σ 2 x ) − (y 2 + z 2 )/(2σ 2 y,z )], one has β = (2π) 3/2 . The number of scattering modes N m can in turn be estimated from the 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 Sec. 5.1), we find N m ≃ 26400. 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 [15] for the back-to-back correlation predicts that its height is given by g (2) BB (0) = 1 + N m /N sc (B.5) 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 back-to-back correlation peak should be approximately given by ∼ 16. This compares favorably with the actual numerical result of 10.2. Similarly, we obtain the back-to-back 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).
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 g 2 − ∆ 2 k = 0 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 δk k r ≃ mg 2 k 2 r = 4πa 00 ρ 0 k 2 r . (C.5) 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 4 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 gt = 1 and gt = 7, corresponding, respectively, to t ≃ 20 µs and t ≃ 140 µs in the present 4 He * example. 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 hours on 7 CPUs running in parallel at 3.6 GHz clock speed.