2207.14070

## Full counting statistics of interacting lattice gases after an expansion: The role of the condensate depletion in the many-body coherence Ga´ etan Herc´ e, 1, ∗ Jan-Philipp Bureik, 1, ∗ Antoine T´ enart, 1 Alain Aspect, 1 Alexandre Dareau, 1 and David Cl´ ement 1 1 Universit´ e Paris-Saclay, Institut d'Optique Graduate School, CNRS, Laboratoire Charles Fabry, 91127, Palaiseau, France (Dated: December 1, 2022) We study the full counting statistics (FCS) of quantum gases in samples of thousands of interacting bosons, detected atom-by-atom after a long free-fall expansion. In this far-field configuration, the FCS reveals the many-body coherence from which we characterize iconic states of interacting lattice bosons by deducing the normalized correlations g ( n ) (0) up to the order n = 6. In Mott insulators, we find a thermal FCS characterized by perfectly-contrasted correlations g ( n ) (0) = n !. In interacting Bose superfluids, we observe small deviations to the Poisson FCS and to the ideal values g ( n ) (0) = 1 expected for a pure condensate. To describe these deviations, we introduce a heuristic model that includes an incoherent contribution attributed to the depletion of the condensate. The predictions of the model agree quantitatively with our measurements over a large range of interaction strengths, that includes the regime where the condensate is strongly depleted by interactions. These results suggest that the condensate component exhibits a full coherence g ( n ) (0) = 1 at any order n up to n = 6 and at arbitrary interaction strengths. The approach demonstrated here is readily extendable to characterize a large variety of interacting quantum states and phase transitions. ## I. INTRODUCTION The dispersion of a physical quantity contains important information, beyond that obtained from its average value. The analysis of quantum and thermal noise is central in various systems, ranging from quantum electronics [1] and quantum optics [2] to quantum gases [3-5]. The ultimate precision on the measurement of noise is given by the full counting statistics (FCS) [6], which is obtained with single-particle-resolved detection methods that provide the number of particles detected in a given time and/or space interval. These methods yield highorder moments of the particle number beyond the variance. Probing high-order moments is a means to study quantum phase transitions [7-9], universality [10, 11], entanglement properties [12] or out-of-equilibrium dynamics [13]. The FCS has indeed successfully characterized various phenomena in mesoscopic conductors [1, 6, 14, 15], Rydberg [16-18] and non-interacting [19, 20] atomic gases. From a Quantum Information perspective, the FCS holds great promises for large ensembles of particles. In contrast to a full-state tomography [21], the FCS is accessible even in large systems as it probes information only about the diagonal part of the n -body density matrices ( i.e. populations). Although it does not contain the total information about the quantum state, the FCS is sufficient to identify many quantum states without resorting to a consuming tomography. A similar idea was introduced by R. Glauber to characterize light sources from photon correlations at any order [22]. For Gaussian states, for which the Wigner function is positive [2], measuring the FCS or the magnitudes of correlation functions ∗ These authors contributed equally to this work. is indeed equivalent. In strongly-correlated quantum states characterized by non-Gaussian Wigner functions, measuring the FCS and many-body correlations is expected to reveal the nontrivial nature of such states [5, 23-25]. Moreover, recent works have shown that applying random unitary transformations before measuring the FCS provides access to non-diagonal correlators [26, 27], further motivating the development of experimental approaches to the FCS in strongly-interacting quantum systems. In this letter, we report the measurement of the full counting statistics in large three-dimensional (3D) ensembles ( ∼ 5 × 10 3 atoms) of interacting lattice bosons after a free-fall expansion (see Fig. 1(a)). This configuration is analogous to the far-field regime of light propagation during which interferences take place, and after which the FCS identifies quantum states through their many-body coherence [28]. In quantum gases, farfield - or momentum - correlations have been measured with single-atom detection in non-interacting and nondegenerate bosonic [29, 30] and fermionic [31, 32] gases and in Bose-Einstein condensates [33]. More recently momentum correlations in interacting lattice bosons [3436] and interacting fermions [37] were studied. However, the FCS has thus far been measured only in 1D noninteracting bosons [30]. Here, we study various regimes of interacting 3D Bose gases across the superfluid-to-Mott phase transition, extending the measurement of manybody correlations to the strongly-interacting regime. This allows us to reveal the role of the depletion of the condensate in the many-body coherence properties of interacting Bose superfluids. We characterize the FCS by measuring the probability distribution of the atom number N Ω falling in a small voxel V Ω (see Fig. 1(a)). As explained below, a crucial asset of our work is the ability to probe many-body coherence into volumes smaller than that occupied by one FIG. 1. (a) Free-fall expansion of interacting quantum gases of metastable Helium-4 atoms from a three-dimensional optical lattice, yielding the 3D positions of individual atoms in the momentum space. The full counting statistics P ( N Ω ) describes the statistics of the atom number N Ω detected in a small voxel of volume V Ω ∼ ( δk ) 3 (red cube). To reveal the many-body coherence properties of the trapped gas of size L , δk is chosen such that δk 2 π/L . (b) Magnitudes g ( n ) ( 0 ) of n -body correlations as a function of the order n , measured in a Mott insulator (black squares) and a superfluid (blue circles). The black solid line is the prediction for thermal states, g ( n ) (0) = n !. (c) g ( n ) ( 0 ) in the superfluid (blue circles) and in a randomized set (orange squares, see main text). A deviation to the prediction for a pure coherent state, g ( n ) (0) = 1, is observed. <details> <summary>Image 1 Details</summary>  ### Visual Description ## Chart/Diagram Type: Multi-Panel Figure ### Overview The image consists of three panels: (a) a schematic diagram of an experimental setup, (b) a log-linear plot showing the correlation function g^(n)(0) versus the order of correlation n for Mott insulator and Superfluid states, and (c) a linear plot showing the correlation function g^(n)(0) versus the order of correlation n for randomized and Superfluid states. ### Components/Axes **Panel (a): Schematic Diagram** * A cloud of atoms is illuminated by a laser (indicated by red sinusoidal waves). * The length of the cloud is labeled as "L" (green). * The direction of gravity is indicated by "g" with an arrow pointing downwards. * The distance between the cloud and the single-atom detector is labeled as "43 cm". * Below the detector, a magnified view shows the momentum space with axes kx, ky, and kz. * A small cube with side length "δk" is shown in the momentum space. * The length 1/L is indicated on the kx axis. **Panel (b): Log-Linear Plot** * Y-axis: g^(n)(0) (logarithmic scale from 10^0 to 10^3) * X-axis: Order of correlation n (linear scale from 1 to 6) * Legend (top-right): * Black line: n! * Gray squares with error bars: Mott insulator * Blue circles with error bars: Superfluid **Panel (c): Linear Plot** * Y-axis: g^(n)(0) (linear scale from 1.0 to 1.1) * X-axis: Order of correlation n (linear scale from 1 to 6) * Legend (top-right): * Yellow squares with error bars: randomized * Blue circles with error bars: Superfluid ### Detailed Analysis or Content Details **Panel (b): Log-Linear Plot** * **n! (Black Line):** The black line represents the factorial of n. It increases exponentially with n. * n=1: g^(n)(0) ≈ 1 * n=2: g^(n)(0) ≈ 2 * n=3: g^(n)(0) ≈ 6 * n=4: g^(n)(0) ≈ 24 * n=5: g^(n)(0) ≈ 120 * n=6: g^(n)(0) ≈ 720 * **Mott insulator (Gray Squares):** The Mott insulator data points increase exponentially with n, closely following the n! curve. * n=1: g^(n)(0) ≈ 1.0 +/- 0.2 * n=2: g^(n)(0) ≈ 2.0 +/- 0.5 * n=3: g^(n)(0) ≈ 7.0 +/- 2.0 * n=4: g^(n)(0) ≈ 25 +/- 5 * n=5: g^(n)(0) ≈ 120 +/- 30 * n=6: g^(n)(0) ≈ 750 +/- 150 * **Superfluid (Blue Circles):** The Superfluid data points remain approximately constant at g^(n)(0) ≈ 1.0 for all values of n. * n=1: g^(n)(0) ≈ 1.0 +/- 0.1 * n=2: g^(n)(0) ≈ 1.0 +/- 0.1 * n=3: g^(n)(0) ≈ 1.0 +/- 0.1 * n=4: g^(n)(0) ≈ 1.0 +/- 0.1 * n=5: g^(n)(0) ≈ 1.0 +/- 0.1 * n=6: g^(n)(0) ≈ 1.0 +/- 0.1 **Panel (c): Linear Plot** * **randomized (Yellow Squares):** The randomized data points remain approximately constant at g^(n)(0) ≈ 1.0 for all values of n. * n=1: g^(n)(0) ≈ 1.0 +/- 0.02 * n=2: g^(n)(0) ≈ 1.0 +/- 0.02 * n=3: g^(n)(0) ≈ 1.0 +/- 0.02 * n=4: g^(n)(0) ≈ 1.0 +/- 0.02 * n=5: g^(n)(0) ≈ 1.0 +/- 0.02 * n=6: g^(n)(0) ≈ 1.1 +/- 0.02 * **Superfluid (Blue Circles):** The Superfluid data points remain approximately constant at g^(n)(0) ≈ 1.0 for all values of n. * n=1: g^(n)(0) ≈ 1.0 +/- 0.02 * n=2: g^(n)(0) ≈ 1.0 +/- 0.02 * n=3: g^(n)(0) ≈ 1.0 +/- 0.02 * n=4: g^(n)(0) ≈ 1.0 +/- 0.02 * n=5: g^(n)(0) ≈ 1.0 +/- 0.02 * n=6: g^(n)(0) ≈ 1.1 +/- 0.02 ### Key Observations * In panel (b), the Mott insulator exhibits factorial scaling of the correlation function with the order of correlation, while the Superfluid remains constant. * In panel (c), both the randomized and Superfluid states exhibit correlation functions close to 1.0. * The error bars are relatively small for the Superfluid and randomized data in panel (c), indicating high precision. * The error bars for the Mott insulator in panel (b) increase with the order of correlation. ### Interpretation The data suggests that the Mott insulator exhibits strong, high-order correlations, as indicated by the factorial scaling of g^(n)(0). This is a characteristic feature of Mott insulators, where particle number fluctuations are suppressed. In contrast, the Superfluid state exhibits weak correlations, with g^(n)(0) remaining close to 1.0, indicating Poissonian statistics. The randomized data in panel (c) also shows weak correlations, similar to the Superfluid. The experimental setup in panel (a) likely involves trapping atoms and measuring their correlations using a single-atom detector. The momentum space representation shows the distribution of atoms in k-space, with δk representing the momentum resolution. The length scale L is related to the size of the atomic cloud. </details> mode in momentum space, i.e. V Ω (2 π/L ) 3 with L the in-trap size of the gas. This possibility is given by the large quantum efficiency of our detector( η = 0 . 53(2)) [36]. Furthermore, we determine the magnitudes g ( n ) (0) of correlation functions (up to n = 6) from the factorial moments of N Ω [38]. As shown in Fig. 1(b), g ( n ) (0) is found to vary by several orders of magnitude between the superfluid and the Mott insulator states. Interestingly, we observe small deviations to the predictions for a pure condensate in Bose superfluids (see Fig. 1(c)), which is in contrast to a previous work [33]. Since ab-initio calculations of many-body correlations with thousands of interacting atoms are beyond state-of-the-art numerics, we interpret our findings from introducing a heuristic model that includes the contribution of the condensate depletion. Its hypotheses rely on an intuitive physical picture in the weakly-interacting regime but, surprisingly, its predictions are in quantitative agreement with our observations even in the strongly-interacting regime. Our observations lead us to conclude that the condensate component exhibits a full coherence, g ( n ) (0) = 1 (at least up to n = 6), at arbitrary strength of interactions. ## II. FULL COUNTING STATISTICS (FCS) OF PURE-STATE BECS AND MOTT INSULATORS Textbook descriptions of Bose-Einstein Condensates (BECs) and Mott insulators are based on pure states. Bose-Einstein condensation is associated with the breaking of phase symmetry [39] whose complex order parameter defines a coherent state describing the BEC. In a grand canonical approach, coherent states have a Poisson counting statistics, P ( N Ω ) = 〈 N Ω 〉 N Ω exp[ -〈 N Ω 〉 ] /N Ω ! where N Ω is the number of detected bosons in the considered volume V Ω , and a full coherence g ( n ) = 1 at any order n of normalized correlations either in position or in momentum space [22]. In our experiment, we determine $$e \text subscript { e- } \quad g ^ { ( n ) } ( 0 ) = g ^ { ( n ) } ( k , k , \dots , k ) = \frac { \langle [ a ^ { \dagger } ( k ) ] ^ { n } [ a ( k ) ] ^ { n } \rangle } { \langle a ^ { \dagger } ( k ) a ( k ) \rangle ^ { n } } , \quad ( 1 )$$ where k is the momentum at which correlations are evaluated, i.e. where the volume V Ω is located. For a pure coherent state we expect g ( n ) (0) = 1 at all orders. In contrast, a 'perfect' Mott insulator - a uniform Mott insulator at zero temperature - is thought of as a Fock state in the (in-trap) position basis. In the momentum basis, which is probed after a long expansion from the trap, it is expected to exhibit thermal statistics [34, 40, 41]. This is because far-field correlations reflect multi-particle interferences from a discrete series of emitters (atoms in the lattice sites) with no coherence between the sites (no tunnelling), a situation analog to that of light emitted by many incoherent sources. Thermal states are characterized by a counting statistics P ( N Ω ) = (1 -q ) q N Ω where q = 〈 N Ω 〉 / 1 + 〈 N Ω 〉 and g ( n ) ( k , k , ...., k ) = n ! [42]. Note that both the Poisson and the thermal FCS are fully determined by a single parameter, the average number of particles 〈 N Ω 〉 , as a result of the Gaussian character of their quantum state [2]. For Gaussian states, a detection efficiency η smaller than one does not affect the measurement of the FCS - nor that of g ( n ) . Whether these many-body coherence properties of pure states describe experiments is not granted. On the one hand, pure states are approximated descriptions of states produced in an experiment because of the coupling to the environment. Determining the level e.g. the order n of correlations -up to which a description in terms of pure states is valid provides a quantitative certification of experimental realizations. In the context of the development of platforms for quantum technologies, such a certification is of interest. On the other hand, the properties of quantum states produced at thermodynamical equilibrium are affected by constraints on macroscopic quantities. In the canonical (or microcanonical) ensemble there are no fluctuations of the total BEC atom number N BEC in a gas of non-interacting bosons at zero temperature [43]: N BEC is fixed to the total atom number, N BEC = N . The statistics of N BEC is therefore not that of a coherent state [44]. To alleviate such global constraints and mimic a grand canonical ensemble, one may probe the FCS in a volume V Ω much smaller than that, V BEC , occupied by the BEC. The volume V BEC -V Ω ∼ V BEC then serves as a reservoir for the sub-volume V Ω where the number of bosons N Ω N may fluctuate [45]. Measuring the FCS in small sub-volumes V Ω is also crucial if correlation functions g ( n ) ( k , k ′ , ...., k ′′ ) are bell shaped with widths l ( n ) c . While this is not an issue for a coherent state, which is fully coherent over the entire volume it occupies, correlation functions of a thermal state must be probed in a volume V Ω much smaller than the coherence volume V ( n ) c = [ l ( n ) c ] 3 , since particles distant by l ( n ) c are essentially uncorrelated. This is well known in the case of the Hanbury-Brown and Twiss effect where the property g (2) ( k , k ′ ) = 2 is expected only if | k -k ′ | is less that the width of the far-field diffraction pattern associated with the source size. In the far-field, the correlation lengths of Mott insulators and of thermal Bose gases are set by the inverse in-trap size, l ( n ) c ∼ 2 π/L [34, 35]. As a result, observing fully-contrasted n -body correlations requires using V Ω (2 π/L ) 3 . This condition also ensures the above-mentioned criterion on probing BECs as the volume occupied by the BEC in the momentum space is set by ∆ k ∼ 1 . 6 /L [46]. With these consider- ations in mind, we compute the magnitudes g ( n ) (0) of correlation functions in a volume V Ω ∼ ( δk ) 3 of the momentum space such that δk × L 2 π (see Fig. 1(a)). As illustrated below, this choice is essential to correctly reveal the full statistical properties. ## III. MEASUREMENT OF THE FCS IN MOTT INSULATORS AND BOSE SUPERFLUIDS Our measurement of the FCS in the far-field exploits the 3D atom-by-atom detection of metastable Helium4 ( 4 He ∗ ) after a long free-fall expansion [47, 48] (see Fig. 1(a)). The detection efficiency is η = 0 . 53(2) per atom, with negligible dark counts. For an individual run, we register the number of atoms in each of the voxels V Ω mentioned above, and we use more than 2000 runs to obtain the probability distribution in each voxel. We apply this approach to probe equilibrium quantum states of 4 He ∗ atoms loaded in the lowest energy-band of a threedimensional (3D) optical lattice [49]. The lattice implements the 3D Bose-Hubbard Hamiltonian whose main parameters are the tunnelling amplitude J and the onsite (repulsive) interaction U . We first investigate the thermal nature of Mott insulators in the momentum space. We realize Mott insulators with N = 6 . 5(6) × 10 3 atoms at U/J = 76, which corresponds to a lattice filling of one atom per site at the trap center [34] and an almost uniform filling of the first Brillouin zone in the momentum space. To compute the counting statistics in this first set of experiments, we divide the first Brillouin zone into cubic voxels V Ω of size δk = 6 × 10 -2 k d and average the probability distributions measured over all these voxels. Here k d = 2 π/d is the momentum associated with the lattice spacing d = 775 nm and the size δk is comparable to that of one mode in momentum space, δk ∼ 2 π/L (see Appendix A). The resulting FCS of the Mott state is shown in Fig. 2(a). It is found to be in excellent agreement with a thermal statistics whose average atom number is that measured in the experiment, 〈 N Ω 〉 = 0 . 46(5). The statistical properties of the Mott state are also revealed through the magnitudes g ( n ) (0) of the normalized correlation functions, as an alternative to the probability distribution P ( N Ω ). We determine g ( n ) (0) from the factorial moments of the detected atom number N Ω in small voxels V Ω with δk 2 π/L (see Appendix B), $$g ^ { ( n ) } ( 0 ) = \frac { \langle N _ { \Omega } ( N _ { \Omega } - 1 ) \dots ( N _ { \Omega } - n + 1 ) \rangle } { \langle N _ { \Omega } \rangle ^ { n } } , \quad ( 2 )$$ transposing a well-known approach in quantum optics [50]. The correlation functions of quantum optics [22] are defined with normal ordering of the destruction operators since a detected photon is destroyed. Even if it is not always emphasized, the correlation functions in quantum optics are indeed calculated from the factorial moments, and it is remarkable that the corresponding quantities in classical statistics are also based on factorial moments FIG. 2. (a) Full Counting statistics P ( N Ω ) to find N Ω atoms in a small volume V Ω of the momentum space when probing a Mott insulator with unity filling (black circles). The predictions for thermal ( resp. Poissonian) statistics is shown as a dashed-dotted black ( resp. dashed blue) line. We use the measured value 〈 N Ω 〉 = 0 . 46(5) for the theoretical predictions (the shaded areas reflect the uncertainty on 〈 N Ω 〉 ). The error bars denote the statistical uncertainty (standard deviation) estimated with the bootstrapping method. (b) Same as (a) in the BEC mode ( k = 0 ) of interacting lattice superfluids (SF) with U/J = 5. The mean atom number is 〈 N Ω 〉 = 5 . 3(2). Error bars are smaller than the dots. <details> <summary>Image 2 Details</summary>  ### Visual Description ## Combined Chart: Atom Number Distribution ### Overview The image presents two plots, (a) and (b), showing the probability distribution P(NΩ) of detected atom numbers NΩ under different conditions. Plot (a) compares a Mott insulator distribution with thermal and Poisson distributions on a semi-log scale. Plot (b) compares a BEC mode distribution with thermal and Poisson distributions on a linear scale. ### Components/Axes **Plot (a):** * **Title:** (a) (top-left) * **X-axis:** Detected atom number NΩ, ranging from 0 to 7. * **Y-axis:** Probability P(NΩ), ranging from 10^-3 to 10^0 (logarithmic scale). * **Legend (top-right):** * Thermal (black dashed-dotted line) * Poisson (blue dashed line) * Mott insulator (gray squares with error bars) **Plot (b):** * **Title:** (b) (top-left) * **X-axis:** Detected atom number NΩ, ranging from 0 to 19. * **Y-axis:** Probability P(NΩ), ranging from 0.00 to 0.15 (linear scale). * **Legend (top-right):** * Thermal (black dashed-dotted line) * Poisson (blue dashed line) * BEC mode (blue circles with error bars) ### Detailed Analysis **Plot (a): Mott Insulator, Thermal, and Poisson Distributions** * **Mott insulator (gray squares):** The data points are approximately: * (0, 0.8) * (1, 0.3) * (2, 0.08) * (3, 0.03) * (4, 0.01) * (5, 0.004) * (6, 0.001) * (7, 0.0003) The Mott insulator distribution shows a decreasing probability with increasing atom number. The gray shaded area around the Mott insulator data points represents the uncertainty. * **Thermal (black dashed-dotted line):** The thermal distribution follows a similar decreasing trend as the Mott insulator, but is slightly higher at larger atom numbers. * **Poisson (blue dashed line):** The Poisson distribution decreases more rapidly than the Mott insulator and thermal distributions. It falls below the Mott insulator data after NΩ = 2. The blue shaded area around the Poisson line represents the uncertainty. **Plot (b): BEC Mode, Thermal, and Poisson Distributions** * **BEC mode (blue circles):** The data points are approximately: * (0, 0.02) * (2, 0.08) * (4, 0.14) * (5, 0.15) * (6, 0.14) * (8, 0.08) * (10, 0.05) * (12, 0.02) * (14, 0.01) * (16, 0.005) * (18, 0.002) The BEC mode distribution shows a peak around NΩ = 5. The blue shaded area around the Poisson line represents the uncertainty. * **Thermal (black dashed-dotted line):** The thermal distribution starts high and decreases monotonically with increasing atom number. * **Poisson (blue dashed line):** The Poisson distribution shows a peak around NΩ = 5, similar to the BEC mode, but the peak is sharper and the distribution decays faster. ### Key Observations * In plot (a), the Mott insulator distribution is broader than the Poisson distribution. * In plot (b), the BEC mode distribution is similar to the Poisson distribution but has a slightly broader peak. * The thermal distribution decreases monotonically in both plots. ### Interpretation The plots illustrate the probability distributions of detected atom numbers for different quantum states. Plot (a) shows that the Mott insulator distribution is broader than the Poisson distribution, indicating a more uniform distribution of atoms across lattice sites. Plot (b) shows that the BEC mode distribution is similar to the Poisson distribution, but with a broader peak, suggesting a less coherent state than a pure Bose-Einstein condensate. The thermal distribution represents a classical, incoherent state, which is characterized by a monotonically decreasing probability with increasing atom number. The shaded regions around the Poisson distribution in both plots, and the Mott insulator data in plot (a), represent the uncertainty in the measurements or theoretical calculations. </details> [51]. In the experiments reported here, each detected He ∗ atom is destroyed ( i.e. it decays to its ground-state) and the results of quantum optics are directly transposed. As explained in section II, fully-contrasted correlation functions are measured only when computed in voxels of small size δk 2 π/L . This requirement is more stringent than the one for measuring the probability distributions (see Appendix B). The results are plotted in Fig. 1(b) and are found in excellent quantitative agreement with the prediction for thermal states, g ( n ) (0) = n !. They represent a significant progress with respect to the literature where 2- [40] and 3-body [34] correlations had only been measured with limited amplitudes g ( n ) (0) < n !. The momentum-space FCS of a Mott insulator is thus identical to that of a statistical mixture of thermal bosons. Their full correlation functions however differ in their sizes l ( n ) c which are determined by their different in-trap sizes and the incompressible ( resp. compressible) character of Mott insulators ( resp. thermal gases) [34]. In a second set of experiments, we address the n -body coherence of the BEC component in lattice superfluids with N = 5(1) × 10 3 at U/J = 5. In the momentum space, the BEC occupies a volume of width ∆ k 0 . 15 k d centered at k = 0 [36]. This volume contains many atoms in each run and the statistics of the atom number falling in a sphere S Ω of a radius δk = 0 . 025 k d ∆ k , centered at k = 0 , is sufficient to extract the counting statistics (see Appendix C). We measure the counting statistics P ( N Ω ) in S Ω over about 2000 experimental runs, with the results plotted in Fig. 2(b). It is compared with Poisson and thermal statistics whose average atom number 〈 N Ω 〉 = 5 . 3(2) is that measured in the experiment. The counting statistics in the BEC mode is close to the Poisson FCS and clearly differs from the thermal FCS. This is confirmed by the measured values of g ( n ) (0) calculated from the normalized factorial moments of N Ω . In Fig. 1(b) we plot g ( n ) (0) as a function of n and find that g ( n ) (0) ∼ 1 at any order n in the BEC mode. These results, predicted by Glauber for a fully coherent sample, are in striking contrast with those of the Mott state, a difference that illustrates the outstanding capabilities of the FCS measured after an expansion to reveal the n -body coherence. Surprisingly, however, our measurements in the BEC mode deviate from the predictions for a coherent state and from a previous observation [33]: the deviation in the FCS (see Fig. 2(b)) is reflected in the fact that g ( n ) (0) > 1, as shown in Fig. 1(c). To verify that the observed deviations are statistically meaningful, we apply our computation of the n -body correlations to a randomized set, with the same numbers of atoms and of runs. This randomized set is obtained by randomly shuffling the detected atoms across the experimental runs. Doing so, atom correlations present within individual runs ( i.e. before shuffling) should vanish, and a Poisson statistics should be observed as a result of the discrete nature of our detection method applied to fully independent events. Indeed we find g ( n ) (0) = 1 . 00(2) at any order n of the randomized ensemble (see the orange squares in Fig. 1(c)), confirming that the deviations in the (nonrandomized) experimental data are significant. The randomization method also yields a Poisson statistics when applied to the Mott insulator data set. Note that the results of the randomization method validate the algorithm used to compute the n -body correlations and provide a means to test the accuracy of the measured statistics (see Appendix D). In the next section, we propose an interpretation of the deviations g ( n ) (0) > 1 revealed by our experiment in the case of lattice superfluids. ## IV. COHERENT FRACTION OF BOSE SUPERFLUIDS IN THE BEC MODE At thermal equilibrium, a fraction of the total atom number is expelled from the BEC by the finite interactions (quantum depletion) and by the finite temperature FIG. 3. (a) -(b) Plots of g ( n ) (0) measured at k = 0 in lattice superfluids with U/J = 5 and U/J = 20. The data in panel (a) is identical to that of shown in Fig. 1(c). The dashed lines ( resp. the shaded areas) are the predictions of the model with the values f coh ( resp. the uncertainties on f coh ) fitted to the data. (c) -(d) Plots of 1D cuts through the momentum densities measured at U/J = 5 and U/J = 20 and normalized to their value at k = 0. The vertical shaded area indicates the volume occupied by the sphere S Ω where the counting statistics is evaluated. The horizontal dashed-dotted lines indicate the fitted values 1 -f coh . The dashed lines are Lorentzian functions fitted to the tails of the densities in the range [0 . 2 k d , 0 . 5 k d ] in order to estimate the density of the depletion at k = 0 (the shaded areas represent the error from the fitted parameters). <details> <summary>Image 3 Details</summary>  ### Visual Description ## Combined Line and Scatter Plots: Correlation Analysis of g(n)(0) and Density ρ(k) ### Overview The image presents four plots arranged in a 2x2 grid. The top row (a and c) displays data for U/J = 5, while the bottom row (b and d) shows data for U/J = 20. Plots (a) and (b) are scatter plots showing the relationship between "Order n" and "g(n)(0)", with error bars and dashed trend lines. Plots (c) and (d) are line plots on a semi-log scale, showing "Density ρ(k)" as a function of "k/kd". Each line plot also includes a dashed horizontal line indicating "1 - fcoh". ### Components/Axes **Plot (a):** * **Type:** Scatter plot * **Title:** (a) - located at the top-left corner of the plot. * **X-axis:** Order n (values: 1 to 6) * **Y-axis:** g(n)(0) (values: 1.00 to 1.15) * **Data Series:** * U/J = 5 (blue circles with error bars) * fcoh = 0.9960(5) (dashed black line) * **Legend:** Located at the top of the plot. **Plot (b):** * **Type:** Scatter plot * **Title:** (b) - located at the top-left corner of the plot. * **X-axis:** Order n (values: 1 to 6) * **Y-axis:** g(n)(0) (values: 1.00 to 2.25) * **Data Series:** * U/J = 20 (gray circles with error bars) * fcoh = 0.971(1) (dashed black line) * **Legend:** Located at the top of the plot. **Plot (c):** * **Type:** Line plot (semi-log scale) * **Title:** (c) - located at the top-left corner of the plot. * **X-axis:** k/kd (values: -0.4 to 0.4) * **Y-axis:** Density ρ(k) (a.u.) (values: 10^-4 to 10^0) * **Data Series:** * U/J = 5 (solid blue line with shaded region) * 1 - fcoh = 0.004 (dashed black line) * **Legend:** Located at the top of the plot. **Plot (d):** * **Type:** Line plot (semi-log scale) * **Title:** (d) - located at the top-left corner of the plot. * **X-axis:** k/kd (values: -0.4 to 0.4) * **Y-axis:** Density ρ(k) (a.u.) (values: 10^-4 to 10^0) * **Data Series:** * U/J = 20 (solid gray line with shaded region) * 1 - fcoh = 0.029 (dashed black line) * **Legend:** Located at the top of the plot. ### Detailed Analysis **Plot (a):** * The blue data points (U/J = 5) show an upward trend. * Order n = 1, g(n)(0) ≈ 1.00 * Order n = 2, g(n)(0) ≈ 1.01 * Order n = 3, g(n)(0) ≈ 1.03 * Order n = 4, g(n)(0) ≈ 1.05 * Order n = 5, g(n)(0) ≈ 1.09 * Order n = 6, g(n)(0) ≈ 1.12 * The dashed black line (fcoh = 0.9960(5)) also shows an upward trend, closely following the data points. **Plot (b):** * The gray data points (U/J = 20) show an upward trend. * Order n = 1, g(n)(0) ≈ 1.00 * Order n = 2, g(n)(0) ≈ 1.08 * Order n = 3, g(n)(0) ≈ 1.22 * Order n = 4, g(n)(0) ≈ 1.42 * Order n = 5, g(n)(0) ≈ 1.73 * Order n = 6, g(n)(0) ≈ 2.15 * The dashed black line (fcoh = 0.971(1)) also shows an upward trend, closely following the data points. **Plot (c):** * The solid blue line (U/J = 5) shows a peak around k/kd = 0. * At k/kd = 0, Density ρ(k) ≈ 10^0 * At k/kd = -0.4 and 0.4, Density ρ(k) ≈ 10^-4 * The dashed black line (1 - fcoh = 0.004) is horizontal. **Plot (d):** * The solid gray line (U/J = 20) shows a peak around k/kd = 0. * At k/kd = 0, Density ρ(k) ≈ 10^0 * At k/kd = -0.4 and 0.4, Density ρ(k) ≈ 10^-4 * The dashed black line (1 - fcoh = 0.029) is horizontal. ### Key Observations * In plots (a) and (b), g(n)(0) increases with Order n for both U/J = 5 and U/J = 20. The rate of increase appears higher for U/J = 20. * In plots (c) and (d), the density ρ(k) peaks at k/kd = 0 for both U/J = 5 and U/J = 20. The peak appears sharper for U/J = 20. * The values of "1 - fcoh" are different for U/J = 5 and U/J = 20. ### Interpretation The plots suggest a correlation between the order 'n' and the function g(n)(0), with a steeper increase observed for U/J = 20 compared to U/J = 5. This indicates that the system's behavior is more sensitive to changes in 'n' when U/J is larger. The density plots show a central peak, indicating a higher probability of finding particles with momentum close to zero. The difference in "1 - fcoh" values suggests that the coherence properties of the system are influenced by the ratio U/J. Specifically, a higher U/J value (20) corresponds to a larger "1 - fcoh" value, which may indicate a decrease in coherence. </details> (thermal depletion). Even if it amounts to a negligible fraction of the atom number N Ω falling in S Ω , the total depletion of the condensate may contribute to the measured statistics. We define the 'coherent fraction' f coh as the fraction of N Ω that belongs to the BEC. Building on these considerations, we introduce a heuristic model that assumes (i) that atoms in the BEC and in the depletion contribute independently to the measured counting statistics in S Ω [44], and (ii) that the BEC is a coherent state while both the thermal and quantum depletion exhibit thermal statistics in S Ω . While a thermal FCS is expected for the thermal depletion, we emphasize that describing the contribution of the quantum depletion with a thermal statistics is an assumption. The statistics of the quantum depletion was shown to be thermal when measured at non-zero momenta, outside the BEC [35]. Whether this is also valid for small momenta within the BEC mode is difficult to assess in a harmonic trap. Finally, there is a priori no reason for our model to be valid in the strongly-interacting regime. In contrast to the weakly-interacting regime where the Bogoliubov approximation holds, neglecting the correlations between the condensate and its depletion at large interaction strengths could be an erroneous assumption. Our measurements however show that the model agrees with the experimental data over a wide range of interaction strengths. With the hypotheses of our model, we obtain an analytical prediction for g ( n ) (0) that depends only on the coherent fraction f coh (see Appendix E), $$g ^ { ( n ) } ( 0 ) - 1 = \sum _ { p = 1 } ^ { n - 1 } \left [ ( n - p ) ! \binom { n } { p } ^ { 2 } - \binom { n } { p } \right ] f _ { c o h } ^ { p } ( 1 - f _ { c o h } ) ^ { n - p } <text><loc_434><loc_205><loc_499><loc_255>(3) While</text> </doctag>$$ While our model straightforwardly predicts the magnitudes g ( n ) (0), this is not the case for the probability distribution P ( N Ω ). The latter results from a complex convolution of the probability distributions of the condensate and of the depletion and it is well established that obtaining P ( X ) from the moments of a random variable X is a difficult problem [52]. In Fig. 3(a), we fit the experimental data shown in Fig. 1(b) with the analytical prediction of Eq. (3), finding a good agreement with the coherent fraction as the only adjustable parameter (found equal to f c = 0 . 9960(5) for the case of U/J = 5). We then repeated our measurements for various condensate fractions. To do so, we change the lattice depth to obtain ratios of on-site interaction U to tunnelling amplitude J ranging from U/J = 2 to U/J = 22. In this range of parameters, the gas remains far from entering the Mott insulator regime expected at the critical ratio ( U/J ) c 25 -30 [49], but it enters the strongly-interacting regime where the condensate is strongly depleted (at U/J = 22 the condensate fraction is f c ∼ 0 . 15). In Fig. 3(b) we plot the magnitude of the n -body correlations for a lattice superfluids with an increased interaction U/J = 20. The deviation from the ideal coherent state is increased, in qualitative agreement with the physical picture at the root of the model. To be quantitative, we fit the experimental data with the analytical prediction of Eq. (3). Firstly, we confirm that Eq. (3) correctly fits the values of g ( n ) (0) with a single adjustable parameter f coh . Secondly, the extracted values of f coh decrease with the interaction strength as intuitively expected. The uncertainty on the values f coh is extremely small, at the ∼ 0 . 1% level. As can be inferred from Fig. 3, the larger the order n of correlations we measure, the smaller the uncertainty on f coh . This illustrates the extreme sensitivity of high-order correlations to probe many-body coherence. A quantitative test of the model would compare the value 1 -f coh to the fraction η D of depleted atoms detected within S Ω . We are not aware of a quantitative analytical prediction for η D in 3D interacting trapped lattice Bose gases. However, an indirect comparison is amenable from measuring the momentum densities. In Fig 3(c)-(d), we plot 1D cuts through the momentum densities measured at U/J = 5 and U/J = 20 and we fit the tails (in the range [0 . 2 k d , 0 . 5 k d ]) with a Lorentzian function to extrapolate the density of the depletion at k = 0. Using a Lorentzian function is an arbitrary choice which happens to correctly fit the tails. This analysis indicates that the values 1 -f coh are compatible with the extrapolated densities, while both quantities vary by one order of magnitude. To assess the validity of the model, we proceed with a quantitative comparison of the measured coherent fraction to the measured condensate fraction [47]. The coherent fraction f coh in S Ω differs from the condensate fraction f c of the entire gas since the radius of S Ω is much smaller than 2 π/L . Moreover, the scaling of f coh with f c is non-linear since f coh reveals the BEC atom number in a small volume S Ω where the BEC contribution is maximum while that of the depletion is small (except when f c 1). The volumes occupied by the BEC ( V c ) and by the depletion ( V d ) set their respective contributions in S Ω . A simple estimate (see Appendix F) leads to $$f _ { c o h } \simeq \frac { f _ { c } } { f _ { c } + ( 1 - f _ { c } ) \mathcal { V } _ { c } / \mathcal { V } _ { d } } . \quad \ \ ( 4 )$$ In Fig. 4 we plot the measured values of f coh and f c , along with the prediction of Eq. (4). Here, f c is obtained similarly to [47] and the ratio V c / V d used to plot Eq. (4) is estimated from the measured density pro- FIG. 4. Coherent fraction f coh as a function of the condensate fraction f c . The dashed-line is the prediction of Eq. (4) where the ratio V c / V d = 0 . 03(1) is evaluated from the measured density profiles (see Appendix F). The shaded area reflects the uncertainty on V c / V d . <details> <summary>Image 4 Details</summary>  ### Visual Description ## Scatter Plot: f_coh vs f_c ### Overview The image is a scatter plot showing the relationship between two variables, f_coh (y-axis) and f_c (x-axis). The plot includes data points with error bars, a dashed orange trend line, and a shaded region around the trend line indicating uncertainty. An inset plot provides a zoomed-in view of the data at higher values of f_c. ### Components/Axes * **X-axis:** Labeled "f_c", with a scale from 0 to 1.0, marked at 0.2, 0.4, 0.6, 0.8, and 1.0. * **Y-axis:** Labeled "f_coh", with a scale from 0.85 to 1.00, marked at 0.85, 0.90, 0.95, and 1.00. * **Data Points:** Orange circles with error bars, indicating the uncertainty in the measurements. * **Trend Line:** A dashed orange line showing the general trend of the data. * **Uncertainty Region:** A shaded orange region around the trend line, representing the confidence interval. * **Inset Plot:** A smaller plot in the lower-right corner, zooming in on the data points with f_c values between approximately 0.6 and 1.0. The inset plot's y-axis ranges from 0.99 to 1.00. ### Detailed Analysis **Main Plot Data Points:** * **Point 1:** f_c ≈ 0.15, f_coh ≈ 0.85. * **Point 2:** f_c ≈ 0.4, f_coh ≈ 0.93. * **Point 3:** f_c ≈ 0.65, f_coh ≈ 0.98. * **Point 4:** f_c ≈ 0.8, f_coh ≈ 0.99. * **Point 5:** f_c ≈ 0.9, f_coh ≈ 0.995. * **Point 6:** f_c ≈ 1.0, f_coh ≈ 0.998. **Inset Plot Data Points:** * **Point 1:** f_c ≈ 0.65, f_coh ≈ 0.991. * **Point 2:** f_c ≈ 0.8, f_coh ≈ 0.997. * **Point 3:** f_c ≈ 0.9, f_coh ≈ 0.998. * **Point 4:** f_c ≈ 1.0, f_coh ≈ 0.999. **Trend Verification:** The trend line in both the main plot and the inset plot slopes upward, indicating a positive correlation between f_c and f_coh. The rate of increase appears to diminish as f_c approaches 1.0. ### Key Observations * There is a positive correlation between f_c and f_coh. As f_c increases, f_coh also increases. * The rate of increase in f_coh diminishes as f_c approaches 1.0, suggesting a saturation effect. * The inset plot provides a more detailed view of the data at higher f_c values, showing that f_coh approaches 1.0 asymptotically. * The error bars indicate some uncertainty in the measurements, but the overall trend is clear. ### Interpretation The data suggests that f_coh is dependent on f_c, with a positive relationship between the two variables. The saturation effect observed at higher f_c values may indicate a physical limit or a diminishing return on increasing f_c. The uncertainty region around the trend line provides a visual representation of the confidence in the relationship between f_c and f_coh. The inset plot confirms the asymptotic behavior of f_coh as f_c approaches 1.0. </details> files (see Appendix F). The quantitative prediction of our model (without any adjustable parameter) correctly reproduces the observed non-linear dependency of f coh with f c . As discussed previously, the agreement in the strongly-interacting regime is rather surprising. Understanding this fact requires a more elaborate theoretical approach than the one presented in this work. We observe that the heuristic model quantitatively captures the deviations to g ( n ) (0) = 1, attributing the latter to the presence of depleted atoms (both from the thermal and the quantum depletion). In turn, this suggests that BEC atoms have the statistics of a fully coherent state, with g ( n ) (0) = 1 at any order n ≤ 6. The role of the depletion in the many-body coherence unveiled in this work was not identified previously [33]. A major difference of our experiment with respect to that of [33] is that the volume V Ω used in that work to compute the statistics is larger than the coherence volume (2 π/L ) 3 , a choice made in light of a smaller detection efficiency η ∼ 0 . 08. Using too large of a volume V Ω entails a convolution of the correlation functions that significantly reduces their amplitudes and, in turn, hides the contribution from the depletion. This is illustrated by the fact that in the thermal case the authors of [33] find g (2) (0) = 1 . 022(2) and g (3) (0) = 1 . 061(6) instead of g (2) (0) = 2 and g (3) (0) = 6 . We conclude that our ability to obtain sufficient statistics in tiny volumes is a major asset to quantitatively probe many-body correlations. ## V. CONCLUSION We have presented measurements of the full counting statistics (FCS) and of perfectly-contrasted n -body correlations in interacting lattice Bose gases. We find that a Mott insulator exhibits a thermal statistics in the mo- mentum space with g ( n ) (0) = n !, while the condensate component is fully coherent with g ( n ) (0) = 1 (at least) up to n = 6. The latter conclusion derives from assuming that the heuristic model we introduce to analyse the FCS in Bose superfluids and to describe the contribution from the depletion of the condensate is valid. We have assessed this validity in the experiment from studying the coherent fraction at increasing interaction strengths. A theoretical validation is beyond the scope of this work. Our results represent the most stringent certifications of the many-body coherence properties of Mott insulators and of BECs, and they validate emblematic pure-state descriptions of n -body coherence up to the order n = 6. In the future, the experimental approach we use is readily extendable to probe the many-body coherence in a large variety of interacting quantum states and phase transitions realized on cold-atom platforms ( e.g. [911, 53]). It relies on (i) a free-fall expansion from the trap and (ii) a single-atom-resolved detection method. The first condition is ensured with lattice or low-dimensional (1D and 2D) gases for which interactions do not affect the expansion, as well as with quantum gases for which interactions are switched off during the expansion [54]. The second condition is met with various detection methods and atomic species [55]. - [1] Y. Blanter and M. B¨ uttiker, Shot noise in mesoscopic conductors, Physics Reports 336 , 1 (2000). - [2] C. Gardiner and P. Zoller, Quantum Noise - A Handbook of Markovian and Non-Markovian Quantum Stochastic Methods with Applications to Quantum Optics (Springer Berlin, Heidelberg, 2004). - [3] E. A. Burt, R. W. Ghrist, C. J. Myatt, M. J. Holland, E. A. Cornell, and C. E. Wieman, Coherence, correlations, and collisions: What one learns about bose-einstein condensates from their decay, Phys. Rev. Lett. 79 , 337 (1997). - [4] E. Altman, E. Demler, and M. D. Lukin, Probing manybody states of ultracold atoms via noise correlations, Phys. Rev. A 70 , 013603 (2004). - [5] T. Schweigler, V. Kasper, S. Erne, I. Mazets, B. Rauer, F. Cataldini, T. Langen, T. Gasenzer, J. Berges, and J. Schmiedmayer, Experimental characterization of a quantum many-body system via higher-order correlations, Nature 545 , 323 (2017). - [6] L. S. Levitov, H. Lee, and G. B. Lesovik, Electron counting statistics and coherent states of electric current, Journal of Mathematical Physics 37 , 4845 (1996). - [7] D. A. Ivanov and A. G. Abanov, Phase transitions in full counting statistics for periodic pumping, EPL (Europhysics Letters) 92 , 37008 (2010). - [8] F. J. G´ omez-Ruiz, J. J. Mayo, and A. del Campo, Full counting statistics of topological defects after crossing a phase transition, Phys. Rev. Lett. 124 , 240602 (2020). - [9] P. Devillard, D. Chevallier, P. Vignolo, and M. Albert, Full counting statistics of the momentum occupation numbers of the tonks-girardeau gas, Phys. Rev. A 101 , 063604 (2020). Another interesting direction consists in realizing recent proposals to access non-trivial n -body correlations [26, 27]. These theoretical works have shown that nondiagonal correlators can be accessed from combining twoparticle unitary transformations e.g. beamsplitters with a measurement of the FCS. Such protocols have many potential applications, from quantifying entanglement to implementing variational algorithms. ## ACKNOWLEDGEMENT We thank S. Butera, A. Browaeys, I. Carusotto, T. Lahaye, T. Roscilde and the members of the Quantum Gas group at Institut d'Optique for insightful discussions. We acknowledge financial support from the R´ egion Ile-de-France in the framework of the DIM SIRTEQ, the 'Fondation d'entreprise iXcore pour la Recherche', the Agence Nationale pour la Recherche (Grant number ANR-17-CE30-0020-01). D.C. acknowledges support from the Institut Universitaire de France. A. A. acknowledges support from the Simons Foundation and from Nokia-Bell labs. - [10] V. Eisler, Universality in the full counting statistics of trapped fermions, Phys. Rev. Lett. 111 , 080402 (2013). - [11] I. Lovas, B. D´ ora, E. Demler, and G. Zar´ and, Full counting statistics of time-of-flight images, Phys. Rev. A 95 , 053621 (2017). - [12] I. Klich and L. Levitov, Quantum noise as an entanglement meter, Phys. Rev. Lett. 102 , 100502 (2009). - [13] M. Esposito, U. Harbola, and S. Mukamel, Nonequilibrium fluctuations, fluctuation theorems, and counting statistics in quantum systems, Rev. Mod. Phys. 81 , 1665 (2009). - [14] S. Gustavsson, R. Leturcq, B. Simoviˇ c, R. Schleser, T. Ihn, P. Studerus, K. Ensslin, D. C. Driscoll, and A. C. Gossard, Counting statistics of single electron transport in a quantum dot, Phys. Rev. Lett. 96 , 076605 (2006). - [15] V. F. Maisi, D. Kambly, C. Flindt, and J. P. Pekola, Full counting statistics of andreev tunneling, Phys. Rev. Lett. 112 , 036801 (2014). - [16] T. C. Liebisch, A. Reinhard, P. R. Berman, and G. Raithel, Atom counting statistics in ensembles of interacting rydberg atoms, Phys. Rev. Lett. 95 , 253002 (2005). - [17] N. Malossi, M. M. Valado, S. Scotto, P. Huillery, P. Pillet, D. Ciampini, E. Arimondo, and O. Morsch, Full counting statistics and phase diagram of a dissipative rydberg gas, Phys. Rev. Lett. 113 , 023006 (2014). - [18] J. Zeiher, R. van Bijnen, P. Schauß, S. Hild, J.-y. Choi, T. Pohl, I. Bloch, and C. Gross, Many-body interferometry of a rydberg-dressed spin lattice, Nature Physics 12 , 1095 (2016). - [19] A. ¨ Ottl, S. Ritter, M. K¨ ohl, and T. Esslinger, Correlations and counting statistics of an atom laser, Phys. Rev. Lett. 95 , 090404 (2005). - [20] M. Perrier, Z. Amodjee, P. Dussarrat, A. Dareau, A. Aspect, M. Cheneau, D. Boiron, and C. I. Westbrook, Thermal counting statistics in an atomic two-mode squeezed vacuum state, SciPost Phys. 7 , 2 (2019). - [21] S. T. Flammia, D. Gross, Y.-K. Liu, and J. Eisert, Quantum tomography via compressed sensing: error bounds, sample complexity and efficient estimators, New Journal of Physics 14 , 095022 (2012). - [22] R. J. Glauber, The quantum theory of optical coherence, Phys. Rev. 130 , 2529 (1963). - [23] P. E. Dolgirev, J. Marino, D. Sels, and E. Demler, Nongaussian correlations imprinted by local dephasing in fermionic wires, Phys. Rev. B 102 , 100301 (2020). - [24] C. Fabre and N. Treps, Modes and states in quantum optics, Rev. Mod. Phys. 92 , 035005 (2020). - [25] M. Walschaers, Non-gaussian quantum states and where to find them, PRX Quantum 2 , 030204 (2021). - [26] E. Brunner, A. Buchleitner, and G. Dufour, Many-body coherence and entanglement from randomized correlation measurements, arXiv 10.48550/arxiv.2107.01686 (2021). - [27] P. Naldesi, A. Elben, A. Minguzzi, D. Cl´ ement, P. Zoller, and B. Vermersch, Fermionic correlation functions from randomized measurements in programmable atomic quantum devices, arXiv 10.48550/arxiv.2205.00981 (2022). - [28] A. Aspect, Hanbury brown and twiss, hong ou and mandel effects and other landmarks in quantum optics: from photons to atoms, in Current Trends in Atomic Physics , edited by A. Browaeys, T. Lahaye, T. Porto, C. S. Adams, M. Weidem¨ uller, and L. F. Cugliandolo (Oxford University Press, 2019). - [29] M. Schellekens, R. Hoppeler, A. Perrin, J. V. Gomes, D. Boiron, A. Aspect, and C. I. Westbrook, Hanbury brown twiss effect for ultracold quantum gases, Science 310 , 648 (2005). - [30] R. G. Dall, A. G. Manning, S. S. Hodgman, W. RuGway, K. V. Kheruntsyan, and A. G. Truscott, Ideal n-body correlations with massive particles, Nature Physics 9 , 341 (2013). - [31] T. Jeltes, J. M. McNamara, W. Hogervorst, W. Vassen, V. Krachmalnicoff, M. Schellekens, A. Perrin, H. Chang, D. Boiron, A. Aspect, and C. I. Westbrook, Comparison of the hanbury brown-twiss effect for bosons and fermions, Nature 445 , 402 (2007). - [32] A. Bergschneider, V. M. Klinkhamer, J. H. Becher, R. Klemt, L. Palm, G. Z¨ urn, S. Jochim, and P. M. Preiss, Experimental characterization of two-particle entanglement through position and momentum correlations, Nature Physics 15 , 640 (2019). - [33] S. S. Hodgman, R. G. Dall, A. G. Manning, K. G. H. Baldwin, and A. G. Truscott, Direct measurement of long-range third-order coherence in bose-einstein condensates, Science 331 , 1046 (2011). - [34] C. Carcy, H. Cayla, A. Tenart, A. Aspect, M. Mancini, and D. Cl´ ement, Momentum-space atom correlations in a mott insulator, Phys. Rev. X 9 , 041028 (2019). - [35] H. Cayla, S. Butera, C. Carcy, A. Tenart, G. Herc´ e, M. Mancini, A. Aspect, I. Carusotto, and D. Cl´ ement, Hanbury brown and twiss bunching of phonons and of the quantum depletion in an interacting bose gas, Phys. Rev. Lett. 125 , 165301 (2020). - [36] A. Tenart, G. Herc´ e, J.-P. Bureik, A. Dareau, and D. Cl´ ement, Observation of pairs of atoms at opposite momenta in an equilibrium interacting bose gas, Nature Physics 17 , 1364 (2021). - [37] M. Holten, L. Bayha, K. Subramanian, S. Brandstetter, C. Heintze, P. Lunt, P. M. Preiss, and S. Jochim, Observation of cooper pairs in a mesoscopic two-dimensional fermi gas, Nature 606 , 287 (2022). - [38] L. Mandel and E. Wolf, Optical Coherence and Quantum Optics (Cambridge University Press, 1995). - [39] L. P. Pitaevskii and S. Stringari, Bose-Einstein Condensation (Clarendon press, Oxford, 2004). - [40] S. F¨ olling, F. Gerbier, A. Widera, O. Mandel, T. Gericke, and I. Bloch, Spatial quantum noise interferometry in expanding ultracold atom clouds, Nature 434 , 481 (2005). - [41] E. Toth, A. M. Rey, and P. B. Blakie, Theory of correlations between ultracold bosons released from an optical lattice, Phys. Rev. A 78 , 013627 (2008). - [42] D. F. Walls and G. J. Milburn, Quantum Optics (Springer Berlin, Heidelberg, 2008). - [43] M. A. Kristensen, M. B. Christensen, M. Gajdacz, M. Iglicki, K. Paw lowski, C. Klempt, J. F. Sherson, K. Rza˙ zewski, A. J. Hilliard, and J. J. Arlt, Observation of atom number fluctuations in a bose-einstein condensate, Phys. Rev. Lett. 122 , 163601 (2019). - [44] Y. Castin and R. Dum, Low-temperature bose-einstein condensates in time-dependent traps: Beyond the u (1) symmetry-breaking approach, Phys. Rev. A 57 , 3008 (1998). - [45] Note that in the experiment the total atom number N is not strictly constant from shot to shot. However, the measurement of the FCS in a tiny volume V Ω is largely immune to small shot-to-shot fluctuations of N . - [46] J. Stenger, S. Inouye, A. P. Chikkatur, D. M. StamperKurn, D. E. Pritchard, and W. Ketterle, Bragg spectroscopy of a bose-einstein condensate, Phys. Rev. Lett. 82 , 4569 (1999). - [47] H. Cayla, C. Carcy, Q. Bouton, R. Chang, G. Carleo, M. Mancini, and D. Cl´ ement, Single-atom-resolved probing of lattice gases in momentum space, Phys. Rev. A 97 , 061609 (2018). - [48] A. Tenart, C. Carcy, H. Cayla, T. Bourdel, M. Mancini, and D. Cl´ ement, Two-body collisions in the time-of-flight dynamics of lattice bose superfluids, Phys. Rev. Research 2 , 013017 (2020). - [49] C. Carcy, G. Herc´ e, A. Tenart, T. Roscilde, and D. Cl´ ement, Certifying the adiabatic preparation of ultracold lattice bosons in the vicinity of the mott transition, Phys. Rev. Lett. 126 , 045301 (2021). - [50] K. Laiho, T. Dirmeier, M. Schmidt, S. Reitzenstein, and C. Marquardt, Measuring higher-order photon correlations of faint quantum light: A short review, Physics Letters A 435 , 128059 (2022). - [51] J. W. Goodman, Statistical optics , edited by W.-I. New York, Vol. 1 (1985). - [52] N. I. Akhiezer, The classical moment problem: and Some Related Questions in Analysis , edited by O. B. L. Edinburgh and London (1965). - [53] D. Contessi, A. Recati, and M. Rizzi, Phase diagram detection via gaussian fitting of number probability distribution, arXiv , 2207.01478 (2022). - [54] I. Bloch, J. Dalibard, and W. Zwerger, Many-body physics with ultracold gases, Rev. Mod. Phys. 80 , 885 (2008). - [55] H. Ott, Single atom detection in ultracold quantum gases: a review of current progress, Reports on Progress in Physics 79 , 054401 (2016). ## Appendix A. PROBABILITY DISTRIBUTION IN MOTT INSULATORS AND MULTIMODE THERMAL STATISTICS We consider the thermal statistics found in the case of Mott insulators. The bell-shaped correlation functions are well-fitted by a Gaussian function [34] and we introduce the two-body correlation length l c as $$g ^ { ( 2 ) } ( k _ { 1 } , k _ { 2 } ) = g ^ { ( 2 ) } ( 0 ) \times \exp \, \left ( \frac { - 2 | k _ { 1 } - k _ { 2 } \, | ^ { 2 } } { l _ { c } ^ { 2 } } \right ) . \quad ( 5 )$$ For the data set of a Mott insulator shown in the main text, the two-body correlation length is l c /k d = 0 . 031(2) and corresponds to the expected value [34], l c /k d ∼ 1 /N site where N site is the number of lattice sites occupied by the Mott insulator in the trap. In the case of a thermal statistics, the two-body correlation length in momentum space is equal to the inverse trap size of the gas, i.e. to the size of a mode in momentum space. We define a mode from its coherence volume V c as a cubic voxel of size 2 l c . When one computes two-body correlations in a cubic volume V Ω of size δk comparable to, or larger than, the correlation length l c , the measured magnitude differs from the fully-contrasted value g (2) (0) due to the spatial averaging of g (2) ( k 1 , k 2 ). For instance, using V Ω = V c , the magnitude of two-body correlations is reduced by a factor 1 / ( √ π/ 8 erf[ √ 2]) 3 ∼ 4 . 7. In contrast, we find that computing the probability distribution P ( N Ω ) is not affected by considering a volume V Ω ∼ V c , as we shall now discuss. In Fig. 5 we plot the probability distributions P ( N Ω ) computed with volumes V Ω equal or larger than V c . Firstly, we observe that the statistics is thermal when V Ω ≤ V c . More specifically, there is no need to ensure the condition V Ω V c to measure the probability distribution of a thermal statistics. Secondly, this statement is confirmed by analysing the measured probability distributions P ( N Ω ) for V Ω > V c . Indeed, the probability distribution computed over a number M of independent modes, each of which exhibiting a thermal statistics with an average occupation 〈 N 〉 , reads [20, 51] $$P _ { M } ( N _ { \Omega } ) = \frac { ( \langle N _ { \Omega } \rangle + M - 1 ) ! } { \langle N _ { \Omega } \rangle ! ( M - 1 ) ! } \frac { ( \langle N _ { \Omega } \rangle / M ) ^ { N _ { \Omega } } } { ( 1 + \langle N \rangle / M ) ^ { N _ { \Omega } + M } } .$$ In Fig. 5 the prediction P M ( N Ω ) for a multi-mode thermal statistics is plotted with a number of modes equal to M = V Ω /V c . These predictions are found to be in excellent agreement with the measured distributions without any adjustable parameters. FIG. 5. Probability distributions P ( N Ω ) measured in a Mott insulator with a varying voxel size V Ω . For V Ω = 0 . 9 V c , the measured distribution (black dots) matches the prediction (black dashed-dotted line) for a thermal distribution. For V Ω > V c , the measured distributions (red squares, orange triangles and blue triangles) agree with the predictions P M ( N Ω ) for a multi-mode thermal distribution (red, orange and blue solid lines). No adjustable parameters are used to plot the predictions P M ( N Ω ) that depend on the measured mean atom number 〈 N Ω 〉 and the chosen ratio M = V Ω /V c . <details> <summary>Image 5 Details</summary>  ### Visual Description ## Chart: Probability Distribution vs. N_Ω for Different V_Ω Values ### Overview The image is a plot showing the probability distribution P(N_Ω) as a function of N_Ω for different values of V_Ω. The plot uses a logarithmic scale for the y-axis (P(N_Ω)) and a linear scale for the x-axis (N_Ω). There are four distinct data series, each representing a different value of V_Ω, and a dashed line. ### Components/Axes * **X-axis:** N_Ω (linear scale, ranging from 0 to 8) * **Y-axis:** P(N_Ω) (logarithmic scale, ranging from 10^-3 to 10^0) * **Legend:** Located in the center-left of the plot. * Black circles: V_Ω = 0.9V_c * Red squares: V_Ω = 2.1V_c * Orange triangles: V_Ω = 4.2V_c * Blue inverted triangles: V_Ω = 7.2V_c ### Detailed Analysis * **V_Ω = 0.9V_c (Black Circles):** The data points for this series show a decreasing trend as N_Ω increases. * N_Ω = 0, P(N_Ω) ≈ 0.7 * N_Ω = 2, P(N_Ω) ≈ 0.07 * N_Ω = 4, P(N_Ω) ≈ 0.007 * **V_Ω = 2.1V_c (Red Squares):** The data points for this series show a decreasing trend as N_Ω increases. * N_Ω = 0, P(N_Ω) ≈ 0.4 * N_Ω = 2, P(N_Ω) ≈ 0.15 * N_Ω = 4, P(N_Ω) ≈ 0.02 * N_Ω = 6, P(N_Ω) ≈ 0.003 * N_Ω = 8, P(N_Ω) ≈ 0.0001 * **V_Ω = 4.2V_c (Orange Triangles):** The data points for this series initially increase and then decrease as N_Ω increases. * N_Ω = 0, P(N_Ω) ≈ 0.2 * N_Ω = 2, P(N_Ω) ≈ 0.15 * N_Ω = 4, P(N_Ω) ≈ 0.05 * N_Ω = 6, P(N_Ω) ≈ 0.01 * N_Ω = 8, P(N_Ω) ≈ 0.005 * **V_Ω = 7.2V_c (Blue Inverted Triangles):** The data points for this series initially increase and then decrease as N_Ω increases. * N_Ω = 0, P(N_Ω) ≈ 0.07 * N_Ω = 2, P(N_Ω) ≈ 0.2 * N_Ω = 4, P(N_Ω) ≈ 0.15 * N_Ω = 6, P(N_Ω) ≈ 0.05 * N_Ω = 8, P(N_Ω) ≈ 0.02 * **Dashed Line (Gray):** This line shows a decreasing trend as N_Ω increases. * N_Ω = 2, P(N_Ω) ≈ 0.1 * N_Ω = 6, P(N_Ω) ≈ 0.001 ### Key Observations * As V_Ω increases, the probability distribution P(N_Ω) shifts towards higher values of N_Ω. * For lower values of V_Ω (0.9V_c and 2.1V_c), the probability distribution decreases monotonically with increasing N_Ω. * For higher values of V_Ω (4.2V_c and 7.2V_c), the probability distribution initially increases and then decreases with increasing N_Ω. * The dashed line has a steeper negative slope than the other lines. ### Interpretation The plot illustrates how the probability distribution of N_Ω changes with varying values of V_Ω. The shift towards higher N_Ω values as V_Ω increases suggests that larger values of V_Ω favor higher occupancy states. The non-monotonic behavior of the probability distribution for higher V_Ω values indicates a more complex relationship between V_Ω and N_Ω, possibly due to the emergence of multiple occupancy states. The dashed line may represent a theoretical prediction or a reference point for comparison. </details> ## Appendix B. MAGNITUDES OF CORRELATION FUNCTIONS IN MOTT INSULATORS ## A. Condition to measure fully-contrasted atom correlations As discussed in section II, obtaining fully-contrasted n -body correlation functions requires to probe the statistics in small voxels V Ω V c . This condition is much stringent than that ( V Ω ≤ V c ) to obtain the probability distribution P ( N Ω ). This poses serious difficulties in the case of a Mott insulator as its momentum-space density is low because a Mott insulator occupies a large volume in momentum-space, extending beyond the first Brillouin zone. The statistics of the atom number N Ω in a single voxel V Ω V c is therefore not sufficient to extract the magnitude g ( n ) (0) of the n -body correlation functions accurately, even when using an average over all the voxels contained in the first Brillouin zone. To circumvent this issue, we first compute the statistics in larger voxels from which the results in small voxels V Ω V c are then extrapolated. This procedure is detailed below. ## B. Magnitudes of correlation functions measured in voxels V Ω ∼ V c We first evaluate the magnitudes of correlation functions in anisotropic voxels of volume V ( δk ⊥ ) = δk × δk 2 ⊥ with δk ⊥ ≥ l c > δk = 0 . 015 k d . The volume V ( δk ⊥ ) ∼ V c is sufficiently large to compute the statistics of N δk ⊥ . To proceed, we divide the first Brillouin zone in voxels of volume V ( δk ⊥ ) and label each voxel with an integer j . The factorial moment in one voxel is denoted 〈 N δk ⊥ ( N δk ⊥ -1) ... ( N δk ⊥ -n + 1) 〉 j and the average of the factorial moments over the first Brillouin zone reads $$\sum _ { j } \, \langle N _ { \delta k _ { \perp } } ( N _ { \delta k _ { \perp } } - 1 ) \dots ( N _ { \delta k _ { \perp } } - n + 1 ) \rangle _ { j } . \quad ( 7 )$$ We obtain the magnitudes of correlation functions from normalizing the average value of the factorial moments, $$g _ { \delta k _ { \perp } } ^ { ( n ) } ( 0 ) = \frac { \sum _ { j } \, \langle N _ { \delta k _ { \perp } } ( N _ { \delta k _ { \perp } } - 1 ) \dots ( N _ { \delta k _ { \perp } } - n + 1 ) \rangle _ { j } } { \sum _ { j } \, \langle N _ { \delta k _ { \perp } } \rangle _ { j } ^ { n } } .$$ In the case of two-body correlations, an alternative to the approach based on the factorial moments is to compute the histogram of atom pairs [34]. We have verified that both approaches yield similar results for the magnitude g (2) δk ⊥ (0). ## C. Magnitudes of correlation functions measured in voxels V Ω V c We repeat the above analysis for voxels of varying volume, by changing the transverse size δk ⊥ , with the results plotted in Fig. 6. These results illustrate the fact that the magnitudes of the correlation functions are reduced when the voxel used to compute the statistics is too large, i.e. of the order or larger than the coherence volume V c . The fully-contrasted magnitude g ( n ) (0), expected in voxels V Ω V c , is the value of g ( n ) δk ⊥ (0) in the limit δk ⊥ → 0. We fit the measured variations of g ( n ) δk ⊥ (0) with δk ⊥ (see dashed-line in Fig. 6) and extrapolate g ( n ) (0) in the limit δk ⊥ → 0. ## Appendix C. PROBABILITY DISTRIBUTION AND MAGNITUDES OF CORRELATION FUNCTIONS IN BOSE SUPERFLUIDS The magnitudes of the normalized correlation functions are obtained from the factorial moments of the atom number N Ω in the small volume V Ω V c . In the case of the BEC mode at k = 0 , the average number of atoms in V Ω is relatively large, 〈 N Ω 〉 ∼ 5. As a consequence, the statistics in V Ω is sufficient to extract accurately the factorial moments. The n -factorial moment is defined as 〈 ( N Ω ) n 〉 = 〈 N Ω ( N Ω -1) ... ( N Ω -n +1) 〉 and g ( n ) (0) are expressed as $$g ^ { ( n ) } ( 0 ) = \frac { \langle N _ { \Omega } ( N _ { \Omega } - 1 ) \dots ( N _ { \Omega } - n + 1 ) \rangle } { \langle N _ { \Omega } \rangle ^ { n } } = \frac { \langle ( N _ { \Omega } ) _ { n } \rangle } { \langle N _ { \Omega } \rangle ^ { n } } \quad t h e r ( 9 ) \quad s t a t i s$$ ⊥ FIG. 6. Plot of the magnitude g ( n ) δk ⊥ (0) as a function of the transverse integration δk ⊥ . The different panels correspond to different orders n of correlation functions, from n = 2 to n = 6. Apart from the data for n = 2 plotted in linear scale, the vertical axes are in log scale. <details> <summary>Image 6 Details</summary>  ### Visual Description ## Scatter Plot Grid: Correlation Functions vs. Momentum Difference ### Overview The image presents a grid of five scatter plots. Each plot displays the correlation function g^(n)_k⊥(0) as a function of the momentum difference δk⊥[kd] for different values of 'n' (n=2, 3, 4, 5, and 6). The y-axis is linear for n=2 and logarithmic for n=3, 4, 5, and 6. Each data point is represented by an orange circle with error bars. A dashed black line is overlaid on each plot, representing a theoretical fit to the data. ### Components/Axes * **Arrangement:** The plots are arranged in a 2x3 grid, with the bottom-right position empty. * **X-axis:** * Label: δk⊥[kd] * Scale: Linear, ranging from 0.0 to approximately 0.3 for n=2, and from 0.0 to approximately 0.6 for n=3, 4, 5, and 6. * **Y-axis:** * Label: g^(n)_k⊥(0) * Scale: * n=2: Linear, ranging from 1.0 to 2.0. * n=3, 4, 5, 6: Logarithmic (base 10). * n=3: Ranges from 10^0 (1) to 10^1 (10). * n=4: Ranges from 10^0 (1) to 10^2 (100). * n=5: Ranges from 10^0 (1) to 10^2 (100). * n=6: Ranges from 10^0 (1) to 10^3 (1000). * **Data Points:** Orange circles with error bars. * **Fit Line:** Dashed black line representing a theoretical fit. * **Titles:** Each plot has a title indicating the value of 'n': n=2 (top-left), n=3 (top-right), n=4 (middle-left), n=5 (middle-right), n=6 (bottom-left). ### Detailed Analysis **Plot for n=2:** * Trend: The correlation function decreases rapidly as δk⊥[kd] increases. * Data Points: * At δk⊥[kd] = 0.0, g^(2)_k⊥(0) ≈ 2.0. * At δk⊥[kd] = 0.1, g^(2)_k⊥(0) ≈ 1.2. * At δk⊥[kd] = 0.2, g^(2)_k⊥(0) ≈ 1.05. * At δk⊥[kd] = 0.3, g^(2)_k⊥(0) ≈ 1.0. **Plot for n=3:** * Trend: The correlation function decreases rapidly as δk⊥[kd] increases. * Data Points: * At δk⊥[kd] = 0.0, g^(3)_k⊥(0) ≈ 8.0. * At δk⊥[kd] = 0.1, g^(3)_k⊥(0) ≈ 1.5. * At δk⊥[kd] = 0.2, g^(3)_k⊥(0) ≈ 0.8. * At δk⊥[kd] = 0.6, g^(3)_k⊥(0) ≈ 0.3. **Plot for n=4:** * Trend: The correlation function decreases rapidly as δk⊥[kd] increases, then plateaus. * Data Points: * At δk⊥[kd] = 0.0, g^(4)_k⊥(0) ≈ 150. * At δk⊥[kd] = 0.1, g^(4)_k⊥(0) ≈ 10. * At δk⊥[kd] = 0.2, g^(4)_k⊥(0) ≈ 2. * For δk⊥[kd] > 0.3, g^(4)_k⊥(0) ≈ 1. **Plot for n=5:** * Trend: The correlation function decreases rapidly as δk⊥[kd] increases, then plateaus. * Data Points: * At δk⊥[kd] = 0.0, g^(5)_k⊥(0) ≈ 120. * At δk⊥[kd] = 0.1, g^(5)_k⊥(0) ≈ 8. * At δk⊥[kd] = 0.2, g^(5)_k⊥(0) ≈ 2. * For δk⊥[kd] > 0.3, g^(5)_k⊥(0) ≈ 1. **Plot for n=6:** * Trend: The correlation function decreases rapidly as δk⊥[kd] increases, then plateaus. * Data Points: * At δk⊥[kd] = 0.0, g^(6)_k⊥(0) ≈ 600. * At δk⊥[kd] = 0.1, g^(6)_k⊥(0) ≈ 30. * At δk⊥[kd] = 0.2, g^(6)_k⊥(0) ≈ 5. * For δk⊥[kd] > 0.4, g^(6)_k⊥(0) ≈ 1. ### Key Observations * As 'n' increases, the initial value of the correlation function g^(n)_k⊥(0) at δk⊥[kd] = 0.0 increases significantly. * For n=4, 5, and 6, the correlation function plateaus at a value close to 1 for larger values of δk⊥[kd]. * The theoretical fit (dashed black line) appears to match the data reasonably well in all plots. ### Interpretation The plots illustrate how the correlation function g^(n)_k⊥(0) changes with the momentum difference δk⊥[kd] for different values of 'n'. The rapid decrease in the correlation function at small δk⊥[kd] suggests a strong correlation at short distances, which diminishes as the distance increases. The plateau observed for n=4, 5, and 6 indicates that the correlation function approaches a constant value at larger momentum differences. The increasing initial value of the correlation function with increasing 'n' suggests that the correlation strength is dependent on the parameter 'n'. The theoretical fit provides a model for the relationship between the correlation function and the momentum difference, which can be used to further analyze the underlying physical processes. </details> We define the error bars on g ( n ) (0) as ∆( N Ω ) n / 〈 N Ω 〉 n , where ∆( N Ω ) n is the standard error of the factorial moment of order n , $$\Delta ( N _ { \Omega } ) _ { n } = \frac { 1 } { \sqrt { N _ { r u n s } } } \sqrt { \langle ( N _ { \Omega } ) _ { n } ^ { 2 } \rangle - \langle ( N _ { \Omega } ) _ { n } \rangle ^ { 2 } } \quad ( 1 0 )$$ ## Appendix D. RANDOMIZED DATA AND UPPER LIMIT OF THE ORDER OF CORRELATIONS We use the randomized data to estimate up to which order n we are able to compute the magnitudes of correlation functions correctly. As explained in the main text, the randomized data sets are built from the actual data sets of the experiment by shuffling the atoms from a given shot into many different shots. By construction, the randomized data has the same number of shots and of atoms per shot than the non-randomized data. Therefore we expect it to provide a means to test the accuracy of the finite data sets we use: because the randomized data should exhibit a perfect Poisson statistics, a systematic deviation to the latter is used as a signal of insufficient statistics. As illustrated in Fig. 7 in the case of Bose superfluids at U/J = 5 (same data set as the one shown in Fig. 1 of the main text), we are able to compute correlations of order larger than n = 6. However, for n > 6 we observe a systematic shift to the expected value g ( n ) (0) = 1 in the randomized data. We attribute this systematic shift to the fact that there are too few shots with an atom number N Ω > 6 to correctly evaluate the probability to find N Ω = n for n > 6. We therefore restrict our analysis to the orders n ≤ 6. A similar analysis for the other data sets shown in this work yields the same result. FIG. 7. Magnitude g ( n ) (0) as a function of the order n of correlations. The blue circles correspond to experimental data measured in Bose superfluids at U/J = 5. The orange squares correspond to the randomized data set. <details> <summary>Image 7 Details</summary>  ### Visual Description ## Chart: Correlation Amplitude vs. Order ### Overview The image contains two line plots showing the correlation amplitude g^(n)(0) as a function of the order 'n'. The top plot displays two data series: one for U/J = 5, which shows an increasing trend, and another for randomized data, which remains relatively constant. The bottom plot shows a zoomed-in view of the randomized data. ### Components/Axes **Top Plot:** * **Y-axis:** Correlation amplitude g^(n)(0), ranging from 1.0 to 1.4. * **X-axis:** Order n, ranging from 1 to 10. * **Legend (top-left):** * Black dashed line: f_coh = 0.9960(5) * Blue circles: U/J = 5 * Orange squares: randomized **Bottom Plot:** * **Y-axis:** g^(n)(0), ranging from 0.98 to 1.00. * **X-axis:** Order of correlation n, ranging from 1 to 10. * **Data:** Orange squares, corresponding to the "randomized" data from the top plot. ### Detailed Analysis **Top Plot:** * **f_coh = 0.9960(5) (Black dashed line):** This line represents a theoretical curve. It starts at approximately 1.0 and increases with 'n'. A shaded blue region surrounds this line, indicating uncertainty. * n=1: g^(n)(0) ≈ 1.0 * n=2: g^(n)(0) ≈ 1.0 * n=4: g^(n)(0) ≈ 1.04 * n=6: g^(n)(0) ≈ 1.1 * n=8: g^(n)(0) ≈ 1.22 * n=10: g^(n)(0) ≈ 1.35 * **U/J = 5 (Blue circles):** This data series shows a clear increasing trend. Error bars are present at each data point. * n=1: g^(n)(0) ≈ 1.0, error bar ±0.01 * n=2: g^(n)(0) ≈ 1.02, error bar ±0.01 * n=3: g^(n)(0) ≈ 1.04, error bar ±0.01 * n=4: g^(n)(0) ≈ 1.09, error bar ±0.01 * n=5: g^(n)(0) ≈ 1.12, error bar ±0.01 * n=6: g^(n)(0) ≈ 1.13, error bar ±0.01 * n=7: g^(n)(0) ≈ 1.22, error bar ±0.02 * n=8: g^(n)(0) ≈ 1.23, error bar ±0.02 * n=9: g^(n)(0) ≈ 1.32, error bar ±0.03 * n=10: g^(n)(0) ≈ 1.35, error bar ±0.05 * **Randomized (Orange squares):** This data series remains relatively constant around 1.0. Error bars are present at each data point. * n=1: g^(n)(0) ≈ 1.0, error bar ±0.01 * n=2: g^(n)(0) ≈ 1.0, error bar ±0.01 * n=3: g^(n)(0) ≈ 1.0, error bar ±0.01 * n=4: g^(n)(0) ≈ 1.0, error bar ±0.01 * n=5: g^(n)(0) ≈ 1.0, error bar ±0.01 * n=6: g^(n)(0) ≈ 1.0, error bar ±0.01 * n=7: g^(n)(0) ≈ 1.0, error bar ±0.01 * n=8: g^(n)(0) ≈ 1.0, error bar ±0.01 * n=9: g^(n)(0) ≈ 1.0, error bar ±0.01 * n=10: g^(n)(0) ≈ 1.0, error bar ±0.01 **Bottom Plot:** * This plot provides a zoomed-in view of the "randomized" data, allowing for a closer examination of its fluctuations around 1.0. The error bars are more visible in this plot. ### Key Observations * The correlation amplitude for U/J = 5 increases significantly with the order 'n'. * The correlation amplitude for the randomized data remains relatively constant around 1.0. * The theoretical curve (f_coh = 0.9960(5)) closely follows the trend of the U/J = 5 data. * The error bars for the randomized data in the bottom plot show the degree of fluctuation around 1.0. ### Interpretation The data suggests that the system with U/J = 5 exhibits increasing correlations as the order 'n' increases, indicating a stronger dependence on higher-order correlations. In contrast, the randomized data shows no such trend, suggesting a lack of significant higher-order correlations. The theoretical curve provides a model that aligns well with the U/J = 5 data, supporting the theoretical framework. The randomized data serves as a baseline, demonstrating the absence of correlation effects in a disordered system. The error bars provide a measure of the uncertainty in the measurements. </details> ## Appendix E. MODEL FOR THE MAGNITUDES OF CORRELATIONS IN BOSE SUPERFLUIDS We introduce a simple model to account for the presence of atoms belonging both to the BEC and to its depletion in the central voxel located at k = 0 . The depletion of the condensate results from the sum of the quantum depletion (induced by interactions) and the thermal depletion (induce by the finite temperature). As we have no mean to distinguish between the quantum and the thermal depletion, the depletion of the condensate we refer to in the following is the sum of both contributions. We assume that the statistical properties of the BEC mode are those of a coherent state, g ( n ) BEC (0) = 1. In contrast, the statistical properties of the depletion are assumed to be those of a thermal state, g ( n ) dep (0) = n !. This latter assumption holds because we consider a small voxel of size dk = 0 . 025 k d l c 0 . 15 k d where the correlation length l c is set by the inverse in-trap size 1 /L of the gas (see the main text). Moreover, we further assume that the BEC and depletion operators are uncorrelated. This is valid because the voxel we consider is much smaller than the volume occupied by the depletion, implying that the correlation between the total number of BEC atoms and the total number of depleted atoms (in the canonical ensemble) can be safely neglected. We are interested in the statistical properties of the atom number ˆ N Ω = a † Ω a Ω falling in the considered voxel, where a Ω = a BEC + a dep is the sum of the operators associated with the BEC and the depletion. Since we assume the operators a BEC and a dep are uncorrelated, one simply obtains $$\langle ( a _ { \Omega } ^ { \dagger } ) ^ { n } a _ { \Omega } ^ { n } \rangle = \sum _ { p = 1 } ^ { n } { n \choose p } ^ { 2 } \langle ( a _ { B E C } ^ { \dagger } ) ^ { p } a _ { B E C } ^ { p } ( a _ { d e p } ^ { \dagger } ) ^ { n - p } a _ { d e p } ^ { n - p } \rangle$$ where by definition a n Ω = ∑ n p =1 ( n p ) a p BEC a n -p dep . This leads to the formula given in the main text, $$g ^ { ( n ) } ( 0 ) - 1 & = \frac { \langle ( a _ { \Omega } ^ { \dagger } ) ^ { n } a _ { \Omega } ^ { n } \rangle } { \langle a _ { \Omega } ^ { \dagger } a _ { \Omega } \rangle ^ { n } } - 1 \\ & = \sum _ { p = 1 } ^ { n - 1 } \left [ ( n - p ) ! \binom { n } { p } ^ { 2 } - \binom { n } { p } \right ] f _ { c o h } ^ { p } ( 1 - f _ { c o h } ) ^ { n - p }$$ ## Appendix F. RELATION f coh VERSUS f c . The coherent fraction is approximatively given by $$f _ { c o h } \simeq \frac { \rho _ { B E C } ( 0 ) } { \rho _ { B E C } ( 0 ) + \rho _ { d e p } ( 0 ) } \quad \ \ ( 1 2 )$$ where ρ BEC (0) ( resp. ρ dep (0)) is the momentum density of the BEC ( resp. of the depletion) at the center of the Brillouin zone, k = 0. The condensate fraction is approximately given by $$f _ { c } \simeq \frac { \rho _ { B E C } ( 0 ) \mathcal { V } _ { c } } { \rho _ { B E C } ( 0 ) \mathcal { V } _ { c } + \rho _ { d e p } ( 0 ) \mathcal { V } _ { d } } \quad ( 1 3 )$$ where V c ( resp. V d ) is the total volume of the momentum-space occupied by the BEC ( resp. by the depletion) normalized to the density at k = 0. Combining these two definitions, one obtains the equation Eq. (4) of the main text. To compare the prediction of Eq. (4) with the measurements, we need to evaluate the volumes V c and V d . We describe the BEC density profile with an isotropic function [46] $$n _ { B E C } ( \mathbf k ) = \rho _ { B E C } ( 0 ) e ^ { - k ^ { 2 } / 2 \sigma _ { b e c } ^ { 2 } } \quad ( 1 4 ) \quad \begin{array} { l l } { { d e n } } \\ { { v / } } \end{array}$$ and the density profile of the depletion with a lorentzian function along each (uncoupled) lattice axis, $$n _ { d e p } ( k ) = \rho _ { d e p } ( 0 ) \, \Pi _ { j = x , y , z } \left ( \frac { \sigma _ { d e p } ^ { 2 } / 4 } { k _ { j } ^ { 2 } + ( \sigma _ { d e p } / 2 ) ^ { 2 } } \right ) . \quad ( 1 5 )$$ We obtain V c = ( √ 2 πσ bec ) 3 and V d = ( σ dep arctan [ k d σ dep ]) 3 . From fitting the measured density profiles with the above dependency, we obtain V c / V d = 0 . 03(1), a value which we use to plot the theoretical prediction in Fig. 4.