#  Suppression of Bogoliubov momentum pairing and emergence of non-Gaussian
correlations in ultracold interacting Bose gases

Jan-Philipp Bureik  These authors contributed equally. Université Paris-
Saclay, Institut d’Optique Graduate School, CNRS, Laboratoire Charles Fabry,
91127, Palaiseau, France  Gaétan Hercé  These authors contributed equally.
Université Paris-Saclay, Institut d’Optique Graduate School, CNRS, Laboratoire
Charles Fabry, 91127, Palaiseau, France  Maxime Allemand  Université Paris-
Saclay, Institut d’Optique Graduate School, CNRS, Laboratoire Charles Fabry,
91127, Palaiseau, France  Antoine Tenart  Université Paris-Saclay, Institut
d’Optique Graduate School, CNRS, Laboratoire Charles Fabry, 91127, Palaiseau,
France  Tommaso Roscilde  Univ Lyon, Ens de Lyon, CNRS, Laboratoire de
Physique, F-69342 Lyon, France  David Clément  Université Paris-Saclay,
Institut d’Optique Graduate School, CNRS, Laboratoire Charles Fabry, 91127,
Palaiseau, France

###### Abstract

Strongly correlated quantum matter Sachdev (2023) – such as interacting
electron systems Morosan  _et al._ (2012) or interacting quantum fluids Adams
_et al._ (2012); Bloch  _et al._ (2022); Chevy and Salomon (2016) – possesses
properties that cannot be understood in terms of linear fluctuations and free
quasi-particles Georges  _et al._ (1996); Metzner  _et al._ (2012). Quantum
fluctuations in these systems are indeed large and generically exhibit non-
Gaussian statistics Georges  _et al._ (2013); Chang  _et al._ (2014) – a
property captured only by inspecting high-order correlations, whose
quantitative reconstruction poses a formidable challenge to both experiments
and theory alike. A prime example of correlated quantum matter is the strongly
interacting Bose fluid, realized by superfluid Helium Pitaevskii and Stringari
(2016); Leggett (2006) and, more recently, ultra-cold atoms Makotyn  _et al._
(2014); Fletcher  _et al._ (2017); Eismann  _et al._ (2016); Klauss  _et al._
(2017); Eigen  _et al._ (2018); Yan  _et al._ (2020). Here, we experimentally
study interacting Bose gases from the weakly to the strongly interacting
regime through single-atom-resolved correlations in momentum space. We observe
that the Bogoliubov pairing among modes of opposite momenta, emblematic of the
weakly interacting regime Bogoliubov (1947); Tenart  _et al._ (2021), is
suppressed as interactions become stronger. This departure from the
predictions of Bogoliubov theory signals the onset of the strongly correlated
regime, as confirmed by numerical simulations that highlight the role of non-
linear quantum fluctuations in our system. Additionally, our measurements
unveil a non-zero four-operator cumulant at even stronger interactions, which
is a direct signature of non-Gaussian correlations. These results shed light
on the emergence and physical origin of non-Gaussian correlations in ensembles
of interacting bosons.

Quantum fluctuations induced by interactions play a prominent role in the
equilibrium properties Sachdev (2023); Paschen and Si (2021) and/or in the
dynamics Polkovnikov  _et al._ (2011) of ensembles of interacting quantum
particles. When fluctuations are weak, they can be captured by a linear theory
which describes the quantum state of the system as a Gaussian state Wang  _et
al._ (2007). In Gaussian states any correlation can be expressed as a
combination of the first two cumulants (mean and variance) of fluctuations. In
contrast, strongly correlated quantum systems exhibit substantial quantum
fluctuations that linear theories fail to account for Georges  _et al._
(1996); Metzner  _et al._ (2012); Chevy and Salomon (2016). The unique nature
of these systems reveals itself in (some) observables exhibiting non-Gaussian
statistics, which are essential to understand the microscopic behavior Armijo
_et al._ (2010); Schweigler  _et al._ (2017); Rispoli  _et al._ (2019) but are
difficult to identify in both experiments and theories.

Bosons with pairwise contact interaction can realize a paradigmatic example of
strongly correlated quantum systems, namely the strongly interacting Bose gas.
In the regime of weak interactions, Bogoliubov theory Bogoliubov (1947) –
leading to the celebrated Lee-Huang-Yang correction to the mean-field equation
of state Lee  _et al._ (1957); Navon  _et al._ (2011); Jørgensen  _et al._
(2018) – treats interaction-induced quantum fluctuations Pitaevskii and
Stringari (1991) as linear fluctuations of a nearly perfect Bose-Einstein
condensate (BEC). The fraction of atoms expelled from the BEC – the quantum
depletion Lopes  _et al._ (2017) – is described in terms of two-mode-squeezed
states for pairs of modes with opposite momenta – a relevant example of a
Gaussian bosonic state. This pairwise population of modes at opposite momenta
(hereafter referred to as Bogoliubov momentum pairing) is the emblematic form
of correlation of the weakly interacting regime, as recently confirmed
experimentally Tenart  _et al._ (2021).

When interactions become strong enough as to push the system away from the
weakly interacting regime, the reconstruction of correlations going beyond
those predicted by Bogoliubov theory becomes a challenge, for both experiments
and theory. On the theory side, the main theoretical tool to treat strong
correlations is offered by numerical calculations Ceperley (1995); Svistunov
_et al._ (2015) – which nonetheless rarely address the correlations probing
non-Gaussian statistics, because of their significant numerical cost. On the
experimental front, ultracold atoms offer the possibility to control the
strength of interactions, giving access to the rich showcase of quantum many-
body phenomena exhibited by strongly interacting bosons Fisher  _et al._
(1989); Bloch  _et al._ (2008); Georgescu  _et al._ (2014). Resonant two-body
interactions tuned with Feshbach resonances realize unitary Bose gases Chevy
and Salomon (2016) hosting few-body Efimov physics Makotyn  _et al._ (2014);
Fletcher  _et al._ (2017); Eismann  _et al._ (2016); Klauss  _et al._ (2017);
Eigen  _et al._ (2018); Yan  _et al._ (2020), which are nonetheless plagued by
large three-body losses Chevy and Salomon (2016) and observing their
equilibrium behavior is challenging. But strongly-interacting Bose gases are
also realized in optical lattices, where three-body losses are negligible,
offering a clear path to study equilibrium many-body physics. Moreover the use
of detection methods with single-atom resolution, as developed in our
experiment using metastable helium-4 (4He∗) atoms Cayla  _et al._ (2018), is
instrumental in reconstructing correlations to all orders Hercé  _et al._
(2023).

Figure 1: Two-body scattering processes from a Bose-Einstein condensate and
momentum correlations. a In Bogoliubov theory, the term $V^{(2)}$ in Eq.(1)
generates two-mode squeezed states at opposite momenta (yellow shaded area). b
Beyond the Bogoliubov approximation, the term $V^{(3)}$ in Eq.(2) involves the
interaction of one atom in the BEC and one atom $\bm{k}\neq\bm{0}$. c The
beyond-Bogoliubov term $V^{(4)}$ in Eq.(3) describes the interaction of two
atoms with $\bm{k},\bm{k}^{\prime}\neq\bm{0}$. At largest order, the terms
$V^{(3)}$ and $V^{(4)}$ belong to series of scattering from the BEC, a
framework from which triplets, quadruplets and $n$-uplets with
${\bm{K}}=\sum_{j=1}^{n}\bm{k}_{j}={\bm{0}}$ emerge (yellow shaded areas). d
Single-shot experimental distribution of individual atoms (blue dots) in
momentum space, from which momentum correlations are computed.

As discussed above, the quantitative description of many-body correlations in
strongly interacting Bose gases is generally a difficult problem. However, the
mechanism by which interactions induce correlations in momentum space is
rather intuitive. In a non-interacting Bose gas, assuming translational
invariance, the populations $N_{\bm{k}}$ of each (quasi-)momentum mode
${\bm{k}}$ are independent of each other (up to the global constraint
$\sum_{\bm{k}}N_{\bm{k}}=N$ with $N$ the total number of particles), and
populations fluctuate incoherently (e.g. thermal fluctuations). When
interactions are present, momentum modes are instead coupled and correlations
between them emerge. The most elementary process of momentum-mode correlations
comes from an elastic interaction between two particles, coupling two incoming
modes $({\bm{k}}_{1},{\bm{k}}_{2})$ to two outgoing ones
$({\bm{k}}_{3},{\bm{k}}_{4})$ which conserve the total momentum,
${\bm{k}}_{1}+{\bm{k}}_{2}={\bm{k}}_{3}+{\bm{k}}_{4}$. As we shall discuss
below, this four-wave mixing process generates a hierarchy of subsets of
correlated momentum modes.

To exploit momentum conservation in the four-wave mixing process, two-body
contact interactions should be expressed in the momentum basis. The contact
interaction term $V_{\rm int}$ between bosons inside a volume $V$ reads
$V_{\rm
int}=\frac{U}{2V}\sum_{{\bm{k}},{\bm{k}}^{\prime},{\bm{k}}^{\prime\prime}}(a^{\dagger}_{{\bm{k}}}a^{\dagger}_{{\bm{k}}^{\prime}}a_{{\bm{k}}^{\prime\prime}}a_{{\bm{k}}+{\bm{k}}^{\prime}-{\bm{k}}^{\prime\prime}}+{\rm
h.c.})$ where $a_{\bm{k}},a_{\bm{k}}^{\dagger}$ are destruction and creation
operators for bosonic momentum mode $\bm{k}$, $U$ is the interaction strength,
and ${\rm h.c.}$ indicates Hermitian conjugates. When interactions are weak
and the temperature is close to zero, the state of the system is a nearly
perfect condensate, namely $N_{\bm{k}=0}=N_{0}\approx N$. Then one can adopt
the Bogoliubov replacement $a_{0},a_{0}^{\dagger}\approx\sqrt{N_{0}}$, and
rewrite the interactions as $V_{\rm int}\approx V^{(2)}+V^{(3)}+V^{(4)}$
Fetter and Walecka (2012) where

$\displaystyle V^{(2)}$ $\displaystyle=$
$\displaystyle\frac{UN_{0}^{2}}{2V}+\frac{UN_{0}}{2V}\sum_{{\bm{k}}\neq{\bm{0}}}\left(a^{\dagger}_{{\bm{k}}}a^{\dagger}_{-{\bm{k}}}+{\rm
h.c.}\right)$ (1) $\displaystyle V^{(3)}$ $\displaystyle=$
$\displaystyle\frac{U\sqrt{N_{0}}}{2V}\sum_{{\bm{k}},{\bm{k}^{\prime}}}\left(\hat{a}^{\dagger}_{{\bm{k}+\bm{k}^{\prime}}}\hat{a}^{\dagger}_{{-\bm{k}^{\prime}}}\hat{a}_{{\bm{k}}}+{\rm
h.c.}\right)$ (2) $\displaystyle V^{(4)}$ $\displaystyle=$
$\displaystyle\frac{U}{2V}{\sum_{{\bm{k}},{\bm{k}}^{\prime},{\bm{k}}^{\prime\prime}}}\left(a^{\dagger}_{{\bm{k}}+{\bm{k}}^{\prime}-{\bm{k}}^{\prime\prime}}a^{\dagger}_{{\bm{k}}^{\prime\prime}}a_{{\bm{k}}^{\prime}}a_{{\bm{k}}}+{\rm
h.c.}\right)~{}~{}.$ (3)

Here the sums are restricted to momenta
${\bm{k}},{\bm{k}^{\prime}},{\bm{k}^{\prime\prime}}$ such that only operators
of modes with non-zero momentum are included. Bogoliubov theory, mentioned in
the introduction, amounts to retaining only the quadratic term Eq. (1) in the
above expression, which corresponds to linearizing quantum fluctuations, and
describing only pairwise correlations among opposite momenta.

Moving beyond the weakly-interacting regime, quantum fluctuations become too
large to be linearized, as the quantum depletion increases monotonically with
the interaction strength. As a consequence the cubic and quartic terms, Eqs.
(2) and (3), can no longer be neglected. Interestingly, momentum correlations
induced by these terms have a clear physical picture when considering
successive scattering events, as illustrated in Fig. 1b-c. The cubic term
$V^{(3)}$ of Eq. (2) correlates triplets of modes whose momenta sum up to
zero, $\bm{K}=\sum_{j=1}^{3}\bm{k}_{j}=\bm{0}$ while the quartic term
$V^{(4)}$ of Eq. (3) correlates (among others) quartets of modes with
$\bm{K}=\sum_{j=1}^{4}\bm{k}_{j}=\bm{0}$.

Figure 2: Momentum density profiles and two-body connected correlations at
opposite momenta. a Normalized 1D cuts of the momentum density
$\rho(k)/\rho(0)$ for varying ratios $u=U/J\in[2;22]$. b Same as panel a in
semi-log scale. The green shaded area represents the integration volume
$\Omega_{k}$ of the depletion over which two-body correlations are computed. c
Condensate fraction $f_{c}$ as a function of $u=U/J$. The dashed line is a
guide to the eye $\propto(1-u/u_{c})^{2\beta}$ with $\beta=0.3485$ and
$u_{c}=26$ Hercé  _et al._ (2021). Error bars reflect one standard deviation.
d Amplitude $\overline{G}_{\rm c}^{(2)}(\bm{K}=\bm{0})$ of the Bogoliubov
momentum pairing at opposite momenta
(${\bm{K}}={\bm{k}}_{1}+{\bm{k}}_{2}={\bm{0}}$) measured in the volume
$\Omega_{k}$ occupied by the depletion of the lattice BEC (orange circles) and
of a harmonically-trapped BEC (orange square). The dashed-line is a guide to
the eye. The vertical error bars correspond to $68\%$-confidence intervals
obtained from a bootstrap method. Inset: $\overline{G}_{\rm
c}^{(2)}(\bm{K}=\bm{0})$ plotted as a fonction of the condensate fraction
$f_{c}$.

The emergence of correlated sets of $n>2$ modes, whose correlations are not
reducible to pairwise ones, clearly brings to a departure from the Gaussian-
state description. A quantitative understanding of this regime, and in
particular of the structure of its correlations, remains challenging (see e.g.
Ref. Carlen  _et al._ (2021) for a recent theoretical effort in continuum
space). To address this problem experimentally, we probe correlated sets of
$n$ momentum modes (depicted in Fig. 1a-c) via connected correlation functions
$G_{\rm c}^{(n)}$ which are $n$-th order multivariate cumulants Ursell (1927).
Here, we probe connected correlations among $n$ momentum populations, $G_{\rm
c}^{(n)}(N_{{\bm{k}}_{1}},N_{{\bm{k}}_{2}},...,N_{{\bm{k}}_{n}})=\langle
N_{{\bm{k}}_{1}}N_{{\bm{k}}_{2}}...N_{{\bm{k}}_{n}}\rangle-G_{\rm
dis}^{(n)}(N_{{\bm{k}}_{1}},N_{{\bm{k}}_{2}},...,N_{{\bm{k}}_{n}})$ where the
disconnected part $G_{\rm dis}^{(n)}$ includes the contributions of all lower-
order correlations $n^{\prime}<n$ to the total $n$-body correlations. The
$G_{\rm c}^{(n)}$ function isolates therefore genuine $n$-mode correlations,
i.e. those not reducible to correlations of lower order. Connected
correlations among zero-sum momentum modes, $G_{\rm c}^{(n)}(\bm{K}=0)$,
represent the most relevant form of quantum correlations induced by
interactions, both at the level of Bogoliubov theory as well as beyond (Fig.
1a-c). On the other hand, incoherent (thermal) fluctuations are not subject to
any momentum-sum constraint (as the momentum sum itself can fluctuate
incoherently). Hence detecting correlations among zero-sum modes is a strong
indication of their quantum origin, caused by processes of the kind depicted
in Fig. 1a-c. In the following, we shall focus our attention on the two-mode,
zero-sum connected correlations $G_{\rm c}^{(2)}(\bm{k},-\bm{k})=\langle
a^{\dagger}_{\bm{k}}a^{\dagger}_{-\bm{k}}a_{\bm{k}}a_{-\bm{k}}\rangle-\langle
N_{\bm{k}}\rangle\langle N_{-\bm{k}}\rangle$. Bogoliubov theory predicts these
correlations to take the form $G^{(2)}_{c}({\bm{k}},-{\bm{k}})=|\langle
a^{\dagger}_{\bm{k}}a^{\dagger}_{-\bm{k}}\rangle|^{2}$ (Supplementary
Information), namely as stemming from momentum pairing in the system.

Our experiment begins with the production of a metastable helium-4
$\left({}^{4}\textrm{He}^{*}\right)$ Bose-Einstein condensate (BEC) in a
crossed optical dipole trap Bouton  _et al._ (2015). The BEC is then
adiabatically transferred to a 3D optical lattice Carcy  _et al._ (2021),
which we use to increase the interactions and the quantum depletion Tenart
_et al._ (2021). In the lattice, the Bose gas is well described by the Bose-
Hubbard Hamiltonian, $H=-J\sum_{\langle
ij\rangle}(b_{i}^{\dagger}b_{j}+b_{j}^{\dagger}b_{i})+V_{\rm int}$ with
tunneling energy $J$ and on-site interaction energy $U$ – here $\langle
ij\rangle$ denotes a pair of nearest neighbors and $b_{j}^{\dagger}$ (resp.
$b_{j}$) the creation (resp. destruction) operator on site $j$. The lattice
wave vector is $k_{d}=2\pi/d$ where the $d=775$ nm is the lattice spacing. We
access momentum correlations in a time-of-flight (TOF) experiment where the
4He∗ atoms fall onto a pair of Micro-Channel Plates (MCPs) located $43~{}$cm
below the optical lattice. The combination of MCPs with delay lines Nogrette
_et al._ (2015) yields the 3D single-atom-resolved distribution of the gas in
momentum space, as illustrated in Fig. 1d. We use a total atom number
$N=5.0(7)\times 10^{3}$ to ensure that the lattice filling is unity in the
center of the trap and that interactions are negligible during the TOF
dynamics. Under these conditions, we are able to measure the in-trap momenta
of individual atoms Cayla  _et al._ (2018); Tenart  _et al._ (2020).

We have recorded atom distributions at various interaction strengths $u=U/J$
in the range $u\in[2,22]$ ($\sim 2300$ distributions for each $u$), far from
the transition to the (finite-temperature) Mott insulator, occurring at the
critical ratio $u_{c}=26(1)$ in our experiment Hercé  _et al._ (2021). This is
illustrated in Figure 2a, where normalized 1D cuts of the momentum density are
shown to exhibit contrasted peaks associated to the presence of large
condensate fractions. The semi-logarithmic scale in Figure 2b highlights the
increase in the depletion density by more than one order of magnitude at the
edge of the first Brillouin zone $k=0.5~{}k_{d}$. The corresponding condensate
fraction $f_{c}$ is plotted as a function of interactions $u=U/J$ in Fig. 2c.
Importantly, the reduced temperature $T/J$ is constant (within error bars) in
the range studied here Carcy  _et al._ (2021). The observed increase of the
depletion in Fig. 2a-c thus results predominantly from the increase of
interactions (i.e. not an effect of heating). In the following we concentrate
our attention on a region of momentum space, denoted by $\Omega_{k}$ (Fig.
2b), which captures the quantum depletion of the condensate. In this work, we
use for $\Omega_{k}$ a cubic corona comprised between two cubes centred on
${\bm{k}}={\bm{0}}$ and of sides $2\times 0.7~{}k_{d}$ and $2\times
0.3~{}k_{d}$ respectively (Fig. 2b and Methods).

From the experimental data, we extract the connected correlations at opposite
momenta, $G_{\rm c}^{(2)}(\bm{k},-\bm{k})$, for momenta $\bm{k}$ belonging to
$\Omega_{k}$. To increase the signal-to-noise ratio, we sum $G_{\rm
c}^{(2)}(\bm{k},-\bm{k})$ over the modes occupying the volume $\Omega_{k}$,
defining the integrated amplitude $\overline{G}_{\rm
c}^{(2)}(\bm{0})=\sum_{\bm{k}\in\Omega_{k}}\ G_{\rm c}^{(2)}(\bm{k},-\bm{k})$.
In Figure 2d, we plot $\overline{G}_{\rm c}^{(2)}(\bm{0})$ as a function of
the interaction strength $u$. The measured values have been multiplied by
$\eta^{-2}$ to account for the finite detection efficiency $\eta=0.53(2)$ of
our detector. At small values of $u$, $\overline{G}_{\rm c}^{(2)}(\bm{0})$
increases with $u$. This observation is compatible with Bogoliubov theory,
which predicts that the two-mode correlated quantum depletion always increases
with $u$ (Fig. 3a) Toth  _et al._ (2008). In contrast, our measurement shows
that $\overline{G}_{\rm c}^{(2)}(\bm{0})$ decreases at larger values of $u$.
The position of the maximum corresponds to a relatively large condensate
fraction of $f_{c}\sim 0.75$ (inset of Fig. 2d), far from the Mott transition.
This strongly suggests that our observations are not specific to lattice
bosons at integer filling, but are generic of strongly interacting Bose gases.

We stress that this non-monotonic variation of $\overline{G}_{\rm
c}^{(2)}(\bm{0})$ results from an adiabatic change of the interaction strength
Carcy  _et al._ (2021), and is not the result of heating. It should also be
contrasted with the variation of the bunching amplitude $g^{(2)}_{\rm
HBT}(0)=\left(\sum_{\bm{k}\in\Omega_{k}}\langle
a^{\dagger}_{\bm{k}}a^{\dagger}_{\bm{k}}a_{\bm{k}}a_{\bm{k}}\rangle\right)/\left(\sum_{\bm{k}\in\Omega_{k}}\langle
N_{\bm{k}}\rangle^{2}\right)$ set by the Bose statistics. Indeed, we find a
constant and perfectly-contrasted value $g^{(2)}_{\rm HBT}(0)=2.00(5)$ in all
regimes of interactions (Supplementary Information). This result, expected at
weak interactions, extends our previous observations Cayla  _et al._ (2020) to
the strongly-interacting regime. Finally, a further observation is the
increase of the fluctuations of $\overline{G}_{\rm c}^{(2)}(\bm{0})$ with
interactions shown by the error bars (see Fig. 2d) obtained from a bootstrap
analysis of the connected correlations measured in each experimental shot.
This analysis suggests that the number of pairs at opposite momenta fluctuates
more as interactions increase. An interesting question for a future work is to
understand whether this behaviour stems from the increasing complexity of the
equilibrium many-body state at strong interactions, made of correlated subsets
with $n>2$ modes.

To investigate the fate of two-mode connected correlations $\overline{G}_{\rm
c}^{(2)}(\bm{0})$ beyond the weakly interacting regime theoretically, we make
use of numerically exact quantum Monte Carlo (QMC) simulations on a 3D quantum
rotor model (Supplementary Information), which captures the physics of the
Bose-Hubbard model at large, integer filling $n_{\rm QR}$, and without
confining potentials Wallin  _et al._ (1994). Calculations on a quantum-rotor
model are not expected to be quantitatively accurate when compared with
experiments conducted at nearly unit filling in a harmonic trap. Yet, unlike
QMC simulations for the Bose-Hubbard model, they give access to arbitrary
correlations in momentum space without the need for specialized update
algorithms Roscilde (2008); Fang  _et al._ (2016); and they do so across all
regimes of temperatures and interactions up to the Mott transition
(Supplementary Information). Fig. 3a shows $\overline{G}_{\rm
c}^{(2)}(\bm{0})$, integrated over quasi-momenta in a cubic corona
$\Omega_{k}$ between two cubes of sides $2\times 0.3k_{d}$ and $2\times
0.5k_{d}$, for a system of $10^{3}$ cubic lattice sites, as a function of the
interaction strength, and for different temperatures. The theoretical data
have been normalized to the atom number $\langle N_{\Omega_{k}}\rangle$
detected experimentally at a given interaction strength $u$ (Methods). The
temperature of the experiment was previously estimated to be around
$k_{B}T/J\approx 2$ Carcy  _et al._ (2021), corresponding to the data for the
quantum rotor model at a temperature $T_{J}=k_{B}T/(2Jn_{\rm QR})=1$.

We observe that the quantum-rotor predictions capture the non-monotonic
behavior of $\overline{G}_{\rm c}^{(2)}(\bm{0})$, as well as the right order
of magnitude of the measured correlations. Our numerical results confirm that
the decay of the connected two-body correlations $\overline{G}_{\rm
c}^{(2)}(\bm{0})$ originate from quantum fluctuations beyond the linearized
Bogoliubov regime. The maximum value of $\overline{G}_{\rm c}^{(2)}(\bm{0})$
can therefore be regarded as a microscopic signature of the Bose gas entering
the strongly correlated regime, whose description requires accounting for
correlations among $n>2$ momentum modes.

Figure 3: Numerical calculations of connected correlations $\overline{G}_{\rm
c}^{(n)}(\bm{0})$ with the quantum rotor model. a Two-mode connected
correlations $\overline{G}_{\rm c}^{(2)}(\bm{0})$ obtained with the quantum
rotor model (solid lines) at different temperatures, ranging from
$T_{J}=k_{B}T/2n_{\rm QR}J=0.2$ (dark blue) to $T_{J}=1.4$ (dark red) in steps
of 0.4. The system size is $10^{3}$. The dotted line is the result from
Bogoliubov theory for the Bose-Hubbard model with $10^{3}$ sites and unity
filling. Inset: log-plot of the data from the main panel. The amplitude of
$\overline{G}_{\rm c}^{(2)}(\bm{0})$ is normalized to the measured atom number
$\langle N_{\Omega_{k}}\rangle_{\rm exp}$ in the experiment (Methods and
Supplementary Information). b Three-mode connected correlations
$\overline{G}_{\rm c}^{(3)}(\bm{0})$ obtained with the quantum rotor model at
the various temperatures. The amplitude of $\overline{G}_{\rm
c}^{(3)}(\bm{0})$ is normalized to the measured atom number (Methods and
Supplementary Information).

To proceed with probing the hierarchy of momentum correlations illustrated in
Fig. 1, we studied the integrated three-mode correlations $\overline{G}_{\rm
c}^{(3)}(\bm{0})=\sum_{\bm{k}_{1},\bm{k}_{2}}G^{(3)}_{\rm
c}(\bm{k}_{1},\bm{k}_{2},-\bm{k}_{1}-\bm{k}_{2})$ that reveal genuine three-
mode correlations with zero momentum sum (Supplementary Information). In the
experiment, the measured values of $\overline{G}_{\rm c}^{(3)}(\bm{0})$ are
consistent with zero for all the interaction strengths we probe. This result
leaves the question open on the magnitude expected for three-mode
correlations. We thus evaluate these correlations theoretically using the
quantum-rotor model; and in Fig. 3b we plot $\overline{G}_{\rm
c}^{(3)}(\bm{0})$, integrated over trios of modes in the cubic corona
$\Omega_{k}$. We observe that genuine three-mode correlations are positive, in
agreement with the picture offered in Fig. 1, and rise along with the two-mode
correlations. This shows that the Bogoliubov picture of the Bose gas with
intermediate interactions is in fact incomplete, and that significant non-
Gaussian correlations are present in the system. Nonetheless the magnitude of
$\overline{G}_{\rm c}^{(3)}(\bm{0})$ is significantly smaller than the one of
$\overline{G}_{\rm c}^{(2)}(\bm{0})$, the temperature and interaction being
the same. While two-mode correlations with total zero (quasi-)momentum
uniquely identify the correlated pairs of modes, three-mode correlations with
zero-sum (quasi-)momentum can couple every mode to $O(N)$ pairs of other
modes, making correlations much weaker within each trio. Correlations are
expected to be even more tenuous upon increasing the number $n$ of modes.
Moreover, revealing triple coincidences of atoms in three modes is made
difficult by the finite-detection efficiency $\eta$ of our detector. These
aspects clearly make the observation of $n$-mode correlations with $n>2$
challenging. Finally, both forms of correlations (two-mode and three-mode)
decrease upon moving to the strongly correlated regime; in the correlated
superfluid regime they should scale linearly with the particle number, but
they both become of $O(1)$ (i.e. not scaling with the system size) in the deep
Mott limit, where they uniquely stem from particle-number conservation
(Supplemental Material).

As we shall now explain, measuring $\overline{G}_{\rm c}^{(3)}(0)\neq 0$ is,
however, not the only means to demonstrate non-Gaussian statistics. In the
spirit of Bogoliubov theory, we relax the assumption of particle-number
conservation, but we consider that all the interaction processes described in
Eqs. (1), (2) and (3) do not lead to the appearance of a finite value of
$\langle a_{\bm{k}}\rangle$ in the depletion. This assumption is compatible
with the perfectly-contrasted bunching, $g^{(2)}_{\rm HBT}(0)=2.00(5)$,
measured in the experiment as a finite value of $\langle a_{\bm{k}}\rangle$
would reduce the bunching amplitude. Hence, the connected two-mode
correlations can then be expressed in terms of the 4-th-order cumulant
$\langle
a^{\dagger}_{\bm{k}}a^{\dagger}_{-\bm{k}}a_{\bm{k}}a_{-\bm{k}}\rangle_{c}$
(Supplemental Material),

$\displaystyle G^{(2)}_{c}({\bm{k}},-{\bm{k}})$ $\displaystyle=$
$\displaystyle|\langle
a^{\dagger}_{\bm{k}}a^{\dagger}_{-\bm{k}}\rangle|^{2}+|\langle
a^{\dagger}_{\bm{k}}a_{-\bm{k}}\rangle|^{2}$ (4) $\displaystyle+$
$\displaystyle\langle
a^{\dagger}_{\bm{k}}a^{\dagger}_{-\bm{k}}a_{\bm{k}}a_{-\bm{k}}\rangle_{c}.$

Eq. (4) shows that negative values of $G^{(2)}_{c}({\bm{k}},-{\bm{k}})$
necessarily imply that $\langle
a^{\dagger}_{\bm{k}}a^{\dagger}_{-\bm{k}}a_{\bm{k}}a_{-\bm{k}}\rangle_{c}\neq
0$, which is a signature of non-Gaussian statistics. Indeed a negative
correlation $\overline{G}_{\rm c}^{(2)}({\bm{k}},-{\bm{k}})<0$ amounts to
anti-bunching, namely the fact that the coincidence counts $\langle
a^{\dagger}_{\bm{k}}a^{\dagger}_{-\bm{k}}a_{\bm{k}}a_{-\bm{k}}\rangle$ is
smaller than what expected in the absence of correlations, $\langle
N_{\bm{k}}\rangle\langle N_{-\bm{k}}\rangle$. This only occurs in the presence
of connected correlations of order $n>2$.

In the experiment, we observe a drastic change in the shape of the integrated
correlation function $\overline{G}_{\rm
c}^{(2)}({\bm{k}},-{\bm{k}}+\delta{\bm{k}})$ at large interactions $u=20$ and
$22$: there is hardly any correlation peak at $\delta{\bm{k}}={\bm{0}}$ but a
constant (negative) offset as a function of $\delta{\bm{k}}$ (Supplemental
Material). In this regime, we thus find $\overline{G}_{\rm c}^{(2)}(\bm{0})<0$
(see Fig. 2d) and $\sum_{{\bm{k}}}\langle
a^{\dagger}_{\bm{k}}a^{\dagger}_{-\bm{k}}a_{\bm{k}}a_{-\bm{k}}\rangle_{c}\neq
0$. This non-zero fourth-order cumulant is an unambiguous signature of non-
Gaussian statistics that does not rely on any theoretical modeling. As such,
it represents a striking evidence of the strongly correlated (non-Gaussian)
nature of Bose gases at large interactions.

In conclusion, by making use of state-of-the-art preparation and detection
techniques for ultracold gases, and of advanced numerical methods, we provide
a physical picture of correlations in an interacting Bose gas beyond the
Bogoliubov regime. Two-mode correlations, at the heart of the Bogoliubov
weakly interacting regime, are suppressed when entering the strongly
interacting regime, in favor of $n$-mode correlations with $n>2$, and non-
Gaussian correlations appear. These results suggest that, upon increasing the
interactions, entanglement in momentum space (that reflects itself in the
correlations we probed in this work) becomes increasingly complex in its
multipartite nature; as well as also increasingly elusive, since each mode
becomes weakly entangled with a myriad (${\cal O}(N^{n-2})$) of $(n-1)$-plets
of other modes to form $n$-mode entangled clusters. Finding a proper
theoretical characterization, as well as experimental signatures, of this
complex entanglement structure represents a fascinating open question for
future studies.

## Acknowledgements

We thank E. Gradova for her contribution to the bootstrap analysis, T.
Chalopin for a careful reading of the manuscript, D. Boiron and the members of
the Quantum Gas group at Institut d’Optique for insightful discussions. We
acknowledge financial support from the Région Ile-de-France in the framework
of the DIM SIRTEQ, the “Fondation d’entreprise iXcore pour la Recherche”, the
French National Research Agency (Grant number ANR-17-CE30-0020-01) and France
2030 programs of the French National Research Agency (Grant number
ANR-22-PETQ-0004). All numerical simulations were made on the PSMN cluster at
the ENS of Lyon.

## Methods

Measurement of connected correlations at opposite momenta. As indicated in the
main text, the volume $\Omega_{k}$ of the depletion in which we compute
correlations is the volume between two cubes of sides $2\times 0.3~{}k_{d}$
and $2\times 0.7~{}k_{d}$. Since $\Omega_{k}$ corresponds to
$1.8\leq|k|\xi\leq 3.0$ at $U/J=2$ and $1.2\leq|k|\xi\leq 2.0$ at $U/J=22$,
with $\xi$ the healing length, it excludes by design the phononic part
$\left(|k|\xi\ll 1\right)$ of the thermal depletion. Two-body correlations at
opposite momenta are computed from (i) calculating the histogram of momentum
sums ${\bm{k}}_{1}+{\bm{k}}_{2}$ for all pairs of atoms detected in a run of
the experiment (from which we subtract the product of average occupations
$\langle N(\bm{k}_{1})\rangle\langle N(\bm{k}_{2})\rangle$), and (ii)
integrating this histogram over a cubic volume centered on
${\bm{k}}_{1}+{\bm{k}}_{2}={\bm{0}}$ and of size $0.18~{}k_{d}$ set by the
width of the two-body correlation function.

Comparison of numerical results from quantum rotor model with experiments.
Connected correlations $\overline{G}_{\rm c}^{(2)}(0)$ and $\overline{G}_{\rm
c}^{(3)}(0)$ are computed numerically with the quantum rotor model over a
cubic corona $\Omega_{k}$ between two cubes of sides $(2\times 0.3~{}k_{d}$
and $2\times 0.5~{}k_{d}$. These calculations are performed with a large,
integer filling $n_{\rm QR}$ of the 3D lattice (Supplemental material), which
leaves us with an adjustable parameter to compare with experiments (performed
in the Bose-Hubbard regime with one atom per site at the trap center). The
numerically-computed population $\langle N_{\Omega_{k}}^{th}\rangle$ of the
modes in momentum space changes with the filling $n_{\rm QR}$. To plot the
numerical results of $\overline{G}_{\rm c}^{(2)}(0)$ (resp. $\overline{G}_{\rm
c}^{(3)}(0)$) in the main text, we have rescaled their amplitude by $(\langle
N_{\Omega_{k}}\rangle/\langle N_{\Omega_{k}}^{th}\rangle)^{2}$ (resp.
$(\langle N_{\Omega_{k}}\rangle/\langle N_{\Omega_{k}}^{th}\rangle)^{3}$),
where $\langle N_{\Omega_{k}}\rangle$ is the measured momentum population
observed in the experiment.

## References

  * Sachdev (2023) Sachdev, S.,  _Quantum Phases of Matter_ , 1st ed. (Cambridge University Press, 2023).
  * Morosan  _et al._ (2012) Morosan, E., Natelson, D., Nevidomskyy, A. H. and Si, Q.,  _Strongly Correlated Materials_ , Advanced Materials 24, 4896 (2012).
  * Adams  _et al._ (2012) Adams, A., Carr, L. D., Schäfer, T., Steinberg, P. and Thomas, J. E.,  _Strongly correlated quantum fluids: ultracold quantum gases, quantum chromodynamic plasmas and holographic duality_ , New Journal of Physics 14, 115009 (2012).
  * Bloch  _et al._ (2022) Bloch, J., Cavalleri, A., Galitski, V., Hafezi, M. and Rubio, A.,  _Strongly correlated electron–photon systems_ , Nature 606, 41 (2022).
  * Chevy and Salomon (2016) Chevy, F. and Salomon, C.,  _Strongly correlated Bose gases_ , Journal of Physics B: Atomic, Molecular and Optical Physics 49, 192001 (2016).
  * Georges  _et al._ (1996) Georges, A., Kotliar, G., Krauth, W. and Rozenberg, M. J.,  _Dynamical mean-field theory of strongly correlated fermion systems and the limit of infinite dimensions_ , Rev. Mod. Phys. 68, 13 (1996).
  * Metzner  _et al._ (2012) Metzner, W., Salmhofer, M., Honerkamp, C., Meden, V. and Schönhammer, K.,  _Functional renormalization group approach to correlated fermion systems_ , Rev. Mod. Phys. 84, 299 (2012).
  * Georges  _et al._ (2013) Georges, A., Medici, L. d. and Mravlje, J.,  _Strong Correlations from Hund’s Coupling_ , Annual Review of Condensed Matter Physics 4, 137 (2013).
  * Chang  _et al._ (2014) Chang, D. E., Vuletić, V. and Lukin, M. D.,  _Quantum nonlinear optics — photon by photon_ , Nature Photonics 8, 685 (2014).
  * Pitaevskii and Stringari (2016) Pitaevskii, L. and Stringari, S.,  _Bose-Einstein Condensation and Superfluidity_ (Oxford University Press, 2016).
  * Leggett (2006) Leggett, A.,  _Quantum Liquids: Bose Condensation and Cooper Pairing in Condensed-matter Systems_, Oxford Graduate Texts (OUP Oxford, 2006).
  * Makotyn  _et al._ (2014) Makotyn, P., Klauss, C. E., Goldberger, D. L., Cornell, E. A. and Jin, D. S.,  _Universal dynamics of a degenerate unitary Bose gas_ , Nature Physics 10, 116 (2014).
  * Fletcher  _et al._ (2017) Fletcher, R. J., Lopes, R., Man, J., Navon, N., Smith, R. P., Zwierlein, M. W. and Hadzibabic, Z.,  _Two- and three-body contacts in the unitary Bose gas_ , Science 355, 377 (2017).
  * Eismann  _et al._ (2016) Eismann, U. _et al._ ,  _Universal Loss Dynamics in a Unitary Bose Gas_ , Phys. Rev. X 6, 021025 (2016).
  * Klauss  _et al._ (2017) Klauss, C. E., Xie, X., Lopez-Abadia, C., D’Incao, J. P., Hadzibabic, Z., Jin, D. S. and Cornell, E. A.,  _Observation of Efimov Molecules Created from a Resonantly Interacting Bose Gas_ , Phys. Rev. Lett. 119, 143401 (2017).
  * Eigen  _et al._ (2018) Eigen, C., Glidden, J. A. P., Lopes, R., Cornell, E. A., Smith, R. P. and Hadzibabic, Z.,  _Universal prethermal dynamics of Bose gases quenched to unitarity_ , Nature 563, 221 (2018).
  * Yan  _et al._ (2020) Yan, Z. Z., Ni, Y., Robens, C. and Zwierlein, M. W.,  _Bose polarons near quantum criticality_ , Science 368, 190 (2020).
  * Bogoliubov (1947) Bogoliubov, N.,  _On the theory of superfluidity_ , Journal of Physics (USSR) 11, 23 (1947).
  * Tenart  _et al._ (2021) Tenart, A., Hercé, G., Bureik, J.-P., Dareau, A. and Clément, D.,  _Observation of pairs of atoms at opposite momenta in an equilibrium interacting Bose gas_ , Nature Physics 17, 1364 (2021).
  * Paschen and Si (2021) Paschen, S. and Si, Q.,  _Quantum phases driven by strong correlations_ , Nature Reviews Physics 3, 9 (2021).
  * Polkovnikov  _et al._ (2011) Polkovnikov, A., Sengupta, K., Silva, A. and Vengalattore, M.,  _Colloquium: Nonequilibrium dynamics of closed interacting quantum systems_ , Rev. Mod. Phys. 83, 863 (2011).
  * Wang  _et al._ (2007) Wang, X.-B., Hiroshima, T., Tomita, A. and Hayashi, M.,  _Quantum information with Gaussian states_ , Physics Reports 448, 1 (2007).
  * Armijo  _et al._ (2010) Armijo, J., Jacqmin, T., Kheruntsyan, K. V. and Bouchoule, I.,  _Probing Three-Body Correlations in a Quantum Gas Using the Measurement of the Third Moment of Density Fluctuations_ , Phys. Rev. Lett. 105, 230402 (2010).
  * Schweigler  _et al._ (2017) Schweigler, T., Kasper, V., Erne, S., Mazets, I., Rauer, B., Cataldini, F., Langen, T., Gasenzer, T., Berges, J. and Schmiedmayer, J.,  _Experimental characterization of a quantum many-body system via higher-order correlations_ , Nature 545, 323 (2017).
  * Rispoli  _et al._ (2019) Rispoli, M., Lukin, A., Schittko, R., Kim, S., Tai, M. E., Léonard, J. and Greiner, M.,  _Quantum critical behaviour at the many-body localization transition_ , Nature 573, 385 (2019).
  * Lee  _et al._ (1957) Lee, T. D., Huang, K. and Yang, C. N.,  _Eigenvalues and Eigenfunctions of a Bose System of Hard Spheres and Its Low-Temperature Properties_ , Phys. Rev. 106, 1135 (1957).
  * Navon  _et al._ (2011) Navon, N., Piatecki, S., Günter, K., Rem, B., Nguyen, T. C., Chevy, F., Krauth, W. and Salomon, C.,  _Dynamics and Thermodynamics of the Low-Temperature Strongly Interacting Bose Gas_ , Phys. Rev. Lett. 107, 135301 (2011).
  * Jørgensen  _et al._ (2018) Jørgensen, N. B., Bruun, G. M. and Arlt, J. J.,  _Dilute Fluid Governed by Quantum Fluctuations_ , Phys. Rev. Lett. 121, 173403 (2018).
  * Pitaevskii and Stringari (1991) Pitaevskii, L. and Stringari, S.,  _Uncertainty principle, quantum fluctuations, and broken symmetries_ , Journal of Low Temperature Physics 85, 377 (1991).
  * Lopes  _et al._ (2017) Lopes, R., Eigen, C., Navon, N., Clément, D., Smith, R. P. and Hadzibabic, Z.,  _Quantum Depletion of a Homogeneous Bose-Einstein Condensate_ , Phys. Rev. Lett. 119, 190404 (2017).
  * Ceperley (1995) Ceperley, D. M.,  _Path integrals in the theory of condensed helium_ , Rev. Mod. Phys. 67, 279 (1995).
  * Svistunov  _et al._ (2015) Svistunov, B. V., Babaev, E. S. and Prokof’ev, N. V.,  _Superfluid States of Matter_ (CRC Press, 2015).
  * Fisher  _et al._ (1989) Fisher, M. P. A., Weichman, P. B., Grinstein, G. and Fisher, D. S.,  _Boson localization and the superfluid-insulator transition_ , Phys. Rev. B 40, 546 (1989).
  * Bloch  _et al._ (2008) Bloch, I., Dalibard, J. and Zwerger, W.,  _Many-body physics with ultracold gases_ , Rev. Mod. Phys. 80, 885 (2008).
  * Georgescu  _et al._ (2014) Georgescu, I. M., Ashhab, S. and Nori, F.,  _Quantum simulation_ , Rev. Mod. Phys. 86, 153 (2014).
  * Cayla  _et al._ (2018) Cayla, H., Carcy, C., Bouton, Q., Chang, R., Carleo, G., Mancini, M. and Clément, D.,  _Single-atom-resolved probing of lattice gases in momentum space_ , Phys. Rev. A 97, 061609 (2018).
  * Hercé  _et al._ (2023) Hercé, G., Bureik, J.-P., Ténart, A., Aspect, A., Dareau, A. and Clément, D.,  _Full counting statistics of interacting lattice gases after an expansion: The role of condensate depletion in many-body coherence_ , Phys. Rev. Res. 5, L012037 (2023).
  * Fetter and Walecka (2012) Fetter, A. L. and Walecka, J. D.,  _Quantum Theory of Many-Particle Systems_ (Dover, 2012).
  * Hercé  _et al._ (2021) Hercé, G., Carcy, C., Tenart, A., Bureik, J.-P., Dareau, A., Clément, D. and Roscilde, T.,  _Studying the low-entropy Mott transition of bosons in a three-dimensional optical lattice by measuring the full momentum-space density_ , Phys. Rev. A 104, L011301 (2021).
  * Carlen  _et al._ (2021) Carlen, E. A., Holzmann, M., Jauslin, I. and Lieb, E. H.,  _Simplified approach to the repulsive Bose gas from low to high densities and its numerical accuracy_ , Phys. Rev. A 103, 053309 (2021).
  * Ursell (1927) Ursell, H. D.,  _The evaluation of Gibbs’ phase-integral for imperfect gases_ , Mathematical Proceedings of the Cambridge Philosophical Society 23, 685–697 (1927).
  * Bouton  _et al._ (2015) Bouton, Q., Chang, R., Hoendervanger, A. L., Nogrette, F., Aspect, A., Westbrook, C. I. and Clément, D.,  _Fast production of Bose-Einstein condensates of metastable helium_ , Phys. Rev. A 91, 061402 (2015).
  * Carcy  _et al._ (2021) Carcy, C., Hercé, G., Tenart, A., Roscilde, T. and Clément, D.,  _Certifying the Adiabatic Preparation of Ultracold Lattice Bosons in the Vicinity of the Mott Transition_ , Phys. Rev. Lett. 126, 045301 (2021).
  * Nogrette  _et al._ (2015) Nogrette, F., Heurteau, D., Chang, R., Bouton, Q., Westbrook, C. I., Sellem, R. and Clément, D.,  _Characterization of a detector chain using a FPGA-based time-to-digital converter to reconstruct the three-dimensional coordinates of single particles at high flux_ , Review of Scientific Instruments 86, 113105 (2015).
  * Tenart  _et al._ (2020) Tenart, A., Carcy, C., Cayla, H., Bourdel, T., Mancini, M. and Clément, D.,  _Two-body collisions in the time-of-flight dynamics of lattice Bose superfluids_ , Phys. Rev. Res. 2, 013017 (2020).
  * Toth  _et al._ (2008) Toth, E., Rey, A. M. and Blakie, P. B.,  _Theory of correlations between ultracold bosons released from an optical lattice_ , Phys. Rev. A 78, 013627 (2008).
  * Cayla  _et al._ (2020) Cayla, H., Butera, S., Carcy, C., Tenart, A., Hercé, G., Mancini, M., Aspect, A., Carusotto, I. and Clément, D.,  _Hanbury Brown and Twiss Bunching of Phonons and of the Quantum Depletion in an Interacting Bose Gas_ , Phys. Rev. Lett. 125, 165301 (2020).
  * Wallin  _et al._ (1994) Wallin, M., So/rensen, E. S., Girvin, S. M. and Young, A. P.,  _Superconductor-insulator transition in two-dimensional dirty boson systems_ , Phys. Rev. B 49, 12115 (1994).
  * Roscilde (2008) Roscilde, T.,  _Bosons in one-dimensional incommensurate superlattices_ , Phys. Rev. A 77, 063605 (2008).
  * Fang  _et al._ (2016) Fang, B., Johnson, A., Roscilde, T. and Bouchoule, I.,  _Momentum-Space Correlations of a One-Dimensional Bose Gas_ , Phys. Rev. Lett. 116, 050402 (2016).
  * Sondhi  _et al._ (1997) Sondhi, S. L., Girvin, S. M., Carini, J. P. and Shahar, D.,  _Continuous quantum phase transitions_ , Rev. Mod. Phys. 69, 315 (1997).
  * Roscilde  _et al._ (2016) Roscilde, T., Faulkner, M. F., Bramwell, S. T. and Holdsworth, P. C. W.,  _From quantum to thermal topological-sector fluctuations of strongly interacting Bosons in a ring lattice_ , New Journal of Physics 18, 075003 (2016).
  * Trotzky  _et al._ (2010) Trotzky, S., Pollet, L., Gerbier, F., Schnorrberger, U., Bloch, I., Prokof’ev, N. V., Svistunov, B. and Troyer, M.,  _Suppression of the critical temperature for superfluidity near the Mott transition_ , Nature Physics 6, 998 (2010).
  * Wessel  _et al._ (2004) Wessel, S., Alet, F., Troyer, M. and Batrouni, G. G.,  _Quantum Monte Carlo simulations of confined bosonic atoms in optical lattices_ , Phys. Rev. A 70, 053615 (2004).

SUPPLEMENTAL MATERIAL – Suppression of Bogoliubov momentum pairing and
emergence of non-Gaussian correlations in ultracold interacting Bose gases

Jan-Philipp Bureik, Gaétan Hercé, Maxime Allemand, Antoine Tenart, Tommaso
Roscilde, David Clément

##  I Normalized two-body correlation functions

We evaluate the normalized two-body correlation function
$g^{(2)}(\bm{k},\bm{k}^{\prime})$ in momentum space in the depletion of
interacting Bose gases (in the volume $\Omega_{k}$ shown in Fig. 2a of the
main text).

In the weakly-interacting regime, the two-body correlation function
$g^{(2)}(\bm{k},\bm{k}^{\prime})$ exhibits two peaks Tenart  _et al._ (2021):
one located at $\bm{k}^{\prime}=\bm{k}$ due to bosonic bunching (also coined
Hanbury-Brown and Twiss (HBT) correlations) and one located at
$\bm{k}^{\prime}=-\bm{k}$ due to interaction-induced quantum fluctuations. We
fit these peaks with Gaussian functions to obtain their amplitudes, $g_{\rm
HBT}^{(2)}(0)$ and $g_{\rm pairs}^{(2)}(0)$, with the results plotted in Fig.
S1.

Figure S1: Magnitudes $g_{\rm HBT}^{(2)}(0)$ of local $k/k$ two-body
correlations (bunching) and $g_{\rm pairs}^{(2)}(0)$ of normalized two-body
correlations at opposite momenta $k/-k$ plotted as a function of the
interaction strength $u$ for data sets shown in Fig. 2 of the main text.
Horizontal and vertical error bars correspond to one standard deviation.

The bunching amplitude is perfectly-contrasted in all regimes of interaction,
$g_{\rm HBT}^{(2)}(0)=2.00(5)$, extending previous observations Cayla  _et
al._ (2020) to the strongly-interacting regime. In contrast, $g_{\rm
pairs}^{(2)}(0)$ is found to decrease with the interaction $u$. This decrease
is compatible with the one expected in Bogoliubov theory: the magnitude
$g_{\rm pairs}^{(2)}(0)-1$ is indeed expected to decrease with the inverse of
the mode population $1/N(\bm{k})$, a property of two-mode squeezed states
Tenart  _et al._ (2021). No qualitative change is expected when the gas enters
the strongly correlated regime where the Bogoliubov approximation fails:
$g_{\rm pairs}^{(2)}(0)$ also decreases when the population of two-mode
squeezed states decreases. This illustrates that normalized correlation
functions are not suited to study the hierarchy of momentum-correlated subsets
introduced in the main text (see Fig. 1a) and it motivates the use of
connected correlations to reveal the correlated subsets of $n$ modes.

##  II Extracting two-body connected correlations

In finite-size systems, correlation functions have a finite extent – the
correlation volume – whose volume is related to that of one mode. In our
experiment, the normalized correlation function
$g^{(2)}(\bm{k},-\bm{k}+\delta\bm{k})$ is expected to be isotropic due to the
presence of an isotropic harmonic trap and is found, empirically, to be well
fitted by an isotropic Gaussian function. We introduce the correlation length
$l_{c}$ as the size of $g^{(2)}(\bm{k},-\bm{k}+\delta\bm{k})$ at $1/e^{2}$:

$g^{(2)}(\bm{k},-\bm{k}+\delta\bm{k})=1+\frac{A\sqrt{2}}{\sqrt{\pi}l_{c}}e^{-\frac{2\delta\bm{k}^{2}}{l_{c}^{2}}}.$

Thus 95% of atom pairs are expected to be found in the spherical volume where
$|\delta\bm{k}|\in[0,l_{c}]$.

The two-body connected correlation $G^{(2)}_{c}({\bm{k}},-{\bm{k}})$ is
obtained from integrating $G^{(2)}_{c}({\bm{k}},-{\bm{k}}+\delta\bm{k})$ over
a cubic volume of size $\Delta k$:

$\begin{split}G^{(2)}_{c}({\bm{k}},-{\bm{k}})=\int_{[-\Delta k/2,\Delta
k/2]^{3}}d\delta\bm{k}\ G^{(2)}_{c}(\bm{k},-\bm{k}+\delta\bm{k})).\end{split}$
(S1)

Since $G^{(2)}_{c}({\bm{k}},-{\bm{k}}+\delta\bm{k})\propto
g^{(2)}(\bm{k},-\bm{k}+\delta\bm{k})-1$ it follows that
$G^{(2)}_{c}({\bm{k}},-{\bm{k}})\propto\left[\mathrm{erf}\left(\frac{\Delta
k}{\sqrt{2}l_{c}}\right)\right]^{3}$. We choose the volume $\Delta k^{3}$ such
that the above integral on $d\delta{\bm{k}}$ saturates: we set $\Delta k=0.18\
k_{d}\simeq 3\times l_{c}$ such that $G^{(2)}_{c}({\bm{k}},-{\bm{k}})$
accounts for 99% of the pairing amplitude, theoretically. At moderate
interaction strengths $U/J<20$, this is verified in the experiment, as
illustrated in Fig. S2 where we plot the integrated version
$\bar{G}^{(2)}_{c}(\bm{0})=\sum_{\bm{k}}G^{(2)}_{c}({\bm{k}},-{\bm{k}})$
normalized by the integration volume $\Delta k^{3}$. We find indeed that
$\bar{G}^{(2)}_{c}(\bm{0})/\Delta k^{3}$ decays to zero over the expected
range $\Delta k\sim 0.18\ k_{d}$, with the variation
$\left[\mathrm{erf}\left(\frac{\Delta
k}{\sqrt{2}l_{c}}\right)\right]^{3}/\Delta k^{3}$ obtained from the bell-
shaped function of the normalized correlations. This is reminiscent of the
fact that we observe a peak in the normalized correlation
$g^{(2)}(\bm{k},-\bm{k}+\delta\bm{k})$ centered on $\delta\bm{k}={\bm{0}}$ for
$U/J<20$.

Figure S2: Two-body connected correlation $\bar{G}^{(2)}_{c}(\bm{0})/\Delta
k^{3}$ as a function of the integration volume $\Delta k^{3}$, at $U/J=5,\
10,\ 15\ \mathrm{and}\ 20$. The error bars correspond to 68%-confidence
intervals obtained from a bootstrap method. The red solid line is the expected
variation $\left[\mathrm{erf}\left(\frac{\Delta
k}{\sqrt{2}l_{c}}\right)\right]^{3}/\Delta k^{3}$ associated with the
correlation length $l_{c}$ of the normalized two-body correlation
$g^{(2)}(\bm{k},-\bm{k}+\delta\bm{k})$. Insets: zoom in the range
$d\delta{\bm{k}}\in[0.14~{}k_{d},0.38~{}k_{d}]$.

In contrast, at $U/J=20$ (and similarly at $U/J=22$),
$\bar{G}^{(2)}_{c}(\bm{0})/\Delta k^{3}$ follows a very different trend. While
it takes slightly positive values at vanishingly small values of $\Delta k$,
it is constant and negative otherwise. This different behaviour is at the
origin of the negative values $\bar{G}^{(2)}_{c}(\bm{0})$ plotted in Fig. 2 of
the main text. Observing $\bar{G}^{(2)}_{c}(\bm{0})<0$ implies that the number
of coincidence counts in modes at opposite momenta is smaller than the
coincidence counts expected in an uncorrelated ensemble of atoms with a
similar momentum density $N({\bm{k}})=\langle
a^{\dagger}_{\bm{k}}a_{\bm{k}}\rangle$, and hence that modes at opposite
momenta are anti-correlated.

##  III Connected correlations $G^{(2)}_{c}({\bf 0})$ and $\langle
a^{\dagger}_{\bf k}a^{\dagger}_{-\bf k}a_{\bf k}a_{-{\bf k}}\rangle_{c}$

We discuss the relation between the connected two-body correlations
$G^{(2)}_{c}({\bm{k}},-{\bm{k}})$ and connected four-operator correlations
$\langle
a^{\dagger}_{\bm{k}}a^{\dagger}_{-\bm{k}}a_{\bm{k}}a_{-\bm{k}}\rangle_{c}$ in
the condensate depletion.

In the condensate depletion and in a number-conserving approach, the connected
correlations of four operators write,

$\displaystyle\langle
a^{\dagger}_{\bm{k}}a^{\dagger}_{-\bm{k}}a_{\bm{k}}a_{-\bm{k}}\rangle_{c}$
$\displaystyle=$ $\displaystyle\langle
a^{\dagger}_{\bm{k}}a^{\dagger}_{-\bm{k}}a_{\bm{k}}a_{-\bm{k}}\rangle
-|\langle a^{\dagger}_{\bm{k}}a^{\dagger}_{-\bm{k}}\rangle|^{2}$
$\displaystyle-$ $\displaystyle|\langle
a^{\dagger}_{\bm{k}}a_{-\bm{k}}\rangle|^{2}-\langle
a^{\dagger}_{\bm{k}}a_{\bm{k}}\rangle\langle
a^{\dagger}_{-\bm{k}}a_{-\bm{k}}\rangle,$

and, in turn, $G^{(2)}_{c}({\bm{k}},-{\bm{k}})$ can be expressed as,

$\displaystyle G^{(2)}_{c}({\bm{k}},-{\bm{k}})$ $\displaystyle=$
$\displaystyle|\langle
a^{\dagger}_{\bm{k}}a^{\dagger}_{-\bm{k}}\rangle|^{2}+|\langle
a^{\dagger}_{\bm{k}}a_{-\bm{k}}\rangle|^{2}$ (S2) $\displaystyle+$
$\displaystyle\langle
a^{\dagger}_{\bm{k}}a^{\dagger}_{-\bm{k}}a_{\bm{k}}a_{-\bm{k}}\rangle_{c}.$

This generic equation is valid for both Gaussian and non-Gaussian states.

Bogoliubov theory. Within the Bogoliubov approximation which considers only
linearized quantum fluctuations, quantum states are Gaussian and connected
correlations of four operators are null, $\langle
a^{\dagger}_{\bm{k}}a^{\dagger}_{-\bm{k}}a_{\bm{k}}a_{-\bm{k}}\rangle_{c}=0$.
In addition, modes at opposite momenta are predicted to be incoherent, i.e.
$|\langle a^{\dagger}_{\bm{k}}a_{-\bm{k}}\rangle|^{2}=0$. Under these
additional assumptions, the above equation Eq. S2 simplifies to

$G^{(2)}_{c}({\bm{k}},-{\bm{k}})_{\rm Bogo}=|\langle
a^{\dagger}_{\bm{k}}a^{\dagger}_{-\bm{k}}\rangle|^{2}.$ (S3)

This illustrates that the connected correlations
$G^{(2)}_{c}({\bm{k}},-{\bm{k}})_{\rm Bogo}$ directly reflect the pairing at
opposite momenta, $\langle a^{\dagger}_{\bm{k}}a^{\dagger}_{-\bm{k}}\rangle$,
under the assumptions of Bogoliubov theory.

##  IV Correlation functions in momentum space and their scaling behavior

###  IV.1 Definitions of the correlation functions

In this work, we consider the following connected correlations of the
(momentum) density operator $N(\bm{k})$ (defined as
$N(\bm{k})=a^{\dagger}_{\bm{k}}a_{\bm{k}}$),

$\displaystyle G_{c}^{(2)}(\bm{k}_{1},\bm{k}_{2})$ $\displaystyle=\langle
N(\bm{k}_{1})N(\bm{k}_{2})\rangle-\langle N(\bm{k}_{1})\rangle\langle
N(\bm{k}_{2})\rangle$ $\displaystyle
G_{c}^{(3)}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})$ $\displaystyle=\langle
N(\bm{k}_{1})N(\bm{k}_{2})N(\bm{k}_{3})\rangle$ $\displaystyle-\langle
N(\bm{k}_{1})\rangle\langle N(\bm{k}_{2})N(\bm{k}_{3})\rangle$
$\displaystyle-\langle N(\bm{k}_{2})\rangle\langle
N(\bm{k}_{1})N(\bm{k}_{3})\rangle$ $\displaystyle-\langle
N(\bm{k}_{3})\rangle\langle N(\bm{k}_{1})N(\bm{k}_{2})\rangle$
$\displaystyle+2\langle N(\bm{k}_{1})\rangle\langle
N(\bm{k}_{2})\rangle\langle N(\bm{k}_{3})\rangle$ (S4)

and their integrals over zero-sum (quasi-)momentum modes

$\displaystyle\overline{G}_{c}^{(2)}(0)=$
$\displaystyle\sum_{\bm{k}\in\Omega_{k}}G_{c}^{(2)}(\bm{k},-\bm{k})$ (S5)
$\displaystyle\overline{G}_{c}^{(3)}(0)=$
$\displaystyle\sum_{\bm{k}_{1},\bm{k}_{2}\in\Omega_{k}}G_{c}^{(3)}(\bm{k}_{1},\bm{k}_{2},-\bm{k}_{1}-\bm{k}_{2})$
(S6)

where the sum in Eq. (S5) is defined over (quasi-) momenta $\bm{k}$ contained
in $\Omega_{k}$, such that ${\bm{k}}+{\bm{Q}}\neq{-\bm{k}}$, where ${\bm{Q}}$
is a reciprocal lattice wavevector; while the sum in Eq. (S6) is defined over
inequivalent pairs of (quasi-) momenta $\bm{k}_{1},\bm{k}_{2}$, namely
${\bm{k}}_{1}\neq{\bm{k}}_{2}+\bm{Q}$, both contained in $\Omega_{k}$ and such
that $\bm{k}_{3}=-\bm{k}_{1}-\bm{k}_{2}$ is also contained in $\Omega_{k}$.

###  IV.2 Comparison between theory and experiment

We define the total particle number in the region $\Omega_{k}$ as $\langle
N_{\Omega_{k}}\rangle=\sum_{\bm{k}\in\Omega_{k}}\langle N(\bm{k})\rangle$, and
the number of modes in the region of interest as
$V_{\Omega_{k}}=\sum_{\bm{k}\in\Omega_{k}}1$. We expect that
$\overline{G}_{c}^{(2)}(0)$ scales as the square of the density of particles
in $\Omega_{k}$, $\rho_{\Omega_{k}}=\langle
N_{\Omega_{k}}\rangle/V_{\Omega_{k}}$, and linearly in the number of modes
$\Omega_{k}$, namely $\overline{G}_{c}^{(2)}(0)\sim\langle
N_{\Omega_{k}}\rangle^{2}/V_{\Omega_{k}}$. Similarly we expect
$\overline{G}_{c}^{(3)}(0)$ to scale as the cube of the density and as the
square of the number of modes, so that $\overline{G}_{c}^{(3)}(0)\sim\langle
N_{\Omega_{k}}\rangle^{3}/V_{\Omega_{k}}$. As a consequence a meaningful
comparison between theory and experiment referring to different particle and
mode numbers requires a renormalization of the theoretical
$\overline{G}_{c}^{(2)}(0)$ value by a factor $\sim(\langle
N_{\Omega_{k}}\rangle_{\rm exp}/\langle N_{\Omega_{k}}\rangle_{\rm
th})^{2}(V_{\Omega_{k},\rm th}/V_{\Omega_{k},\rm exp})$; and of the
theoretical $\overline{G}_{c}^{(3)}(0)$ value by a factor $\sim(\langle
N_{\Omega_{k}}\rangle_{\rm exp}/\langle N_{\Omega_{k}}\rangle_{\rm
th})^{3}(V_{\Omega_{k},\rm th}/V_{\Omega_{k},\rm exp})$. The ratio
$r_{V}=V_{\Omega_{k},\rm th}/V_{\Omega_{k},\rm exp}$ is purely geometrical and
independent of the Hamiltonian parameters. Yet estimating it unambiguously is
not obvious. Only the core of the atomic cloud is expected to be superfluid at
the temperature of the experiment and to contribute to the measured momentum-
space correlations; and it is this core which in principle dictates the size
of the momentum modes in momentum space. On the other hand, the ratio
$r_{N}=\langle N_{\Omega_{k}}\rangle_{\rm exp}/\langle
N_{\Omega_{k}}\rangle_{\rm th}$ depends significantly on the interaction
strength and can be unambiguously determined. We have therefore decided to
rescale the theoretical data for $\overline{G}_{c}^{(2)}(0)$ and
$\overline{G}_{c}^{(3)}(0)$ by the factors $r_{N}^{2}$ and $r_{N}^{3}$,
respectively, assuming that $r_{V}\sim O(1)$. Indeed we observe that the
rescaled theoretical curve $r_{N}^{2}\overline{G}_{c}^{(2)}(0)$ matches the
experimental one in order of magnitude.

Figure S3: a $\overline{G}^{(2)}(0)$ integrated correlations for the quantum-
rotor model with size $V=10^{3}$ for the $\Omega_{k}$ region defined as in the
main text, and for a temperature $T_{J}=T/(2Jn_{\rm QR})=1.0$. The solid line
indicates the deep normal-gas limit. b $\overline{G}^{(3)}(0)$ integrated
correlations for the quantum-rotor model with size $V=6^{3}$ for an
$\Omega_{k}$ region corresponding to a cubic corona contained within two cubes
with sides $k=2k_{d}\pi/3$ and $k=k_{d}$, and at a temperature
$T_{J}=T/(2Jn_{\rm QR})=1.0$. The solid line indicates again the deep normal-
gas limit.

###  IV.3 Scaling of the correlation functions in the superfluid vs. in the
Mott (or normal gas) limit

We discuss here the different scaling regimes of the momentum correlations
predicted by the theory on quantum rotors. We restrict our attention to a
system in the _canonical_ ensemble at fixed particle number $N$, and at fixed
density in the lattice $n_{\rm QR}=N/V$ where $V$ is the system volume, so
that particle number and volume are proportional. In the superfluid regime we
expect the two-mode correlations $G_{c}^{(2)}(\bm{k},-\bm{k})$ to be $\sim
O(1)$, as predicted by Bogoliubov theory, and therefore we expect
$\overline{G}_{c}^{(2)}(0)\sim O(N)$ since $V_{\Omega_{k}}\sim V$. This
scaling is indeed observed in our quantum-rotor calculations. On the other
hand we do not know of a simple theoretical prediction for
$G_{c}^{(3)}(\bm{k}_{1},\bm{k}_{2},-\bm{k}_{1}-\bm{k}_{2})$ – indeed
Bogoliubov theory predicts this type of correlation to be vanishing. Hence we
resort to numerics in order to estimate the scaling of
$\overline{G}_{c}^{(3)}(0)$, and we find from our quantum-rotor calculations
that $\overline{G}_{c}^{(3)}(0)\sim O(N)$. This implies that three-mode
correlations are much more tenuous than two-mode ones, namely
$G_{c}^{(3)}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})\sim O(N^{-1})$ so that their
integral over $O(N^{2})$ trios of modes is $O(N)$.

The scaling of the integrated correlations in the superfluid regime should be
contrasted with that of a deep Mott-insulating or normal-gas regime, in which
all quasi-momentum modes become statistically independent, modulo the
conservation of the particle number. The population of each mode is Poisson-
distributed, with average occupation $\langle N(\bm{k})\rangle=n$ where $n$ is
the lattice filling; moreover one has $\langle N^{2}(\bm{k})\rangle=n(n+1)$
and $\langle N^{3}(\bm{k})\rangle=n^{3}+3n^{2}+n)$, according to Poisson
statistics. A lengthy but straightforward calculation leads to the result

$\displaystyle G_{c}^{(2)}(\bm{k}_{1},\bm{k}_{2})$
$\displaystyle\to-\frac{n}{V-1}$ $\displaystyle
G_{c}^{(3)}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})$
$\displaystyle\to\frac{2n}{(V-1)(V-2)}$ (S7)

(for inequivalent quasi-momenta $\bm{k}_{1}$, $\bm{k}_{2}$ and $\bm{k}_{3}$).
This implies that

$\displaystyle\overline{G}_{c}^{(2)}(0)$
$\displaystyle\to-\frac{nV_{\Omega_{k}}}{V-1}$
$\displaystyle\overline{G}_{c}^{(3)}(0)$ $\displaystyle\to\frac{2n{\cal
N}_{3,\Omega_{k}}}{(V-1)(V-2)}$ (S8)

where ${\cal N}_{3,\Omega_{k}}$ is the number of momentum trios $\bm{k}_{1}$,
$\bm{k}_{2}$ and $-\bm{k}_{1}-\bm{k}_{2}$ contained in $\Omega_{k}$. Since
$V_{\Omega_{k}}\sim V$ and ${\cal N}_{3,\Omega_{k}}\sim V^{2}$, both
quantities take limiting values which do _not_ scale with the particle number,
marking a strong suppression of correlations with respect to the superfluid
regime.

Fig. S3 shows results for the $\overline{G}_{c}^{(2)}(0)$ and
$\overline{G}_{c}^{(3)}(0)$ correlations of the quantum-rotor model without
any further renormalization factor. We observe that the asymptotic behavior
predicted in the deep normal-gas limit at large interaction strength is
correctly reproduced by our results. Interestingly, the integrated
$\overline{G}_{c}^{(3)}(0)$ does not reach this limit monotonically from
above, but it rather shows an intermediate minimum, corresponding roughly to
the superfluid/normal-gas transition.

##  V Finite-temperature Bogoliubov theory

The results from Bogoliubov theory presented in the main text have been
obtained by a simple self-consistent generalization of the theory to access
finite temperature. Bogoliubov theory is based on the following approximate
grand-canonical Hamiltonian on a 3D cubic lattice

$\displaystyle H-\mu N$ $\displaystyle\approx$
$\displaystyle\frac{UN_{0}^{2}}{2V}-\mu N_{0}+\sum_{\bm{k}\neq
0}(e_{\bm{k}}+2Un_{0}-\mu)a^{\dagger}_{\bm{k}}a_{\bm{k}}$ (S9)
$\displaystyle+$ $\displaystyle\frac{Un_{0}}{2}\sum_{\bm{k}\neq
0}\left(a_{\bm{k}}a_{-\bm{k}}+{\rm h.c.}\right)$

where $e_{\bm{k}}=-2J\left[\cos(k_{x})+\cos(k_{y})+\cos(k_{z})\right]$, and
$n_{0}=N_{0}/V$ is the condensate density. The chemical potential $\mu$ is
fixed by requiring the condensate to satisfy the Gross-Pitaevskii equation,
which leads to the condition $\mu=Un_{0}$.

The Bogoliubov Hamiltonian is then diagonalized by the transformation
$b_{\bm{k}}=u_{\bm{k}}a_{\bm{k}}+v_{\bm{k}}a^{\dagger}_{-\bm{k}}$ to the form

$H-\mu N\approx\sum_{\bm{k}\neq
0}\epsilon_{\bm{k}}b_{\bm{k}}^{\dagger}b_{\bm{k}}+{\rm const.}$ (S10)

where
$\epsilon_{\bm{k}}=\sqrt{A_{\bm{k}}^{2}-B_{\bm{k}}^{2}}=\sqrt{e_{\bm{k}}(e_{\bm{k}}+2Un_{0})}$
with $A_{\bm{k}}=e_{\bm{k}}+Un_{0}$, $B_{\bm{k}}=Un_{0}$,
$u_{\bm{k}}=\sqrt{\frac{1}{2}(A_{\bm{k}}/E_{\bm{k}}+1)}$, and
$v_{\bm{k}}=\sqrt{\frac{1}{2}(A_{\bm{k}}/E_{\bm{k}}-1)}$. The population of
finite-momentum modes at finite temperature reads then

$\langle
a_{\bm{k}}^{\dagger}a_{\bm{k}}\rangle=v_{\bm{k}}^{2}+n_{B}(\epsilon_{\bm{k}},T)(u_{\bm{k}}^{2}+v_{\bm{k}}^{2})$
(S11)

where $n_{B}(\epsilon_{\bm{k}},T)=(e^{\epsilon_{\bm{k}}/(k_{B}T)}-1)^{-1}$ is
the Bose distribution for Bogoliubov quasiparticles. This expression depends
parametrically on the condensate density $n_{0}$ which enters into the
dispersion relation $\epsilon_{\bm{k}}$ as well as in the $u_{\bm{k}}$ and
$v_{\bm{k}}$ coefficients. Fixing the total density at temperature $T$

$n=n_{0}+\frac{1}{V}\sum_{\bm{k}\neq 0}\langle
a_{\bm{k}}^{\dagger}a_{\bm{k}}\rangle$ (S12)

amounts then to imposing a non-linear equation on the condensate density
$n_{0}$ (or equivalently on the chemical potential), which can be solved
numerically. In the example given in the main text we have chosen $n=1$.

Once the condensate density satisfying Eq. (S12) is found, one can evaluate
the pairing correlation

$\langle
a^{\dagger}_{\bm{k}}a^{\dagger}_{-\bm{k}}\rangle=-u_{\bm{k}}v_{\bm{k}}\left[1+2n_{B}(\epsilon_{\bm{k}},T)\right]$
(S13)

leading to the connected correlations of Eq. (S3), shown in the main text.

##  VI Quantum-rotor model and Quantum Monte Carlo

The quantum-rotor model is an approximation to the Bose-Hubbard model valid at
large, integer filling $n_{\rm QR}$ Wallin  _et al._ (1994). Its starting
point is the amplitude-phase decomposition $a_{i}=e^{i\phi}\sqrt{n_{i}}$ of
the bosonic field, with $[\phi,n_{i}]=1$. Rewriting the density as
$n_{i}=n_{\rm QR}+\delta n_{i}$, and assuming that $\sqrt{\langle(\delta
n_{i})^{2}\rangle}\ll n_{\rm QR}$ (since $n_{\rm QR}\gg 1$) in all physical
regimes of interest, then one can write

$a_{i}^{\dagger}a_{j}+{\rm h.c.}\approx 2n_{\rm QR}\cos(\phi_{i}-\phi_{j})$
(S14)

neglecting density fluctuations. Density-density interactions are nonetheless
introducing quantum fluctuations of the phase, since

$n_{i}^{2}=-\frac{\partial^{2}}{\partial\phi_{i}^{2}}~{}.$ (S15)

Working at integer filling eliminates all linear terms in the density from the
Hamiltonian, leading to the quantum rotor Hamiltonian

$H_{\rm QR}=-2Jn_{\rm QR}\sum_{\langle
ij\rangle}\cos(\phi_{i}-\phi_{j})-\frac{U}{2}\sum_{i}\frac{\partial^{2}}{\partial\phi_{i}^{2}}~{}.$
(S16)

Using a coherent-state path-integral approach Wallin  _et al._ (1994); Sondhi
_et al._ (1997), the quantum statistical physics of the quantum-rotor
Hamiltonian in $d$ dimensions can be mapped onto an effective classical
problem of rotors i.e. planar spins, or angles $\phi_{i,p}$, where $i$ is the
spatial position and $p$ is the position in the extra (Trotter) dimension. A
quantum Monte Carlo simulation of the quantum rotor model amounts to sampling
the configurations of this $(d+1)$-dimensional classical model Wallin  _et
al._ (1994); Roscilde  _et al._ (2016). Given the approximation of Eq. (S14),
to each configuration of the classical rotors $\\{\phi_{i,p}\\}$ at a given
Trotter step $p$ one can uniquely associate a momentum distribution

$N_{p}(\bm{k})=\frac{n_{\rm
QR}}{V}\sum_{ij}e^{i\bm{k}\cdot\left({\bm{r}}_{i}-{\bm{r}}_{j}\right)}e^{-i\left(\phi_{i,p}-\phi_{j,p}\right)}~{}.$
(S17)

Hence sampling the partition function of the classical rotors allows for the
reconstruction of arbitrary correlation functions of the momentum populations
$N(\bm{k})$. This aspect is at the core of our use of the quantum rotor model
to probe the fluctuations of the momentum distribution.

The original Bose-Hubbard model of interest can also be studied in its
equilibrium thermodynamics via different quantum Monte Carlo approaches
Trotzky  _et al._ (2010); Wessel  _et al._ (2004). Yet all these approaches
use Fock states in real space as computational basis, and therefore the
momentum distribution is an off-diagonal observable. Estimating off-diagonal
observables requires a dedicated update scheme (a so-called ”worm algorithm”).
Complex worm-like algorithms have been developed in the past by some of us to
probe the second moment of the fluctuations of the momentum distribution in
the one-dimensional Bose-Hubbard model Roscilde (2008); Fang  _et al._ (2016),
but third moments have never been probed to the best of our knowledge. An
important direction for the future will then be to develop such dedicated
algorithms for the study of the fluctuations of the momentum distribution of
the three-dimensional Hubbard model.