2207.14070v2

## 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 ## Diagram: Quantum Correlation Measurement Setup ### Overview The image depicts a quantum correlation measurement system with three components: (a) a schematic of the experimental setup, (b) a log-scale graph of correlation strength vs. order of correlation, and (c) a comparison of correlation strength between different states. ### Components/Axes #### (a) Experimental Setup Diagram - **Key Elements**: - Cone-shaped beam labeled with length **L = 43 cm** - Beam splitter (green horizontal line at top) - Single-atom detector (gray circle at base) - Momentum space inset showing: - Axes: **k_x**, **k_y**, **k_z** - Red box labeled **δk** (momentum uncertainty) - Green line labeled **1/L** (spatial confinement) #### (b) Log-Scale Graph - **Axes**: - **x-axis**: Order of correlation **n** (1–6, integer values) - **y-axis**: Correlation strength **g^(n)(0)** (log scale, 10⁰–10³) - **Legend**: - **n!** (black solid line) - **Mott insulator** (black squares) - **Superfluid** (blue circles) #### (c) Randomized Data Comparison - **Axes**: - Same as (b): **n** (1–6) and **g^(n)(0)** - **Legend**: - **Randomized** (orange squares) ### Detailed Analysis #### (b) Log-Scale Graph Trends - **n! (black line)**: Exponential growth (e.g., 1.0 → 720 for n=6) - **Mott insulator (black squares)**: - Follows **n!** trend closely - Data points: - n=1: ~1.0 - n=2: ~2.0 - n=3: ~6.0 - n=4: ~24.0 - n=5: ~120.0 - n=6: ~720.0 - Error bars: ±5–10% of measured values - **Superfluid (blue circles)**: - Flat line at **g^(n)(0) ≈ 1.0** for all n - Error bars: ±0.1–0.2 #### (c) Randomized Data Comparison - **Randomized (orange squares)**: - Slightly above 1.0 for all n - Values: - n=1: ~1.05 - n=2: ~1.08 - n=3: ~1.10 - n=4: ~1.12 - n=5: ~1.15 - n=6: ~1.20 - **Superfluid (blue circles)**: - Consistent with (b): **g^(n)(0) ≈ 1.0** - **Mott insulator (black squares)**: - Consistent with (b): Exponential growth ### Key Observations 1. **Mott insulator** exhibits strong quantum correlations, with **g^(n)(0)** matching **n!** exactly. 2. **Superfluid** shows no correlation decay, maintaining **g^(n)(0) ≈ 1.0**. 3. **Randomized data** demonstrates weak, linear growth in correlation strength, suggesting residual correlations from classical noise. ### Interpretation The data suggests: - **Mott insulator** systems preserve quantum correlations across all orders (n), behaving like an ideal quantum state. - **Superfluid** lacks long-range correlations, consistent with its thermodynamic equilibrium state. - **Randomized data** indicates that classical noise introduces weak, artificial correlations, highlighting the importance of quantum state preparation fidelity. The experimental setup (a) uses a 43 cm beam path to spatially confine atoms, enabling momentum-resolved correlation measurements. The **δk** parameter in the momentum space inset likely represents the momentum resolution critical for distinguishing correlation orders. </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 ## Line Graphs: Probability Distributions of Detected Atom Numbers ### Overview The image contains two line graphs (a) and (b) comparing probability distributions of detected atom numbers (N_Ω) under different theoretical models (Thermal, Poisson) and experimental modes (Mott insulator, BEC mode). Both graphs use logarithmic and linear scales for probability (P(N_Ω)) and linear scales for N_Ω, respectively. --- ### Components/Axes #### Graph (a) - **Y-axis**: Probability P(N_Ω) (logarithmic scale: 10⁻³ to 10⁰) - **X-axis**: Detected atom number N_Ω (linear scale: 0 to 6) - **Legend**: - Dashed black line: Thermal distribution - Dashed blue line: Poisson distribution - Gray squares with error bars: Mott insulator experimental data - **Legend Position**: Top-right corner #### Graph (b) - **Y-axis**: Probability P(N_Ω) (linear scale: 0 to 0.15) - **X-axis**: Detected atom number N_Ω (linear scale: 0 to 15) - **Legend**: - Dashed black line: Thermal distribution - Dashed blue line: Poisson distribution - Blue circles with error bars: BEC mode experimental data - **Legend Position**: Top-right corner --- ### Detailed Analysis #### Graph (a) - **Thermal Distribution**: - Dashed black line decreases monotonically from ~10⁰ at N_Ω=0 to ~10⁻³ at N_Ω=6. - No shaded region. - **Poisson Distribution**: - Dashed blue line decreases similarly to Thermal but with a shaded blue region (likely confidence interval) around it. - Shaded region widens at lower N_Ω (0–2), indicating higher uncertainty. - **Mott Insulator Data**: - Gray squares with error bars align closely with the Poisson line. - Error bars increase in size at lower N_Ω (e.g., ±0.1 at N_Ω=0, ±0.05 at N_Ω=6). #### Graph (b) - **Thermal Distribution**: - Dashed black line decreases monotonically from ~0.15 at N_Ω=0 to ~0.005 at N_Ω=15. - **Poisson Distribution**: - Dashed blue line peaks at N_Ω=5 (~0.12) and decreases afterward. - No shaded region. - **BEC Mode Data**: - Blue circles with error bars align with the Poisson peak at N_Ω=5. - Error bars are smaller at higher N_Ω (e.g., ±0.01 at N_Ω=10, ±0.02 at N_Ω=15). --- ### Key Observations 1. **Graph (a)**: - Mott insulator data closely follows the Poisson distribution, suggesting Poissonian statistics dominate at low N_Ω. - Thermal distribution deviates slightly at higher N_Ω but remains close to Poisson. - Shaded region around Poisson in (a) implies greater uncertainty at low N_Ω. 2. **Graph (b)**: - BEC mode data peaks at N_Ω=5, matching the Poisson distribution’s maximum. - Thermal distribution declines steadily without a peak, contrasting with Poisson/BEC mode. - Error bars in BEC mode are smaller at higher N_Ω, indicating more precise measurements. --- ### Interpretation - **Poisson vs. Thermal**: The Poisson distribution (dashed blue) consistently aligns with experimental data (Mott insulator in (a), BEC mode in (b)), suggesting quantum systems exhibit Poissonian statistics, likely due to coherent state behavior or shot noise. - **Thermal Distribution**: The dashed black line represents a classical thermal model, which diverges from experimental results, indicating quantum effects dominate over thermal noise. - **Shaded Regions**: In (a), the shaded area around Poisson highlights uncertainty in low-N_Ω measurements, possibly due to experimental limitations or quantum fluctuations. - **BEC Mode Peak**: The sharp peak in BEC mode data at N_Ω=5 (graph b) mirrors the Poisson distribution’s maximum, reinforcing the role of quantum coherence in atom number statistics. --- ### Notable Trends - **Low-N_Ω Uncertainty**: Larger error bars at N_Ω=0–2 in graph (a) suggest challenges in detecting small atom numbers. - **High-N_Ω Precision**: Smaller error bars at N_Ω>10 in graph (b) indicate improved measurement accuracy for larger atom numbers. - **Model Agreement**: Experimental data (Mott/BEC) closely follow Poisson predictions, validating the theoretical model for these systems. --- ### Conclusion The graphs demonstrate that quantum systems (Mott insulator, BEC mode) exhibit Poissonian statistics for detected atom numbers, contrasting with classical thermal models. The shaded regions and error bars emphasize experimental uncertainties, particularly at low N_Ω. These results highlight the dominance of quantum coherence over thermal noise in these systems. </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 ## 2x2 Grid of Plots: Coherence and Density Analysis ### Overview The image contains four plots arranged in a 2x2 grid, analyzing coherence (f_coh) and density distributions (ρ(k)) for two values of U/J (5 and 20). Each plot includes error bars, shaded regions, and dashed reference lines. --- ### Components/Axes #### Plot (a) & (b): Coherence vs. Order - **X-axis**: "Order n" (integer values 1–6) - **Y-axis**: "g^(n)(0)" (unitless, range 1.00–2.25) - **Legend**: - Dashed line: f_coh = 0.9960(5) (U/J = 5) or 0.971(1) (U/J = 20) - Markers: Blue circles (U/J = 5) or black circles (U/J = 20) - **Shaded region**: ±σ uncertainty around the dashed line #### Plot (c) & (d): Density Distributions - **X-axis**: "k/k_d" (unitless, range -0.4 to 0.4) - **Y-axis**: "Density ρ(k)" (a.u., log scale 10⁻⁴ to 10⁰) - **Legend**: - Solid line: Density distribution for U/J = 5 (blue) or 20 (gray) - Dashed line: 1 - f_coh = 0.004 (U/J = 5) or 0.029 (U/J = 20) - **Shaded region**: ±σ uncertainty around the solid line --- ### Detailed Analysis #### Plot (a): U/J = 5 - **Trend**: g^(n)(0) increases monotonically with n (slope ≈ 0.02 per order). - **Data points**: - n=1: 1.00 ± 0.01 - n=2: 1.02 ± 0.01 - n=3: 1.04 ± 0.01 - n=4: 1.06 ± 0.01 - n=5: 1.09 ± 0.01 - n=6: 1.12 ± 0.01 - **Dashed line**: f_coh = 0.9960 ± 0.0005 (high coherence) #### Plot (b): U/J = 20 - **Trend**: g^(n)(0) increases more steeply (slope ≈ 0.03 per order). - **Data points**: - n=1: 1.00 ± 0.01 - n=2: 1.03 ± 0.01 - n=3: 1.07 ± 0.01 - n=4: 1.12 ± 0.01 - n=5: 1.17 ± 0.01 - n=6: 1.22 ± 0.01 - **Dashed line**: f_coh = 0.971 ± 0.001 (reduced coherence) #### Plot (c): U/J = 5 - **Density peak**: Sharp peak at k/k_d = 0 (height ≈ 10⁻¹). - **Dashed line**: 1 - f_coh = 0.004 (low disorder). - **Shaded region**: ±σ uncertainty spans ±0.05 in k/k_d. #### Plot (d): U/J = 20 - **Density peak**: Broader peak at k/k_d = 0 (height ≈ 10⁻²). - **Dashed line**: 1 - f_coh = 0.029 (higher disorder). - **Shaded region**: ±σ uncertainty spans ±0.1 in k/k_d. --- ### Key Observations 1. **Coherence decay**: Higher U/J (20) reduces f_coh by ~2.5% compared to U/J = 5. 2. **Density broadening**: U/J = 20 shows a wider density distribution (peak height reduced by ~10× and width increased by ~2×). 3. **Error margins**: Uncertainty in g^(n)(0) is consistent (~0.01), but density distributions show larger relative uncertainty for U/J = 20. --- ### Interpretation - **Coherence vs. Disorder**: The inverse relationship between f_coh and U/J suggests that stronger interactions (higher U/J) disrupt coherence, increasing disorder (1 - f_coh). - **Density distributions**: The broadening of ρ(k) for U/J = 20 indicates enhanced fluctuations or thermal disorder, consistent with reduced coherence. - **Error analysis**: The shaded regions and error bars highlight measurement precision limits, critical for validating the trends. The dashed lines (1 - f_coh) directly quantify disorder, linking coherence to system stability. This analysis demonstrates how interaction strength (U/J) governs coherence and disorder, with implications for understanding phase transitions or critical phenomena in the system. </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 ## Line Chart: Relationship Between f_c and f_coh ### Overview The image depicts a line chart with a shaded confidence interval and experimental data points. The chart explores the relationship between two variables: f_c (x-axis) and f_coh (y-axis). A dashed theoretical model curve is plotted alongside experimental data points with error bars, and an inset provides a zoomed-in view of the high-f_c region. ### Components/Axes - **X-axis (f_c)**: Ranges from 0.2 to 1.0 in increments of 0.2. Labeled "f_c" with gridlines. - **Y-axis (f_coh)**: Ranges from 0.85 to 1.00 in increments of 0.05. Labeled "f_coh" with gridlines. - **Legend**: - Dashed orange line: "Theoretical model" - Shaded orange region: "Confidence interval (±1σ)" - Orange data points with error bars: "Experimental measurements" - **Inset**: A smaller chart inset in the bottom-right corner, focusing on f_c values from 0.6 to 1.0. ### Detailed Analysis #### Main Chart Data Points | f_c | f_coh | Error Bars (approx.) | |-----|-------|----------------------| | 0.2 | 0.85 | ±0.02 | | 0.4 | 0.92 | ±0.03 | | 0.6 | 0.95 | ±0.02 | | 0.8 | 0.98 | ±0.01 | | 1.0 | 1.00 | ±0.01 | #### Inset Chart Data Points | f_c | f_coh | Error Bars (approx.) | |-----|-------|----------------------| | 0.6 | 0.95 | ±0.02 | | 0.8 | 0.98 | ±0.01 | | 1.0 | 1.00 | ±0.01 | #### Theoretical Model Curve - The dashed orange line follows a sigmoidal trend, starting at (0.2, 0.85) and asymptotically approaching 1.00 as f_c increases. - The curve's slope steepens between f_c = 0.4 and 0.8, then plateaus near f_c = 1.0. ### Key Observations 1. **Positive Correlation**: f_coh increases monotonically with f_c, suggesting a direct relationship. 2. **Model Fit**: Experimental data points align closely with the theoretical model, with error bars within the shaded confidence interval. 3. **Inset Insight**: The zoomed-in view reveals a steeper slope in the high-f_c region (0.6–1.0), indicating a nonlinear relationship. 4. **Uncertainty**: Error bars decrease as f_c increases, implying higher precision in measurements at higher f_c values. ### Interpretation The chart demonstrates that f_coh is strongly dependent on f_c, with the theoretical model accurately predicting experimental outcomes. The confidence interval (±1σ) suggests the model's predictions are reliable within this uncertainty range. The inset highlights that the relationship becomes more sensitive to changes in f_c at higher values, which could imply a threshold effect or saturation behavior near f_c = 1.0. The decreasing error bars at higher f_c values may reflect improved measurement techniques or reduced variability in the system under study. This data could be critical for optimizing systems where f_c and f_coh are interdependent parameters, such as in signal processing or material science applications. </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 ## Log-Log Plot: Probability P(N_Ω) vs. N_Ω for Different V_Ω Values ### Overview The image is a log-log plot showing the relationship between probability \( P(N_\Omega) \) (y-axis) and \( N_\Omega \) (x-axis) for four distinct values of \( V_\Omega \). The data is represented by four distinct symbols and colors, with a dashed reference line. The plot demonstrates how \( P(N_\Omega) \) decreases as \( N_\Omega \) increases, with varying trends depending on \( V_\Omega \). --- ### Components/Axes - **X-axis (N_Ω)**: Labeled \( N_\Omega \), ranging from 0 to 8 in integer increments. Logarithmic scale. - **Y-axis (P(N_Ω))**: Labeled \( P(N_\Omega) \), ranging from \( 10^{-3} \) to \( 10^0 \). Logarithmic scale. - **Legend**: Located in the top-left corner, mapping symbols/colors to \( V_\Omega \) values: - **Black circles**: \( V_\Omega = 0.9V_c \) - **Red squares**: \( V_\Omega = 2.1V_c \) - **Orange triangles**: \( V_\Omega = 4.2V_c \) - **Blue inverted triangles**: \( V_\Omega = 7.2V_c \) - **Dashed reference line**: Diagonal line from top-left to bottom-right, likely representing a theoretical prediction or baseline. --- ### Detailed Analysis 1. **Black Circles (\( V_\Omega = 0.9V_c \))**: - Starts at \( P(N_\Omega) \approx 10^{-1} \) when \( N_\Omega = 0 \). - Declines steeply, crossing the red squares (\( V_\Omega = 2.1V_c \)) near \( N_\Omega = 2 \). - Ends at \( P(N_\Omega) \approx 10^{-3} \) for \( N_\Omega = 8 \). 2. **Red Squares (\( V_\Omega = 2.1V_c \))**: - Begins at \( P(N_\Omega) \approx 10^{-1} \) for \( N_\Omega = 0 \). - Declines gradually, intersecting the black circles at \( N_\Omega \approx 2 \). - Reaches \( P(N_\Omega) \approx 10^{-3} \) at \( N_\Omega = 8 \). 3. **Orange Triangles (\( V_\Omega = 4.2V_c \))**: - Starts at \( P(N_\Omega) \approx 10^{-1} \) for \( N_\Omega = 0 \). - Declines more gradually than red squares, ending at \( P(N_\Omega) \approx 10^{-3} \) for \( N_\Omega = 8 \). 4. **Blue Inverted Triangles (\( V_\Omega = 7.2V_c \))**: - Begins at \( P(N_\Omega) \approx 10^{-1} \) for \( N_\Omega = 0 \). - Declines the slowest among all series, ending at \( P(N_\Omega) \approx 10^{-3} \) for \( N_\Omega = 8 \). 5. **Dashed Reference Line**: - Follows a steep negative slope, suggesting a theoretical \( P(N_\Omega) \propto N_\Omega^{-k} \) relationship. - All data series lie above this line, indicating deviations from the predicted trend. --- ### Key Observations - **Universal Decline**: All \( V_\Omega \) values show a decreasing trend of \( P(N_\Omega) \) with increasing \( N_\Omega \). - **Crossover Point**: The black circles (\( V_\Omega = 0.9V_c \)) and red squares (\( V_\Omega = 2.1V_c \)) intersect near \( N_\Omega = 2 \), suggesting a threshold effect. - **Divergence from Theory**: The dashed line’s steep slope contrasts with the gradual decline of the data series, implying the model underestimates \( P(N_\Omega) \) at higher \( N_\Omega \). - **V_Ω Dependency**: Higher \( V_\Omega \) values (e.g., \( 7.2V_c \)) exhibit slower decay in \( P(N_\Omega) \), indicating stronger resilience or stability. --- ### Interpretation The plot demonstrates that \( P(N_\Omega) \) decreases with increasing \( N_\Omega \), but the rate of decay depends critically on \( V_\Omega \). Higher \( V_\Omega \) values (e.g., \( 7.2V_c \)) maintain higher probabilities at larger \( N_\Omega \), suggesting they represent more robust or stable configurations. The crossover between \( V_\Omega = 0.9V_c \) and \( 2.1V_c \) implies a critical threshold where system behavior shifts. The dashed reference line’s divergence highlights a discrepancy between the observed data and the theoretical model, potentially pointing to unaccounted variables (e.g., noise, external perturbations) or nonlinear effects in the system. This analysis could inform optimization strategies for systems governed by \( V_\Omega \)-dependent probabilities. </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 ## Line Graphs: Relationship Between δk⊥ and g^(n)_δk⊥(0) for n=2 to 6 ### Overview The image contains five line graphs arranged in a 2x2 grid with one additional graph at the bottom right. Each graph represents a different value of the parameter `n` (2–6) and plots the relationship between the perpendicular wavevector deviation `δk⊥` (normalized by `k_d`) and the function `g^(n)_δk⊥(0)` on a logarithmic scale. Data points are represented by orange circles with error bars, and a dashed black reference line is present in all graphs. --- ### Components/Axes 1. **X-Axis**: - Label: `δk⊥ [k_d]` - Range: 0.0 to 0.6 (linear scale) - Position: Bottom of all graphs 2. **Y-Axis**: - Label: `g^(n)_δk⊥(0)` - Scales: - `n=2`: Linear (1.0–2.0) - `n=3`: Logarithmic (10⁰–10¹) - `n=4`: Logarithmic (10⁰–10¹) - `n=5`: Logarithmic (10⁰–10²) - `n=6`: Logarithmic (10⁰–10³) - Position: Left of all graphs 3. **Legend**: - Located in the top-right corner of each graph. - Format: Orange circle with `n=` value (e.g., `n=2`, `n=3`, etc.). 4. **Data Points**: - Orange circles with vertical error bars. - Error bars are small relative to the data points. 5. **Reference Line**: - Dashed black line in all graphs. - Position: Overlays data points, suggesting a theoretical or fitted trend. --- ### Detailed Analysis #### Graph 1: `n=2` - **Y-Axis Scale**: Linear (1.0–2.0). - **Trend**: Data points decrease smoothly from ~2.0 at `δk⊥=0.0` to ~1.0 at `δk⊥=0.3`. - **Error Bars**: Minimal (~±0.05). #### Graph 2: `n=3` - **Y-Axis Scale**: Logarithmic (10⁰–10¹). - **Trend**: Data points drop from ~10¹ at `δk⊥=0.0` to ~10⁰ at `δk⊥=0.6`. - **Error Bars**: ~±0.1 (log scale). #### Graph 3: `n=4` - **Y-Axis Scale**: Logarithmic (10⁰–10¹). - **Trend**: Similar to `n=3`, but steeper decline. - **Error Bars**: ~±0.1 (log scale). #### Graph 4: `n=5` - **Y-Axis Scale**: Logarithmic (10⁰–10²). - **Trend**: Data points start at ~10¹ at `δk⊥=0.0` and drop to ~10⁰ at `δk⊥=0.6`. - **Error Bars**: ~±0.1 (log scale). #### Graph 5: `n=6` - **Y-Axis Scale**: Logarithmic (10⁰–10³). - **Trend**: Sharp initial drop from ~10² at `δk⊥=0.0` to ~10⁰ at `δk⊥=0.4`, then plateaus. - **Error Bars**: ~±0.1 (log scale). --- ### Key Observations 1. **Universal Decline**: All graphs show a decreasing trend of `g^(n)_δk⊥(0)` with increasing `δk⊥`. 2. **Scale Dependence**: Higher `n` values require logarithmic scales to capture the full range of `g^(n)_δk⊥(0)`. 3. **Convergence**: For `n≥4`, the data points align closely with the dashed reference line, suggesting a consistent theoretical model. 4. **Error Consistency**: Error bars are smallest at `δk⊥=0.0` and grow slightly with increasing `δk⊥`. --- ### Interpretation The graphs demonstrate that `g^(n)_δk⊥(0)` decreases monotonically with increasing perpendicular wavevector deviation `δk⊥`. The use of logarithmic scales for `n≥3` indicates that the magnitude of `g^(n)_δk⊥(0)` spans several orders of magnitude, particularly for `n=6`. The dashed reference line likely represents a theoretical prediction (e.g., exponential decay or power-law behavior), which the data closely follows. The error bars suggest high precision in measurements, especially at low `δk⊥`. The parameter `n` may correspond to a system property (e.g., dimensionality, coupling strength, or disorder level). The steeper decline for higher `n` implies that larger `n` values amplify the sensitivity of `g^(n)_δk⊥(0)` to `δk⊥`. This could reflect phenomena such as increased dissipation, localization, or interference effects in a physical system. The spatial arrangement of graphs emphasizes the systematic variation of `n`, allowing direct comparison of trends across different parameter regimes. The absence of outliers or anomalies suggests robustness in the observed behavior. </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 ## Line Chart with Scatter Plots: Correlation Amplitude vs. Order n ### Overview The image contains two vertically stacked plots. The top plot shows a line graph with two scatter series and a dashed reference line, while the bottom plot displays a single scatter series with error bars. Both plots share the same x-axis ("Order n") but have distinct y-axis labels. The data suggests a comparison between theoretical predictions (dashed line), numerical simulations (blue circles), and randomized data (orange squares). ### Components/Axes **Top Plot:** - **X-axis (Order n):** Integer values from 2 to 10, labeled "Order n". - **Y-axis (Correlation amplitude g^(n)(0)):** Range from 0.98 to 1.4, labeled "Correlation amplitude g^(n)(0)". - **Legend (top-left):** - Dashed black line: `f_coh = 0.9960(5)` (theoretical prediction). - Blue circles: `U/J = 5` (numerical simulation). - Orange squares: "randomized" (control data). **Bottom Plot:** - **X-axis (Order n):** Same as top plot (2–10). - **Y-axis (g^(n)(0)):** Range from 0.98 to 1.02, labeled "g^(n)(0)". - **Legend (top-left):** Same as top plot, but only orange squares ("randomized") are plotted. ### Detailed Analysis **Top Plot Trends:** 1. **Dashed Line (`f_coh`):** - Starts at ~1.0 for n=2, increases steadily to ~1.4 at n=10. - Uncertainty band (shaded gray) widens with increasing n, suggesting growing error margins. 2. **Blue Circles (`U/J = 5`):** - Follows the dashed line trend but with slight scatter (e.g., n=4: 1.05, n=6: 1.11). - Error bars (not explicitly shown but implied by shaded area) align with the dashed line's uncertainty. 3. **Orange Squares (`randomized`):** - Flat at ~1.0 for all n (2–10), with minimal variation (±0.01). **Bottom Plot Trends:** - **Orange Squares (`randomized`):** - All data points cluster tightly around 1.0 (e.g., n=2: 1.00, n=10: 0.99). - Error bars vary slightly (e.g., ±0.01 at n=2, ±0.02 at n=10). ### Key Observations 1. **Theoretical vs. Simulated Data:** - The `f_coh` dashed line and `U/J = 5` blue circles exhibit nearly identical upward trends, confirming strong agreement between theory and simulation. 2. **Randomized Data Consistency:** - The "randomized" orange squares remain flat at 1.0 in both plots, indicating no correlation amplitude dependence on order n under randomization. 3. **Uncertainty Growth:** - The widening shaded band around `f_coh` suggests increasing uncertainty in theoretical predictions at higher orders n. ### Interpretation The data demonstrates that: - **Higher-order correlations (n)** amplify the correlation amplitude for both the theoretical (`f_coh`) and simulated (`U/J = 5`) cases, with the effect becoming less certain at larger n. - **Randomization** (orange squares) suppresses order-dependent correlations, maintaining a baseline amplitude of 1.0. This implies that the observed trends in `f_coh` and `U/J = 5` are not artifacts of uncontrolled variables. - The **uncertainty in `f_coh`** grows with n, highlighting potential limitations in the theoretical model or numerical precision at higher orders. The results suggest that the correlation amplitude is sensitive to system parameters (e.g., `U/J = 5`) but robust to randomization, providing insights into the role of order in correlation dynamics. </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.