Solving the quantum many-body problem via correlations measured with a momentum microscope

TL;DR

This paper demonstrates measuring multi-particle momentum correlations in ultracold helium atoms using a quantum momentum microscope, providing operational insights into the quantum many-body problem and revealing classical inequality violations.

Contribution

It introduces a method to measure higher-order momentum correlations in a quantum many-body system, advancing understanding of complex quantum phenomena.

Findings

Measured multi-particle momentum correlations up to third order.

Extracted the pairing field amplitude as a key correlation component.

Observed a record violation of the classical Cauchy-Schwarz inequality.

Abstract

In quantum many-body theory, all physical observables are described in terms of correlation functions between particle creation/annihilation operators. Measurement of such correlation functions can therefore be regarded as an operational solution to the quantum many-body problem. Here we demonstrate this paradigm by measuring multi-particle momentum correlations up to third order between ultracold helium atoms in an s-wave scattering halo of colliding Bose-Einstein condensates, using a quantum many-body momentum microscope. Our measurements allow us to extract a key building block of all higher-order correlations in this system|the pairing field amplitude. In addition, we demonstrate a record violation of the classical Cauchy-Schwarz inequality for correlated atom pairs and triples. Measuring multi-particle momentum correlations could provide new insights into effects such as unconventional superconductivity and many-body localisation.

Figures

Tables

Table 1. Table 1: Complete set of 4th order correlation functions in terms of normal and anomalous lattice mode occupancies, n 𝐤 i subscript 𝑛 subscript 𝐤 𝑖 n_{\mathbf{k}_{i}} and m 𝐤 i subscript 𝑚 subscript 𝐤 𝑖 m_{\mathbf{k}_{i}} . All momenta 𝐤 i subscript 𝐤 𝑖 \mathbf{k}_{i} ( i = 1 , 2 , 3 , 4 𝑖 1 2 3 4 i=1,2,3,4 ) are drawn from the halo peak centered radially at | 𝐤 i | = k r subscript 𝐤 𝑖 subscript 𝑘 𝑟 |\mathbf{k}_{i}|=k_{r} , and for notational simplicity we have omitted the momentum index (i.e. n 𝐤 i ≡ n subscript 𝑛 subscript 𝐤 𝑖 𝑛 n_{\mathbf{k}_{i}}\equiv n and m 𝐤 i ≡ m subscript 𝑚 subscript 𝐤 𝑖 𝑚 m_{\mathbf{k}_{i}}\equiv m ). The unnormalised correlation function is denoted via G ( 4 ) ​ ( 𝐤 1 , 𝐤 2 , 𝐤 3 , 𝐤 4 ) superscript 𝐺 4 subscript 𝐤 1 subscript 𝐤 2 subscript 𝐤 3 subscript 𝐤 4 G^{(4)}(\mathbf{k}_{1},\mathbf{k}_{2},\mathbf{k}_{3},\mathbf{k}_{4}) and refers to the numerator in Eq. ( S2 ).

Case	$G^{\left(4\right)}({𝐤}_1,{𝐤}_2,{𝐤}_3,{𝐤}_4)$	$g^{\left(4\right)}({𝐤}_1,{𝐤}_2,{𝐤}_3,{𝐤}_4)$	Decays to case:
(1)	$n^4$	$1$	N/A
(2)	$4{\left|m\right|}^4+16n^2{\left|m\right|}^2+4n^4$	$24+24/n+4/n^2$	(3)
(3)	$2n^2\left(2{\left|m\right|}^2+n^2\right)$	$6+4/n$	(6)
(4)	$2n^2\left(2{\left|m\right|}^2+n^2\right)$	$6+4/n$	(5)
(5)	$2n^4$	$2$	(1)
(6)	$n^2\left(n^2+{\left|m\right|}^2\right)$	$2+1/n$	(1)
(7)	${\left(n^2+{\left|m\right|}^2\right)}^2$	$4+4/n+1/n^2$	(6)
(8)	$6n^2\left(n^2+3{\left|m\right|}^2\right)$	$24+18/n$	(3)
(9)	$6n^2\left(n^2+3{\left|m\right|}^2\right)$	$24+18/n$	(10)
(10)	$6n^4$	$6$	(5)
(11)	$4n^4$	$4$	(5)
(12)	$2n^2\left(n^2+{\left|m\right|}^2\right)$	$6+4/n$	(6)
(13)	$2n^2\left(n^2+{\left|m\right|}^2\right)$	$6+4/n$	(5)
(14)	$24n^4$	$24$	(10)

Equations

$${\bar{g}}_{BB}^{(2)}(0)=(n_{{\mathbf{k}}_0}^2+|m_{{\mathbf{k}}_0}{|}^2)/n_{{\mathbf{k}}_0}^2\simeq 2+1/n.$$

$${\bar{g}}_{BB}^{(3)}(0,0)=(2n_{{\mathbf{k}}_0}^3+4n_{{\mathbf{k}}_0}|m_{{\mathbf{k}}_0}{|}^2)/n_{{\mathbf{k}}_0}^3\simeq 6+4/n.$$

$$\begin{gathered} g^{(N)}({\mathbf{x}}_1,{\mathbf{x}}_2,...,{\mathbf{x}}_N)=\frac{⟨:\hat{n}({\mathbf{x}}_1)\hat{n}({\mathbf{x}}_2)...\hat{n}({\mathbf{x}}_N):⟩}{⟨\hat{n}({\mathbf{x}}_1)⟩⟨\hat{n}({\mathbf{x}}_2)⟩...⟨\hat{n}({\mathbf{x}}_N)⟩} \\ =\frac{⟨{\hat{\Psi }}^{†}({\mathbf{x}}_1){\hat{\Psi }}^{†}({\mathbf{x}}_2)...{\hat{\Psi }}^{†}({\mathbf{x}}_N)\hat{\Psi }({\mathbf{x}}_N)...\hat{\Psi }({\mathbf{x}}_2)\hat{\Psi }({\mathbf{x}}_1)⟩}{⟨\hat{n}({\mathbf{x}}_1)⟩⟨\hat{n}({\mathbf{x}}_2)⟩...⟨\hat{n}({\mathbf{x}}_N)⟩}. \end{gathered}$$

$$g^{(N)}({\mathbf{k}}_1,{\mathbf{k}}_2,...,{\mathbf{k}}_N)$$

$$\frac{d{\hat{a}}_{\mathbf{k}}(t)}{dt}$$

$$\frac{d{\hat{a}}_{-\mathbf{k}}^{†}(t)}{dt}$$

$${\hat{a}}_{\mathbf{k}}(t)$$

$${\hat{a}}_{-\mathbf{k}}^{†}(t)$$

$${\alpha }_{\mathbf{k}}(t)$$

$${\beta }_{\mathbf{k}}(t)$$

$$n_{\mathbf{k}}(t)$$

$$m_{\mathbf{k}}(t)$$

$$|m_{\mathbf{k}}{|}^2=n_{\mathbf{k}}\left(n_{\mathbf{k}}+1\right),$$

$$g_{CL}^{(2)}(\mathbf{k})$$

$$g_{BB}^{(2)}(\mathbf{k})$$

$$n=\frac{1}{M}\sum _{\mathbf{k}\in V_S}⟨{\hat{n}}_{\mathbf{k}}⟩=\frac{1}{M}\sum _{\mathbf{k}\in V_S}{n}_{\mathbf{k}}\simeq n_{{\mathbf{k}}_0},$$

$$|m|=\frac{1}{M}\sum _{\mathbf{k}\in V_S}|m_{\mathbf{k}}|\simeq |m_{{\mathbf{k}}_0}|.$$

$$V_S\simeq 4\pi \sqrt{2\pi }{k}_r^{2}w$$

$$V_M\simeq (2\pi {)}^{3/2}({\sigma }_k{)}^3,$$

$$g_{BB}^{(2)}(\Delta \mathbf{k})=\frac{{\sum }_{\mathbf{k}\in V}⟨:{\hat{n}}_{\mathbf{k}}{\hat{n}}_{-\mathbf{k}+\Delta \mathbf{k}}:⟩}{{\sum }_{\mathbf{k}\in V}⟨{\hat{n}}_{\mathbf{k}}⟩⟨{\hat{n}}_{-\mathbf{k}+\Delta \mathbf{k}}⟩},$$

$$g_{BB}^{(3)}(\Delta {\mathbf{k}}_1,\Delta {\mathbf{k}}_2)=\frac{{\sum }_{{\mathbf{k}}_3\in V}⟨:{\hat{n}}_{{\mathbf{k}}_3}{\hat{n}}_{{-\mathbf{k}}_3+\Delta {\mathbf{k}}_1}{\hat{n}}_{{-\mathbf{k}}_3+\Delta {\mathbf{k}}_2}:⟩}{{\sum }_{{\mathbf{k}}_3\in V}⟨{\hat{n}}_{{\mathbf{k}}_3}⟩⟨{\hat{n}}_{-{\mathbf{k}}_3+\Delta {\mathbf{k}}_1}⟩⟨{\hat{n}}_{-{\mathbf{k}}_3+\Delta {\mathbf{k}}_2}⟩}.$$

$$g_{CL}^{(2)}(\Delta \mathbf{k})=\frac{{\sum }_{\mathbf{k}\in V}⟨:{\hat{n}}_{\mathbf{k}}{\hat{n}}_{\mathbf{k}+\Delta \mathbf{k}}:⟩}{{\sum }_{\mathbf{k}\in V}⟨{\hat{n}}_{\mathbf{k}}⟩⟨{\hat{n}}_{\mathbf{k}+\Delta \mathbf{k}}⟩},$$

$$g_{CL}^{(3)}(\Delta {\mathbf{k}}_1,\Delta {\mathbf{k}}_2)=\frac{{\sum }_{{\mathbf{k}}_3\in V}⟨:{\hat{n}}_{{\mathbf{k}}_3}{\hat{n}}_{{\mathbf{k}}_3+\Delta {\mathbf{k}}_1}{\hat{n}}_{{\mathbf{k}}_3+\Delta {\mathbf{k}}_2}:⟩}{{\sum }_{{\mathbf{k}}_3\in V}⟨{\hat{n}}_{{\mathbf{k}}_3}⟩⟨{\hat{n}}_{{\mathbf{k}}_3+\Delta {\mathbf{k}}_1}⟩⟨{\hat{n}}_{{\mathbf{k}}_3+\Delta {\mathbf{k}}_2}⟩},$$

$${\bar{g}}_{BB}^{(2)}(\Delta k)=A{e}^{-\Delta k^2/2{\sigma }_{BB}^2}+O,$$

$$V_{HOM}=1-\frac{1}{1+C_2}.$$

Full text

In quantum physics, fully understanding and characterising complex systems, comprising a large (often macroscopic) number of interacting particles, is an extremely challenging problem. Solutions within the standard framework of (first-quantised) quantum mechanics generally require the knowledge of the full quantum many-body wavefunction. This necessitates an exponentially large amount of information to be encoded and simulated using the many-body Schrödinger equation. In an equivalent (second-quantised) quantum field theory formulation, the fundamental understanding of quantum many-body systems comes through the description of all physical observables via correlation functions between particle creation and annihilation operators. Here, the exponential complexity of the quantum many-body problem is converted into the need to know all possible multi-particle correlation functions, starting from two-, three-, and increasing to arbitrary $N$-particle (or higher-order) correlations.

From an experimental viewpoint, an operational solution to the quantum many-body problem is therefore equivalent to measuring all multi-particle correlations. In certain cases, however, knowing only a specific set of (few-body or lower-order) correlations is sufficient to allow a solution of the many-body problem to be constructed. This was recently shown for phase correlations between two coupled one-dimensional (1D) Bose gases Schweigler et al. . Apart from facilitating the description of physical observables, characterising multi-particle correlations is important for introducing controlled approximations in many-body physics, such as the virial- and related cluster-expansion approaches that rely on truncation of the Bogolyubov-Born-Green-Kirkwood-Yvon hierarchy  Shavitt and Bartlett (2009); Kira and Koch (2011). Momentum correlations up to 6th order R. G. Dall et al. (2013) and phase correlations up to 8th Schweigler et al. and 10th order Langen et al. (2015) have so far been measured in ultracold atomic gases. More generally, multi-particle correlation functions have been used to experimentally characterise the fundamental properties of various systems, such as thermal Bose and Fermi gases Jeltes et al. (2007), weakly and strongly interacting 1D Bose gases Kinoshita et al. (2005); Armijo et al. (2010); Fang et al. (2016), tunnel-coupled 1D tubes Langen et al. (2015); Schweigler et al. , collision halos Perrin et al. (2007); Kheruntsyan et al. (2012); R. I. Khakimov et al. (2016), and phenomena such as prethermalisation Gring et al. (2012) and transverse condensation RuGway et al. (2013). Correlations between multiple photons are also routinely used in numerous quantum optics experiments including ghost imaging Chan et al. (2009); Chen et al. (2010), defining criteria for non-classicality Vogel (2008); Ding et al. (2015), analysing entangled states generated by parametric down conversion Bussières et al. (2008) and characterising single photon sources U’Ren et al. (2005).

Here, we demonstrate an experimental solution of the many-body problem as outlined above by measuring second- and third-order correlations between momentum-correlated atoms in a collisional halo between two Bose-Einstein condensates (BECs). The halo is generated by spontaneous $s$-wave scattering of two colliding BECs Perrin et al. (2007); R. I. Khakimov et al. (2016); SOM , creating a spherical shell of pair-correlated atoms (see Fig. 1). After a time-of-flight expansion, we detect the positions of individual atoms, which are mapped back to the initial momenta of the atoms directly after the collision bal . This means that we reconstruct momentum correlation functions from the momenta of individual atoms with full 3D resolution. Thus our detector setup can be regarded as a quantum many-body momentum microscope, complementary to the quantum gas in situ microscopes created using optical lattices Bakr et al. (2009); Sherson et al. (2010); Cheneau et al. (2012); Cheuk et al. (2015); Haller et al. (2015); Omran et al. (2015) or arrays of optical tweezers Lester et al. (2015). We characterise and compare all possible back-to-back and collinear atomic correlation functions for two and three atoms, showing the relationship between the different correlation functions and demonstrating a record violation of the classical Cauchy-Schwarz inequality between the peak values of correlation functions.

The experiments start with a BEC of $\sim 10^{6}$ 4He∗ atoms magnetically trapped in the $m_{J}=+1$ sublevel of the long lived metastable (${2}^{3}$S1) state Vassen et al. (2012). The $s$-wave scattering halos of correlated pairs are produced via a two-step process. First, we transfer $\sim$ 95 $\%$ of the BEC atoms to the untrapped $m_{J}=0$ sublevel with a Raman pulse, giving the atoms a downward (i.e. along the $\mathbf{\hat{z}}$ direction) momentum of $\mathbf{K}=-\sqrt{2}k_{0}\mathbf{\hat{z}}$ SOM in wave-number units, where $k_{0}=2\pi/\lambda$ and $\lambda=1083.2$ nm is the wavelength of a diffraction photon. The untrapped BEC is then diffracted using a second pulse into two or more diffraction orders, using either Bragg or Kapitza-Dirac diffraction SOM . Adjacent pairs of diffracted condensates then collide, producing spherical halos of spontaneously scattered atom pairs via $s$-wave collisions Perrin et al. (2007). Each halo has a radius in momentum space of $k_{r}\approx k_{0}/\sqrt{2}$ and a radial Gaussian width of $w\approx 0.03k_{r}$. The average mode occupancy in each halo ranges from $n=0.0017(17)$ to $n=0.44(2)$. As the scattering in our experiment is always in the spontaneous pair-production regime SOM , the scattering halo can be approximated by an overall quantum many-body state that is the product of independent two-mode squeezed vacuum states analogous to those produced by parametric down-conversion in quantum optics.

The expanding halos then fall $\sim\!850$​ mm (time of flight [TOF] $\sim\!416\!$ ms) onto a multi-channel plate and delay-line detector. Due to the $19.8$ eV internal energy of the ${2}^{3}$S1 state, the individual positions of atoms can be reconstructed in 3D, with a spatial resolution of $\sim\!\leavevmode\nobreak\ \!120$ $\mu$m in $x,y$ (momentum resolution $\sim 0.0044k_{r}$), and a temporal resolution along $z$ of $\sim\!2$ ns ($\equiv 8$ nm or $3\times 10^{-7}k_{r}$). As we are interested in correlations between atoms in different momentum modes, we convert position and time to momentum centered on each halo SOM .

From the reconstructed momentum for each atom, we construct various momentum correlation functions from coincidence counts between atoms within each experimental run that are averaged over all runs SOM . First, we look at momentum correlations between three atoms with momenta $\mathbf{k}_{3}$, $\mathbf{k}_{1}=-\mathbf{k}_{3}+\Delta\mathbf{k}_{1}$ and $\mathbf{k}_{2}=-\mathbf{k}_{3}+\Delta\mathbf{k}_{2}$, i.e. the second two atoms on the opposite side of the halo to the first [see Fig. 2(a-c) for illustration]. We define the relevant correlation function as $\bar{g}_{BB}^{(3)}(\Delta k_{1},\Delta k_{2})$ and refer to it as back-to-back (BB), which is averaged with respect to $\mathbf{k}_{3}$ over the halo and spherically integrated with respect to the directions of vectors $\Delta\mathbf{k}_{1}$ and $\Delta\mathbf{k}_{2}$. Thus it is a function of the absolute values $\Delta k_{1}=|\Delta\mathbf{k}_{1}|$ and $\Delta k_{2}=|\Delta\mathbf{k}_{2}|$ SOM .

Fig. 2 (g) shows a typical surface plot of $\bar{g}_{BB}^{(3)}(\Delta k_{1},\Delta k_{2})$ for the $s$-wave halo generated by Kapitza-Dirac orders $l\!=\!(-2,-3)$. This surface plot also contains other many-body correlation functions within it, shown schematically in Figs. 2 (a)-(c). When $\Delta k_{2}\!=\!0$, we can plot $\bar{g}_{BB}^{(3)}(\Delta k_{1},0)$ (red line in Fig. 2 (g)), which will asymptotically approach $\bar{g}_{BB}^{(2)}(0)$—the two-particle correlation function with one atom on each side of the halo—for $\Delta k_{1}\!\gg\!\sigma_{BB}$, where $\sigma_{BB}$ is the two-particle back-to-back correlation length. Taking $\Delta k_{2}\!\gg\!\sigma_{BB}$, we can also plot $\bar{g}_{BB}^{(3)}(\Delta k_{1},\Delta k_{2}\!\gg\!\sigma_{BB})$ (blue line), which is equivalent to $\bar{g}_{BB}^{(2)}(\Delta k_{1})$ and approaches the uncorrelated case of $\bar{g}_{BB}^{(2)}(\Delta k_{1}\!\gg\!\sigma_{BB})=1$ for large values of $\Delta k_{1}$ (see SOM for a full discussion of the relationship between various correlation functions).

We also measure the collinear (CL) three-atom correlation function [shown in Fig. 2 (d-f)], defined analogously as $\bar{g}_{CL}^{(3)}(\Delta k_{1},\Delta k_{2})$, where now $\Delta k_{1}\!=\!|\mathbf{k}_{3}-\mathbf{k}_{1}|$ and $\Delta k_{2}\!=\!|\mathbf{k}_{3}-\mathbf{k}_{2}|$. A surface plot of this function, measured for the Bragg halo with maximum mode occupancy, is shown in Fig. 2 (h). Like $\bar{g}_{BB}^{(3)}(\Delta k_{1},\Delta k_{2})$, this full correlation function also contains other many-body correlations: for example, $\bar{g}_{CL}^{(3)}(\Delta k_{1},0)$ [red line in Fig. 2 (h)] will asymptotically approach $\bar{g}_{CL}^{(2)}(0)$ (the two-atom collinear correlation function), while $\bar{g}_{CL}^{(3)}(\Delta k_{1},\Delta k_{2}\!\gg\!\sigma_{CL})$ (blue line) will yield $\bar{g}_{CL}^{(2)}(\Delta k_{1})$ SOM . Figs. 2 (g) and (h) therefore show a full characterisation of the hierarchy of all three-body and two-body correlation functions present in our system.

Our collisional halo is an example of a quantum many-body system which, in the spontaneous scattering regime, satisfies Wick’s factorisation scheme SOM . This requires knowledge of both the normal and anomalous second-order operator moments in momentum space, $n_{\mathbf{k},\mathbf{k}+\Delta\mathbf{k}}\!=\!\langle\hat{a}^{\dagger}_{\mathbf{k}}\hat{a}_{\mathbf{k}+\Delta\mathbf{k}}\rangle$ and $m_{\mathbf{k},-\mathbf{k}+\Delta\mathbf{k}}\!=\!\langle\hat{a}_{\mathbf{k}}\hat{a}_{\mathbf{-k}+\Delta\mathbf{k}}\rangle$, with $\hat{a}^{\dagger}_{\mathbf{k}}$ and $\hat{a}_{\mathbf{k}}$ being the respective mode creation and annihilation operators, and the diagonal element of $n_{\mathbf{k},\mathbf{k}+\Delta\mathbf{k}}$ giving the average mode occupancy $n_{\mathbf{k}}\equiv n_{\mathbf{k},\mathbf{k}}$. Knowledge of these quantities is sufficient to reconstruct all higher-order correlation functions and thus completely solve the many-body problem for our system. Here, the anomalous occupancy $m_{\mathbf{k}}\equiv m_{\mathbf{k},-\mathbf{k}}$ (related to the anomalous Green’s function in quantum field theory) describes the pairing field amplitude between atoms with equal but opposite momenta, and is similar to the expectation value of the Cooper pair operator in the Bardeen-Cooper-Schrieffer theory of superconductivity, although in our case the pairing is between two identical bosons.

To examine these factorization properties further, we analyse the dependence of peak correlation amplitudes on the peak halo mode occupancy $n_{\mathbf{k}_{0}}$ and compare them with theoretical predictions. The theory relies on the relationship between the peak anomalous occupancy $|m_{\mathbf{k}_{0}}|$ and $n_{\mathbf{k}_{0}}$: $|m_{\mathbf{k}_{0}}|^{2}=n_{\mathbf{k}_{0}}(n_{\mathbf{k}_{0}}+1)$ SOM . In Fig. 3 (a), we plot the measured peak back-to-back correlation amplitude between two-atoms $\bar{g}_{BB}^{(2)}(0)$, for values of average mode occupancy $n\simeq n_{\mathbf{k}_{0}}$ that span more than two orders of magnitude. $\bar{g}_{BB}^{(2)}(0)$ is extracted by fitting $\bar{g}_{BB}^{(2)}(\Delta k)$ with a Gaussian (for details and plots, see SOM ). From analytic theory we expect $\bar{g}_{BB}^{(2)}(0)$ to scale with $n$ as SOM :

[TABLE]

This relation is plotted as the dashed line in Fig. 3 (a), which matches the data well considering that it is a no free parameters fit. For comparison we also plot the peak collinear correlation between two-atoms, $\bar{g}_{CL}^{(2)}(0)$, shown by squares in Fig. 3 (a) and extracted from $\bar{g}_{CL}^{(2)}(\Delta k)$ in the same way as $\bar{g}_{BB}^{(2)}(0)$. We see values of $\bar{g}_{CL}^{(2)}(0)\simeq 1.5$, seemingly independent of the mode occupancy. This trend is expected theoretically, although in the limit of perfect resolution we would expect $\bar{g}_{CL}^{(2)}(0)=2$ (as in the Hanbury Brown–Twiss effect Perrin et al. (2007, 2008)) for all values of $n$.

From the measured $\bar{g}_{BB}^{(2)}(0)$ and $\bar{g}_{CL}^{(2)}(0)$ at each $n$ we are able to extract the key nontrivial component of all higher-order correlations in the scattering halo—the absolute value of the average anomalous occupancy $|m|\simeq|m_{\mathbf{k}_{0}}|$. This is found through the relation $|m_{\mathbf{k}_{0}}|^{2}/n_{\mathbf{k}_{0}}^{2}\!=\!2C_{2}-1$, where $C_{2}\!\equiv\!\bar{g}_{BB}^{(2)}(0)/\bar{g}_{CL}^{(2)}(0)$. Using the data of Fig. 3 (a) to calculate $C_{2}$, we plot $(2C_{2}-1)^{1/2}\simeq|m|/n$ in Fig. 3 (b). A value of $|m|>n$ is necessary for any system to exhibit non-classical (quantum) behaviour, such as two-mode quadrature squeezing, Einstein-Podolsky-Rosen quadrature-entanglement Kheruntsyan et al. (2005), and Bell inequality violation Lewis-Swan and Kheruntsyan (2015). The fact that we measure values of $|m|/n\!>\!1$ for all $n$ (with $|m|/n\gg 1$ for smallest $n$) is a further demonstration of the strong quantum nature of our system. Since all order correlation functions for this system can be expressed as a function of $n$ and $|m|$, measuring these parameters is essentially equivalent to solving the many-body problem for the collisional halo.

Following a similar analysis for the peak three-atom back-to-back correlation amplitude SOM , extracting $\bar{g}_{BB}^{(3)}(0,0)$ from Gaussian fits to $\bar{g}_{BB}^{(3)}(\Delta k,\Delta k)$, we plot these peak values as a function of $n$ in Fig. 3 (c). Theoretically, we expect $\bar{g}_{BB}^{(3)}(0,0)$ to scale with $|m|$ and $n$ as SOM

[TABLE]

This reflects the enhancement in the correlation amplitude due to both the back-to-back and collinear correlations SOM . In Fig. 3 (c) we plot Eq. (2) as a dashed green line, which agrees quite well with the experimental data.

Additionally, we can construct $\bar{g}_{BB}^{(3)}(0,0)$ from our measured values of $C_{2}$, through the relation $\bar{g}_{BB}^{(3)}(0,0)=8C_{2}-2$. We plot these values in Fig. 3 (c), which match the theory well. This is a direct demonstration of how lower-order correlation functions can be used to construct higher-order correlation functions, showing that measuring a finite number of correlation functions can be operationally equivalent to solving the many-body problem.

The low probability associated with four or more atom coincidence events means that we are unable to perform a full, quantitative analysis of the hierarchy of fourth- and higher-order correlation functions. However, we are able to measure the back-to-back correlation function $\bar{g}_{BB}^{(4)}(\Delta k_{1},\Delta k_{2},\Delta k_{3})$ for four atoms, two on each opposite side of the halo SOM , for $n\!=\!0.31(12)$. This yields $\bar{g}_{BB}^{(4)}(0,0,0)\!=\!70(40)$, compared to the theoretically expected value of $\bar{g}_{BB}^{(4)}(0,0,0)\!\simeq\!24+24/n+4/n^{2}\!\simeq\!143$ for this mode occupancy SOM .

An important feature of our BEC collision experiments compared to previous work Perrin et al. (2007); Kheruntsyan et al. (2012) is that we are able to explore a much larger parameter space, including relatively low values of $n$ and small correlation lengths SOM . Because of this, the values of $\bar{g}_{BB}^{(2)}(0)$ that we measure greatly exceed the maximum possible collinear correlation value of $\bar{g}_{CL}^{(2)}(0)\!=\!2$. Thus our results are the first measurements in the regime $\bar{g}_{BB}^{(2)}(0)\!\gg\!\bar{g}_{CL}^{(2)}(0)$. This is a violation of the simplest formulation of the Cauchy-Schwarz inequality Kheruntsyan et al. (2012) for our system, which dictates that classically we would be restricted to $\bar{g}_{BB}^{(2)}(0)\leq\bar{g}_{CL}^{(2)}(0)$. All previous similar measurements with ultracold atoms were limited to peak correlation amplitudes $\bar{g}_{BB}^{(2)}(0)\!\simeq\!\bar{g}_{CL}^{(2)}(0)$ Perrin et al. (2007); Kheruntsyan et al. (2012). This meant that they were only able to show a violation of the Cauchy-Schwarz inequality for volume-integrated atom numbers, rather than bare peak correlations SOM . Therefore our measurement of $\bar{g}_{BB}^{(2)}(0)\!\gg\!\bar{g}_{CL}^{(2)}(0)$, with $C_{2}\!=\!\bar{g}_{BB}^{(2)}(0)/\bar{g}_{CL}^{(2)}(0)\!>\!100$, represents a more straightforward and much stronger violation of the Cauchy-Schwarz inequality (cf. the maximum value of the corresponding correlation coefficient $C_{2}\!\simeq\!1.2$ measured in Ref. Kheruntsyan et al. (2012)). In fact, to the best of our knowledge even for optical experiments the largest value measured is $C_{2}\!\simeq\!58$ Lee et al. (2006); SOM , meaning that our result of $C_{2}\!>\!100$ is a record for any source.

The Cauchy-Schwarz inequality can also be formulated for higher-order correlation functions. For three-atom correlations in our system it states $\bar{g}_{BB}^{(3)}(0,0)\leq\left(\bar{g}_{CL}^{(2)}(0)\right)^{3/2}$ Ding et al. (2015). Again, we violate this inequality for all data in Fig. 3, with a maximum violation of $\simeq\!100$.

To summarise, we have used a quantum many-body momentum microscope to analyse the spontaneous $s$-wave scattering halos of correlated atom pairs with a range of halo mode occupancies $n$ spanning over two orders of magnitude. We measured the third-order correlation functions $\bar{g}^{(3)}_{BB}$ and $\bar{g}^{(3)}_{CL}$ and confirmed the non-trivial many-body nature of the correlations present. Unlike previous similar measurements, we were able to extract the absolute value of the anomalous occupancy $|m|$ as a function of $n$. $|m|$ and $n$ are all that is required for understanding and predicting all higher-order correlation functions in this system, hence solving the quantum many-body problem in this case. We have also demonstrated a high degree of violation of the classical Cauchy-Schwarz inequality for both two and three atom correlations. This is the first measurement for three atoms, while our two atom result beats the only previous experiment with atoms Kheruntsyan et al. (2012) by nearly two orders of magnitude.

This demonstrated ability to measure higher-order quantum correlations in a complex many-body system (an $s$-wave scattering halo) means that a momentum microscope will be a valuable tool for probing other many-body effects in quantum simulators that possess non-trivial correlations (although this may require additional considerations SOM ). Such effects include many-body localisation and glassy dynamics He et al. (2012), unconventional superconductivity Rey et al. (2009), universal three-body recombination and Efimov resonances Kunitski et al. (2015). Other possible applications include the use of such a microscope as a direct dynamical probe of non-equilibrium many-body effects in TOF expansion.

We thank I. Bloch, C. Regal and A.-M. Rey for helpful comments and B. Henson and D. Shin for technical assistance. This work was supported through Australian Research Council (ARC) Discovery Project grants DP120101390, DP140101763 and DP160102337. S.S.H. is supported by ARC Discovery Early Career Research Award DE150100315. A.G.T. is supported by ARC Future Fellowship grant FT100100468.

Supplemental Material

1. Theoretical background and connection to experimental quantities

In many-body theory, multi-particle correlations are characterised by the $N$-th order ($N$-point) correlation function which, in the normalized and normally-ordered form, can be expressed in terms of the creation and annihilation quantum field operators, $\hat{\Psi}^{\dagger}(\mathbf{x})$ and $\hat{\Psi}(\mathbf{x})$, as

[TABLE]

Here, we have considered equal-time density [$\hat{n}(\mathbf{x}_{j})\!\equiv\!\hat{\Psi}^{\dagger}(\mathbf{x}_{j})\hat{\Psi}(\mathbf{x}_{j})$, $j\!=\!1,2,...,N$] correlations in position space for simplicity, but the definition and the context can be extended to, e.g., correlations in the reciprocal momentum space. The physical meaning of the above correlation function is that it represents the joint probability of detecting $N$ particles in their respective positions $\{\mathbf{x}_{1},\mathbf{x}_{2},...,\mathbf{x}_{N}\}$, normalized to the product of single-particle detection probabilities.

Operationally, sufficient knowledge of higher-order correlation functions is equivalent to solving the quantum many-body problem. For example, consider the case of a thermal state, which is an example of a broad class of Gaussian states. Here, all higher-order correlation functions factorize (as per Wick’s theorem) into a sum of terms involving only products of normal densities $\langle\hat{\Psi}^{\dagger}(\mathbf{x}_{i})\hat{\Psi}(\mathbf{x}_{j})\rangle$ ($i,j=1,2,...,N$). Hence their knowledge alone is sufficient to predict all these higher-order correlations. For the same-point correlation function $g^{(N)}(\mathbf{x},\mathbf{x},...,\mathbf{x})$, this property leads to the simple result that $g^{(N)}(\mathbf{x},\mathbf{x},...,\mathbf{x})=N!$, which for $N=2$ implies the Hanbury Brown–Twiss bunching effect Jeltes et al. (2007). More generally, Wick’s factorisation for Gaussian states requires inclusion of terms also containing products of anomalous densities, $\langle\hat{\Psi}(\mathbf{x}_{i})\hat{\Psi}(\mathbf{x}_{j})\rangle$, which do not have to be zero in general (unlike for thermal states). The anomalous density (also known as the anomalous Green’s function in quantum field theory) involves a product of two annihilation operators and characterises the amplitude of nonlocal pair correlations present in the system.

Momentum-space correlations studied in this work can be defined analogously to Eq. (S1) as

[TABLE]

with $\hat{a}^{\dagger}_{\mathbf{k}}$ and $\hat{a}_{\mathbf{k}}$ being the momentum mode creation and annihilation operators in plane-wave basis, satisfying equal-time bosonic commutation relation $\left[\hat{a}_{\mathbf{k}},\hat{a}^{\dagger}_{\mathbf{k}^{\prime}}\right]=\delta_{\mathbf{k},\mathbf{k}^{\prime}}$, and $n_{\mathbf{k}}\equiv\langle\hat{n}_{\mathbf{k}}\rangle=\langle\hat{a}^{\dagger}_{\mathbf{k}}\hat{a}_{\mathbf{k}}\rangle$ giving the distribution of halo mode occupancies. We have omitted the explicit time-dependence of the operators for notational simplicity.

The scattering problem under investigation can be described within the undepleted and constant (in time) source BEC approximation, valid when the total number of scattered atoms constitutes only a small fraction ($<5\%$) of the initial number of atoms in the source condensates and when the mean kinetic energy of the colliding atoms is much larger than the mean interaction energy per atom Ziń et al. (2005, 2006); Ögren and Kheruntsyan (2009). Furthermore, given that the initial trapping frequencies in all three dimensions are rather weak in our experiments and that the scattering occurs predominantly from the condensate central (high-density) region where the real-space density is nearly constant, we can approximate the source condensate as a uniform system, in which case the scattered atoms can be described by the following pairwise coupled Heisenberg equations of motion Ögren and Kheruntsyan (2009):

[TABLE]

Here, $\chi=U\bar{\rho}/\hbar$ is the effective coupling, $U=4\pi\hbar^{2}a_{s}/m$ is the $s$-wave interaction strength characterised by the scattering length $a_{s}$, and $\bar{\rho}$ is the density of the source condensate prior to its splitting into two colliding (counterpropagating) halves. In addition, $\Delta_{k}\equiv\hbar(k^{2}-k_{r}^{2})/\left(2m\right)$ is the effective detuning from the scattering resonance, with $k^{2}=|\mathbf{k}|^{2}$ and $k_{r}$ the collision momentum, i.e., the momentum kick imparted onto each colliding condensate in the center-of-mass frame.

The explicit solutions to these equations, which are analogous to the equations of motion for producing a two-mode squeezed vacuum state in quantum optics by means of parametric downconversion Walls and Milburn (2007), are Kheruntsyan and Drummond (2002); Savage et al. (2006); Lewis-Swan and Kheruntsyan (2014):

[TABLE]

where the coefficients $\alpha_{\mathbf{k}}(t)$ and $\beta_{\mathbf{k}}(t)$ are given by

[TABLE]

and satisfy $|\beta_{\mathbf{k}}|^{2}-|\alpha_{\mathbf{k}}|^{2}=1$.

From these solutions, and for vacuum initial conditions for all scattering modes, one finds that the only nonzero second-order moments of creation and annihilation operators in this model are the normal and anomalous mode occupancies given, respectively, by:

[TABLE]

An important relationship that follows immediately from these solutions is that the amplitude of the anomalous occupancy $|m_{\mathbf{k}}|$ for any mode is related to the normal occupancy $n_{\mathbf{k}}$ by

[TABLE]

where $n_{\mathbf{k}}(t)=\frac{\chi^{2}}{\chi^{2}-\Delta^{2}_{k}}\sinh^{2}\left(\sqrt{\chi^{2}-\Delta^{2}_{k}}\,t\right)$. This follows from the bosonic commutation relation and the conservation of particle number difference in equal but opposite momentum modes in the $s$-wave pair-production process. Note that for any resonant (or peak) halo mode $\mathbf{k}_{0}$, for which $\Delta_{k}\!=\!0$ and therefore $|\mathbf{k}_{0}|\!=\!k_{r}$, the mode occupation acquires the well-known form for the two-mode squeezed vacuum state, $n_{\mathbf{k_{0}}}(t)\!=\!\sinh^{2}\left(\chi\,t\right)$ Walls and Milburn (2007); Kheruntsyan and Drummond (2002); Savage et al. (2006).

The linearity of Eqs. (S3)–(S4) ensures that the evolution will result in a Gaussian many-body state and hence Wick’s factorisation scheme for all higher-order correlation functions will apply. This allows one to express any $N$-point correlation as a sum of terms containing all possible distinct products of the nonzero second-order moments. Taking the second-order correlation function, i.e., for $N\!=\!2$ in Eq. (S2), between two collinear (CL) atoms with the same momenta as an example, this means that it can be evaluated as

[TABLE]

Similarly, the back-to-back (BB) second-order correlation function between two atoms with equal but opposite momenta will be given by

[TABLE]

The third- and any higher-order peak correlations and their expressions in terms of the mode occupancy $n_{\mathbf{k}}$ can be evaluated similarly; the complete set of third-order correlations is listed in Fig. 2 of the main text, while we list fourteen nontrivial fourth-order correlations in Table 1.

Quantitatively, the best agreement between the results based on the uniform source-condensate model and those that can be obtained numerically for nonuniform condensates is achieved when the size of the quantisation box in the uniform model is matched with the characteristic size of the source Savage et al. (2006); Ögren and Kheruntsyan (2009, 2010). In this case, the spacing between the plane-wave modes $\Delta k=2\pi/L$ (which we assume is for a cubic box of side $L$, for simplicity) takes the role of the coherence or correlation length beyond which any higher-order correlation function containing pairs of nearly equal or nearly opposite momenta decays (stepwise) from the peak value to the respective uncorrelated value. For nonuniform systems, on the other hand, the physical correlation length may span over several plane-wave modes depending on the quantization volume adopted for the problem. In this case, the results of the uniform model, while being applicable to peak amplitudes of correlation functions for nonuniform systems, cannot address the quantitative details of the decay of these correlations with distance. Such details can be addressed numerically using stochastic approaches in phase-space Chwedeńczuk et al. (2006); Deuar and Drummond (2007); Krachmalnicoff et al. (2010); Deuar et al. (2014); Lewis-Swan and Kheruntsyan (2014); Lewis-Swan (2016) or analytic approaches based on perturbation theory Bach et al. (2002); Ziń et al. (2006); Chwedeńczuk et al. (2008); Ögren and Kheruntsyan (2009); Lewis-Swan (2016), which are, however, beyond the scope of this work.

Thus, in order to use the simple analytic results of this section for comparison with experimentally measured quantities (see below), we will restrict ourselves to quantitative results that refer to the peak amplitudes of correlation functions. Furthermore, we note that in order to improve the counting statistics the experimental data is averaged over a certain spherical shell (excluding the regions occupied by the source BECs), which is justified due to the spherical symmetry of the scattering halo and the fact that the halo mode size in our experiments is comparable with the radial thickness of the shell. This means that the average halo-mode occupancy ($n$), used in the analysis of experimental results, is identified with

[TABLE]

in our analytic theory; similarly, the absolute value of the average anomalous occupancy $|m|$ refers to

[TABLE]

Here, $V_{S}$ is the averaging volume chosen to be a spherical shell centered (radially) at the halo peak $|\mathbf{k}_{0}|=k_{r}$, whereas $M$ is the number of modes contained in $V_{S}$ (see the next section for details).

2. Experimental details

The experiment begins with a He∗ BEC trapped in a bi-planar quadrupole Ioffe configuration magnetic trap Dall and Truscott (2007). The trap has harmonic frequencies of $\{\omega_{x},\omega_{y},\omega_{z}\}/2\pi\!\simeq\!\{15,25,25\}$ Hz and a bias magnetic field $B_{0}\!=\!1.31(1)$ G along the $\hat{\mathbf{x}}$-axis. Note that the near uniformity of the trapping frequencies justifies the angular integration performed to convert $g_{BB}^{(3)}(\Delta\mathbf{k}_{1},\Delta\mathbf{k}_{2})$ to $\bar{g}_{BB}^{(3)}(\Delta k_{1},\Delta k_{2})$ described in the main text. Our detector consists of a pair of $80$ mm diameter multichannel plates and a delay line, located $\sim\!850$ mm below the trap center ($416$ ms fall time), which has a quantum efficiency of $\approx$10$\%$.

2.1 Halo generation

Similarly to our previous work Manning et al. (2015) we employ the same laser beams for both Raman and Bragg pulses, only changing the relative frequency detuning of the waveforms, which is set by the bias $B_{0}$ and geometrical angle between the beams (90∘).

Bragg diffraction

We use scattering off a grating formed by the same two laser beams in either the Raman, Bragg and Kapitza-Dirac regimes to generate the halos, similar to our previous work Manning et al. (2015); R. I. Khakimov et al. (2016). The laser beams are blue detuned by $\approx$ 2 GHz from the $2^{3}$S${}_{1}\rightarrow 2^{3}$P0 transition. For all experiments we start with a Raman pulse to transfer $\approx 95\%$ of the BEC atoms in the $\mathbf{k}=0$ momentum mode of the $m_{J}=+1$ sublevel into the magnetically insensitive $m_{J}=0$ sublevel. These untrapped atoms are then diffracted using a second pulse into a number of diffraction orders, depending on the desired mode occupancy of the final $s$-wave halo as follows. For the halos with the five highest mode occupancies, the second diffraction pulse uses Bragg scattering to transfer between $\sim\!10\%$ and $\sim\!\leavevmode\nobreak\ 50\%$ of the atoms back into the $\mathbf{k}=0$ momentum mode. The two momentum modes then collide, producing a sphere of scattered atom pairs through $s$-wave collisions Perrin et al. (2007). The sphere has a radius in momentum space approximately equal to $k_{r}\simeq k_{0}/\sqrt{2}$ (note that there is a small deviation due to the halo not being quite spherical Krachmalnicoff et al. (2010); Deuar et al. (2014) and the angle between the diffraction beams not being exactly 90*∘*). This allows us to produce halos with average mode occupancies ranging from $n=0.104(8)$ to $n=0.44(2)$. Here $n=N/M$, where $N$ is the average number of atoms per halo and $M$ is the average number of modes for that halo (see below). The halos for this subset of the data are produced using five different population transfer fractions, with $7$,$305$ shots for each transfer fraction.

Kapitza-Dirac diffraction

To generate halos with smaller mode occupancies, we use higher-order Kapitza-Dirac scattering Kapitza and Dirac (1933); Gould et al. (1986); Ovchinnikov et al. (1999) to produce multiple BECs that each act as distinct sources of $s$-wave scattering halos in an individual experimental run R. I. Khakimov et al. (2016). This yields 12 separate diffraction orders $l=-6...+5$, with each adjacent pair then colliding to produce 11 different halos (see Fig. S1). Each halo still has a radius of $k_{r}$ in momentum space, centred in between the diffracted condensates. As the diffraction orders are populated with varying numbers of atoms, according to the Kapitza-Dirac effect, a consequence is that the average mode occupancy in the resulting halos ($n^{l,l+1}$) ranges from $n^{0,+1}=0.077(10)$ down to $n^{+5,+6}=0.0017(17)$. By creating and measuring $11$ different halo mode occupancies in a single experimental run, this increases the data acquisition rate by over an order of magnitude, which is crucial to achieve good signal to noise for halos with such a low number of counts per shot. The Kapitza-Dirac data comes from $54$,$473$ separate individual experimental runs. This means the total number of different individual halos presented in this manuscript is $635$,$728$.

For both methods of generating $s$-wave scattering halos, it is convenient to operate in momentum space co-ordinates with the origin at the centre of each halo as follows. Knowing the time-of-flight to the detector $T_{f}$, we convert measured atom spatial coordinates $\mathbf{r}$ in the detector plane to momenta $\mathbf{k}$ via $\mathbf{r}=T_{f}\mathbf{v}$, using $\mathbf{v}=\hbar\mathbf{k}/m_{He}$, where $m_{He}$ is the mass of a 4He* atom.

Stimulated versus spontaneous scattering

The combination of the two halo generation methods produces a dataset of halos where the mode occupancy spans more than two orders of magnitude. An important parameter that characterises whether the scattering is in the spontaneous or stimulated regime is the ratio of two characteristic timescales $\gamma=t_{col}/t_{int}$, the collision duration $t_{col}$ and the interaction timescale $t_{int}$. The characteristic collision duration itself is determined by the smallest of two timescales, $t_{col}=\min\{t_{exp},t_{sep}\}$: the characteristic timescale for expansion of the colliding BECs $t_{exp}\simeq 1/\omega_{y,z}$ (which in our experiment is determined by the largest trapping frequency) and the characteristic timescale for geometric separation of the colliding condensates $t_{sep}\simeq 2R_{TF}/3v_{r}$, with $R_{TF}$ being the Thomas-Fermi radius of the BEC in the collision direction ($\sim 70$ $\mu$m) and $v_{r}=\hbar k_{r}/m_{He}$ the collision velocity. The numerical factor of $2/3$ in $t_{sep}$ is introduced to account for the fact that, as the colliding condensates separate in space, the density-dependent scattering rate Ziń et al. (2005); Perrin et al. (2008) becomes significantly smaller well before the condensate overlap region vanishes completely. With these definitions and for our experimental parameters, we estimate that $t_{exp}\sim 8$ ms, whereas $t_{sep}\sim 0.7$ ms, and therefore the characteristic collision duration is $t_{col}\sim 0.7$ ms.

The characteristic interaction timescale, on the other hand, is given by $t_{int}\simeq 1/\chi_{0}=\hbar/U\rho_{0}$ and defines the intrinsic pair production rate should the scattering proceed at a constant effective coupling strength $\chi_{0}$, in the simple undepleted and uniform condensate model of Section 1, Eqs. (S5)-(S6). Within this model, where the halo peak mode occupancy grows as $n_{\mathbf{k}}(t)=\sinh^{2}(\chi_{0}t)$, the spontaneous scattering regime corresponds to $\chi_{0}t\ll 1$ where $n_{\mathbf{k}}(t)\simeq(\chi_{0}t)^{2}$ grows nearly quadratically with time, whereas the stimulated regime corresponds to $\chi_{0}t>1$ where the growth of $n_{\mathbf{k}}(t)$ becomes near exponential. The predictions of this model will remain approximately valid as long as the condensate density in the overlap region is not significantly diminished due to their expansion in free space or their geometric spatial separation. Accordingly, if the stimulated regime is never reached due to $t_{int}$ remaining smaller than the actual collision duration $t_{col}$, the pair production process will remain in the spontaneous regime. Estimating $t_{int}$ for our experiment, with $\chi_{0}=U\bar{\rho}/\hbar$, $U=4\pi\hbar^{2}a_{s}/m_{He}$, $a_{s}=7.51$ nm, and taking $\bar{\rho}$ equal to the average BEC density within the typical parabolic Thomas-Fermi profile, we obtain (for our highest density samples with peak density of $\rho(0)\simeq 1.7\times 10^{18}$ m*-3*) $t_{int}\sim 1$ ms. This gives $t_{col}/t_{int}\simeq 0.7$.

Thus, for the halos with maximum mode occupancy, the largest value for our experiment is $\gamma\simeq 0.7$, meaning it should always be in the spontaneous regime of pair production. Increasing $\gamma>1$ would cause the system to enter the stimulated regime, where Bose-enhancement leads to a few modes dominating the scattering and a range of other effects, such as the formation of phase grains, become apparent Deuar et al. (2013). However, as long as the condensates are not significantly depleted, the Wick’s factorisation scheme described in Section 1 still holds and we would expect to observe similar correlation properties.

2.2 Halo mode occupancy

The mode occupancy in the halo is given by $n=N/M$, where $N$ is the total number of atoms in the halo and $M$ is the number of modes in the halo. Following Ref. Perrin et al. (2008), we take $M=V_{S}/V_{M}$, where

[TABLE]

is the momentum-space volume of the scattering shell (assuming a Gaussian profile radially), with $w$ the rms width of the shell and $k_{r}$ the shell radius. The mode volume, $V_{M}$, is given by

[TABLE]

where $\sigma_{k}$ is the rms width of the momentum distribution of the source BEC, assuming that it can be approximated by an isotropic Gaussian distribution. For a simple Thomas-Fermi parabolic density profile of the source BEC (which our BECs can be well approximated by), the momentum distribution is given by a Bessel function Zambelli et al. (2000); Ögren and Kheruntsyan (2009); the bulk of the Bessel function can be fitted by a Gaussian, which is then used to define the rms width $\sigma_{k}$ in Eq. (S17).

Previous analytic and numerical calculations of the second-order back-to-back correlation function have demonstrated the relation $\sigma_{k}\approx\sigma_{BB}/1.1$ Ögren and Kheruntsyan (2009), which motivates us to use the experimentally measured values of $\sigma_{BB}$ (see the discussion in the following section) to extract $\sigma_{k}$. Note that although $\sigma_{k}$ will be slightly different along each dimension corresponding to the magnetic trap axes, the spherically averaged correlation function measured experimentally in this paper yields only a single, spherically averaged correlation length $\sigma_{BB}$. The assumption of an isotropic correlation length, i.e. that it is the same along all directions, is an approximation. However, there is only a small asymmetry in our trap frequencies, which in our previous work R. I. Khakimov et al. (2016) led to a measured variation of $<15\%$ in the correlation lengths between different axes and thus the approximation should be valid.

From the average Gaussian halo thickness (measured across all halos) of $w_{av}\approx 0.031k_{r}$ and the fitted back-to-back correlation length we find $M$ ranging from 2,900(200)–32,000(21,000).

2.3 Correlation functions

Back-to-back correlations

From the aggregate (over many repeated runs of the experiment) detected atoms within each experimental run, we can reconstruct the normalised two-atom back-to-back correlation function,

[TABLE]

between atom pairs with nearly equal but opposite momenta, $\mathbf{k}$ and $-\mathbf{k}+\Delta\mathbf{k}$. Due to the spherically symmetric nature of the scattering halo, the summation over momenta $\mathbf{k}$ lying within a certain spherical shell of volume $V$ corresponds to first averaging the unnormalised correlation function over the bulk of the halo (to improve the overall statistics) and then normalising the result with respect to the uncorrelated counterpart of the same quantity. Thus, the normalisation and averaging in Eq. (S18) ensures that $g_{BB}^{(2)}(\Delta\mathbf{k})=1$ for uncorrelated momenta, while any dependence on $\mathbf{k}$ within $V$ is lost due to the averaging.

The integration volume $V$ in Eq. (S18), as in all correlation functions described in this paper, is performed over a thin shell in momentum space encompassing the halo but excluding the condensates and any background atoms which lie outside the halo. Radially, this volume includes all atoms with $0.55k_{r}\!<\!|\mathbf{k}|\!<\!1.28k_{r}$, which means all atoms in the halo (gaussian width $\approx\!0.03k_{r}$) will be counted. The volume around the condensates is windowed out by taking $-0.44k_{r}\!<k_{z}<\!0.44k_{r}$, which is somewhat conservative but avoids potential saturation problems that we observe on the delay-line detector when the flux becomes too high. Note that although the experimental integration volume $V$ is somewhat thicker radially than $V_{S}$, very few correlated pairs are measured outside $V_{S}$. Therefore the correlation amplitudes are dominated by atoms in $V_{S}$ and thus the use of Eq. (S14) in theoretical predictions will yield a reasonable approximation for comparison with the experimental results.

Furthermore, for consistency across all correlation functions and to allow easy visualisation of higher-order functions, we additionally integrate over the angles of the displacement vector $\Delta\mathbf{k}$ in spherical coordinates, both for the top and bottom lines of Eq. (S18). We denote the resulting correlation function as $\bar{g}_{BB}^{(2)}(\Delta k)$, which we plot (see Fig. S2 for examples) as a function of the scalar distance between atoms $\Delta k=|\Delta\mathbf{k}|$.

Similarly, the three-atom back-to-back correlation function is given by

[TABLE]

Again, to allow clear visualisation of the correlation function, as well as improve the signal to noise, it is more practical to plot the correlations as a function of the scalar distances $\Delta k_{1}=|\Delta\mathbf{k}_{1}|$ and $\Delta k_{2}=|\Delta\mathbf{k}_{2}|$, which involves spherically integrating over all angles for both the top and bottom terms of Eq. (S19). We denote the corresponding correlation function $\bar{g}_{BB}^{(3)}(\Delta k_{1},\Delta k_{2})$, and this function is plotted in Fig. 2 (g) of the main text.

As a consistency check, we can use the inter-dependency of the different correlation functions to reconstruct $\bar{g}_{BB}^{(2)}(\Delta k)$ from the asymptotes of the full three-atom correlation function $\bar{g}_{BB}^{(3)}(\Delta k_{1},\Delta k_{2})$. This is highlighted by the blue lines $\bar{g}_{BB}^{(3)}(\Delta k_{1},\Delta k_{2}\gg\sigma_{BB})$ and $\bar{g}_{BB}^{(3)}(\Delta k_{1}\gg\sigma_{BB},\Delta k_{2})$ in Fig. 2 (g). We average the two lines and fit a Gaussian to extract $\bar{g}_{BB}^{(2)}(0)$. These extracted values are shown in Fig. S3 along with the theory curve Eq. (1). The agreement between experiment and theory is comparable to the values of $\bar{g}_{BB}^{(2)}(0)$ measured directly [Fig. 3 (a)].

Collinear correlations

Similarly to Eq. (S18), the collinear two-particle correlation function is defined as

[TABLE]

with the same spherical integration as described above performed to transform from $g_{CL}^{(2)}(\Delta\mathbf{k})$ to $\bar{g}_{CL}^{(2)}(\Delta k)$, where $\Delta k=|\Delta\mathbf{k}|$. Due to the low number of counts for some values of $n$ the fits have a large uncertainty and in fact for two halos the statistics are so poor that no meaningful fit is possible. However, theoretically there is no expected trend with $n$, as in the limit of small detector resolution and small correlation function bins we expect $\bar{g}_{CL}^{(2)}(0)\simeq 2$ for all $n$ Perrin et al. (2008); Ögren and Kheruntsyan (2009).

The correlation function between three collinear atoms is then analogously defined as

[TABLE]

where now $\mathbf{k}_{1}=\mathbf{k}_{3}+\Delta\mathbf{k}_{1}$ and $\mathbf{k}_{2}=\mathbf{k}_{3}+\Delta\mathbf{k}_{2}$. As with Eq. (S19), the actual correlation function we plot is integrated over all angles and denoted $\bar{g}_{CL}^{(3)}(\Delta k_{1},\Delta k_{2})$, which is shown in Fig. 2  (h) of the main text. This surface plot is equivalent to previous measurements of $\bar{g}^{(3)}(\Delta k_{1},\Delta k_{2})$ for cold thermal Bose gases Hodgman et al. (2011); R. G. Dall et al. (2013), although note that in the present experiment our BEC source would be much colder and without the collision would not display such correlations Hodgman et al. (2011).

Higher-order correlations

The poor statistics associated with three-atom collinear correlations in the regimes investigated in this manuscript mean that not enough values of $n$ yield $\bar{g}_{CL}^{(3)}(\Delta k_{1},\Delta k_{2})$ with sufficient signal to extract the dependence of $\bar{g}_{CL}^{(3)}(0,0)$ on $n$. However, theoretically we expect $\bar{g}_{CL}^{(3)}(0,0)=6$ for all values of $n$. Even for three-atom back-to-back correlations $\bar{g}_{BB}^{(3)}(\Delta k_{1},\Delta k_{2})$, reliable correlation functions can only be extracted for 11 different values of $n$, hence only 11 points are shown in Fig. 3 (c) of the main text.

The lack of statistics in this dataset also prevented us from extracting the complete 4th and higher-order correlation functions. However, by combining the 5 highest count rate data into a single dataset (with $n=0.31(12)$) we are able to reconstruct the maximally correlated case of two atoms on each side of the halo $\bar{g}_{BB}^{(4)}(\Delta k_{1},\Delta k_{2},\Delta k_{3})$. The four dimensional nature of this correlation function and the poor signal to noise make a clear visual representation problematic, so for simplicity we only plot the diagonal case of $\Delta k_{1}=\Delta k_{2}=\Delta k_{3}=\Delta k$, shown in Fig. S4.

Correlation function fits

To extract the maximum correlation amplitudes [$\bar{g}^{(2)}_{BB}(0),\bar{g}^{(3)}_{BB}(0,0)$ etc.], the individual correlation functions are fitted with a 1D Gaussian fit function Hodgman et al. (2011); R. G. Dall et al. (2013). For $\bar{g}^{(2)}_{BB}(\Delta k)$, this is given by

[TABLE]

where the fit parameters $A$, $\sigma_{BB}$ and $O$ correspond to the correlation amplitude, two-atom correlation length and offset respectively. Note that this means that $\bar{g}^{(2)}_{BB}(0)=A+O$. Similar fit functions are defined for the other correlation functions.

For the three-atom correlation functions, we fit the 1D correlation functions to $\Delta k_{1}=\Delta k_{2}$ [see Fig. 2 (g) and (h)], while the $\bar{g}^{(2)}_{BB}(0)$ values shown in Fig. 3 (c) are taken from the average of fits $\bar{g}^{(3)}_{BB}(0,\Delta k_{2}\gg\sigma_{BB})$ and $\bar{g}^{(3)}_{BB}(\Delta k_{1}\gg\sigma_{BB},0)$ [see Fig. 2 (g) and (h) of the main text]. Due to poor statistics, the uncorrelated values $\bar{g}^{(3)}_{CL}(\Delta k_{1}\gg\sigma_{BB},\Delta k_{2}\gg\sigma_{BB})$ $\bar{g}^{(4)}_{BB}(\Delta k_{1}\gg\sigma_{BB},\Delta k_{2}\gg\sigma_{BB},\Delta k_{3}\gg\sigma_{BB})$ are also slightly larger than one.

Due to the angular integration inherent in our definition of the correlation functions (see, e.g., Eq. (S20)), smaller values of $\Delta k$ will have relatively fewer pairs in both the numerator and denominator of the correlation function compared to larger values of $\Delta k$ (the same is also true of $\Delta k_{1}$ and $\Delta k_{2}$). This means the relative noise will be larger for small values of $\Delta k$. To account for this, when fitting correlation functions to fits such as Eq. (S22), we perform three fits: one to all data, one excluding the smallest value of $\Delta k$ and one excluding the two smallest values of $\Delta k$. In the data shown in Fig. S2, the average of each fit parameter is then used, while the spread is taken as the fit uncertainty, which can be quite large for some data.

$g^{(2)}_{BB}(\Delta k)$ correlation length

From the Gaussian fits shown in Fig. S2, we can extract the two-particle back-to-back correlation length $\sigma_{BB}$ from Eq. (S22), which is plotted against average mode occupancy in Fig. S5. The correlation length is seen to increase with higher mode occupancy. Theoretically, we would expect Ögren and Kheruntsyan (2009) the length scale of $\sigma_{BB}$ to be set by the Thomas-Fermi radius $R_{\mathrm{TF}}$ of the source condensate that is being split into two equal halves in the Bragg diffraction regime. Even in the Kapitza-Dirac regime of splitting into multiple pairs of colliding condensates, the characteristic size of any particular pair is set by $R_{\mathrm{TF}}$ of the source condensate prior to any diffraction, and therefore $\sigma_{BB}$ should be approximately constant for halos of different average mode occupancy $n$. On the other hand, $\sigma_{BB}$ and $n$ have both previously been predicted (using numerical simulations Ögren and Kheruntsyan (2009)) to increase with increasing collision times, which is consistent with our observed result of increasing $\sigma_{BB}$ with larger $n$ (Fig. 6). However, due to subtle differences in our experimental regime (e.g., in our setup $n$ is the final mode occupation, which is varied due to the different total number of atoms in a given pair of colliding condensates, rather than by changing the collision duration), a direct comparison is not possible. A complete understanding of the relationship between $\sigma_{BB}$ and $n$ for our experiment would require detailed numerical simulations, which are beyond the scope of this work.

Note that the correlation lengths we measure are smaller than the experimental values for previously published work Perrin et al. (2007, 2008). Small correlation lengths are desirable in this system, as the correlation length determines the halo mode volume, and thus smaller correlation lengths enable lower mode occupancies to be accessed. The only caveat to this is that the correlation lengths should be larger than the detector resolution, as otherwise it will be impossible to differentiate between adjacent modes. In our system we ensured that $\sigma_{BB}$ is always at least 3 times larger than our worst detector resolution.

2.4 General applicability of momentum microscopes to other many-body systems

As mentioned in the main text, additional considerations may be necessary when looking to apply a similar momentum microscope to study other many-body systems. In some strongly interacting systems the TOF expansion may not initially be ballistic, so in these cases alternative experimental steps or further minimal input from theory might be required, such as a treatment within short-time expansion or linear response theory, to relate the experimentally measured quantities to the in-trap momentum or quasi-momentum distributions. An example of the latter case is the expansion of a 1D Lieb-Liniger gas, for which the long-time asymptotic density distribution can be mapped to the in-trap distribution of quasi-momenta (also known as Bethe rapidities) Campbell et al. (2015).

2.5 Cauchy-Schwarz violation

There are a number of subtleties which need to be considered when comparing our Cauchy-Schwarz violation to experiments with photons. Firstly, in a quantum optics experiment the source is usually parametric down conversion in the low-gain regime, which produces a state approximating a twin-photon Fock state: $|\psi\rangle\sim|0,0\rangle+\alpha|1,1\rangle$, with $\alpha$ a constant. Although there are higher occupancy states ($|22\rangle,|33\rangle$…) generated, in practice they occur with such low probability that any direct experimental measurement of $g^{(2)}_{CL}(0)$ would not result in any statistically significant signal. This means that $g^{(2)}_{CL}(0)$, which enters into the definition of the correlation coefficient $C_{2}$ that quantifies the degree of Cauchy-Schwarz violation, can only accurately be inferred, for example from the visibility $V_{HOM}$ of a Hong-Ou-Mandel dip Hong et al. (1987); Lewis-Swan and Kheruntsyan (2014), where

[TABLE]

To the best of our knowledge, the largest $V_{HOM}$ reported is $V_{HOM}\approx 0.98$ ($98\%$ visibility) Lee et al. (2006), which from Eq. (S23) would imply a value of $C_{2}\approx 58$, hence the claim that our Cauchy-Schwarz violation is the best on record. Note that for another experiment to beat our value of $C_{2}>100$ it would need to measure $V_{HOM}>0.99$ ($>99\%$ visibility). The high degree of violation that we observe clearly demonstrates the quantum nature of our halo.

All previous similar collision halo experiments with ultracold atoms were only able to measure peak correlation amplitudes $\bar{g}_{BB}^{(2)}(0)\!\simeq\!\bar{g}_{CL}^{(2)}(0)$ Perrin et al. (2007); Kheruntsyan et al. (2012). Therefore they were only able to demonstrate a violation of the Cauchy-Schwarz inequality using volume-integrated atom numbers, rather than bare peak correlations. This was made possible by the larger correlation volume of $g^{(2)}_{BB}$ in their system, which was due to the fact that the colliding sources were phase-fluctuating quasicondensates, rather than nearly pure condensates as is the case in our experiment.