{'fraction_non_alphanumeric': 0.05084523223371082, 'fraction_numerical': 0.016231741342524208, 'mean_word_length': 4.394510845506861, 'pattern_counts': {'":': 0, '<': 15, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
First-principles quantum simulations of dissociation of molecular condensates: Atom correlations in momentum space

We investigate the quantum many-body dynamics of dissociation of a Bose-Einstein condensate of molecular dimers into pairs of constituent bosonic atoms and analyze the resulting atom-atom correlations. The quantum fields of both the molecules and atoms are simulated from first principles in three dimensions using the positive-P representation method. This allows us to provide an exact treatment of the molecular field depletion and s-wave scattering interactions between the particles, as well as to extend the analysis to nonuniform systems. In the simplest uniform case, we find that the major source of atom-atom decorrelation is atom-atom recombination which produces molecules outside the initially occupied condensate mode. The unwanted molecules are formed from dissociated atom pairs with non-opposite momenta. The net effect of this process -- which becomes increasingly significant for dissociation durations corresponding to more than about 40% conversion -- is to reduce the atom-atom correlations. In addition, for nonuniform systems we find that mode-mixing due to inhomogeneity can result in further degradation of the correlation signal. We characterize the correlation strength via the degree of squeezing of particle number-difference fluctuations in a certain momentum-space volume and show that the correlation strength can be increased if the signals are binned into larger counting volumes.


I. INTRODUCTION
Dissociation of a diatomic molecule produces two quantum mechanically entangled atoms with equal and opposite momenta in the molecule's rest frame. These atoms have Einstein-Podolsky-Rosen (EPR) type correlations in position and momentum [1], and hence are of fundamental interest [2,3]. Experimental advances in coherent manipulation of quantum gases are now enabling the production of ultracold molecular gases [4] and even molecular Bose-Einstein condensates (BECs) [5] from atomic condensates and degenerate Fermi gases. These molecules can in turn be dissociated [6,7,8] into strongly correlated ensembles of bosonic [9,10,11] or fermionic [12,13] atoms, thus extending possible fundamental tests of quantum mechanics into macroscopic regimes [14]. This has close analogies with experiments in quantum optics using continuous-variable quadrature correlations generated via parametric down-conversion [15,16].
An important recent development in the study of ultracold quantum gases is the direct measurement of atomatom correlation functions and quantum statistics via the analysis of the noise in absorption images [8,36,37,38], atom counting using microchannel plate detectors [39,40], and fluorescence imaging combined with high-finesse optical cavities [41,42]. These experiments have close parallels with the pioneering photon correlation measurements of Hanbury Brown and Twiss (HBT) [43] and initiated further theoretical studies of HBT interferometry as a sensitive tool for probing quantum many-body states of ultracold gases [44,45,46,47,48]. All these correlation measurements are ultimately related to the quantum-atom optics counterpart of Glauber's secondorder, or density-density, correlation function [49]. While giving access to the simplest higher-order correlations, these techniques may in the future lead to more sophisticated correlation measurements, such as between matterwave quadratures required to demonstrate the version of the EPR paradox proposed in Ref. [14].
In correlation measurements via the noise in absorption images of ballistically expanded clouds [8,36], the extracted correlation signal reflects the momentum correlations before expansion [37]. This is particularly suitable for studies of molecular dissociation, since the correlations arise from momentum conservation and are therefore most fundamentally between atoms with equal and opposite momenta, k and − k.
In this paper, we study the quantum dynamics of dis-sociation of a molecular BEC and analyze the atomatom correlations in momentum space. Spatial correlations after expansion will be studied in a subsequent paper [50]. Our analysis is based on first-principles quantum simulations of the coupled atomic-molecular system in three dimensions (3D) using the positive Prepresentation method [51,52,53,54,55]. This allows us to go beyond the analytical results known for the case of a spatially uniform condensate and undepleted, classical molecular field [10,13], as well as beyond the pairing mean-field method used in Ref. [12]. In previous related work our simulations were limited to one spatial dimension and to the calculation of the number difference squeezing between all atoms with positive momenta and all with negative momenta [10]. In this paper we simulate three spatial dimensions and calculate squeezing localized in momentum space.
In the present treatment, we are able to extend these results and analyze the molecular depletion without approximations. In addition, the positive-P method allows us to treat nonuniform condensates, as well as s-wave scattering interactions between the species. The main limitation of the method is that the positive-P simulations eventually fail as the simulation time increases, particularly when two-body scattering is included [56]. For this reason, the simulations involving s-wave interactions are performed with slightly reduced values of the scattering lengths, or else are limited to short dissociation times. Nevertheless, the results allow us to gain quantitative insights into the role of these interactions in atom-atom correlations, which have not been analyzed before.
More importantly, we find that mode mixing due to inhomogeneity of the condensates can potentially have a quite disruptive effect on the atom-atom correlations. We give quantitative estimates of this effect. Similarly, we quantify the disruptive role of atom-atom recombination which becomes increasingly important past the initial spontaneous dissociation regime.
We emphasize that for the correct interpretation of the results on correlations, especially when comparing the uniform and nonuniform results, it is important to operate with observables that have a well-defined operational meaning. We find that the correlation strength is most conveniently quantified via the normalized numberdifference fluctuations. This is defined with respect to the number of atoms in a certain momentum-space volume around pairs of the carrier momenta of interest. This measure of the correlation strength can be defined with respect to different counting volumes, which corresponds to the procedure of binning often employed in experiments. We find that the strength of the correlation signal increases for larger bin sizes, which is in agreement with the results of analogous experimental measurements of Ref. [8] performed in position space. This paper is organized as follows. Section II describes the system we study including relevant approximations. The positive-P simulation method and the full stochastic equations we simulate are presented. Section III reviews the idealized analytic model presented by Kherutsyan [13]. This is the reference against which we can determine the effect of the various realistic processes included in our simulations. Section IV is the heart of the paper and reports numerical simulations that quantify the degradation of quantum correlations in the presence of: molecular depletion (Sec. IV B), scattering (Sec. IV C), and spatial nonuniformity (Sec. IV D). The effect of binning data is briefly discussed in Sec. IV D 3. The paper ends with a summary, followed by appendices on Fourier transforming the fields.

II. THE MODEL
The quantum field theory effective Hamiltonian describing the system we shall study is given by [10,57]: The molecular and atomic fields are, respectively, described by the bosonic operatorsΨ 0 (x, t) andΨ 1 (x, t) satisfying the following equal-time commutation relation: The first term in the Hamiltonian (1) describes the kinetic energy, where m 1 and m 0 = 2m 1 are the atomic and molecular masses. The trapping potentials, including internal energies E i , are given by V i (x). This determines the detuning 2∆ which corresponds to the overall energy mismatch 2 ∆ = [2V 1 (0) − V 0 (0)] = 2E 1 − E 0 between the free two-atom state at the dissociation threshold and the bound molecular state.
The coupling term χ is responsible for coherent conversion of molecules into atom pairs, e.g., via Raman transitions or a Feshbach resonance (see, for example, Refs. [57,58,59,60,61,62] and Ref. [63] for a recent review), while U ij are the strengths of the twobody s-wave interactions describing atom-atom, atommolecule, and molecule-molecule scattering. The diagonal terms are given by U ii = 4π a ii /m i , where a 11 and a 00 are the atom-atom and molecule-molecule s-wave scattering lengths. The off-diagonal terms U 01 = U 10 are given by U 01 = 2π a 01 /µ 01 = 3π a 10 /m 1 , where a 01 is the atom-molecule s-wave scattering length and µ 01 = m 0 m 1 /(m 0 + m 1 ) = 2m 1 /3 is the reduced mass.
The Hamiltonian (1) with delta-function interactions implicitly assumes a momentum cutoff k max 1/a ij [64], and we associate the coupling terms with the observed scattering lengths. In the numerical simulations on a finite computational lattice, where there is always a maximum cutoff, this description is valid as long as the lattice spacing in coordinate space is larger than a ij . Otherwise, a more careful renormalization procedure [62] is needed.
Starting from a stable (E 0 < 2E 1 ) molecular BEC, dissociation into the constituent atoms may be achieved by a rapid Feshbach sweep into the atomic side of the resonance or by a coherent Raman transition [6,8,10]. If the transition has a negative detuning 2 ∆ < 0 after the sweep (E 0 > 2E 1 ), the molecules become unstable against (spontaneous) dissociation into free atom pairs. The energy level configuration after the sweep is the initial condition for our simulations. We assume that the molecular BEC is initially in a coherent state, whereas the atoms are in the vacuum state. In addition, we assume that once the dissociation coupling χ is switched on, the trapping potentials are simultaneously switched off, so that the evolution of the atomic and molecular fields takes place in free space.
For molecules at rest, the excess of potential energy 2 |∆| is converted into kinetic energy, 2 |∆| → 2 2 k 2 /(2m 1 ), of atom pairs with equal but opposite momenta around ±k 0 , where k 0 = |k 0 | = 2m 1 |∆|/ . This is the physical origin of the expected correlations between the atoms, which we will study below in greater detail and in the context of many-body field theory.
To investigate the quantum dynamics we use stochastic differential equations corresponding to the positive-P representation of the density matrix [51,52,53,54,55]. The essence of this method is the mapping of the operator equations of motion into c number stochastic differential equations that can be solved numerically. This requires two complex stochastic fields Ψ i (x, t) and Φ i (x, t) , corresponding to the oper-atorsΨ i (x, t) andΨ † i (x, t), respectively. Ensemble averages of the stochastic fields over a large number of trajectories . . . st correspond to quantum mechanical ensemble averages of normally ordered operator moments. For example The stochastic differential equations governing the quantum dynamical evolution under the Hamiltonian (1) are, in a rotating frame: Here, the noise terms ζ j (j = 1, ..., 10) are real, independent Gaussian noises with ζ j (x, t) st = 0 and the following nonzero correlations We note that the Hamiltonian (1) conserves the total number of atomic particleŝ whether they are counted as bound pairs in the molecular state or as free atoms,N i (t) = dxΨ † i (x, t)Ψ i (x, t) (i = 0, 1). Since we are interested in dissociation of a pure molecular BEC, with no atoms present initially, thenN is a known constant operator given byN = 2N 0 (0), wherê N 0 (0) is the total initial number of molecules. Expectation values of these operators will be denoted by removing the hats, e.g. N i = N i .

III. ANALYTIC SOLUTIONS IN THE UNIFORM SYSTEM
The simplest treatment of dissociation resulting in analytic solutions is achieved under an undepleted molecular field approximation in a uniform system [10,13]. The approximation is valid for short dissociation times during which the total number of atoms produced is only a small fraction ( 10%) of the initial number of molecules.
In this treatment we consider a uniform molecular BEC in a coherent state with a constant amplitude Ψ 0 = √ ρ 0 (where ρ 0 is the uniform density), which we assume is real without loss of generality. The condensate is contained within a cubic box of side L (volume = L 3 ) with periodic boundary conditions, extending from −L/2 to L/2 in each dimension. The dissociation coupling χ is turned on suddenly, and subsequently assumed to be constant.
The Heisenberg equations for the atomic field opera-tors are then where g ≡ χ √ ρ 0 is the coupling and ∆ eff ≡ ∆ + (U 01 − U 00 /2)ρ 0 is the effective detuning that takes into account the initial mean-field energy shifts due to atom-molecule and molecule-molecule interactions. We note that in our simulations we will be interested in cases with large bare detuning |∆| ≫ |U 01 − U 00 /2|ρ 0 so that the mean-field energy shifts are negligible.
Introducing creation and annihilation operatorsâ k (t) andâ † k (t) for the plane-wave atomic modes with momentum k and commutation relations [â k (t),â † k ′ (t)] = δ k,k ′ (see Appendix A), the Heisenberg equations (5) transform into These have the following solution where the coefficients are and we have introduced ∆ k ≡ k 2 /(2m 1 ) + ∆ eff and g k ≡ [g 2 − ∆ 2 k ] 1/2 , where k = |k|. The coefficients A k (t) and B k (t) satisfy |A k (t)| 2 − B 2 k (t) = 1. Given the initial vacuum state of the atomic field we can calculate any expectation values of the field operators. In particular, we find that the only nonzero secondorder moments are the particle number per mode and the pairing field per mode: wheren k (t) =â † k (t)â k (t) andm k (t) =â k (t)â −k (t) are the respective operators. All other second-order moments are equal to zero. Higher-order moments factorize according to Wick's theorem and can be expressed in terms of the above second-order moments. From |A k (t)| 2 − B 2 k (t) = 1, we find in addition that n k (t) and m k (t) are related by Using these solutions, the total number of atoms produced is The treatment of an infinite system in free space is achieved as usual by taking the limit of L → ∞ (∆k → 0) and transforming the mode creation and annihilation operators into continuous Fourier transformŝ a(k) andâ † (k) that satisfy the commutation relation [â(k),â † (k ′ )] = δ(k − k ′ ) (see appendix A). In this case, the solutions for the normal and anomalous moments are A. Atom-atom correlations We are interested in the atomic correlations resulting from momentum conservation in the molecular dissociation. These may be quantified using a number of different, but related, density-density correlation functions. Furthermore, these measures may correlate fields at points in position space or momentum space. Since the fundamental correlation is between momenta, the momentum correlation functions are simpler to interpret. However, experimentally it is easier to measure spatial correlations after time-of-flight expansion of the cloud. In the present paper we will concentrate on analyzing atomic correlations in momentum space. Spatial correlations after expansion will be studied in a future work [50].
Starting from the atomic vacuum, and assuming no initial momentum spread in the molecular BEC and no scattering, molecular dissociation produces pairs of atoms with equal and opposite momenta. In the short time limit this is well approximated by the analytic solutions of the previous subsection. The strength of the atomic correlations may be quantified via Glauber's second-order correlation function [49] g (2) It is defined in terms of normally ordered operator products and is normalized so that it is dimensionless. The pair correlation describes the ratio of the probability of joint detection of pairs of atoms with k and k ′ to the product of probabilities of independent atom detection events at k and k ′ . For example, g (2) = 1 for uncorrelated atoms, g (2) = 2 for thermally bunched atoms, and g (2) > 2 indicates super-thermal bunching. Note that the conditions for an expansion using Wick's theorem are met by the solutions for the uniform system under the undepleted molecular field approximation, Eq. (6). In this case, the pair correlation function g (2) (k, k ′ , t) can be factorized and be expressed in terms of products of the second-order normal and anomalous moments â † k (t)â k ′ (t) and â k (t)â k ′ (t) . For example, in the case of k ′ = −k, the pair correlation function is where we have taken into account that n k (t) = n −k (t) , â † k (t)â −k (t) = 0 and used Eq. (11). For other pairs of momenta the correlation function is obtained in a similar way, with the overall result given by As we see, the correlation function for atom pairs with equal but opposite momenta is the strongest [except for the special case of k = k ′ = 0]. Quantitatively, strong superbunching in g (2) (k, −k, t) is achieved for low mode populations, n k (t) ≪ 1. In this regime, the numerator in Eq. (15) is approximately proportional to the mode population n k (t), while the denominator is the product of populations, and therefore g (2) As the mode occupancies grow with time, the pair correlation approaches the level of thermal bunching Thus, at high densities Glauber pair correlation becomes a less sensitive measure of the correlation strength. This is despite the fact that the correlation between the atoms with equal and opposite momenta is still maximal at any given density, in this analytically soluble model. There is perfect squeezing of particle number-difference fluctuations below the shot-noise level. This is quantified via the normalized variance where ∆Â =Â − Â is the fluctuation and the colons :: indicate normal ordering of the creation operators before the annihilation operators. This definition uses the conventional normalization with respect to the shot-noise level resulting from a Poissonian number probability distribution, such as for a coherent state. This is given by the sum of the mean occupation numbers, n k (t) + n k ′ (t) . Variance smaller than one, V (k, k ′ , t) < 1, implies reduction of fluctuations below the shot-noise level and is due to correlation between particle number fluctuations in the k and k ′ -modes, while V (k, k ′ , t) = 1 for uncorrelated modes.
For the solutions of Eq. (6) the number-difference variance for equal but opposite momenta is given by V (k, −k, t) = 0 at all times, i.e., for all occupation numbers n k (t) = n −k (t). This is the ideal case, implying perfect correlation between fluctuations inn k (t) andn −k (t) and corresponding to 100% squeezing below the shotnoise level.
We note that for two equally occupied modes, the variance V (k, −k, t) and the second-order correlation function g (2) where we have used g (2) (k, k, t) = g (2) (−k, −k, t) due to the spherical symmetry of the problem. Therefore, perfect noise reduction, V (k, −k, t) = 0, implies that the maximum degree of correlation is g Using the thermal level of autocorrelation that occurs in the analytic solutions, Eq. (17), which coincides exactly with the analytic correlation, Eq. (17).
In a similar way, one can show that the maximum (ideal) degree of correlation Eq. (20) corresponds to the maximum anomalous moment m k (t) which again coincides with the solution, Eq. (11), of the present analytically soluble model. This result follows immediately from the Wick-factorized expression for g (2) (k, −k, t), Eq. (16), assuming â † k (t)â −k (t) = 0. More generally, i.e., for nonideal cases, the numberdifference variance V (k, −k, t) would be larger than zero and can serve as a sensitive measure of the correlation strength between particle number fluctuations in different modes. In the nonideal cases, V (k, −k, t) > 0 also implies that the following inequalities hold: provided that one still has g (2) (k, k, t) = 2 and â † k (t)â −k (t) = 0. Quantifying the strength of correlations via the number-difference variance V (k, −k, t) is especially useful at large mode occupancies n k (t) ≫ 1, when the pair correlation approaches g (2) (k, −k, t) → 2 and therefore becomes a less sensitive measure of the correlation strength.

IV. NUMERICAL RESULTS AND DISCUSSION
Our goal in the analysis of atom correlations in molecular dissociation is to go beyond the undepleted molecular field approximation and to treat the quantum dynamics for realistic nonuniform systems. Analytic solutions for these cases are no longer available, and we will resort to numerical solution of the positive-P equations (3), which describe all quantum effects from first principles.
Using the positive-P method, we are able to extend the analytic results of the previous section in three ways, to include: (a) the depletion of the molecular condensate; (b) s-wave scattering interactions; and (c) nonuniform condensates.
In the first two cases, we will consider uniform systems, while the nonuniform case will be simulated within the undepleted molecular field approximation. The reason for doing this is to give a precise quantitative description of each of these effects and to assess their relative importance in altering the atom-atom correlations. In addition, the results obtained in each case can serve as a benchmark for other approximate numerical methods. An example is the pairing mean-field method of Ref. [12] which treats the molecular field depletion at the level of a coherent state, while the dynamics of the atomic field is treated via the normal and anomalous densities. The range of validity of this method will be examined below using the exact results of the present treatment.
The positive-P numerical method is implemented on a uniform spatial lattice with periodic boundary conditions, and the continuous fields are therefore represented by discrete mode amplitudes. The momentum space lattice is reciprocal to the spatial lattice, with ∆k = 2π/L being the lattice spacing and L the length of the spatial domain in each dimension.
Since the numerical method is able to treat both uniform systems in a finite box as well as nonuniform systems in free space, the momentum-space field amplitudes are commonly treated via the lattice-discretized momentum componentsâ(k, t) andâ † (k, t), which correspond to the continuous Fourier transforms (see Appendix A) in the limit ∆k → 0 (L → ∞). (The correspondences between the operator Fourier components and those of the stochastic fields are outlined in Appendix B). Accordingly, for any finite computational lattice one can formally identify a set of momentum modes described by creation and annihilation operatorsâ Similarly, the (dimensionless) mode populations n k (t) and the pairing fields per mode m k (t) are related to the normal and anomalous densities n(k, t) = n(k, t) = â † (k, t)â(k, t) and m(k, t) = m(k, In the continuous limit, these should be understood as corresponding to the expectation values of the respective operators defined as integrals over the momentum-space where ∆k is small enough so thatn(k ′ , t) andm k (k ′ , t) under the integrals do not vary much within the integration volume and can be replaced byn(k, t) andm(k, t).
The total number of atoms is The normalized second-order correlation function is now defined via which is equivalent to Eq. (15), provided thatâ(k, t) and a k (t) are defined on the same momentum-space lattice. Finally, we define the normalized variance of the particle number differencen k (t) −n k ′ (t), where the particle number operators are to be understood according to Eq. (27). The subscript v in Eq. (31) signifies the fact that in the continuous model describing an infinite system, the variance is defined for the counting volume v = (∆k) 3 , and therefore Eq. (31) is equivalent to Eq. (18). Accordingly, the results for V v (k, k ′ , t) depend on this volume; unlike the pair correlation function g (2) (k, k ′ , t), the scaling of the numerator and of the denominator with respect to (∆k) 3 does not lead to cancellation. The explicit dependence on (∆k) 3 becomes evident if we rewrite Eq. (31) via the densitiesn(k, t), which are independent of (∆k) 3 : n(k, t) + n(k ′ , t) .
(32) As we see, for a given negative value of the normally ordered variance : [∆ (n(k, t) −n(k ′ , t))] 2 : , the degree of squeezing below the shot noise level degrades (V v becomes closer to one) as the counting volume is decreased.
should not be confused with the requirement that physical observables should not depend on the choice of the lattice spacing ∆k. Rather, the calculated variance for a given computational grid corresponds physically to fluctuations in the difference of particle numbers in the lattice volume element (∆k) 3 . For a given density, the number of particles in a smaller volume is lower and the correlations between number-difference fluctuations becomes weaker. Stronger correlation signals can be obtained via the procedure of binning (see below), in which case the counting volume is chosen to be larger than (∆k) 3 .
To further clarify this point we note that the unnormalized variance of the continuous density-difference fluctuations, for k = k ′ , includes (after rewriting it in terms of the normally ordered operator products) a term proportional to a deltafunction δ (3) (0), which is the shot-noise level. While strong quantum correlation and the reduction of densitydifference fluctuations below the shot-noise level is generically due to negative values of the normally ordered variance, the degree of squeezing below the delta-function shot-noise is not well defined via Eq. (33). On the other hand, operating with the variance of fluctuations between the particle numbers, as in Eq. (31), makes the definition of the shot-noise unambiguous and therefore the degree of squeezing is well defined. Albeit, the squeezing now depends on the counting volume (∆k) 3 .

A. Parameter values
Taking a 87 Rb (m 1 = 1.43 × 10 −25 kg, m 0 = 2m 1 ) experiment as our specific example, we have performed simulations for the following set of parameters. The initial molecular BEC density was chosen to be ρ 0 = 5 × 10 19 m −3 for the uniform system. For the nonuniform systems, the same value was chosen to be the molecular BEC peak density. The size of the uniform system L was the same in all three spatial dimensions, with L = 1.38×10 −5 m, and the lattice grid contained 64 3 points. Thus, the lattice spacing in momentum space was ∆k = 2π/L = 4.56 × 10 5 m −1 , while the maximum cutoff momentum was k max = 1.46 × 10 7 m −1 , which is about twice the resonant momentum k 0 = 2m 1 |∆ eff |/ and exceeds the physically relevant range of momenta. The total initial number of molecules was N 0 (0) = ρ 0 L 3 = 1. With a spatially uniform molecular field and no depletion or two-body interactions, we validated our code by comparing the numerical results with those of the analytically soluble model, Eq. (5). The numerical code for solving the positive-P stochastic differential equations was implemented using the XMDS software package [65]. Within the sampling errors, the simulated momentum space mode populations n k (t) = n(k, t) (∆k) 3 and the correlation functions were in excellent agreement with the analytic results given by Eqs. (9), (10) and (17). The subsections below extend these results to include the molecular depletion, s-wave scattering interactions, and the treatment of nonuniform systems.
B. Role of the molecular field depletion

Total particle numbers and atom-atom recombination
Here, we simulate Eqs. (3) for a uniform system and all s-wave scattering interactions set to zero (a ij = 0, and hence ∆ = ∆ eff ). Including the molecular field in the dynamics leads to its depletion and a corresponding reduction in the number of atoms produced compared to the analytically soluble case.
In Fig. 1 we plot the total fractional number of molecules and atoms, N 0 (t)/N 0 (0) and N 1 (t)/2N 0 (0), as a function of time. The dashed line shows the total atom number in the undepleted molecular field approximation, Eq. (12), for comparison. We see that the undepleted molecular-field result agrees well with the exact quantum result for time durations resulting in the conversion of less than 10% molecules. In the present example, 10% conversion corresponds to 2.6 × 10 4 atoms produced [N 1 (t)/2N 0 (0) ≃ 0.1] and occurs at t/t 0 ≃ 1. At this time the discrepancy between the exact numerical and the approximate analytic results is ∼ 9%. For longer times the discrepancy increases and the approximate analytic result eventually produces an exponentially growing output. This unphysical behavior is an obvious consequence of the fact that the undepleted molecular field approximation is no longer valid.
Physically, the total number of atoms produced saturates due to the depletion of the molecular condensate containing a finite number of molecules to start with. At the same time, the dynamics of dissociation -after the initial spontaneous regime -is affected by the reverse process of recombination of atom pairs into molecules. The dashed-dotted line in Fig. 1 shows the total number of molecules outside the condensate mode k = 0. Since the dissociation in this uniform system starts with all molecules being initially in the condensate mode, the population of the k = 0 modes can only occur (in the absence of elastic collisions, a ij = 0) due to atom-atom recombination. Note that the recombination can in principle involve atom pairs with any momenta, and not necessarily those with equal and opposite momenta. This has the effect of reducing atom-atom correlations in the long time limit (see below), hence the term "rogue association" which we will use.
The population of the molecular noncondensate modes implies effective heating of the molecular gas and can be regarded as the reverse of the rogue dissociation known to occur in the opposite process of association of an atomic BEC into a molecular BEC [62,66]. In the example of Fig. 1, 18% of the remaining molecules at t/t 0 ≃ 2.5 are outside the condensate mode. At the same time, the total number of molecules remaining is about 45% of the initial number.

Atomic momentum distribution
In Fig. 2(a) we plot a slice through the origin of the atomic momentum density distribution n(k, t) at t/t 0 = 1. We see a clear ring structure corresponding to the distribution of atoms around the surface of a sphere of radius k 0 . Similar ring structures have been observed in the time-of-flight spatial column densities of atoms produced in dissociation experiments using 87 Rb 2 and 40 K 2 molecules [6,8].
To simplify the presentation and comparison of the results we make use of the spherical symmetry of the problem and also plot angle-averaged distributions n(k, t) ϕ,θ and n k (t) ϕ,θ [see Fig. 2(b)] as a function of the absolute momentum |k|. Here, the brackets . . . ϕ,θ refer to the procedure of averaging of the quantum me- chanical expectation values over the angles ϕ and θ in a spherical coordinate system. Comparison of the exact numerical result for the angleaveraged population distribution n k (t) ϕ,θ and the respective analytic result, Eq. (9), shows that the discrepancy at t/t 0 = 1 is less than 8%, with the peak value of n k (t) ϕ,θ at |k|/k 0 = 1 being reduced from the undepleted result of sinh 2 (1) = 1.38 to 1.27. For longer time durations the discrepancy increases. The small oscillatory structure seen in the wings of the distribution functions (for both the exact numerical and analytic curves) is additional evidence that the sampling error due to stochastic averaging is small. This oscillatory behavior is characteristic of the highly detuned modes outside the spherical shell as these modes do not experience an exponential growth but rather oscillate at the spontaneous level. In terms of the approximate analytic solutions, the oscillations take place for the modes with k = |k| for which the gain coefficient g k ≡ [g 2 − ∆ 2 k ] 1/2 is pure imaginary. This occurs when ∆ 2 k = [ k 2 /(2m 1 ) + ∆ eff ] 2 > g 2 ,  (2) (k, −k, t) at t = t0. The corresponding density profile is shown in Fig. 2(a). The sampling errors due to stochastic averaging are small only within the region of 0.65 < |k|/k0 < 1.25, where the correlation function ranges between ∼ 2.8 and 6. (b) Angle-averaged correlation functions g (2) (k, k ′ , t) ϕ,θ (for k ′ = −k, k, k90) at t = t0 as a function of the absolute momentum |k| within the spherical shell 0.65 < |k|/k0 < 1.25. The solid line with the dots (right scale) is the respective angle-averaged variance Vv(k, −k, t) ϕ,θ .
in which case the sinh 2 (g k t) = sinh 2 (i|g k |t) term in Eq. (9) turns into sin 2 (|g k |t) thus producing the oscillations.

Atom-atom correlations
The pair correlation function for the atoms with opposite momenta, g (2) (k, −k, t), is shown in Fig. 3(a). Here, we plot a slice through the origin of the 3D correlation function at t/t 0 = 1. Due to the finite number of stochastic trajectories and the normalization of the correlation function to the product of mode occupancies, the stochastic ensemble average gives highly noisy results in the regions of kspace where the mode populations are small (n k 0.3). This is seen within the central (|k|/k 0 < 0.65) and outer (|k|/k 0 > 1.25) region of the slice plot where the sampling error is large [67]. Within the spherical shell of 0.65 < |k|/k 0 < 1.25, on the other hand, the mode populations are higher and the sampling errors are negligible. Within this region the results for g (2) (k, −k, t) scale according to the colormap shown on the right of Fig. 3(a). Figure 3(b) shows the angle-averaged correlation functions g (2) (k, k ′ , t) ϕ,θ within the spherical shell 0.65 < |k|/k 0 < 1.25, at t/t 0 = 1. As expected, for k ′ = −k the correlations are super bunched, g (2) (k, −k, t) ϕ,θ > 2, while for k ′ = k we see the thermal level of bunching g (2) (k, k, t) ϕ,θ = 2 characteristic of Gaussian statistics. The line showing an uncorrelated level of g (2) (k, k 90 , t) ϕ,θ = 1 is for the pair of momenta k and k ′ = k 90 , where k 90 corresponds to a 90 • degree rotation of the density distribution in the (k z, k y ) plane around the k z axis.
Apart for the fact that the pump depletion results in slightly decreased mode populations relative to the analytic result, the correlations follow closely Eqs. (17), except that n k (t) is now the actual, decreased, mode population. This observation is valid for at least the simulated time window of t/t 0 = 2.5, corresponding to about 55% conversion. This is expected since the correlations are due to momentum conservation, which is unaffected by the decreasing rate of conversion, which is the major effect of the molecular depletion. Thus, we find that the pair correlation for equal and opposite momenta is still well approximated by Since n k (t) ϕ,θ is slightly less than in the undepleted molecular field approximation, the pair correlation Eq. (34) is accordingly higher. In Fig. 3(b) the result of Eq. (34) is plotted by the dotted line and is almost invisible as it follows closely the actual numerical value calculated using Eq. (30). The fact that the correlations can be approximated by Eq. (34) shows that Wick factorization continues to approximately hold for time durations corresponding to ∼ 55% conversion. We recall that the ingredients of this result are the vanishing off-diagonal normal moments, â † k (t)â −k (t) = 0, and the relationship of Eq. (11) between the normal and anomalous moments. These ingredients continue to hold in the present example. Figure 4 shows the temporal behavior of the angleaveraged pair correlation functions g (2) (k, k ′ , t) ϕ,θ at |k|/k 0 = 1 corresponding to the peak momentum in the density distribution. The approximate result for g (2) (k, −k, t) ϕ,θ , given by Eq. (34), is plotted with a dotted line, but is coincident with the actual numerical value shown by the solid line.

Number-difference variance
As we mentioned earlier, the pair correlation function becomes a less sensitive measure of the correlation strength when the mode occupation numbers increase. In particular, as n k (t) increases with time the pair correlation approaches g (2) (k, −k, t) → 2 which coincides with the level of thermal bunching for the autocorrelation function. In this regime, the variation in g (2) (k, −k, t) is no longer significant. As a result, distinguishing between strong and weak correlations becomes problematic. A more sensitive measure of the correlation strength is provided by the variance of the particle numberdifference fluctuations between the correlated modes, Eq. (18). The angle-averaged number-difference variance for the modes with equal but opposite momenta, Figs. 3(b) and 4 by solid lines with dots. As we see from Fig. 3(b), the variance is suppressed well below the shot-noise level V v (k, −k, t) ϕ,θ < 1 for the entire spectral range.
Similarly, the variance V v (k, −k, t) ϕ,θ shows subshot-noise fluctuations in the entire simulated time window (see Fig. 4). It increases above the ideal result of V v (k, −k, t) ϕ,θ = 0 towards the end of the simulated time window. In particular, at t/t 0 = 2.5 it is given by V v (k, −k, t) ϕ,θ ≃ 0.22 which corresponds to 78% squeezing below the shot-noise level. The reduction of squeezing is attributed to atom-atom recombination into molecules with nonzero momenta, k = 0 (see Fig. 1). Since this process involves pairs of atoms with nonopposite momenta, a certain fraction of the atoms is left without their correlated partners of opposite momentum. As a result, the overall strength of the correlation signal for the atoms with equal and opposite momenta is reduced.
The reduction of the correlation strength due to atomatom recombination becomes a sizeable effect after dissociation times (t/t 0 2) corresponding to more than ∼ 35% conversion. At very short times (less than 10% conversion) our exact numerical results confirm that the predictions of the analytically soluble model with the undepleted molecular field are a good approximation. At longer times the correct treatment of the system requires the inclusion of the molecular field depletion.
As shown in Ref. [12], the molecular depletion can be taken into account using a pairing mean-field method which is simpler to implement numerically than the present (exact) positive-P method. In the mean-field aij (nm) Uij (×10 −17 m 3 /s) tsim/t0 (t0 = 0.  method, the molecules are described as a coherent field in the condensate mode b k=0 (t) while the atoms are treated at the level of diagonal normal and anomalous populations, â † k (t)â k (t) and â k (t)â −k (t) . While imposing Wick's factorization scheme and assuming perfect correlation between the atoms with opposite momenta at all times, this treatment does not allow for atom-atom recombination into molecules outside the condensate mode. Accordingly, the predictions of the theory are a good approximation only for dissociation durations resulting in no more than ∼ 35% conversion, as can be seen from the comparison with the present exact results.
The pairing mean-field method performs better in the long time limit if the only observables of interest are the total particle numbers. Its performance can be expected to improve for calculating the correlation functions if it is expanded to include the molecular noncondensate modes and the off-diagonal terms in the atomic normal and anomalous densities.
Before turning to the analysis of atom correlations in the presence of s-wave scattering interactions, we note that the positive-P simulations without the two-body collisions fail for durations longer than t/t 0 ≃ 2.5. At this time about 55% of molecules have been converted into atoms. The failure of the simulations for longer durations is due to the generic boundary term problem that occurs in the positive-P method [55,56]. Despite this, our results are a significant improvement and extend beyond the regime of validity of the simple undepleted molecular field approximation or the pairing mean-field theory. A possible route to increase the simulation time is to use stochastic gauges [55,68], which is, however, beyond the scope of the present paper.

C. Role of the s-wave scattering interactions
Here, in addition to simulating the molecular field depletion we include the s-wave scattering interaction terms. Table I shows the different cases considered, together with the simulated time window t sim . The simulations are shorter than in the absence of the s-wave scattering since the positive-P boundary term problem [55,56] becomes more severe for quartic interactions.
The simulation times were chosen to be less than those for which the sampling error was observed to grow rapidly [55].
The shortest t sim is for the case of a 11 = 5.3 nm which is the natural background scattering length for 87 Rb atoms. Longer t sim are achieved with smaller values of a ij . The scattering length can in principle be tuned during dissociation by applying a magnetic field in the vicinity of a Feshbach resonance. In this case the dissociation itself would have to be invoked independently, using optical Raman transitions. While this is unlikely to simultaneously result in the specific values of a ij quoted in Table I, our main goal is the understanding of the role of each physical process separately -by quantitatively estimating their relative importance and the potential disruptive effect on atom-atom correlations.
The results of simulations of the cases of Table I show that the momentum distribution and the pair correlation for atoms with opposite momenta remain (within the simulated time windows) very close to the values obtained with no s-wave interactions. For example, in the case of a 11 = 2 nm, the peak occupancy reduces from 1.70 to 1.62 at t/t 0 = 1.25, while the corresponding pair correlation still follows Eq. (34). The reduction in the mode population is attributed to the mean-field energy shift due to atom-atom interactions, which dynamically changes the phase matching condition and hence results in a drift of the effective detuning ∆ eff as the populations evolve. As a result, the resonant momentum also drifts causing less efficient conversion into the initially resonant mode. Nevertheless, since the bare detuning |∆| in our examples is still much larger than all mean-field phase shifts U ij Ψ † j (x, t)Ψ j (x, t) , the resonant momentum is still well approximated by k 0 ≃ 2m 1 |∆|/ .
The expected reduction of the actual correlation strength due to s-wave collisions and the resulting redistribution of the atomic momenta is best revealed through the variance of the number-difference fluctuations, V v (k, −k, t) ϕ,θ . For example, in the case of a ij = 2 nm (which is the worst case, among the simulations surviving up to durations of t/t 0 = 1.25), the variance at the resonant momentum k 0 increases from the collisionless result of 0.024 to 0.032 at t/t 0 = 1.25. This implies the reduction of number-difference squeezing from 97.6% to 96.8%. The overall effect is, however, smaller than decorrelation due to atom-atom recombination discussed earlier, which in this example is responsible for the reduction of squeezing from 100% down to 97.6%.

D. Nonuniform systems and mode mixing
We now turn to the analysis of atom-atom correlations in nonuniform systems corresponding to dissociation of realistic trapped BECs. The net effect of inhomogeneity is mode mixing, which can dramatically affect the correlations. The quantitative details are best understood in the undepleted molecular field approximation. In this case, the atom-atom recombination and s-wave scattering interactions are absent, and any reduction in atom-atom correlation strength -compared to the uniform case -is due to mode mixing.
Hence, the molecular field amplitude Ψ 0 (x, 0) = ρ 0 (x) [where ρ 0 (x) is the density] can be absorbed into an effective coupling g(x) ≡ χ ρ 0 (x). This leads to the same equations of motion for the atomic field operator Ψ 1 (x, t) as Eq. (5), except that the coupling constant g is now a function of x. Introducing the Fourier transformg(k) = (2π) −3/2 dxg(x) exp(−ik · x), the equation of motion for the atomic Fourier componentâ(k, t) is As we see,â(k, t) couples to the range of momentum components k ′ in the conjugate fieldâ † (k ′ , t), which we refer to as mode mixing. This is in contrast to the uniform case, whereg(k + k ′ ) is the delta functioñ g(k + k ′ ) = (2π) 3/2 gδ(k + k ′ ), so thatâ(k, t) couples only to the conjugate field at the opposite momentum a † (−k, t).
The overall effect of the mode mixing is to correlate the atoms with momentum k not only with the opposite momentum −k, but also with the atoms distributed in a range of momenta around −k. As a result, for a given strength ofg(k) the pair correlation between k and −k is expected to be reduced compared to the uniform case.
To characterize mode mixing quantitatively we simulate the positive-P stochastic differential equations that are equivalent to Eqs. (35). The initial molecular BEC is assumed to be formed in a spherically symmetric harmonic trap with equal trap oscillation frequencies ω ≡ ω x,y,z . It is assumed to be in a coherent state, with the density distribution given by the Thomas-Fermi parabola ρ 0 (x) = ρ 0 (0)(1−|x| 2 /R 2 T F ) for |x| < R T F , and ρ 0 (x) = 0 for |x| ≥R T F . Here, ρ 0 (0) is the peak density, while R T F = 2 U 00 ρ 0 (0)/(m 0 ω 2 ) is the Thomas-Fermi (TF) radius.

Case 1: weak inhomogeneity
The choice of parameter values in this set of simulations is specifically targeted to give a situation which is directly comparable with the results of the previous uniform system. The uniform results for a cubic box of side L u can be regarded as an approximation to a realistic nonuniform system if L u is matched with the characteristic size ∼ 2R T F of the initial molecular BEC, while the uniform density ρ 0 is matched with the peak density ρ 0 (0). More specifically, we choose R T F = (8π/15) −1/3 L u and ρ 0 (0) = ρ 0 , which gives the same initial total number of molecules in both cases, N 0 (0) = (8π/15)ρ 0 (0)R 3 T F = ρ 0 L 3 u .  Thus, in the present example (with the results shown in Fig. 5, column 1) the peak density is ρ 0 (0) = 5 × 10 19 m −3 , while the TF radius is R T F = 1.16 × 10 −5 m, which would correspond to a molecule-molecule scattering length of a 00 = 2 nm in a relatively weak trap with oscillation frequency ω/2π = 8 Hz. This gives the total initial number of molecules N 0 (0) = 1.3 × 10 5 . The simulations are performed on the spatial domain of L = 2.76 × 10 −5 m in each dimension, which is twice the size of the earlier uniform system L u . The lattice grid contained 128 3 points. Thus, the lattice spacing in momentum space is now ∆k = 2π/L = 2.28 × 10 5 m −1 (half that in the uniform system), while the maximum cutoff momentum is the same k max = 1.46 × 10 7 m −1 .
From Fig. 5 (column 1) we see that the atomic density in momentum space (b) and the total number of atoms (c) follow closely the predictions of the size-matched uniform system shown in Figs. 2(b) and 1(a). At the same time, the mode populations are about eight times smaller, which is simply due to the fact that the simulation of the nonuniform system is performed within a spatial domain which is twice as large, L = 2L u . Accordingly, the respective counting volumes (∆k) 3 = (2π/L) 3 and (∆k u ) 3 = (2π/L u ) 3 are different by a factor of eight, so that the approximately equal densities n(k, t) ≃ n u (k, t) give different mode populations n k (t) = n(k, t)(∆k) 3 and n (u) k (t) = n u (k, t)(∆k u ) 3 . Scaling the uniform result for n u (k, t) with respect to the counting volume used in the nonuniform system gives n u (k, t)(∆k) 3 that approximates well the actual calculated value of n k (t) shown in Fig. 5.
While the results for the momentum-space density and the atom number are in good quantitative agreement with the respective results of the size-matched uniform system, the same is not true for the correlation function g (2) (k, −k, t). The angle-averaged pair correlation g (2) (k, −k, t) ϕ,θ at |k| = k 0 in the present nonuniform system is shown in Fig. 5(e) by the full line. Before comparing this result with that of the uniform system [given by the idealized result of g (2) max (k, −k, t) = 2 + 1/n k (t)], we note that the correct interpretation of the normalized pair correlation is in the excess of the probability of joint detection of pairs of atoms at k and −k relative to that probability in an uncorrelated state. Operationally, this corresponds to determining the number of particles in a "detection" volume (∆k) 3 . Therefore, the pair correlation g (2) (k, −k, t) ϕ,θ for the nonuniform system can be compared with the idealized result g (2) max (k, −k, t) = 2 + 1/n k (t) provided that n k (t) is evaluated for the same counting volume (∆k) 3 as in the nonuniform simulation.
The idealized result g (2) max (k, −k, t) obtained in this way is shown in Fig. 5(e) by the dotted line. We immediately see that the strength of the pair correlation for the nonuniform system g (2) (k, −k, t) ϕ,θ is substantially lower than g (2) max (k, −k, t). For comparison, the same quantity calculated numerically in the uniform system was well approximated by the idealized result [see discussions of Figs. 3(b) and 4], where it was almost invisible behind the solid line for g (2) (k, −k, t) ϕ,θ .
We note that the correlation between the atoms with opposite momenta is still stronger than in an uncorrelated state, g (2) (k, −k, t) ϕ,θ > 1. We have also calculated the pair correlation functions g (2) (k, k, t) ϕ,θ and g (2) (k, k 90 , t) ϕ,θ . As in the uniform case, these were given by 2 and 1, respectively, implying that the mode mixing has negligible effect on these correlations.
Another measure of the reduction of the correlation strength between the atoms with opposite momenta can be obtained using the variance of the number-difference fluctuations V v (k, −k, t), Eq. (31). The angle-averaged variance V v (k, −k, t) ϕ,θ at |k| = k 0 as a function of time t is shown in Fig. 5(f), and is given approximately 3 is the elementary volume element of the computational grid. This corresponds to 28% squeezing of fluctuations below the shot-noise level, which is a significant degradation compared to the ideal case of 100% squeezing obtained in the uniform systems in the absence of molecular depletion or s-wave scattering interactions. The reduction of the correlation strength due to mode mixing is a much stronger effect than decorrelation due to the atom-atom recombination and s-wave scattering, analyzed in Secs. IV B and IV C.
Finally, in Fig. 5(g) we plot the correlation coefficient as a function of time t, at |k| = k 0 . According to the inequality (24), this gives a useful measure of the departure [C k (t) < 1] of the anomalous mode population from the maximally correlated, ideal result of Eq. (21) corresponding to C k (t) = 1. In addition, the calculated anomalous population m k (t) gives us a direct verification of the validity of Wick's factorization, according to which the pair correlation g (2) (k, −k, t) in the present example can be evaluated as in Eq. (22). We have checked that the offdiagonal normal moment â † k (t)â −k (t) gives a negligible contribution, | â † k (t)â −k (t) | 2 /n 2 k (t) ≪ 1.

Cases 2 and 3: strong inhomogeneity
Columns 2 and 3 in Fig. 5 show the results of simulations of a nonuniform system in a tighter trap, with ω/2π = 19.3 Hz and R T F = 4.82 × 10 −6 m. The total initial number of molecules is N 0 (0) = 9.4 × 10 3 , while all other parameters are as in the previous example. Thus, in the present case the inhomogeneity in position space is stronger. Conversely, the effective couplingg(k) is broader in momentum space and therefore the mode mixing has a stronger effect on the reduction of correlations between the atom pairs with opposite momenta.
The results presented in column 2 of Fig. 5 are obtained on the same computational grid as in column 1, i.e., using L = 2.76 × 10 −5 m and 128 3 lattice points. Therefore the quantities depending on the counting volume v = (∆k) 3 -through the mode occupation number n k (t) = n(k, t)(∆k) 3 -are directly comparable.
Apart from the obvious reduction in the momentum space density and total atom number produced as a function of time, the fraction of atoms relative to the initial number of molecules is almost the same as in the example of column 1. The density-density correlation function does not reveal significant quantitative change either. However, as we explained in the previous subsection, the quantitative aspect of the reduction of the correlation strength should be assessed relative to the idealized result for g (2) max (k, −k, t) = 2 + 1/n k (t). In the present example, g (2) max (k, −k, t) is off the scale of Fig. 5(e), signaling a much more dramatic reduction in the correlation strength than before. This is further seen through the number-difference variance, V v (k, −k, t) ϕ,θ ≃ 0.98, implying only 2% of squeezing below the shot noise-level. Similarly, the correlation coefficient C k (t) is much lower than the ideal result of C k (t) = 1.
In column 3, Fig. 5, we show the results of simulation of the same physical system as in column 2 except that the simulation is performed on a smaller (64 3 ) computational grid, with half the length L = 1.38 × 10 −5 m. Thus, the elementary volume element v = (∆k) 3 , where ∆k = 2π/L, is eight times larger and therefore the results depending on the counting volume v are scaled by a factor of eight. For example, the number-difference variance is now given by V v (k, −k, t) ϕ,θ ≃ 0.84 (16% squeezing), demonstrating that the larger counting volume results in a stronger correlation signal between the particle number-difference fluctuations.

Binning
The comparison between cases 2 and 3 as a function of the counting volume (∆k) 3 is similar to the procedure of binning except that it is done implicitly -via the variation of the length L of the simulated spatial domain in each dimension. Since the minimum acceptable length must be larger than the characteristic size of the actual physical system, the relationship ∆k = 2π/L puts an upper bound on the volume of the elementary bin that can be treated in this implicit way.
To correlate the signals in larger counting volumes, the binning must be done explicitly. We introduce the binned number operator which is defined on a {K}-sublattice corresponding to the central momenta of the bins, where v(K) is the bin where the integration over the bin volume must precede the ensemble averaging. In Fig. 6 we present the results of simulations of the same physical system as in Fig. 5, column 1, except that the results are binned. The bin counting volume for calculating the occupation numbers n K (t) = n K (t) and the number-difference variance 3 . This is eight times larger than the elementary volume element of the computational grid (∆k) 3 .
As expected, the binned variance shows a stronger correlation signal. The degree of squeezing increases to ∼ 60% ( V v (K, −K, t) ϕ,θ ≃ 0.4, for |K| = k 0 ), compared to 28% squeezing in the respective unbinnned variance of Fig. 5(f), column 1. Similar results are obtained for the cases represented in columns 2 and 3 of Fig. 5.
In practice, similar issues arise in time-of-flight images due to the finite resolution of the imaging system, and binning was employed in the recent spatial correlation measurements of Ref. [8]. The correlation signal between fluctuations of the binned particle numbers was indeed found to be stronger for larger bins.

V. SUMMARY
To summarize, we have performed first-principles quantum dynamical simulations of dissociation of a BEC of molecular dimers into correlated atom pairs in three dimensions. The simulations are done using the positive-P representation method which allows us to obtain exact results for the atomic and molecular populations and atom-atom correlations in momentum space. We have analyzed the effects of molecular depletion, s-wave scattering interactions, and mode mixing in nonuniform condensates.
We find that the most useful measure of the strength of atom-atom correlations is given by the variance of atom number-difference fluctuations in a certain range of momenta around the central momenta of interest k and k ′ . The strongest correlation signal resulting in squeezing of particle number-difference fluctuations below the shot noise level is obtained for pairs of equal but opposite momenta, k ′ = −k. The degree of squeezing depends on the geometry of the system and is stronger for BEC samples that are spatially larger and more uniform.
The major source of atom-atom decorrelation compared to the ideal result (achievable in uniform systems and in the absence of molecular depletion) is mode mixing due to inhomogeneity of the molecular BEC. This suggests that a preferred experimental strategy to obtain a strong correlation signal is to start with a large, low density sample. The actual correlation strength depends also on the counting volume around the central momenta and can be increased by binning the signals into larger bins.
The next important source of decorrelation is the atom-atom recombination, which produces increasingly large number of molecules in the initially unpopulated noncondensate modes as the dissociation proceeds. In our example the fraction of these molecules is about 20% of the remaining total number once the overall conversion reaches 50%. At this stage, the reduction in numberdifference squeezing due to recombination is about 20%.
Finally, the simulated s-wave scattering interactions are found to have a less severe effect on atom-atom correlations at least for time durations corresponding to less than 10% conversion. For longer time durations, the simulations fail due to the limitations of positive-P method, unless the atom-atom scattering length is reduced below the typical background value.
To treat the uniform system in a cubic box of side L with periodic boundary conditions, we expand the field operatorsΨ 1 (x, t) in terms of the plane-wave momentum modesâ k (t):Ψ The inverse transforms arê a k (t) = 1 The mode creation and annihilation operators satisfy the usual commutation relation [â k (t),â † k ′ (t)] = δ k,k ′ , where δ k,k ′ is the Kroenecker delta function, k = (k x , k y , k z ) is the momentum (in wave-number units), and k i = (2π/L)n i (i = x, y, z, n i = 0, ±1, ±2, . . .).
The treatment of the infinite system in free space corresponds to taking the limit of L → ∞ (∆k = 2π/L → 0), together with converting the sums k ≡ n = nx ny nz into integrals according to ni → dn i = (L/2π) dk i . The corresponding continuous Fourier transforms are defined according tô with the inverse transforms of a(k, t) = 1 (2π) 3/2 dxΨ 1 (x, t)e −ik·x , (A7) The Fourier componentsâ(k) andâ † (k, t) satisfy the commutation relation [â(k),â † (k ′ )] = δ(k − k ′ ). When represented on a discrete computational lattice, the continuous Fourier components are approximated bŷ a(k, t) =â k (t)/ (∆k) and the continuous delta function by where ∆k = 2π/L is the mode spacing in each dimension. The Fourier componentsâ(k) (unlike the mode annihilation operatorsâ k ) have units; n(k, t) = â † (k, t)â(k, t) corresponds to the particle number density in momentum space, while n k (t) = â † k (t)â k (t) is the number of particles in the mode k.

APPENDIX B: FOURIER TRANSFORMS OF THE STOCHASTIC FIELDS
Here, we outline the correspondences between the Fourier transforms of the positive-P stochastic fields and those for the field operators. This is important for correct implementation of the routines aimed at the calculation of various operator moments in Fourier space.
We recall that the two independent stochastic (c number) fields Ψ 1 (x, t) and Φ 1 (x, t) are chosen here to correspond to the annihilation and creation operatorsΨ 1 (x, t) andΨ † 1 (x, t) according to the following operator correspondences:Ψ with similar relationships valid for the molecular field. The Fourier transforms of the stochastic fields give: Thus, the Fourier component β(k, t) corresponds tô a † (−k, t), rather than toâ † (k, t) as might have been expected. Therefore, care must be taken in assigning the operator correspondences in Fourier space.