Temporal and spatio-temporal correlation functions for trapped Bose gases

Density correlations unambiguously reveal the quantum nature of matter. Here, we study correlations between measurements of density in cold-atom clouds at different times at one position, and also at two separated positions. We take into account the effects of finite-size and -duration measurements made by light beams passing through the atom cloud. We specialise to the case of Bose gases in harmonic traps above critical temperature, for weakly-perturbative measurements. For overlapping measurement regions, shot-noise correlations revive after a trap oscillation period. For non-overlapping regions, bosonic correlations dominate at long times, and propagate at finite speeds. Finally, we give a realistic measurement protocol for performing such experiments.


I. INTRODUCTION
Correlations between density measurements in quantum gases are related to the underlying correlation functions of the quantum field. In cold-atom experiments which detect individual atoms extracted from a cloud or a beam, the correlations appear in the arrival times and positions of the detected atoms [1][2][3][4][5][6][7][8]. Atoms are more (bosons) or less (fermions) likely to be coincident than classical statistics would suggest.
By imaging a cloud after expansion and inspecting the fluctuations in the detected light, it is possible to infer correlations in momentum space [9]. This technique has been used to probe quantum many-body states such as the bosonic Mott insulator in a lattice [10] and fermion pairing [11]. Destructive in situ imaging has revealed spatial correlation information about one-dimensional quasicondensates (through both pair [12] and triple [13] correlations) and Fermi gases [14,15]. More sophisticated optical detection schemes are used to measure correlations in spin polarisation [16,17] rather than number density. Spatial density correlations can be inferred from fluctuations about the mean density [12] . Fluctuations are related to temperature and to thermodynamic susceptibilities like compressibility via fluctuation-dissipation theorems [18]. As susceptibilities tend to diverge close to phase transitions (e.g. Bose-Einstein condensation), so the measurable fluctuations become larger. However, fluctuation-dissipation theorems give information about the magnitude of fluctuations but not their length or time scales [19]. Naraschewski and Glauber [20] showed how the underlying first-order correlation function relates straightforwardly to observable density-density [21] * Current address: QUEST Institute, Physikalisch-Technische Bundesanstalt and Leibniz Universität Hannover, Bundesallee 100, 38116 Braunschweig, Germany † Correspondence should be addressed to r.nyman@imperial.ac.uk (second-order) correlations under the strong assumption that the system being measured is in the grand canonical ensemble at thermal equilibrium. In this paper we adapt their theory to include temporal correlations.
Many of the measurements of correlation functions in atomic gases were done by detecting atoms which had already been extracted from an expanding cloud or beam. Other experiments measured correlations between spatially separated regions of a cloud at one given instant or fluctuations from one experimental run to another. In the recent experiments of Guarrera et al. [5,7] atoms were extracted one by one from a single cloud and the temporal correlations of their arrival times measured. There, the extraction of the first atom and its detection had no direct effect on the subsequent behaviour of the cloud: the correlations were present independently of the measurement, and the measurement itself was in some sense weakly perturbative [22].
In contrast to previous experiments, we consider here density-density (not atom-atom) correlation measurements using light beams which pass through a single atom cloud at different times. The most general version of the measurement includes sampling the density at differing positions and times. These retarded correlations have not yet been observed in cold atoms.
Correlations persist over thermal timescales like /k B T at temperature T . A typical length scale is the thermal de Broglie wavelength: λ T = 2π 2 /M k B T for atoms of mass M . Taking a cloud of 87 Rb at 200 nK the typical correlation times and lengths (without timeof-flight expansion) are 40 µs and 0.4 µm. Measurements can happen faster than 10 µs [23], with spatial resolution down to about 1 µm, limited by photon collection and diffraction respectively. We therefore expect to be able to resolve correlations more easily in the time domain than in space.
In this article we first state our method for calculating correlation functions for thermal bosons and then show how optical measurements of atomic density relate to the correlation functions. We then specialise our discussion to the example of non-condensed Bose gases in harmonic traps. Finally, we discuss plausible experimental implementations using weakly-destructive measurements.

II. SECOND-ORDER CORRELATION MEASURES IN THE GRAND CANONICAL ENSEMBLE
Following the formalism of Ref. [20], but retaining time dependence wherever possible, we start by defining the atomic-field annihilation operator for an atom at position r and time t asΨ(r, t). Its Hermitian conjugate is the creation operator, and they follow the commutation relations for bosonic operators The function D depends on the Hamiltonian governing the system (see Appendix A). Density is an observable which is the expectation of the number density operator, n(r, t) = n(r, t) wherê n(r, t) =Ψ † (r, t)Ψ(r, t). The quantum expectation value is denoted with angle brackets, and for readability we will use both notations n(r, t) and n r,t for expectation of the density at a specified position and time.
We define normally-ordered, time-ordered (t ′ ≥ t), first-and second-order atom-atom correlation functions [20,21]: The first-order correlation function (with equal time arguments) G (1) is the spatial representation of the oneparticle density operator. Also, the density and correlation function are simply related: n(r, t) = G (1) (r, r, t, t).
The correlation between pairs of atoms is found in G (2) . An obvious local, second-order, density correlation measure, as distinct to an atom-atom correlation function, would be n r,tnr ′ ,t ′ − n r,t n r ′ ,t ′ . However, the operator in the first term is not Hermitian, and so this measure is not a good observable. We conjecture that the appropriate observable is simply the real part.
We suppose that the measurements to be made are in some sense very weak, so that it does not matter in what order the measurements are made. The observable should then be symmetric under exchange of coordinates, { r, t } and { r ′ , t ′ }. We note that when the variables are exchanged, the commutator in Eqn. (2) is complex conjugated, and thus the real part is a correct, symmetric operator. Therefore, the correlation measure that we will use is: C r,r ′ ,t,t ′ = 1 2 n r,tnr ′ ,t ′ +n r ′ ,t ′n r,t − n r,t n r ′ ,t ′ (5) The purpose of defining this density correlation measure is to account for measurements over finite volumes and times, as will be shown in section Section II B. Applying the definitions of G (1) and G (2) , and the commutation relations, we note that ensemble average product of non-Hermitian operators can then be expressed in terms of the correlation functions: Since G (2) is real: Re G For a non-interacting gas in a confining potential (with localised eigenfunctions) at thermal equilibrium, above the critical temperature for bosons, in the grand canonical ensemble, it can be shown that [20] where the second term is due to bosonic exchange symmetry. The correlation measure can be expressed in terms of the first-order correlation function G (1) only: The first term corresponds to counting (shot noise) number fluctuations. The second term comes from quantum exchange correlations.
A. The first-order correlation function G (1) for non-interacting, trapped, thermal bosons G (1) is conveniently expressed in a basis of the eigenstates of the trapping potential. First we write the field creation and annihilation operators in this basis set: where the state labelled i has energy ǫ i , spatial variation u i (r) and annihilation operatorâ i . The sum runs over all states. The expectation â † iâ j = δ ij N i , with N i being the thermal occupation number of the state, is given by the Bose-Einstein distribution, is the inverse temperature, and ζ = e βµ the fugacity with µ being the chemical potential. It follows that the first-order correlation function is: (11) For bosons, it turns out to be helpful to expand the thermal occupation factor in Eqn. (11) as a power series of the fugacity, of the form x/(1 − x) = k=1 x k . Combining exponentiated terms, we find: The inner sum is in the same form as the representation of the propagator, or Green's function [24], for the singleparticle HamiltonianĤ: Using this propagator with a complex time argument, the correlation function can be evaluated: At this point, we note that the expansion is effectively over increasing thermal occupation numbers, and that correlations with larger particle numbers will decay rapidly for non-condensed systems.
The fugacity can be found by integrating the density over all space to give the total atom number: If the atom number is to be fixed in simulation, then the fugacity must be found so that this equation is satisfied.

B. Correlation measurements over finite volumes
Imaging experiments measure the change of amplitude, phase or polarisation of a probe light beam due to the atom-light interaction [21]. In the limit of small optical density, the change to the light beam is proportional to the average atomic density seen by the beam as it passes through the cloud, necessarily integrating along line of sight. There will also be a finite resolution transverse to the beam due to the optical geometry, ultimately limited by diffraction. Therefore, any optical measurement of density covers a finite volume. Real measurements also take finite times to count sufficient photons. This finite detection volume and time will tend to smear out and reduce measured correlations.
The expectation of the time-averaged number of atoms in the probe beam during a measurement is [25]: with R a label for the central position of the probe beam for measurements around time τ . The probe beam has a spatial and temporal intensity profile I R,τ (r, t) which is normalised such that dt I R,τ (R, t) = 1 (dimensionless) at its spatial maximum [26]. Photons arriving at detectors will be used to infer this atom number. A general second-order correlation measurement includes a series of paired number measurements at two times, τ and τ ′ and with probe beams in modes centred at different positions R and R ′ , as shown schematically in Fig. 1. This measurement corresponds to measuring the product of densities, which is the only kind of correlation that can be measured using near-resonant light [21]. A good estimator for the correlations, the covariance, noted with the tilde to remind us that it's spatially integrated and averaged over a series of realisations, is: This finite-volume correlation estimator is the double integral of the local correlation measure, Eqn. (5), weighted by the probe beam intensities. This can in turn be expressed in terms of G (1) , as in Eqn. (9). We find: It is common to work with correlation estimators which are normalised by the standard deviations of each of the component measurements. In this case, we choose to approximate the standard deviations of measured numbers with Poissonian statistics in the limit of large number: Var (N ) = E (N ), where Var (X) is the sample variance.
The normalised correlation estimator is then:

III. LOCALISED CORRELATIONS IN BOSE GASES IN HARMONIC TRAPS
We now turn to the specific case of thermal bosons in a harmonic trap. Our task is to evaluate the local correlation estimator, C r,r ′ ,t,t ′ . The Hamiltonian for an atom of mass M in a trap of frequencies Ω x , Ω y and Ω z The eigenstates u lmn (r) are the Hermite polynomials with energies ǫ lmn = (lΩ x + mΩ y + nΩ z ), with zero-point energy subtracted for ease of notation. Since the Hamiltonian is a sum of three terms, one for each dimension, both the eigenstates and the propagator are separable: The propagator along the x direction, for the kth term in the expansion of the correlation function, is given by [24]: M With a real argument this propagator oscillates, so there will be damped oscillations in the correlation function from the complex argument. The correlation function is then a sum over products of partial propagators: Note that this does not mean that the correlation function is separable. By numerical evaluation of the fugacity constraint Eqn. (15): it is possible to choose the atom number, and then to iteratively estimate the fugacity ζ to make this equation consistent. 500 nK in a cigar-shaped trap of 1 kHz × 1 kHz × 20 Hz, along the centre of the trap at time zero. In Fig. 2, we see that the series Eqn. (12) converges rapidly for temperatures significantly above critical. Each term in ζ k corresponds to a thermal occupation number of the ground state by k particles. Since the temperature is above critical for this example, the probability of thermal occupation of large numbers of particles in the ground state is very small. The correlation function is significant over length scales around the de Broglie wavelength. For further discussion, we define the normalised correlation functions: which are independent of the number of particles in the trap and start at g (1) = 1 and g (2) = 2 for t = t ′ = 0 and r = r ′ .

B. Two-time correlations at one position
Turning our attention to the normalised, two-time correlation function, Fig. 3 shows a decay of the two-particle correlations on the time scale τ c = β (inset). The correlations revive and decay as particles oscillate in the confining potential, returning partially after the half the shortest oscillation period, i.e 500 µs. The correlations return completely on the longest time scale in the system, half the axial oscillation period (25 ms in this case). Any two particles will separate and then return to their original position after half an oscillation in the trap, so the two-particle correlation function revives. The real part of g (1) is negative for most of the time between the beginning and the first revival. This phase shift could in principle be observable in an two-path interferometer based on extracting atoms from the cloud, and then relatively delaying one of the paths. Beats between the different trap frequencies appear in the correlation function.  (2) for the same atom cloud as Fig. 2, but offset from the trap centre along the short direction by a distance x0. We see that the revival goes from twice per oscillation period to just once as position is moved away from the trap centre.
Inspecting g (2) at a fixed position but off-centre of the trap, Fig. 4 shows that the revival changes from occurring twice per oscillation period to just once. The distance off-axis required for this effect is approximately the de Broglie wavelength.   5 shows the temporal evolution of g (2) when the measurement at the second time τ ′ is done away from the trap centre. The correlation function peaks at later times at points further away from the centre. We can therefore define a speed at which the peak of the correlation function travels. It can be calculated by combining the typical correlation time scale τ c = β with the de Broglie wavelength λ T to make the correlation velocity v c = 2πk B T /M . We see from figure

IV. MEASUREMENTS OVER FINITE VOLUMES IN OBLATE HARMONIC TRAPS
We now apply the general form of Section II B to the specific case of harmonically trapped Bose gases probed by Gaussian-shaped beams. The co-ordinate system is defined in Fig. 1. We make the approximation that the Rayleigh range of the probe beam is much longer than the thickness of the cloud in the direction of propagation. This approximation makes it possible to separate the detection response function, into components I R,τ (r, t) = I X (x)I Y (y)I Z (z)I τ (t).
Measurements can be performed quickly compared to the typical correlation time scales, β, so we will approximate the temporal response function of the detection with an infinitesimal response time: I R,τ (r, t) = I R (r)δ(t − τ ). The separated factors of the spatial detector response functions are I X (x) = 1, I Y (y) = e −2(y−Y ) 2 /w 2 0 , and I Z (z) = e −2(z−Z) 2 /w 2 0 , where R = (X, Y, Z) and w 0 is the waist of the probe beam.
We consider a pancake shaped atomic cloud, which is thin along the line of sight in order to minimise the loss of coherence through integration along the probe axis; the other axes are kept loose. The trap contains 10 5 atoms at 100 nK and has trap frequencies of 413 × 20 × 9 Hz, much as in Ref. [5].
The full expression for the correlation estimator with a finite detection volume, Eqn. (18), requires a numerically intractable 8-dimensional integral. Since the detector response contains a δ-function we can execute the time integrals. The separability of both propagators and detector response can be used to express the 6-dimensional spatial integral as a product of two-dimensional integrals, reducing the numerical complexity.
Making use of the result of Eqn. (A4), we can now numerically evaluate Eqn. (18) as: The first term is due to quantisation of the matter field into atoms (shot noise). This term is dominant around the peaks of D(r, . This equation gives an upper bound for the shot noise contributions. The second term in Eqn. (25) is due to bosonic fluctuations and will be labelledC bs . It decays more slowly than the first term when increasing τ ′ − τ .

A. Two-position correlations at one time
The correlation measure taken for two short (compared to β) measurements of the number of atoms in the probe beam, at nearly equal times, but two positions is shown in Fig. 6. The number of atoms in the beam estimated from equation Eqn. (16) is 37 atoms. The measures are taken 10 µs apart since the numerical integrals do not converge for exactly equal times. The bosonic-enhancement of fluctuations is up to 5 atoms. The bosonic correlations decay on the length scale of the measurement probe beam size, in this case 1.5 µm. The bosonic-correlation length scale intrinsic to the atom cloud is λ T = 0.59 µm here, much smaller than beam size. The shot-noise fluctuation at the origin is 24 atoms and therefore less than would be expected from Poissonian statistics, where variance equals the mean. The total mean atom number in the probe region is made of many atoms in the wings of the Gaussian probe beams, which are only partially detected, and hence fluctuations are reduced. Un-normalised fluctuationsC, for a cloud of 10 5 atoms at 100 nK in a pancake-shape trap of 413×20×9 Hz as in Ref. [5], for two measurements at nearly equal times (just 10 µs apart). The phase-space density is 0.82. The scale over which correlations decay is the size of the probe beam used for detection, which is 1.5 µm in this case. For comparison, the de Broglie wavelength is 0.59 µm. The maximum shot noise fluctuation calculated with Eqn. (26) is 24 atoms. The mean atom number in the probe beam is 37, indicating that the fluctuations are effectively suppressed by the presence of atoms in the Gaussian wings of the detection volume.

B. Correlations at two times
We now turn to normalised correlation estimators at two times as defined in Eqn. (19). Fig. 7 shows how the one-particle autocorrelation [first part of Eqn. (25)] decays faster than the bosonic correlations. Furthermore, the shot noise demonstrates anticorrelations as particles which start in the middle move away. When they come back to the centre, a revival can be seen in the correlation function, roughly one full oscillation period along the tightest trap axis. We also note that fluctuations become weaker as the temperature becomes much higher than the critical temperature for Bose-Einstein condensation. Furthermore, the time scale on which the correlations decrease is dictated by the time it takes atoms with the mean thermal speed to cross the detection region (again 1.5 µm here). Since this speed increases with temperature, so the correlations die away faster at higher temperatures, with a typical timescale proportional to 1/ √ T and not the thermal coherence time scale, ∝ 1/T . When the two measurements at different times are also done at different positions, we see again that correlations take a finite time to move from the first point of measurement to the second. As shown in Fig. 8, when the probe beams are well separated (Z ′ > w 0 ), the maximum of correlation is found at finite (not zero) time. This is an observable manifestation of the propagation of correlations at finite speed. For well-separated beams, the revival in correlations is not detectable.

V. REALISTIC MEASUREMENT PROTOCOLS AND UNCERTAINTIES
The determination of atom number will come with some unavoidable, fundamental uncertainties. We will now estimate how many measurements it will take to distinguish the bosonic part of the correlations from the shot noise.
The uncertainty in the correlation measure, Var C , assuming approximately Gaussian statistics, is [27]: where Var (N R,τ ) is the uncertainty in the measured number of atoms in the probe beam at position R and time τ , and N runs is the number of experimental runs. Optical detection of atoms always comes with incoherent photon scattering [28], no matter what property of the light beam is measured. The more photons spontaneously emitted, the more precise the measurement can be. However, spontaneous emission destroys the secondorder coherence that we want to measure. We therefore have to limit the precision of each measurement, in favour of more, independent, experimental runs to reduce the statistical uncertainties. The measurement uncertainty per spontaneously emitted photon can be decreased by using small probe beams, so it pays to use high numerical aperture optics. In a shot-noise limited measurement with a Gaussian beam of size w 0 the uncertainty in the inferred atom number is [26]: where n sc is the number of spontaneous emission events per atom and σ 0 is the resonant atom-light scattering cross section, which can be up to 3λ 2 /2π for λ the resonant optical wavelength. Note that this formula is only valid if the transverse size of the probe beam is significantly smaller than the atom cloud.
To detect the bosonic enhancement of the correlations, the variance on the correlations measure has to be less than the bosonic part of the correlation measure, the second term in Eqn. (25): Using Eqn. (27) and Eqn. (28) we get an expression for the number of measurements required to measure the bosonic part of the correlation measure with a signalto-noise ratio of one: Taking the example of Ref. [5] for experimental parameters, using a 1.5 µm Gaussian-shape probe laser beam and limiting the number of photons scattered to 0.2 per atom, we find that a signal-to-noise ratio of 1 is achieved for about 750 experimental runs, for measurements at nearly equal positions and times. The cloud extends over a few times the size of the detection region, so it is possible to measure correlations for several pairs of probe beams in parallel, pairs being separated by much more than intra-pair spacing. Several measurements can be made on one cloud, although atom number will reduce and temperature will rise.
Only certain kinds of optical measurements are suitable. The optical densities (exponents in the Beer-Lambert law) of the Bose gas clouds discussed are typically larger than 1, so the detector response is not linear in atom number and Eqn. (16) is not valid. Phase shift measurements (far off resonance) have a linear response, but will induce a net, localised, mean phase shift on the atoms due to the AC Stark effect. Such a phase shift will lead to mechanical effects on the atom cloud, such as heating or sound waves, which might disguise the atom-atom correlation signals. We have already developed a suitable two-frequency measurement technique which cancels this mean phase shift [23].

VI. CONCLUSIONS
We have shown that the results of temporal and spatiotemporal correlation measurements in non-interacting thermal Bose gases reflect the underlying atom-atom correlation functions. We have explained how to measure these quantities optically, for Bose gases in harmonic traps in thermal equilibrium. Correlations unsurprisingly take time to travel finite distances. The revivals in the correlation functions are in practice probably too weak to be visible in finite-volume correlation measurements, especially in the presence of non-zero atom-atom scattering [29].
The two-time correlation functions can, in principle, be calculated in other quantum gas systems such as Bose-Einstein condensates, interacting Bose gases, Fermi gases, and for other kinds of potentials, e.g. box potentials or optical lattices. Correlation measures for these systems would provide unambiguous evidence of exotic phase transitions.
One can expect that strong, projective measurements will take the ensemble out of thermal equilibrium. Extending the theory beyond thermal equilibrium is possible using, for example, stochastic time-evolution techniques [30]. We would be very interested to read about the results of any such calculations. the propagator, Eqn. (20). This shows that the commutator is non-zero even for unequal positions. A physical picture, for bound states, is that an initial deltafunction wavepacket spreads out over time. Starting with two such wavepackets, with different positions at different times, there will be some overlap between them due to this spreading. Hence, their field operators do not commute, except for special cases like revivals in a harmonic oscillator. In those cases, the commutator tends to a Dirac delta function function of position in the limit where either the time arguments are equal or separated a multiple of the characteristic oscillation period.
In general, a measurement done at time t will modify the state of the system and therefore change D(r, r ′ , t − t ′ ). For the purpose of this article, however, we made the strong assumption that the measurement is only weakly perturbative and can therefore be neglected.
We also note that the numerical implementation of Eqn. (25) has two particular difficulties. First, the integration of complex numbers is not well handled by some low-level integration routines. Secondly, and more severely, the propagators tend to Dirac delta functions of position in the limit where the time arguments are equal (or separated by an integer number of harmonic oscillator periods). That poses no problem for the second term due to the thermal, imaginary-time component, but means that numerical integrals of the first part can be non-convergent, or even worse, biased. In Fig. 8, we have not presented numerical results close to the divergences of the commutator for that reason. We note that evaluation of Fig. 8 took about 80 hours of processor time.