Extracting density-density correlations from in situ images of atomic quantum gases

We present a complete recipe to extract the density-density correlations and the static structure factor of a two-dimensional (2D) atomic quantum gas from in situ imaging. Using images of non-interacting thermal gases, we characterize and remove the systematic contributions of imaging aberrations to the measured density-density correlations of atomic samples. We determine the static structure factor and report results on weakly interacting 2D Bose gases, as well as strongly interacting gases in a 2D optical lattice. In the strongly interacting regime, we observe a strong suppression of the static structure factor at long wavelengths.


Introduction
Fluctuations and correlations result from the transient dynamics of a many-body system deviating away from its equilibrium state. Generally, fluctuations are stronger at higher temperatures and when the system is more susceptible to the external forces (as governed by the fluctuation-dissipation theorem, see [1,2]). Local fluctuations and their correlations can thus be a powerful tool to probe thermodynamics, and to identify phase transition of a many-body system due to the sudden change of the susceptibility to the thermodynamic forces.
Measurement of fluctuations and correlations on degenerate atomic gases can reveal much information about their quantum nature [3]. Experiment examples include the quantum statistics of the atoms [4,5,6,7,8], pairing correlations [9] and quantum phases in reduced dimensions [10,11]. In these experiments, images of the sample are taken after the time-of-flight expansion in free space, from which the momentum-space correlations are extracted.
In situ imaging provides a new and powerful tool to examine the density fluctuations in real space [12,13,14,15,16,17,18,19,20], offering a complimentary description of the quantum state. This new tool has been used to resolve spatially separated thermodynamic phases in inhomogeneous samples. From in situ measurements, both Mott insulator density plateaus and a reduction of local density fluctuations were observed [13,17,18,21]. Furthermore, a universal scaling behavior was observed in the density fluctuations of 2D Bose gases [22].
Precise measurements of spatial correlations, however, present significant technical challenges. In in situ imaging, one typically divides the density images into small unit cells or pixels and then evaluates the statistical correlation of the signals in the cells. If both the dimension of the cell and the imaging resolution are much smaller than the correlation length of the sample, the interpretation of the result is straightforward. In practice, because the correlation length of quantum gases is typically on the order of 1 µm, comparable to the optical wavelength that limits the image resolution, interpreting experimental data is often more difficult. Finite image resolution, due to either diffraction, aberrations or both, contributes to systematic errors and uncertainties in the fluctuation and correlation measurements.
In this paper, we present a general method to determine density-density correlations and static structure factors of quantum gases by carefully investigating and removing systematics due to imaging imperfections. In Section 2, we review the static structure factor and its relation to the real space density fluctuations. In Section 3, we describe how the density fluctuation power spectrum of a non-interacting thermal gas can be used to calibrate systematics in an imperfect imaging system, and show that the measurement can be explained by aberration theory. In Section 4, we present measurements of density fluctuations in weakly interacting 2D Bose gases and strongly interacting gases in a 2D optical lattice, and extract their static structure factors from the density-density correlations.

The density-density correlation function and the static structure factor
We start by considering a 2D, homogeneous sample at a mean densityn. The densitydensity correlation depends on the separation r 1 − r 2 between two points, and the static correlation function ν(r) is defined as [23] ν(r 1 − r 2 ) =n −1 δn(r 1 )δn(r 2 ) where ... denotes the ensemble average, and δn(r) = n(r) −n is the local density fluctuation around its mean valuen. The Dirac delta function δ(r 1 − r 2 ) represents the autocorrelation of individual atoms, and Ψ † (r 1 )Ψ † (r 2 )Ψ(r 1 )Ψ(r 2 ) = G (2) (r 1 − r 2 ) is the second-order correlation function [24]. When the sample is completely uncorrelated, only atomic shot noise is present and ν(r 1 − r 2 ) = δ(r 1 − r 2 ). At sufficiently high phase space density, when the inter-particle separation becomes comparable to the thermal de Broglie wavelength λ dB or the healing length, density-density correlation becomes nonzero near this characteristic correlation length scale and ν(r) deviates from the simple shot noise behavior. The static structure factor is the Fourier transform of the static correlation function [23,25] S(k) = ν(r)e −ik·r dr, where k is the spatial frequency wave vector. We can rewrite the static structure factor in terms of the density fluctuation in the reciprocal space as [25] where δn(k) = δn(r)e −ik·r dr, and N is the total particle number. Here, δn(−k) = δn * (k) since the density fluctuation δn(r) is real. The static structure factor is therefore equal to the density fluctuation power spectrum, normalized to the total particle number N . A non-correlated gas possesses a structureless, flat spectrum S(k) = 1 while a correlated gas shows a non-trivial S(k) for k near or smaller than inverse of the correlation length ξ −1 . The static structure factor reveals essential information on the collective and the statistical behavior of thermodynamic phases [2,25,26,27]. It has been shown that, through the generalized fluctuation-dissipation theorem [28], the static structure factor of a Bose condensate is directly related to the elementary excitation energy (k) as [25,29] where m is the atomic mass, T is the temperature,h is the Planck constant h divided by 2π, and k B is the Boltzmann constant. See references [25,26,27] on the static structure factor of a general system with complex dynamic density response in the frequency domain. Previous experimental determinations of the static structure factor in the zerotemperature limit, based on two-photon Bragg spectroscopy, have been reported for weakly interacting Bose gases [30,31] and strongly interacting Fermi gases [32]. Here, we show that S(k) at finite temperatures can be directly determined from in situ density fluctuation and correlation measurements.
Experimental determination of S(k) from density fluctuations is complicated by artificial length scales introduced by the measurement, including, for example, finite image resolution and size of the pixels in the charge-coupled device (CCD). Figure 1 shows a comparison between the measurement length scales (the resolution limited spot size and CCD pixel size), the correlation length ξ, and the thermal de Broglie wavelength λ dB . Ideally, a density measurement should count the atom number inside a detection cell (pixel) with sufficiently high image resolution, and the dimension of the cell should be small compared to the atomic correlation length. In our experiment, the image resolution is determined by the imaging beam wavelength λ = 852 nm, the numerical aperture N.A. = 0.28, and the aberrations of the imaging system. The image of a single atom on the CCD chip would form an Airy-disk like pattern with a radius comparable to or larger than λ dB or ξ. The imaging magnification was chosen such that the CCD pixel size √ A = 0.66 µm in the object plane is small compared to the diffraction limited spot radius ∼ 1.8 µm. The atom number N j recorded on the j-th CCD pixel is related to the atom number density drn(r)P(r j − r), assuming the point spread function P(r) is approximately flat over the length scale of a single pixel, where r j is the center position of the j-th pixel in the object plane, n(r) is the atom number density distribution, and the integration runs over the entire x − y coordinate space. The atom number fluctuation measured at the j-th pixel is related to the densitydensity correlation as where δN j = N j − N j is the atom number fluctuation around its mean value N j . In the Fourier space, Eq. (5) can be written as where δn exp (k l ) ≡ j δN j e −ik l ·r j is the discrete Fourier transform of δN j , approximating the continuous Fourier transform. Here, k l = 2π L (l x , l y ), L is the linear size of the image, l x and l y are integer indices in k-space. From Eq. (3) and (7), the power spectrum of the density fluctuation is related to the static structure factor as where the modulation transfer function M(k) = |P(k)| accounts for the imaging system's sensitivity at a given spatial frequency k, and is determined by the point spread function. Also, from Eq. (6), the pixel-wise atom number fluctuation is related to the weighted static structure factor integrated over the k-space, Generalization of the above calculations to arbitrary image resolution and detection cell size is straightforward. In addition to convolving with the point spread function, the measured atom number density must also be convolved with the detection cell geometry. Equation (5) can therefore be written as N j /A = dkn(k)P(k)H(k)e ik·r j , where H(k) = A e ik·r dr/A and the integration goes over the area A of the detection cell. This suggests simply replacing M 2 (k) by M 2 (k)H 2 (k) to generalize Eq. (8) and (9). Finally, in all cases, the discrete Fourier transform defined in Eq. (7) should well approximate the continuous Fourier transform for spatial frequencies smaller than the sampling frequency 1/ √ A. We view the factor M 2 (k)H 2 (k) as the general imaging response function, describing how the imaging apparatus responds to density fluctuations occurring at various spatial frequencies. To extract the static structure factor from in situ density correlation measurements, one therefore needs to characterize the imaging response function at all spatial frequencies to high precision. Since our pixel-size is much smaller than the diffraction and aberration limited spot size, we will from here forward assume 3. Measuring the imaging response function M 2 (k) In this section, we show how to use density fluctuations of thermal atomic gases to determine the imaging response function M 2 (k). Other approaches based on imaging individual atoms can be found, for example, in Refs. [33,34].

Experiment
Measuring density fluctuations in low density thermal gases provides an easy way to precisely determine the imaging response function. An ideal thermal gas at low phasespace density has an almost constant static structure factor up to k = λ −1 dB [24] which, in our case, is larger than the sampling frequency 1/ √ A. Therefore, the density fluctuation power spectrum derived from an ideal thermal gas reveals the square of the modulation transfer function, as indicated by Eq. (8).
We prepare a 2D thermal gas by first loading a three-dimensional 133 Cs Bose-Einstein condensate with 2×10 4 atoms into a 2D pancake-like optical potential with trap frequencies ω z = 2π × 2000 Hz (vertical) and ω r = 2π × 10 Hz (horizontal) [21,22]. We then heat the sample by applying magnetic field pulses near a Feshbach resonance. After sufficient thermalization time, we ramp the magnetic field to 17 G where the scattering length is nearly zero. The resulting thermal gas is non-interacting at a temperature T = 90 nK and its density distribution is then recorded through in situ absorption imaging [22].
We evaluate the density fluctuation power spectrum M 2 exp (k l ) = |δn exp (k l )| 2 using 60 thermal gas images (size: 256 × 256 pixels). Figure 2(a) shows a sample of the noise recorded in the images. Outside the cloud (whose boundary roughly follows the dashed line), the noise is dominated by the optical shot noise, and is therefore uncorrelated and independent of spatial frequency. In the presence of thermal atoms, we observe excess noise due to fluctuations in the thermal atom density. The noise power spectrum is shown as an image in Fig. 2(b), with line-cuts shown in Fig. 2(d). We note that the power spectrum acquires a flat offset extending to the highest spatial frequency, due to the photon shot noise in the imaging beam. Above the offset, the contribution from atomic density fluctuations is non-uniform and has a hard edge corresponding to the finite range of the imaging response function. Close to the edge, ripples in the noise power spectrum appears because of aberrations of the imaging optics, discussed in later paragraphs. Finally, the bright peak at the center corresponds to the large scale density variation due to the finite extent of the trapped atoms, and is masked out in our following analysis.
To fully understand the imaging response function with imaging imperfections, we compare our result with calculations based on Fraunhofer diffraction and aberration theory [35] as described in the following paragraphs.

Point spread function in absorption imaging
We consider a single atom illuminated by an imaging beam, the latter is assumed to be a plane wave with a constant phase across both the object and the image planes. The atom, driven by the imaging electric field, scatters a spherical wave (dark field) which interferes destructively with the incident plane wave [36]. The dark field is clipped by the limiting aperture of the imaging optics and is distorted by the imaging aberrations. An exit pupil function p is used to describe the aberrated dark field at the exit of the imaging optics [35], and its Fourier transform p(k) with k ∝ R describes the dark field distribution on the CCD chip, where R is the position in the image plane. The image of a single atom is then an absorptive feature formed by the interference between the dark field and the incident plane wave in the image plane.
We extend this to absorption imaging of many atoms with density n(r) = i δ(r − r i ), where r i is the location of the i-th atom in the object plane. The total scattered field in the image plane is ‡ ∆E = i p(k−k i ), with each atom contributing a dark field amplitude , and k relates to the position r in the object plane through k = r/ad, where a is the radius of the limiting aperture and d = λ/ (2πN.A.). The dark field ∝ e iδs E 0 is proportional to the incident field E 0 , and carries with a phase shift δ s associated with the laser beam detuning from atomic resonance. For a thin sample illuminated by an incident beam with intensity I 0 , the beam transmission is t 2 = |E 0 + ∆E| 2 /|E 0 | 2 ≈ 1 + 2 [∆E/E 0 ] and the atomic density n exp ∝ − ln(t 2 ) + (1 − t 2 )I 0 /I sat [37,22] leads to n exp ∝ −2(1 Here, [.] refers to the real part and I sat is the saturation intensity for the imaging transition. Comparing the above expression to Eq. (5), we derive the point spread function as P(r) ∝ [e iδs p(k)]| k=r/ad , ‡ We consider the density n of the 2D gas in a range that the photon scattering cross section remains density-independent [38].
in contrast to the form |p(k)| 2 in the case of fluorescence or incoherent imaging.
When the dark field passes through aberrated optics, neither the amplitude nor the phase at the exit pupil is uniform, but is distorted by imperfections of the imaging system. To account for attenuation and phase distortion, We can modify the exit pupil function as p(r p , θ p ) = U (r p /a, θ p )e iΘ(rp/a,θp) , where r p and θ p are polar coordinates on the exit pupil, U (ρ, θ) is the transmittance function, and Θ(ρ, θ) is the wavefront aberration function. We assume U to be azimuthally symmetric and model it as is the Heaviside step function setting a sharp cutoff when r p reaches the radius a of the limiting aperture. The factor e −ρ 2 /τ 2 , with 1/e radius r p = aτ , is used to model the weaker transmittance at large incident angle due to, e.g., finite acceptance angle of optical coatings. For the commercial objective used in the experiment, we need only to include a few terms in the wavefront aberration function where the parameters used to quantify the aberrations are: S 0 for spherical aberration, α for astigmatism (with φ the azimuthal angle of the misaligned optical axis), and β for defocusing due to atoms not in or leaving the focal plane during the imaging. Using the exit pupil function in Eq. (10), we can evaluate the point spread function via P(r) ∝ [e iδs p(k)]| k=r/ad with proper normalization. We can also calculate the modulation transfer function M(k) = |P(k)| . In fact, determination of any one of p(r p ), P(r), or M(k) leads to a complete characterization of the imaging system including its imperfections.

Modeling the imaging response function
We fit the exit pupil function p in the form of Eq. (10), using a discrete Fourier transform, by comparing to the thermal gas noise power spectrum M 2 exp shown in Fig. 2(b). Here, FT (.) denotes the discrete Fourier transform. Figure 2(c) shows the best fit to the measurement, which captures most of the relevant features in the experiment data. Sample line-cuts with uniform angular spacing are shown in Fig. 2(d). This experimental method can in principle be applied to general coherent imaging systems, provided the signal-to-noise ratio of the power spectrum image is sufficiently good to resolve all features contributed by the aberrations. Moreover, one can obtain analytic expressions for the point spread function and the modulation transfer function once the exit pupil function is known (see Appendix A).
Having determined the imaging response function, one can remove systematic contributions from imaging imperfections to the static structure factor as extracted from the power spectrum of atomic density fluctuations, see Eq. (8).

Measuring density-density correlations and static structure factors of interacting 2D Bose gases
We measure the density-density correlations of interacting 2D Bose gases based on the method presented in the previous sections. This study is partially motivated by a finding in our earlier work that the local density fluctuation of a 2D Bose gas is suppressed when it enters the Berezinskii-Kosterlitz-Thouless (BKT) fluctuation and the superfluid regions [22]. We attributed this phenomenon to the emergence of long densitydensity correlation length exceeding the size of the imaging cell and the resolution. This results in a smaller pixel-wise fluctuation δN 2 /A than the simple product of the thermal energy k B T and the compressibility κ, as is expected from the classical fluctuationdissipation theorem (FDT) [1]. Below, we present a careful analysis of the densitydensity correlations of interacting 2D Bose gases and discuss the role of correlations in the FDT. To extract local properties from a trapped sample, we limit our analysis to a small central area of the sample where the density is nearly flat. In addition, the area is chosen to be large enough to offer sufficient resolution in the Fourier space. We choose the patch size to be 32×32 pixels. Figure 3 shows a typical image and the density fluctuations inside the patch.
To ensure that we obtain an accurate static structure factor using the small patch, we perform a measurement on a non-interacting 2D thermal gas at a phase space density nλ 2 dB = 0.5 and compare the measured static structure factor to the theory prediction [24]. We first calculate the imaging response function M 2 (k) for a patch size of 32 × 32 pixels and divide the thermal gas noise power spectrum by M 2 (k). The resulting spectrum should represent the static structure factor of an ideal 2D thermal gas. In Fig. 4, we plot the azimuthally averaged static structure factor with data points uniformly spaced in k, up to the resolution limited spatial frequency k = 2πN.A./λ. The measured static structure factor is flat and agrees with the expected value of S(k) ≈ 1.3 for k < λ −1 dB = 2 µm −1 §. Applying the same analysis to interacting 2D Bose gases, we observe very different strengths and length scales for the density fluctuations. In Fig. 4(a-c), we present the single-shot image noise of samples prepared under three different conditions: weakly interacting gases at the temperature T = 40 nK (below the BKT critical point), with dimensionless interaction strength † g = 0.05 and 0.26 (phase space density nλ 2 dB = 9 and 7); and a strongly interacting 2D gas at the temperature T = 8 nK, prepared in a 2D optical lattice at a mean site occupancy number of 2.6, and a depth of 7 E R , where E R = h × 1.3 kHz is the recoil energy. Due to the tight confinement, the sample in the optical lattice has a high effective interaction strength † g eff = 1.0 [39]. Details on the preparation of the 2D gases in the bulk and in an optical lattice can be found in Refs. [22] and [39], respectively.
The difference in the density fluctuations shown in Fig. 4(a-c) can be characterized in their static structure factors shown in Fig. 4(d). We observe positive correlations above the shot noise level S(k) = 1 in the two weakly interacting samples. The one at g = 0.05 shows stronger density correlations at small k than does the sample at g = 0.26. The enhanced density correlations S(k) > 1 at low momenta are expected since the thermally induced phonon excitations can populate states with length scale 1/k longer than the healing length ξ = 1/ √ ng. For gases with stronger interactions, excitations cost more energy and the excited states are less populated. At smaller g, the correlation length is longer and, therefore, the static structure factor decays at a smaller k.
The most intriguing observation is the negative correlations S(k) < 1 in the strongly interacting gas with g eff = 1.0. We observe a below-shot-noise spectrum at low momentum k, showing that long wave-length excitations are strongly suppressed due to a stronger interaction energy nh 2 g eff /m = k B × 34 nK compared to the thermal energy k B ×8 nK. As the momentum k increases, the excitation populations gradually return to the shot noise level. Our observation is consistent with the prediction in Ref. [29] that when the thermal energy drops below the interaction energy, global density fluctuations in a superfluid are suppressed.
Finally, we discuss the contribution of finite density-density correlations in the FDT. Including correlations, we can write the FDT as k B T κ(r) = δn(r)δn(r ) dr = n(r)S(0) [29]. We compare the measured static structure factor, extrapolated to zero-k, to the value of k B T κ/n, where κ = ∂n/∂µ is the experimentally determined compressibility [22], and indeed find that nS(0) equals to k B T κ to within our § Following the calculation in Ref. [24], we find the static correlation function of an ideal 2D thermal gas is ν(r) = δ(r) + |g 1 (z, e −πr 2 /λ 2 dB )| 2 /g 1 (z, 1)λ 2 dB , where z = e µ/k B T is the local fugacity and g γ (x, y) = ∞ k=1 x k y 1/k /k γ is the generalized Bose function. Fourier transforming ν(r) to obtain the static structure factor S(k), we find S(k) ≈ 1.3 remains flat for kλ dB < 1. † The dimensionless interaction strength of a weakly interacting 2D Bose gas is g = √ 8πa s /l z , where a s is the atomic scattering length and l z = h/mω z is the vertical harmonic oscillator length. For a 2D gas in a 2D optical lattice, the effective interaction strength is g eff = mU l 2 /h 2 [39], where U is the on-site interaction and l is the lattice constant. experimental uncertainties of 10 ∼ 20% for all three interacting samples. This agreement shows that the measured correlations and thus the static structure factor can be linked to the thermodynamic quantities via the FDT. Our ability to determine S(0) and κ from in situ images also suggests a new scheme to determine temperature of the sample from local observables as k B T = nS(0)/κ.

Conclusion
We demonstrated the extraction of density-density correlations and static structure factors from in situ images of 2D Bose gases. Careful analysis and modeling of the imaging response function allow us to fully eliminate the systematic effect of imaging imperfections on our measurements of density-density correlations. For thermal gases, our measurement of the static structure factor agrees well with theory. For interacting 2D gases below the BKT critical temperature, intriguingly, we observe positive densitydensity correlations in weakly interacting samples (g 1) and negative correlations in the strong interaction regime (g eff = 1.0). For all interacting gases, our static structure factor measurements agree with the prediction of the FDT as S(0) = k B T κ/n. Extension of our 2D measurement can further test the prediction of anomalous density fluctuations in a condensate [23,40,41,42] and strong correlations in the quantum critical region [43,44]. Finally, our analysis can be applied to perform precise local thermometry [45] and can potentially be used to extract the local excitation energy spectrum through the application of the generalized fluctuation-dissipation theorem [25,28]. Here, we describe our approach to characterize imaging imperfections using extended Nijboer-Zernike diffraction theory [46]. To obtain the point spread function from the exit pupil function p(r p , θ p ), it is convenient to first decompose the exit pupil function using a complete set of orthogonal functions on the unit disk in the polar coordinates p(r p , θ p ) = where J n (z) is the n-th order Bessel function of the first kind. The point spread function P(r, θ) is the real part of e iδs p(k, θ)| k=r/ad with proper normalization P(r, θ) = 1 N where N = 2πd 2 [e iδs p(r p )]| rp=0 = 2πd 2 cos δ s is the normalizing factor such that P(r)d 2 r = 1. For a unaberrated system, only β 0 0 is non-zero and the above equation reduces to the form J 1 (z)/z, as expected.

Acknowledgement
Using the above equations, we can derive the point spread function from the fitted exit pupil function shown in Fig. A1(a). We calculate the expansion coefficients β m n and evaluate the corresponding point spread function, see Fig. A1(b-c).

Appendix A.2. Modulation transfer function
It is straightforward to evaluate the modulation transfer function M(k) = |P(k)| and the imaging response function M 2 (k) = |P(k)| 2 directly from the exit pupil function p(r p , θ p ). Since the point spread function can be written as P(r) = [e iδs p(k) + e −iδs p * (k)]/4πa 2 N | k=r/ad , its Fourier transform is P(k) = πd 2 N [e iδs p(r p , θ + π) + e −iδs p * (r p , θ)]| rp=kad , (A.5) where k = |k| is the spatial frequency and θ is the polar angle of k in the image plane. From Eq. (A.5), the imaging response function is M 2 (k) ∝ |p(kad, θ + π) + e −2iδs p * (kad, θ)| 2 . This result shows that the phase shift δ s is important since M 2 (k) depends on the interference between p(kad, θ + π) and p * (kad, θ). The transmittance U , defined in the exit pupil function Eq. (10), accounts for the radial envelope in M 2 (k), leading to the sharp edge at k = d −1 = 2πN.A./λ. Either the continuous function