# Absence of sizable superconductivity in hydrogen boride: A first principles
study

Antonella Meninno1,2<EMAIL_ADDRESS>Ion Errea1,2,3 1Centro de
Física de Materiales (CSIC-UPV/EHU), Manuel de Lardizabal Pasealekua 5, 20018
Donostia/San Sebastián, Spain  2Fisika Aplikatua 1 Saila, Gipuzkoako
Ingeniaritza Eskola, University of the Basque Country (UPV/EHU), Europa Plaza
1, 20018 Donostia/San Sebastián, Spain  3Donostia International Physics Center
(DIPC), Manuel de Lardizabal Pasealekua 4, 20018 Donostia/San Sebastián, Spain

###### Abstract

The recently synthesized hydrogen boride monolayer in the $Cmmm$ phase is a
promising superconductor due to its similarity to MgB2 and the large hydrogen
content in its structure. Making use of first-principles calculations based on
density functional theory, we study its electronic, vibrational, and
superconducting properties and conclude that, despite the expectations,
hydrogen boride does not have a sizable superconducting critical temperature.
The presence of hydrogen in the system alters the boron-boron bonding,
weakening the electron-phonon interaction. We have studied the effect of
enhancing the critical temperature by doping the system, but the inclusion of
electrons or holes reveals ineffective. We attribute the small critical
temperature of this system to the vanishing hydrogen character of the states
at the Fermi level, which are dominated by boron $p$ states. Our results
determine that a large proportion of hydrogen-like states are needed at the
Fermi level to attain a large superconducting critical temperature in
hydrogenated monolayers.

##  I Introduction

The observation of high temperature superconductivity in hydrogen-based
superconductors, with critical temperatures ($T_{c}$’s) surpassing 200 K, is
one of the most astonishing results in physics of the last years. Some
examples of such “superhydrides” are H3S [1], with a maximum $T_{c}$ of 203 K
at 155 GPa; LaH10 [2, 3], with $T_{c}=250$ K at 150 GPa; YH9 [4, 5], reaching
a $T_{c}$ of around 250 K at approximately 200 GPa; and YH6 [5, 6], with
$T_{c}=224$ K at 166 GPa.

The actual challenge is to reduce the pressure at which superconductivity
occurs in superhydrides, with the ultimate goal of understanding if
high-$T_{c}$ superconductivity is possible at ambient pressure for these type
of compounds. Interestingly, the enhancement of superconductivity by hydrogen
absorption is a well known phenomenon, observed for example in palladium and
thorium hydrides at ambient pressure [7, 8]. In both cases, the
superconducting critical temperature of the original element (Pd or Th)
increases from around 1 K up to approximately 10 K in the stoichiometric
hydride. The role of hydrogen in enhancing the $T_{c}$ has been further
studied in the context of monolayers at ambient pressures. In previous
theoretical works based on first principles calculations, it has been shown
that a MgB2 monolayer, which has a critical temperature of 39 K in the bulk
[9], can reach a $T_{c}$ of 69 K after hydrogenation [10], and that $p$-doped
graphane, a hydrogenated graphene monolayer, reaches a critical temperature of
around $80$ K even if graphene itself is not superconducting [11].
Hydrogenated monolayers offer thus a promising platform to synthesize
high-$T_{c}$ hydrogen-based compounds at ambient pressures, overcoming the
limits of the high-pressure superhydrides.

The recently synthesized stoichiometric hydrogen boride monolayer (HB) [12] is
a promising superconductor [13]. HB was obtained through a procedure of wet
chemical exfoliation at ambient pressure from MgB2 through ion-exchange
treatment, and it is believed to adopt a quasi-hexagonal structure with a
boron honeycomb layer reminiscent of the structure of MgB2, with H atoms out
of plane (see Fig. 2) [12, 14, 13]. The structure of HB observed
experimentally seems to adopt a $Cmmm$ space group as anticipated by _ab
initio_ crystal structure prediction calculations [15]. More recently,
different polymorphs of hydrogen and boron have been synthesized in the two-
dimensional limit [16], showing the potential of “borophane” to form different
stable compounds at the nanoscale and broaden its applications in nanodevices.

Despite its semimetallic character [13], $Cmmm$ HB is a good superconducting
candidate for several reasons. Mainly because several parent “borophene”
monolayers have been predicted to have $T_{c}$’s of around 20 K [17] and
hydrogenation, as mentioned above, is expected to enhance it. Also, because
the bonding nature of boron in $Cmmm$ HB is reminiscent of the $sp^{2}$
bonding in MgB2, which is crucial for its large $T_{c}$, even if hydrogen has
a large impact on the band structure of the compound [13]. The goal of the
present work is to determine the potential superconductivity of $Cmmm$ HB
through first principles calculations based on density-functional theory
(DFT). We have also analyzed the possibility of enhancing the critical
temperature of the system by doping.

The paper is organized as follows. In section 3 we overview the theoretical
framework of the calculations performed, in section III we present the results
of our calculations, in section IV we analyze and interpret the results
further, and in section V we summarize the results.

Figure 1: Different views of the crystal structure of the hydrogen boride (HB)
monolayer in the $Cmmm$ phase. The lattice parameters are indicated, together
with the $\gamma$ angle

Figure 2: Electron localization function (ELF) for the non doped HB monolayer.
In the top panel the $z=0$ plane is represented, while in the bottom panel the
$y=0$ plane (putting the origin of coordinates at the position where the
lattice vectors in Fig. 2 start).

##  II Computational details

Figure 3: Electronic band structure (left panel) and density of states (DOS)
(right panel) for $Cmmm$ hydrogen boride. The bands and DOS are also plotted
with electron and hole doping. The projection of the DOS onto atomic orbitals
is also included.

All computations have been performed using DFT within the Quantum ESPRESSO
package [18, 19], making use of the Perdew-Burke-Ernzerhof (PBE)
parametrization of the exchange-correlation functional [20]. The electron-ion
interaction has been modeled making use of ultrasoft pseudopotentials,
including $2s^{2}$ and $2p^{1}$ electrons of boron in the valence. The self-
consistent DFT computation has been performed with a grid of
42$\times$42$\times$1 k points for the integrals over the Brillouin zone and a
Methfessel-Paxton first-order spreading smearing of 0.02 Ry [21]. A kinetic
energy cutoff of $50$ Ry has been used for the plane-wave basis and a $500$ Ry
cutoff for the charge density. In order to avoid spurious interactions, we
have left a vacuum between monolayer replicas of 35 $a_{0}$. Phonon
frequencies have been computed in the harmonic approximation within density-
functional perturbation theory (DFPT) [22, 23]. Dynamical matrices have been
calculated explicitly in an 8$\times$8$\times$1 grid of q points, and we have
used Fourier interpolation to obtain the phonon spectrum.

The electron-phonon interaction has been computed within DFPT. The electron-
phonon contribution to the phonon linewidth [24] can be approximated at low
temperatures as

$\displaystyle\gamma_{\nu}({\bf q})$ $\displaystyle=$
$\displaystyle\frac{2\pi\omega_{\nu}({\bf q})}{N_{{\bf
k}}}\sum_{nn^{\prime}{\bf k}}|g^{\nu}_{n^{\prime}{\bf k+q},n\bf k}|^{2}$ (1)
$\displaystyle\times$ $\displaystyle\delta(\epsilon_{n^{\prime}{\bf k+\bf
q}}-\epsilon_{F})\delta(\epsilon_{n\bf k}-\epsilon_{F}),$

where $\epsilon_{n\bf k}$ is the energy of an electron in the band $n$ and
wave number $\bf k$, $\omega_{\nu}({\bf q})$ the frequency of a phonon in the
mode $\nu$ and wave number $\bf q$, $\epsilon_{F}$ is the Fermi energy,
$N_{{\bf k}}$ the number of ${\bf k}$ points in the sum, and
$g^{\nu}_{n^{\prime}{\bf k+q},n\bf k}$ are the electron-phonon matrix elements
associated to the scattering of electrons with band energies $\epsilon_{n\bf
k}$ and $\epsilon_{n^{\prime}{\bf k}+{\bf q}}$ with a phonon of frequency
$\omega_{\nu}({\bf q})$. Superconducting properties can be then calculated
with the Eliashberg function

$\alpha^{2}F(\omega)=\frac{1}{2\pi N(\epsilon_{F})N_{{\bf q}}}\sum_{\nu{\bf
q}}\frac{\gamma_{\nu}({\bf q})}{\omega_{\nu}({\bf
q})}\delta(\omega-\omega_{\nu}({\bf q})),$ (2)

where $N(\epsilon_{F})$ is the density of states at the Fermi energy and
$N_{{\bf q}}$ the number of $\bf q$ points in the sum. For instance, the
electron-phonon coupling constant $\lambda$ is the $\omega\to\infty$ limit of
the so-called integrated electron-phonon coupling constant

$\displaystyle\lambda(\omega)=2\int_{0}^{\omega}\frac{\alpha^{2}F(\omega^{\prime})}{\omega^{\prime}}\mathrm{d}\omega^{\prime}.$
(3)

We have estimated $T_{c}$ by solving the Allen-Dynes equations [25, 26, 27],
choosing the effective Coulomb potential $\mu^{*}$ between $0.08$ and $0.12$.

The calculation of the electron-phonon properties has been done by making use
in Eq. (1) of a 60$\times$60$\times$1 grid of $\bf k$ points and an
8$\times$8$\times$1 grid of q points in Eq. (2). To approximate the delta
functions in Eq. (1), we have used Gaussians with a broadening of $0.002$ Ry.
The $\bf k$-point grid yielded converged $\lambda$ values even with such a
small broadening of the Delta funcitons.

We have also studied the possibility of enhancing superconductivity by doping
the system. We have doped the system by directly inserting or removing
electrons and solving the Kohn-Sham equations, or including an external
electric field in the calculations as it corresponds to a field effect
transistor (FET) doping setup [28, 29]. After doping, the crystal structure
has been relaxed again until the forces and the strain tensor vanish. We have
repeated the procedure subtracting electrons from the system. We have
recomputed the electron and phonon bands, as well as the electron-phonon
coupling with the inclusion of extra electrons or holes. The inclusion of the
electric field in the calculations gave electronic band structures in
agreement with those obtained with direct inclusion of electrons, validating
the simpler approach. Therefore, here we report exclusively the results
obtained by direct addition and subtraction of electrons.

##  III Results

###  a Non doped case

####  a.1 Structure and bonding

We describe the centered rectangular lattice of $Cmmm$ hydrogen boride in the
primitive cell as shown in Fig. 2, which includes two boron and two hydrogen
atoms. In this description the two basis vectors $\bf a$ and $\bf b$ have the
same length, which according to our DFT relaxations has a value of
$a=b=5.739a_{0}$. The angle between them is $\gamma=59.33^{\circ}$. All boron
atoms form a quasihexagonal honeycomb-like lattice and the hydrogen atoms sit
at a distance $d$ from the boron plane. More precisely the boron atoms are at
$\alpha(\boldsymbol{a}+\boldsymbol{b})$ and
$(1-\alpha)(\boldsymbol{a}+\boldsymbol{b})$, while the hydrogen atoms sit at
$1/2(\boldsymbol{a}+\boldsymbol{b})\pm d\hat{\boldsymbol{k}}$, where
$\hat{\boldsymbol{k}}$ is the unitary vector in the out-of-plane direction.
According to our calculations $\alpha=0.328$ and $d=1.828a_{0}$. All the
parameters that determine the structure are summarized in Table 1. The
obtained structural parameters are in agreement with previous calculations
[14, 13].

As illustrated by the calculated electron localization function (ELF)
presented in Fig. 2, the presence of hydrogen affects the B-B bonding
significantly. In a system like MgB2, where the B layer is perfectly
hexagonal, all boron atoms are covalently bonded at a distance of 3.34 $a_{0}$
[13]. The presence of hydrogen, however, changes the bonding pattern by
shortening the B-B covalent that does not have a hydrogen atom above and
below, while increasing the B-B length in the other bond due to the creation
of a bridging B-H-B bond [30]. In fact, the presence of hydrogen destroys the
covalent bond between the boron atoms by creating this bridging bonding. The
B-B distances in the covalent bond are 3.25 $a_{0}$ while 3.44 $a_{0}$ between
those affected by the B-H-B bonding.

|  $a$ $(a_{0})$ | $\alpha$ | $\gamma$ |  $d$ $(a_{0})$
---|---|---|---|---
Non doped | 5.739 | 0.328 | 59.33∘ | 1.828
Doped electrons | 5.739 | 0.326 | 59.36∘ | 1.810
Doped holes | 5.739 | 0.330 | 59.27∘ | 1.838
Table 1: Structural parameters of $Cmmm$ HB in the non doped and doped cases.

####  a.2 Electronic bands and Fermi surface

The band structure and the density of states (DOS) projected onto atomic
orbitals are plotted in Fig. 3. The bands, which are in agreement with
previous DFT and tight-binding results [13, 14], show a semimetallic
character, with a hole pocket at $\Gamma$ and an electron pocket at Y, with
two distinct bands crossing the Fermi level. The character of these two bands
is mainly associated to the boron $p$ orbitals. The hole pocket at $\Gamma$ is
coming from in-plane $p_{x}$ and $p_{y}$ orbitals, while the $Y$ pocket from
boron $p_{z}$ [13].

The contribution of hydrogen to these bands that cross the Fermi level is
scarce. Looking at the DOS projected onto atomic orbitals, we can see that
hydrogen contributes to the DOS at the Fermi level only a 1.6%. The DOS at the
Femi level is thus dominated by boron $p$ orbitals. Bands above approximately
2 eV above the Fermi level at $\Gamma$ do have, on the contrary, some
considerable H character.

The hole pocket centered at $\Gamma$ and the electron pocket centered at Y are
evident in the Fermi surface of Fig. 4. Both pockets are ellipsoidal, but, due
to the semimetallic character of HB, the area enclosed by the pockets is
small. Indeed, the system is not far from an insulating state, as it has been
suggested that a gap can be opened with a small strain [31].

####  a.3 Phonons and electron-phonon coupling

Figure 4: Fermi surfaces for the non doped system (left), the one doped with
electrons (center), and the one with holes (right) for $Cmmm$ HB. High-
symmetry points are marked in green. Reciprocal lattice vectors ${\bf a}^{*}$
and ${\bf b}^{*}$ are included as well. The path used for the calculation of
the electronic bands and phonon spectra is illustrated.

We have calculated the phonon bands and the phonon DOS (PDOS), as well as the
projection of the PDOS onto boron and hydrogen atoms (see Fig. 5). As expected
for a 2D material [32], the flexural out-of-plane acoustic mode has a
quadratic dispersion due to rotational symmetry. The phonon spectrum has three
distinctive regions. Below $\sim$700 cm-1 there are three acoustic and another
three optical modes. These modes show, in general, a mixed character between
hydrogen and boron. The hydrogen character of these low-energy modes comes
from the vibrations of H atoms along the $y$ direction. At higher energies,
between approximately 750 cm-1 and 1100 cm-1, there are two isolated phonon
bands that describe in-plane boron vibrations, reminiscent of the $E_{2g}$
modes that are strongly coupled to the electrons in the similar MgB2 [33, 34,
35]. For comparison, the energy of the $E_{2g}$ mode in MgB2 is around 600
cm-1 [36], smaller than the analogous modes in HB. The remaining four phonon
modes appear after a large energy gap above 1500 cm-1 and have a dominant
hydrogen character, with vibrations along the $x$ and $z$ directions.

As we can observe in Fig. 6, where the Eliashberg function
$\alpha^{2}F(\omega)$ and the integrated electron-phonon coupling constant
$\lambda(\omega)$ are shown, phonon modes with energies between $500$ and
$1000$ cm-1 contribute more to the electron-phonon coupling. The modes mainly
contributing to the peaks in this energy range involve boron and hydrogen
displacements, as well as the in-plane pure boron modes, analogous to the
$E_{2g}$ of MgB2. There is another significant peak at around $1500$ cm-1,
related to hydrogen-character modes, but, as expected, the lower frequency
modes are the ones that contribute most to the electron-phonon coupling.
Unfortunately, the integrated electron-phonon coupling $\lambda$ reaches a
very low value of approximately $0.2$, which yields a very low superconducting
$T_{c}$ of only 11 mK.

Figure 5: Phonon spectrum (left panel) for the non doped $Cmmm$ HB, together
with the doped cases with electrons and holes. The phonon density of states
(PDOS) (right panel) is plotted, including also its decomposition onto
different atoms, only for the non doped case. The PDOS decomposition is
similar in the doped cases.

###  b Doped cases

The DOS of 2D materials is energy independent in the free electron-like limit
and, thus, since $\lambda$ is proportional to $N(\epsilon_{F})$, doping is
expected to not affect the electron-phonon coupling constant. However, there
are examples in the literature where electrostatic doping has been used to
induce superconductivity in two-dimensional materials such as amorphous
bismuth [37] or the LaIO3/SrTiO3 interface [38], in analogy to the common
behavior of the cuprates [39]. Also, this FET doping mechanism has been used
to induce and tune supercconductivity in transition metal dichalcogenides
[40], for instance in MoS2, where doping induces first an insulator to metal
transition, and a superconducting state at further electron doping, with a
$T_{c}$ that can reach a value as high as 10 K [41].

Motivated by these experimental results and the possibility of FET doping two-
dimensional materials in the laboratory, we have studied the effect of both
electron and hole doping in the structural, electronic, phononic, and
superconducting properties of $Cmmm$ HB, with the hope that its low $T_{c}$
may be enhanced. In the following we report the results obtained with electron
and hole doping, which will be referred as the electron/holes doped cases
respectively. We dope the system with $0.075$ electrons(holes) per unit cell,
which corresponds to a $9.39\times 10^{13}$ electrons(holes)/cm2.

####  b.1 Structure and bonding

By relaxing the structure before with DFT in the doped system, we can
determine how the structure and the bonding are affected when extra electrons
or holes are included in the system. The modified structural parameters are
summarized in Table 1. Doping leaves the lattice parameter $a$ unvaried, but
increases the $\gamma$ angle when electrons are added and decreases it if,
instead, extra holes are present. Similarly, the distance of hydrogen from the
plane decreases in the electron doped case, while increases in the hole doped
case. These results suggest that electron doping strengthens the B-B and B-H-B
covalent bonds.

####  b.2 Electronic bands and Fermi surface

We report the electronic bands, together with the DOS, in Fig. 3. We observe
that, as expected, the conduction bands of the electron doped case tend to go
down with respect to the non doped case, while the bands of the holes doping
tend to go up with respect to the non doped case. The effect is that in the
electron doped case the hole pocket at $\Gamma$ decreases and the electron
pocket at Y increases. The opposite happens in the hole doped case. This
effect is clear in the Fermi surface in Fig. 4. The character of the bands is
barely affected by the doping and the percentage of H states in the DOS at the
Fermi energy remains minimal, as in the non doped case. Indeed, with the
doping values assumed, $N(\epsilon_{F})$ remains practically unchanged, which
is reasonable given the constant DOS in the vicinity of the Fermi level.

####  b.3 Phonons and electron-phonon coupling

As we can see in Fig. 5, electron or hole doping barely affects the lowest six
phonon modes. The differences become more appreciable for the isolated two
boron in-plane bands in the frequency range between 750 and 1100 cm-1.
Electron doping increases the splitting between these two modes, while whole
doping decreases it. Finally, the hydrogen-dominated modes at high energies
are softened in the electron doped case and hardened in the hole doped case.

Taking a look at the values of the Eliashberg function $\alpha^{2}F(\omega)$
and the integrated electron-phonon coupling constant $\lambda(\omega)$,
plotted in Fig. 6, we can see that in the doped cases still the modes in the
500-1000 cm-1 range are those with the most important contribution to
$\lambda$. However, specially the low-energy acoustic modes, are more coupled
to the electrons in the whole doped case. Also the modes in the 500-1000 cm-1
range seem to be more prone to the electrons in the doped cases. We attribute
these differences to the change in the electron-phonon coupling matrix
elements and not to the shift of the phonon frequencies. The slight changes in
the bonding pattern can explain a subtle change of the electron-phonon matrix
elements themselves. Nevertheless, doping, either hole or electron, is not
able to enhance $\lambda$ sufficiently to reach a sizable critical temperature
in HB. The value of $T_{c}$ in the hole doped case is 21 mK. Thus, despite the
opposite is the case in other systems, doping does not seem to be an efficient
strategy to increase the critical temperatures in hydrogen boride, at least,
in the $Cmmm$ phase synthesized experimentally.

Figure 6: Eliashberg function $\alpha^{2}F(\omega)$ and integrated electron-
phonon coupling constant $\lambda(\omega)$ for the non doped system, as well
as the doped systems with electrons and holes.

##  IV Discussion

Neither the non doped nor the doped HB present a sizable superconductivity in
the $Cmmm$ phase. This can be traced to a variety of reasons. It has been
recently argue [42] that hydrogen-based compounds can reach a high $T_{c}$
when an electronic bonding network between localized units is created. The
_networking value_ $\phi$, which was defined as the highest value of the ELF
that creates an isosurface spanning through the whole crystal in all three
Cartesian directions, measures precisely that network. Interestingly, $\phi$
correlates better than any other descriptor with $T_{c}$. By multiplying the
networking value with the fraction of hydrogen in the compound $H_{f}$ (0.5
for HB) and with the third root of the hydrogen fraction of the total DOS at
the Fermi energy $H_{\mathrm{DOS}}$, i.e. defining $\Phi_{\mathrm{DOS}}=\phi
H_{f}\sqrt[3]{H_{\mathrm{DOS}}}$, a rather good correlation was obtained
between $\Phi_{\mathrm{DOS}}$ and $T_{c}$, capable of estimating the critical
temperature with the empirical [42]

$T_{c}=(750\Phi_{\mathrm{DOS}}-85)\mathrm{K}$ (4)

equation . Here we assume that the networking value in a 2D system should be
determined by the highest value of the ELF that creates an isosurface that
expands in the whole crystal in the plane. With this adapted definition we
obtain $\phi=0.53$, a rather large value. However, the fact that the hydrogen
fraction of the DOS at the Fermi level is so low in this system suppresses all
the $T_{c}$ according to Eq. (4), both in the undoped and doped cases.

Furthermore, the presence of hydrogen alters the bonding of the boron layer
creating the B-H-B bridges. This fact suppresses the electron-phonon coupling
of the boron in-plane modes due to the large impact that the presence of
hydrogen has in the band structure [13, 14]. Therefore, the presence of
hydrogen also suppresses the electron-phonon coupling in the boron layer
itself. Hydrogen doping of borophene is not thus favoring the emergence of
high-temperature superconductivity in hydrogen boride, it has the opposite
effect at least in the $Cmmm$ phase.

##  V Conclusions

From the results of our analysis we can conclude that, despite the large
expectations in the literature [13], hydrogen boride is not a strongly coupled
superconductor mainly because of the weak hydrogen character of the density of
states at the Fermi level and the effect of the presence of hydrogen in the
bonding of the boron layer. The $T_{c}$ values are not sizable, below 1 K,
even if the system is doped. Our results suggest that hydrogenated monolayers
will be high-temperature superconductors, comparable in their $T_{c}$’s to
high-pressure superhydrides, only if the electronic states at the Fermi level
have a large hydrogen character.

## Acknowledgements

We have received funding from the European Research Council (ERC) under the
European Union’s Horizon 2020 research and innovation programme (grant
agreement No 802533).

## References

  * Drozdov  _et al._ [2015] A. Drozdov, M. Eremets, I. Troyan, V. Ksenofontov, and S. I. Shylin, Nature 525, 73 (2015).
  * Somayazulu  _et al._ [2019] M. Somayazulu, M. Ahart, A. K. Mishra, Z. M. Geballe, M. Baldini, Y. Meng, V. V. Struzhkin, and R. J. Hemley, Physical Review Letters 122, 027001 (2019).
  * Drozdov  _et al._ [2019] A. Drozdov, P. Kong, V. Minkov, S. Besedin, M. Kuzovnikov, S. Mozaffari, L. Balicas, F. Balakirev, D. Graf, V. Prakapenka, _et al._ , Nature 569, 528 (2019).
  * Snider  _et al._ [2021] E. Snider, N. Dasenbrock-Gammon, R. McBride, X. Wang, N. Meyers, K. V. Lawler, E. Zurek, A. Salamat, and R. P. Dias, Phys. Rev. Lett. 126, 117003 (2021).
  * Kong  _et al._ [2021] P. Kong, V. S. Minkov, M. A. Kuzovnikov, A. P. Drozdov, S. P. Besedin, S. Mozaffari, L. Balicas, F. F. Balakirev, V. B. Prakapenka, S. Chariton, D. A. Knyazev, E. Greenberg, and M. I. Eremets, Nature Communications 12, 5075 (2021).
  * Semenok  _et al._ [2020] D. V. Semenok, I. A. Troyan, A. G. Kvashnin, A. G. Ivanova, M. Hanfland, A. V. Sadakov, O. A. Sobolevskiy, K. S. Pervakov, A. G. Gavriliuk, I. S. Lyubutin, _et al._ , arXiv preprint arXiv:2012.04787 (2020).
  * Skoskiewicz [1972] T. Skoskiewicz, physica status solidi (a) 11, K123 (1972).
  * Satterthwaite and Toepke [1970] C. Satterthwaite and I. Toepke, Physical Review Letters 25, 741 (1970).
  * Nagamatsu  _et al._ [2001] J. Nagamatsu, N. Nakagawa, T. Muranaka, Y. Zenitani, and J. Akimitsu, Nature 410, 63 (2001).
  * Bekaert  _et al._ [2019] J. Bekaert, M. Petrov, A. Aperis, P. M. Oppeneer, and M. Milošević, Physical review letters 123, 077001 (2019).
  * Savini  _et al._ [2010] G. Savini, A. Ferrari, and F. Giustino, Physical review letters 105, 037002 (2010).
  * Nishino  _et al._ [2017] H. Nishino, T. Fujita, N. T. Cuong, S. Tominaka, M. Miyauchi, S. Iimura, A. Hirata, N. Umezawa, S. Okada, E. Nishibori, A. Fujino, T. Fujimori, S.-i. Ito, J. Nakamura, H. Hosono, and T. Kondo, Journal of the American Chemical Society 139, 13761 (2017), pMID: 28926230, https://doi.org/10.1021/jacs.7b06153 .
  * Tateishi  _et al._ [2019] I. Tateishi, N. Cuong, C. Moura, M. Cameau, R. Ishibiki, A. Fujino, S. Okada, A. Yamamoto, M. Araki, S. Ito, _et al._ , Physical Review Materials 3, 024004 (2019).
  * Kawamura  _et al._ [2019] R. Kawamura, N. T. Cuong, T. Fujita, R. Ishibiki, T. Hirabayashi, A. Yamaguchi, I. Matsuda, S. Okada, T. Kondo, and M. Miyauchi, Nature communications 10, 1 (2019).
  * Jiao  _et al._ [2016] Y. Jiao, F. Ma, J. Bell, A. Bilic, and A. Du, Angewandte Chemie International Edition 55, 10292 (2016), https://onlinelibrary.wiley.com/doi/pdf/10.1002/anie.201604369 .
  * Li  _et al._ [2021] Q. Li, V. S. C. Kolluru, M. S. Rahn, E. Schwenker, S. Li, R. G. Hennig, P. Darancet, M. K. Y. Chan, and M. C. Hersam, Science 371, 1143 (2021), https://www.science.org/doi/pdf/10.1126/science.abg1874 .
  * Gao  _et al._ [2017] M. Gao, Q.-Z. Li, X.-W. Yan, and J. Wang, Phys. Rev. B 95, 024505 (2017).
  * Giannozzi  _et al._ [2009] P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, _et al._ , Journal of physics: Condensed matter 21, 395502 (2009).
  * Giannozzi  _et al._ [2017] P. Giannozzi, O. Andreussi, T. Brumme, O. Bunau, M. B. Nardelli, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, M. Cococcioni, _et al._ , Journal of Physics: Condensed Matter 29, 465901 (2017).
  * Perdew  _et al._ [1996] J. P. Perdew, K. Burke, and M. Ernzerhof, Physical review letters 77, 3865 (1996).
  * Methfessel and Paxton [1989] M. Methfessel and A. T. Paxton, Phys. Rev. B 40, 3616 (1989).
  * Baroni  _et al._ [1987] S. Baroni, P. Giannozzi, and A. Testa, Physical review letters 58, 1861 (1987).
  * Baroni  _et al._ [2001] S. Baroni, S. de Gironcoli, A. Dal Corso, and P. Giannozzi, Rev. Mod. Phys. 73, 515 (2001).
  * Allen [1972] P. B. Allen, Physical Review B 6, 2577 (1972).
  * Eliashberg [1960] G. Eliashberg, Sov. Phys. JETP 11, 696 (1960).
  * Allen and Mitrović [1983] P. B. Allen and B. Mitrović, in  _Solid state physics_ , Vol. 37 (Elsevier, 1983) pp. 1–92.
  * Allen and Dynes [1975] P. B. Allen and R. Dynes, Physical Review B 12, 905 (1975).
  * Brumme  _et al._ [2014] T. Brumme, M. Calandra, and F. Mauri, Physical Review B 89, 245406 (2014).
  * Brumme  _et al._ [2015] T. Brumme, M. Calandra, and F. Mauri, Physical Review B 91, 155436 (2015).
  * Tominaka  _et al._ [2020] S. Tominaka, R. Ishibiki, A. Fujino, K. Kawakami, K. Ohara, T. Masuda, I. Matsuda, H. Hosono, and T. Kondo, Chem 6, 406 (2020).
  * Mortazavi  _et al._ [2018] B. Mortazavi, M. Makaremi, M. Shahrokhi, M. Raeisi, C. V. Singh, T. Rabczuk, and L. F. C. Pereira, Nanoscale 10, 3759 (2018).
  * Katsnelson and Fasolino [2013] M. I. Katsnelson and A. Fasolino, Accounts of chemical research 46, 97 (2013).
  * Choi  _et al._ [2002] H. J. Choi, D. Roundy, H. Sun, M. L. Cohen, and S. G. Louie, Nature 418, 758 (2002).
  * An and Pickett [2001] J. M. An and W. E. Pickett, Phys. Rev. Lett. 86, 4366 (2001).
  * Yildirim  _et al._ [2001] T. Yildirim, O. Gülseren, J. W. Lynn, C. M. Brown, T. J. Udovic, Q. Huang, N. Rogado, K. A. Regan, M. A. Hayward, J. S. Slusky, T. He, M. K. Haas, P. Khalifah, K. Inumaru, and R. J. Cava, Phys. Rev. Lett. 87, 037001 (2001).
  * d’Astuto  _et al._ [2007] M. d’Astuto, M. Calandra, S. Reich, A. Shukla, M. Lazzeri, F. Mauri, J. Karpinski, N. D. Zhigadlo, A. Bossak, and M. Krisch, Phys. Rev. B 75, 174508 (2007).
  * Parendo  _et al._ [2006] K. A. Parendo, K. S. B. Tan, and A. Goldman, Physical Review B 73, 174527 (2006).
  * Caviglia  _et al._ [2008] A. Caviglia, S. Gariglio, N. Reyren, D. Jaccard, T. Schneider, M. Gabay, S. Thiel, G. Hammerl, J. Mannhart, and J.-M. Triscone, Nature 456, 624 (2008).
  * Ahn  _et al._ [1999] C. Ahn, S. Gariglio, P. Paruch, T. Tybell, L. Antognazza, and J.-M. Triscone, Science 284, 1152 (1999).
  * Saito  _et al._ [2016] Y. Saito, T. Nojima, and Y. Iwasa, Nature Reviews Materials 2, 1 (2016).
  * Saito  _et al._ [2015] Y. Saito, Y. Kasahara, J. Ye, Y. Iwasa, and T. Nojima, Science 350, 409 (2015), https://www.science.org/doi/pdf/10.1126/science.1259440 .
  * Belli  _et al._ [2021] F. Belli, T. Novoa, J. Contreras-García, and I. Errea, Nature Communications 12, 1 (2021).


# Full counting statistics of interacting lattice gases after an expansion:
quantifying many-body coherence from atom correlations

Gaétan Hercé  These authors contributed equally to this work. Université
Paris-Saclay, Institut d’Optique Graduate School, CNRS, Laboratoire Charles
Fabry, 91127, Palaiseau, France  Jan-Philipp Bureik  These authors contributed
equally to this work. Université Paris-Saclay, Institut d’Optique Graduate
School, CNRS, Laboratoire Charles Fabry, 91127, Palaiseau, France  Antoine
Ténart  Université Paris-Saclay, Institut d’Optique Graduate School, CNRS,
Laboratoire Charles Fabry, 91127, Palaiseau, France  Alain Aspect  Université
Paris-Saclay, Institut d’Optique Graduate School, CNRS, Laboratoire Charles
Fabry, 91127, Palaiseau, France  Alexandre Dareau  Université Paris-Saclay,
Institut d’Optique Graduate School, CNRS, Laboratoire Charles Fabry, 91127,
Palaiseau, France  David Clément  Université Paris-Saclay, Institut d’Optique
Graduate School, CNRS, Laboratoire Charles Fabry, 91127, Palaiseau, France

###### Abstract

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, suggesting that the condensate component exhibits a full coherence
$g^{(n)}(0)=1$ at any order $n$ up to $n=6$. The approach demonstrated here is
readily extendable to characterize a large variety of interacting quantum
states and phase transitions.

††preprint:

##  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, 4, 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 high-order moments of the particle number, beyond the variance.
Probing high-order moments is a means to study quantum phase transitions [7,
8, 9], universality [10, 11], entanglement properties [12] or out-of-
equilibrium dynamics [13]. The FCS has indeed successfully characterized
various phenomena in mesoscopic conductors [6, 1, 14, 15], Rydberg [16, 17,
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 system, the FCS is
sufficient to identify many (but not all) 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 is indeed equivalent. Moreover,
recent works have shown that applying random unitary transformations before
measuring the FCS provides access to non-diagonal correlators [23, 24],
further motivating the development of experimental approaches to the FCS, in
particular in interacting systems.

Figure 1:  (Left) Free-fall expansion of interacting quantum gases of
metastable Helium-4 from a three-dimensional optical lattice. (Right) Sketch
of the 3D positions of individual atoms detected after an expansion. The full
counting statistics $P(N_{\Omega})$ describes the statistics of the atom
number $N_{\Omega}$ detected in a small voxel of volume $V_{\Omega}\sim(\delta
k)^{3}$ (red cube). To reveal the many-body coherence properties of the
trapped gas of size $L$, $\delta k$ is chosen inferior to $1/L$, $\delta
k\times L\ll 1$.

In this letter, we report the measurement of the full counting statistics in
large three-dimensional (3D) ensembles ($\sim 5\times 10^{3}$ atoms) of
interacting lattice bosons after a free-fall expansion (see Fig. 1). 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 [25]. In quantum gases, far-field –
or momentum – correlations have been measured with single-atom detection in
non-interacting and non-degenerate bosonic [26, 27] and fermionic [28, 29]
gases and in Bose-Einstein condensates [30]. More recently interacting lattice
bosons [31, 32, 33] and interacting fermions [34] were probed. However, the
FCS and high-order correlations have thus far been measured only in 1D thermal
gases of non-interacting bosons [27]. Here, we study various regimes of
interacting 3D Bose gases across the superfluid-to-Mott phase transition, from
weakly-interacting Bose superfluids to Mott insulators.

We characterize the FCS by measuring the probability distribution of the atom
number $N_{\Omega}$ falling in a small voxel $V_{\Omega}$ (see Fig. 1) and we
determine the magnitudes $g^{(n)}(0)$ of correlation functions (up to $n=6$)
based on the factorial moments of $N_{\Omega}$. We compare our findings to the
many-body coherence properties expected for pure state descriptions of both
the superfluid and the Mott phase. Although we find an almost exact agreement,
we observe small deviations to the predictions for a pure condensate in the
case of Bose superfluids. Since, to our knowledge, ab-initio calculations of
the measured quantities with thousands of interacting atoms are beyond state-
of-the-art numerics, we introduce a heuristic model that includes the
contribution of the condensate depletion. Its predictions are in quantitative
agreement with our observations, which leads us to conclude that the
condensate component exhibits a full coherence, $g^{(n)}(0)=1$ at least up to
$n=6$.

##  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 [35] 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_{\Omega})=\langle
N_{\Omega}\rangle^{N_{\Omega}}\exp[-\langle N_{\Omega}\rangle]/N_{\Omega}!$
where $N_{\Omega}$ is the number of detected bosons in the considered volume
$V_{\Omega}$, 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 $g^{(n)}(0)=g^{(n)}({\bm{k}},{\bm{k}},....,{\bm{k}})$, where
${\bm{k}}$ is the momentum at which correlations are evaluated, i.e. where the
volume $V_{\Omega}$ is located, and 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
[36, 37, 31]. 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 a thermal state without site-to-site phase coherence. Thermal states
are characterized by a counting statistics $P(N_{\Omega})=(1-q)q^{N_{\Omega}}$
where $q=\langle N_{\Omega}\rangle/1+\langle N_{\Omega}\rangle$ and
$g^{(n)}({\bm{k}},{\bm{k}},....,{\bm{k}})=n!$ [38]. Note that both the Poisson
and the thermal FCS are fully determined by a single parameter, the average
number of particles $\langle N_{\Omega}\rangle$, as a result of the Gaussian
character of their quantum state [2]. For Gaussian states, a detection
efficiency $\eta$ smaller than one does not affect the measurement of the FCS
– nor that of $g^{(n)}$ – which is of practical interest since $\eta<1$ in
experiments.

Whether these many-body coherence properties associated to pure states
describe the outcome of experiments is not granted. On the one hand, pure
states are only 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 certainly 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 micro-canonical)
ensemble there are no fluctuations of the _total_ BEC atom number $N_{\rm
BEC}$ in a gas of non-interacting bosons at zero temperature [39]: $N_{\rm
BEC}$ is fixed to the total atom number, $N_{\rm BEC}=N$. The statistics of
$N_{\rm BEC}$ is therefore not that of a coherent state [40]. To alleviate
such global constraints and mimic a grand canonical ensemble, one may probe
the FCS within a volume $V_{\Omega}$ much smaller than that, $V_{\rm BEC}$,
occupied by the BEC. The volume $V_{\rm BEC}-V_{\Omega}\sim V_{\rm BEC}$ then
serves as a reservoir for the sub-volume $V_{\Omega}$ where the number of
bosons $N_{\Omega}\ll N$ may fluctuate [41].

Measuring the FCS in small sub-volumes $V_{\Omega}$ is also crucial if
correlation functions
$g^{(n)}({\bm{k}},{\bm{k}^{\prime}},....,{\bm{k}^{\prime\prime}})$ are bell
shaped with widths $l_{c}^{(n)}$. While this is not an issue for a coherent
state which is fully coherent over the entire volume it occupies, the
statistics of a thermal state must be probed over an elementary cell of the
phase space. This is well known in the case of the Hanbury-Brown and Twiss
effect where the property $g^{(2)}({\bm{k}},{\bm{k}^{\prime}})=2$ is expected
only if $|{\bm{k}}-{\bm{k}^{\prime}}|$ is less that the width of the far-field
diffraction pattern associated with the source size. Here, it means that
fully-contrasted correlations are expected only when computed in a volume
smaller than $V_{c}^{(n)}=\large[l_{c}^{(n)}\large]^{3}$ since particles
distant by $l_{c}^{(n)}$ are essentially uncorrelated. In the far-field, the
correlation lengths of Mott insulators and of thermal Bose gases are set by
the inverse in-trap size, $l_{c}^{(n)}\sim 1/L$ [31, 32]. As a result,
observing fully-contrasted FCS and $n$-body correlations requires using
$V_{\Omega}\ll 1/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 $\Delta k\sim 1.6/L$ [42]. With these considerations in mind,
in the following we compute the counting statistics in a volume
$V_{\Omega}\sim(\delta k)^{3}$ of the momentum space such that $\delta k\times
L\ll 1$, where $L$ is the in-trap size of the gas (see Fig. 1). 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

Figure 2:  (a) Full Counting statistics $P(N_{\Omega})$ to find $N_{\Omega}$
atoms in a small volume $V_{\Omega}$ 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 $\langle N_{\Omega}\rangle=0.46(5)$ for
the theoretical predictions (the shaded areas reflect the uncertainty on
$\langle N_{\Omega}\rangle$). The error bars denote the statistical
uncertainty (standard deviation) estimated with the bootstrapping method. (b)
Same as (a) in the BEC mode (${\bf k}={\bf 0}$) of interacting lattice
superfluids (SF) with $U/J=5$. The mean atom number is $\langle
N_{\Omega}\rangle=5.3(2)$. Error bars are smaller than the dots. (c) Magnitude
$g^{(n)}({\bm{0}})$ of $n$-body correlations plotted as a function of the
order $n$ of the correlations. The black squares correspond to a Mott
insulator (panel a) and the blue circles to a Bose superfluid (panel b). The
prediction for thermal Gaussian states, $g^{(n)}(0)=n!$, is shown as a black
solid line. (d) $g^{(n)}({\bm{0}})$ in the SF (blue circles, same data as in
(c)) and in the randomized set (orange squares).

Our measurement of the FCS in the far-field exploits 3D atom-by-atom detection
of metastable Helium-4 (4He∗) after a long free-fall expansion [43, 44], with
a detection efficiency of $\eta=0.53(2)$ per atom and negligible dark counts
(see Fig. 1). For an individual run, we register the number of atoms in each
of the voxels $V_{\Omega}$ 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 4He∗ atoms loaded in the lowest energy-
band of a three-dimensional (3D) optical lattice [45]. The lattice implements
the 3D Bose-Hubbard Hamiltonian whose main parameters are the tunnelling
amplitude $J$ and the on-site (repulsive) interaction $U$.

We first investigate the thermal nature of Mott insulators in the momentum
space. We realize Mott insulators with $N=10(3)\times 10^{3}$ atoms at
$U/J=95$, which corresponds to a lattice filling of one atom per site at the
trap center [31] 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_{\Omega}$
of size $\delta k=4\times 10^{-2}k_{d}\ll 1/L$ and average the probability
distributions measured over all these voxels. Here $k_{d}=2\pi/d$ is the
momentum associated with the lattice spacing $d=775~{}$nm. 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, $\langle N_{\Omega}\rangle=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_{\Omega})$. We determine
$g^{(n)}(0)$ from the factorial moments of the detected atom number,
transposing a well-known approach in quantum optics [46]. The correlation
functions of quantum optics [22] are defined with normal ordering of the
destruction operators, since a detected photon is destroyed and there is one
photon less in the volume of measurement. 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 [47]. In the
experiments reported here, each detected He∗ atom is destroyed and the results
of quantum optics can be directly transposed. Finally, we use a similar
approach as that used for $P(N_{\Omega})$ – averaging over many runs and then
over many voxels within the first Brillouin zone.

We first evaluate the factorial moments of the atom number $N_{\delta
k_{\perp}}$ detected in voxels of adjustable volume $V(\delta
k_{\perp})=\delta k\delta k_{\perp}^{2}$ and normalize them with the average
atom number in the same voxels, to obtain

$g^{(n)}_{\delta k_{\perp}}(0)=\frac{\sum_{j} \langle N_{\delta
k_{\perp}}(N_{\delta k_{\perp}}-1)...(N_{\delta
k_{\perp}}-n+1)\rangle_{j}}{\sum_{j} \langle N_{\delta
k_{\perp}}\rangle_{j}^{n}}$ (1)

where $\sum_{j}$ refers to the sum over all the voxels (labelled by the index
$j$) belonging to the first Brillouin zone. The magnitudes $g^{(n)}(0)$ of the
$n$-body correlation functions are obtained from extrapolating
$g^{(n)}_{\delta k_{\perp}}(0)$ at $\delta k_{\perp}\rightarrow\delta k$ (see
Appendix A). The results are plotted in Fig. 2(c), along with the prediction
for thermal states, $g^{(n)}(0)=n!$. The excellent quantitative agreement
demonstrates that the $n$-body coherence of a unity-filling Mott insulator is
thermal [36, 37, 31]. 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_{c}^{(n)}$ which are
determined by their different in-trap sizes and the incompressible (resp.
compressible) character of Mott insulators (resp. thermal gases) [31].

Figure 3:  (a)-(b) Plots of $g^{(n)}(0)$ measured at ${\bm{k}}={\bm{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. 2(c). The dashed lines (resp. the shaded
areas) are the predictions of the model with the values $f_{\rm coh}$ (resp.
the uncertainties on $f_{\rm 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_{\Omega}$ where the counting statistics is
evaluated. The horizontal dashed-dotted lines indicate the fitted values
$1-f_{\rm coh}$. The dashed lines are Lorentzian functions fitted to the tails
of the densities in the range $[0.2k_{d},0.5k_{d}]$ in order to estimate the
density of the depletion at $k=0$ (the shaded areas represent the error from
the fitted parameters).

In a second set of experiments, we address the $n$-body coherence of the BEC
component in lattice superfluids with $N=5(1)\times 10^{3}$ at $U/J=5$. In the
momentum space, the BEC occupies a volume of width $\Delta k\simeq 0.15k_{d}$
centered at ${\bm{k}}={\bm{0}}$ [33]. This volume contains many atoms in each
run and the statistics of the atom number falling in a sphere $S_{\Omega}$ of
a radius $\delta k=0.025k_{d}\ll\Delta k$, centered at ${\bm{k}}={\bm{0}}$, is
sufficient to extract the counting statistics. We measure the counting
statistics $P(N_{\Omega})$ in $S_{\Omega}$ 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 $\langle
N_{\Omega}\rangle=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_{\Omega}$,

$g^{(n)}(0)=\frac{\langle
N_{\Omega}(N_{\Omega}-1)...(N_{\Omega}-n+1)\rangle}{\langle
N_{\Omega}\rangle^{n}}$ (2)

In Fig. 2c we plot $g^{(n)}(0)$ as a function of $n$ and find that
$g^{(n)}(0)\sim 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 [30]: the
deviation in the FCS (see Fig. 2(b)) is reflected in the fact that
$g^{(n)}(0)>1$, as shown in Fig. 2(d). 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. 2(d)), confirming that the deviations in the (non-
randomised) experimental data are significant. Note that the randomization
method also yields a Poisson statistics when applied to the Mott insulator
data set. In the next section, we propose an interpretation of the tiny
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

A Poisson statistics is expected when only atoms from the BEC, which is
thought of as a coherent state, contribute to the measured statistics in a
volume $S_{\Omega}$ that is small compared to the coherence volume. At a
macroscopic level though, a fraction of the total atom number is expelled from
the BEC by the finite interactions (quantum depletion) and by the finite
temperature (thermal depletion). Even if it amounts to a negligible fraction
of the atom number $N_{\Omega}$ falling in $S_{\Omega}$, the total depletion
of the condensate may contribute to the measured statistics. 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_{\Omega}$ [40], and (ii) that the BEC is a coherent state
while the depletion exhibits thermal statistics in $S_{\Omega}$. We define the
“coherent fraction” $f_{\rm coh}$ as the fraction of the BEC atoms falling
within $S_{\Omega}$. We emphasize that describing the contribution of the
quantum depletion in $S_{\Omega}$ with a thermal statistics is an assumption.
We have shown that the statistics of the quantum depletion is thermal when
measured at non-zero momenta, outside the BEC [32]. However, whether this is
also valid for small momenta belonging to the BEC mode is difficult to assess
in a harmonic trap. In contrast to homogeneous systems for which the momentum
is a good quantum number, in the harmonically-trapped case the quantum
depletion contributes to the FCS measured at ${\bf k}={\bf 0}$. This
contribution is non-trivial since the ground-state of weakly-interacting
bosons is a coherent superposition of the BEC and the quantum depletion.

With the hypotheses of our model, we obtain an analytical prediction for
$g^{(n)}(0)$ that depends only on the coherent fraction $f_{\rm coh}$ (see
Appendix C),

$g^{(n)}(0)-1=\sum_{p=1}^{n-1}\left[(n-p)!\binom{n}{p}^{2}-\binom{n}{p}\right]f_{\rm
coh}^{p}(1-f_{\rm coh})^{n-p}$ (3)

In Fig. 3(a), we fit the experimental data shown in Fig. 2(d) 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}\simeq 25-30$
[45].

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), see Fig. 3. Firstly, we
confirm that Eq. (3) correctly fits the values of $g^{(n)}(0)$ with a single
adjustable parameter $f_{\rm coh}$. Secondly, the extracted values of $f_{\rm
coh}$ decrease with the interaction strength as intuitively expected. The
uncertainty on the values $f_{\rm coh}$ is small, at the $\sim 0.1\%$ level.
As can be inferred from Fig. 3(a)-(b), the larger the order $n$ of
correlations we measure, the smaller the uncertainty on $f_{\rm 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_{\rm coh}$ to
the fraction $\eta_{D}$ of depleted atoms detected within $S_{\rm BEC}$. We
are not aware of a quantitative analytical prediction for $\eta_{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.2k_{d},0.5k_{d}]$) with a Lorentzian function
to extrapolate the density of the depletion at $k=0$. Note that using a
Lorentzian function is an arbitrary choice which happens to correctly fit the
tails. This analysis indicates that the values $1-f_{\rm coh}$ are compatible
with the extrapolated densities, while both quantities vary by one order of
magnitude.

Figure 4:  Coherent fraction $f_{\rm coh}$ as a function of the condensate
fraction $f_{c}$. The dashed-line is the prediction of Eq. (4) where the ratio
$\mathcal{V}_{c}/\mathcal{V}_{d}=0.03(1)$ is evaluated from the measured
density profiles (see Appendix D). The shaded area reflects the uncertainty on
$\mathcal{V}_{c}/\mathcal{V}_{d}$.

To further assess the validity of the model, we proceed with a quantitative
comparison of the measured coherent fraction to the measured condensate
fraction [43]. The coherent fraction $f_{\rm coh}$ in $S_{\Omega}$ differs
from the condensate fraction $f_{c}$ of the entire gas since the radius of
$S_{\Omega}$ is much smaller than $1/L$. Moreover, the scaling of $f_{\rm
coh}$ with $f_{c}$ is non-linear since $f_{\rm coh}$ reveals the BEC atom
number in a small volume $S_{\Omega}$ where the BEC contribution is maximum
while that of the depletion is small (except when $f_{c}\ll 1$). The volumes
occupied by the BEC ($\mathcal{V}_{c}$) and by the depletion
($\mathcal{V}_{d}$) set their respective contributions. A simple estimate (see
Appendix D) leads to

$f_{\rm
coh}\simeq\frac{f_{c}}{f_{c}+(1-f_{c})\mathcal{V}_{c}/\mathcal{V}_{d}}.$ (4)

In Fig. 4 we plot the measured values of $f_{\rm coh}$ and $f_{c}$, along with
the prediction of Eq. (4). Here, $f_{c}$ is obtained similarly to [43] and the
ratio $\mathcal{V}_{c}/\mathcal{V}_{d}$ used to plot Eq. (4) is estimated from
the measured density profiles (see Appendix D). The quantitative prediction of
our model (without any adjustable parameter) correctly reproduces the observed
non-linear dependency of $f_{\rm coh}$ with $f_{c}$. We thus conclude that the
heuristic model we introduced quantitatively captures the observed deviations
to $g^{(n)}(0)=1$ at ${\bm{k}}={\bm{0}}$, attributing the latter to the
presence of depleted atoms. In turns, this suggests that atoms in the BEC have
the statistics of a fully coherent state, with $g^{(n)}(0)=1$ at any order
$n\leq 6$.

Note that in their evaluation of $g^{(3)}(0)$ in a BEC [30], Hodgman et al. do
not observe the small deviation reported here. A major difference with our
experiment is that their observation volume $V_{\Omega}$ is quite large, as
illustrated by the fact that in the thermal case it entails a convolution
reducing the expected values of $g^{(2)}(0)$ or $g^{(3)}(0)$ from 2 and 6 to
$1.022(2)$ and $1.061(6)$, respectively. We conclude that our ability to
obtain good statistics in tiny volumes is a major asset of our experimental
device.

##  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 have investigated two iconic states of interacting lattice bosons, Mott
insulators and lattice superfluids. We find that a Mott insulator exhibits a
thermal statistics in the momentum 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 from studying the coherent fraction at increasing interaction
strengths. The results of this work 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. [11, 10, 9, 48]). 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 [49]. The second condition is met with various
detection methods and atomic species [50].

Another interesting direction consists in realizing recent proposals to access
non-trivial $n$-body correlations [23, 24]. These theoretical works have shown
that non-diagonal correlators can be accessed from combining two-particle
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égion 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.

## References

  * Blanter and Büttiker [2000] Blanter, Y. and Büttiker, M.,  _Shot noise in mesoscopic conductors_ , Physics Reports 336, 1 (2000).
  * Gardiner and Zoller [2004] Gardiner, C. and Zoller, P.,  _Quantum Noise - A Handbook of Markovian and Non-Markovian Quantum Stochastic Methods with Applications to Quantum Optics_ (Springer Berlin, Heidelberg, 2004).
  * Burt  _et al._ [1997] Burt, E. A., Ghrist, R. W., Myatt, C. J., Holland, M. J., Cornell, E. A. and Wieman, C. E.,  _Coherence, Correlations, and Collisions: What One Learns about Bose-Einstein Condensates from Their Decay_ , Phys. Rev. Lett. 79, 337 (1997).
  * Altman  _et al._ [2004] Altman, E., Demler, E. and Lukin, M. D.,  _Probing many-body states of ultracold atoms via noise correlations_ , Phys. Rev. A 70, 013603 (2004).
  * 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).
  * Levitov  _et al._ [1996] Levitov, L. S., Lee, H. and Lesovik, G. B.,  _Electron counting statistics and coherent states of electric current_ , Journal of Mathematical Physics 37, 4845 (1996).
  * Ivanov and Abanov [2010] Ivanov, D. A. and Abanov, A. G.,  _Phase transitions in full counting statistics for periodic pumping_ , EPL (Europhysics Letters) 92, 37008 (2010).
  * Gómez-Ruiz  _et al._ [2020] Gómez-Ruiz, F. J., Mayo, J. J. and del Campo, A.,  _Full Counting Statistics of Topological Defects after Crossing a Phase Transition_ , Phys. Rev. Lett. 124, 240602 (2020).
  * Devillard  _et al._ [2020] Devillard, P., Chevallier, D., Vignolo, P. and Albert, M.,  _Full counting statistics of the momentum occupation numbers of the Tonks-Girardeau gas_ , Phys. Rev. A 101, 063604 (2020).
  * Eisler [2013] Eisler, V.,  _Universality in the Full Counting Statistics of Trapped Fermions_ , Phys. Rev. Lett. 111, 080402 (2013).
  * Lovas  _et al._ [2017] Lovas, I., Dóra, B., Demler, E. and Zaránd, G.,  _Full counting statistics of time-of-flight images_ , Phys. Rev. A 95, 053621 (2017).
  * Klich and Levitov [2009] Klich, I. and Levitov, L.,  _Quantum Noise as an Entanglement Meter_ , Phys. Rev. Lett. 102, 100502 (2009).
  * Esposito  _et al._ [2009] Esposito, M., Harbola, U. and Mukamel, S.,  _Nonequilibrium fluctuations, fluctuation theorems, and counting statistics in quantum systems_ , Rev. Mod. Phys. 81, 1665 (2009).
  * Gustavsson  _et al._ [2006] Gustavsson, S., Leturcq, R., Simovič, B., Schleser, R., Ihn, T., Studerus, P., Ensslin, K., Driscoll, D. C. and Gossard, A. C.,  _Counting Statistics of Single Electron Transport in a Quantum Dot_ , Phys. Rev. Lett. 96, 076605 (2006).
  * Maisi  _et al._ [2014] Maisi, V. F., Kambly, D., Flindt, C. and Pekola, J. P.,  _Full Counting Statistics of Andreev Tunneling_ , Phys. Rev. Lett. 112, 036801 (2014).
  * Liebisch  _et al._ [2005] Liebisch, T. C., Reinhard, A., Berman, P. R. and Raithel, G.,  _Atom Counting Statistics in Ensembles of Interacting Rydberg Atoms_ , Phys. Rev. Lett. 95, 253002 (2005).
  * Malossi  _et al._ [2014] Malossi, N., Valado, M. M., Scotto, S., Huillery, P., Pillet, P., Ciampini, D., Arimondo, E. and Morsch, O.,  _Full Counting Statistics and Phase Diagram of a Dissipative Rydberg Gas_ , Phys. Rev. Lett. 113, 023006 (2014).
  * Zeiher  _et al._ [2016] Zeiher, J., van Bijnen, R., Schauß, P., Hild, S., Choi, J.-y., Pohl, T., Bloch, I. and Gross, C.,  _Many-body interferometry of a Rydberg-dressed spin lattice_ , Nature Physics 12, 1095 (2016).
  * Öttl  _et al._ [2005] Öttl, A., Ritter, S., Köhl, M. and Esslinger, T.,  _Correlations and Counting Statistics of an Atom Laser_ , Phys. Rev. Lett. 95, 090404 (2005).
  * Perrier  _et al._ [2019] Perrier, M., Amodjee, Z., Dussarrat, P., Dareau, A., Aspect, A., Cheneau, M., Boiron, D. and Westbrook, C. I.,  _Thermal counting statistics in an atomic two-mode squeezed vacuum state_ , SciPost Phys. 7, 2 (2019).
  * Flammia  _et al._ [2012] Flammia, S. T., Gross, D., Liu, Y.-K. and Eisert, J.,  _Quantum tomography via compressed sensing: error bounds, sample complexity and efficient estimators_ , New Journal of Physics 14, 095022 (2012).
  * Glauber [1963] Glauber, R. J.,  _The Quantum Theory of Optical Coherence_ , Phys. Rev. 130, 2529 (1963).
  * Brunner  _et al._ [2021] Brunner, E., Buchleitner, A. and Dufour, G.,  _Many-body coherence and entanglement from randomized correlation measurements_ , arXiv (2021), 10.48550/arxiv.2107.01686.
  * Naldesi  _et al._ [2022] Naldesi, P., Elben, A., Minguzzi, A., Clément, D., Zoller, P. and Vermersch, B.,  _Fermionic correlation functions from randomized measurements in programmable atomic quantum devices_ , arXiv (2022), 10.48550/arxiv.2205.00981.
  * Aspect [2019] Aspect, A.,  _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 Browaeys, A., Lahaye, T., Porto, T., Adams, C. S., Weidemüller, M. and Cugliandolo, L. F. (Oxford University Press, 2019).
  * Schellekens  _et al._ [2005] Schellekens, M., Hoppeler, R., Perrin, A., Gomes, J. V., Boiron, D., Aspect, A. and Westbrook, C. I.,  _Hanbury Brown Twiss Effect for Ultracold Quantum Gases_ , Science 310, 648 (2005).
  * Dall  _et al._ [2013] Dall, R. G., Manning, A. G., Hodgman, S. S., RuGway, W., Kheruntsyan, K. V. and Truscott, A. G.,  _Ideal n-body correlations with massive particles_ , Nature Physics 9, 341 (2013).
  * Jeltes  _et al._ [2007] Jeltes, T. _et al._ ,  _Comparison of the Hanbury Brown–Twiss effect for bosons and fermions_ , Nature 445, 402 (2007).
  * Bergschneider  _et al._ [2019] Bergschneider, A., Klinkhamer, V. M., Becher, J. H., Klemt, R., Palm, L., Zürn, G., Jochim, S. and Preiss, P. M.,  _Experimental characterization of two-particle entanglement through position and momentum correlations_ , Nature Physics 15, 640 (2019).
  * Hodgman  _et al._ [2011] Hodgman, S. S., Dall, R. G., Manning, A. G., Baldwin, K. G. H. and Truscott, A. G.,  _Direct Measurement of Long-Range Third-Order Coherence in Bose-Einstein Condensates_ , Science 331, 1046 (2011).
  * Carcy  _et al._ [2019] Carcy, C., Cayla, H., Tenart, A., Aspect, A., Mancini, M. and Clément, D.,  _Momentum-Space Atom Correlations in a Mott Insulator_ , Phys. Rev. X 9, 041028 (2019).
  * 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).
  * 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).
  * Holten  _et al._ [2022] Holten, M., Bayha, L., Subramanian, K., Brandstetter, S., Heintze, C., Lunt, P., Preiss, P. M. and Jochim, S.,  _Observation of Cooper pairs in a mesoscopic two-dimensional Fermi gas_ , Nature 606, 287 (2022).
  * Pitaevskii and Stringari [2004] Pitaevskii, L. P. and Stringari, S.,  _Bose-Einstein Condensation_ (Clarendon press, Oxford, 2004).
  * Fölling  _et al._ [2005] Fölling, S., Gerbier, F., Widera, A., Mandel, O., Gericke, T. and Bloch, I.,  _Spatial quantum noise interferometry in expanding ultracold atom clouds_ , Nature 434, 481 (2005).
  * 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).
  * Walls and Milburn [2008] Walls, D. F. and Milburn, G. J.,  _Quantum Optics_ (Springer Berlin, Heidelberg, 2008).
  * Kristensen  _et al._ [2019] Kristensen, M. A., Christensen, M. B., Gajdacz, M., Iglicki, M., Pawłowski, K., Klempt, C., Sherson, J. F., Rzażewski, K., Hilliard, A. J. and Arlt, J. J.,  _Observation of Atom Number Fluctuations in a Bose-Einstein Condensate_ , Phys. Rev. Lett. 122, 163601 (2019).
  * Castin and Dum [1998] Castin, Y. and Dum, R.,  _Low-temperature Bose-Einstein condensates in time-dependent traps: Beyond the $U(1)$ symmetry-breaking approach_, Phys. Rev. A 57, 3008 (1998).
  * [41] 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_{\Omega}$ is largely immune to small shot-to-shot fluctuations of $N$.
  * Stenger  _et al._ [1999] Stenger, J., Inouye, S., Chikkatur, A. P., Stamper-Kurn, D. M., Pritchard, D. E. and Ketterle, W.,  _Bragg Spectroscopy of a Bose-Einstein Condensate_ , Phys. Rev. Lett. 82, 4569 (1999).
  * 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).
  * 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. Research 2, 013017 (2020).
  * 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).
  * Laiho  _et al._ [2022] Laiho, K., Dirmeier, T., Schmidt, M., Reitzenstein, S. and Marquardt, C.,  _Measuring higher-order photon correlations of faint quantum light: A short review_ , Physics Letters A 435, 128059 (2022).
  * Goodman [1985] Goodman, J. W.,  _Statistical optics_, edited by New York, W.-I., Vol. 1 (1985).
  * Contessi  _et al._ [2022] Contessi, D., Recati, A. and Rizzi, M.,  _Phase Diagram Detection via Gaussian Fitting of Number Probability Distribution_ , arXiv , 2207.01478 (2022).
  * Bloch  _et al._ [2008] Bloch, I., Dalibard, J. and Zwerger, W.,  _Many-body physics with ultracold gases_ , Rev. Mod. Phys. 80, 885 (2008).
  * Ott [2016] Ott, H.,  _Single atom detection in ultracold quantum gases: a review of current progress_ , Reports on Progress in Physics 79, 054401 (2016).

##  Appendix A Data analysis of correlations in Mott insulators

The Mott insulator is spread over the entire Brillouin zone in the momentum
space, with a small density. As a consequence, the number of atoms falling in
voxels of tiny volumes is small. In particular, the statistics of the atom
number in a voxel smaller than the coherence volume is not sufficient to
extract the magnitude $g^{(n)}(0)$ of the $n$-body correlation functions. To
circumvent this issue, we proceed as follows.

Figure 5: Plot of the magnitude $g^{(n)}_{\delta k_{\perp}}(0)$ as a function
of the transverse integration $\delta k_{\perp}$. 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.

We first consider anisotropic voxels of volume $V(\delta k_{\perp})=\delta
k\times\delta k_{\perp}^{2}$ with $\delta k_{\perp}>l_{c}\gg\delta k$. The
transverse integration over $\delta k_{\perp}$ raises the statistics of the
detected atom number $N_{\delta k_{\perp}}$, allowing us to compute the
factorial moments of $N_{\delta k_{\perp}}$ in each voxel labeled by the index
$j$, $\langle N_{\delta k_{\perp}}(N_{\delta k_{\perp}}-1)...(N_{\delta
k_{\perp}}-n+1)\rangle_{j}$. To further raise the statistics, we average the
factorial moments over many such voxels, $\sum_{j} \langle N_{\delta
k_{\perp}}(N_{\delta k_{\perp}}-1)...(N_{\delta k_{\perp}}-n+1)\rangle_{j}$
where $\sum_{j}$ refers to the sum over all the voxels belonging to the first
Brillouin zone.

The normalized correlation function $g^{(n)}_{\delta k_{\perp}}(\delta k)$ in
voxels of volume $V(\delta k_{\perp})$ is obtained from normalizing the
factorial moments,

$g^{(n)}_{\delta k_{\perp}}(\delta k)=\frac{\sum_{j} \langle N_{\delta
k_{\perp}}(N_{\delta k_{\perp}}-1)...(N_{\delta
k_{\perp}}-n+1)\rangle_{j}}{\sum_{j} \langle N_{\delta
k_{\perp}}\rangle_{j}^{n}}$ (5)

We then proceed by fitting the dependency of $g^{(n)}_{\delta
k_{\perp}}(\delta k)$ with $\delta k$ with a Gaussian function, from which we
extract its magnitude $g^{(n)}_{\delta k_{\perp}}(0)$ at the considered
transverse integration $\delta k_{\perp}$. At large transverse integrations,
$\delta k_{\perp}\gg l_{c}$, $g^{(n)}_{\delta k_{\perp}}(0)$ is reduced with
respect to the value $g^{(n)}(0)$ expected at vanishingly small integration,
$\delta k_{\perp}\rightarrow 0$. To obtain $g^{(n)}(0)$, we repeat the
calculation of $g^{(n)}_{\delta k_{\perp}}(0)$ for many values of $\delta
k_{\perp}$, with the results shown in Fig. 5. From the plots in Fig. 5, we
extrapolate the value $g^{(n)}(0)$ at $\delta k_{\perp}\rightarrow 0$.

##  Appendix B Data analysis of correlations in Bose superfluids

The magnitudes of the normalized correlation functions are obtained from the
factorial moments of the atom number $N_{\Omega}$ in the small volume
$V_{\Omega}$. The $n$-factorial moment is defined as
$\langle(N_{\Omega})_{n}\rangle=\langle
N_{\Omega}(N_{\Omega}-1)...(N_{\Omega}-n+1)\rangle$ and $g^{(n)}(0)$ are
expressed as

$g^{(n)}(0)=\frac{\langle
N_{\Omega}(N_{\Omega}-1)...(N_{\Omega}-n+1)\rangle}{\langle
N_{\Omega}\rangle^{n}}=\frac{\langle(N_{\Omega})_{n}\rangle}{\langle
N_{\Omega}\rangle^{n}}$ (6)

We define the error bars on $g^{(n)}(0)$ as $\Delta(N_{\Omega})_{n}/\langle
N_{\Omega}\rangle^{n}$, where $\Delta(N_{\Omega})_{n}$ is the standard error
of the factorial moment of order $n$,

$\Delta(N_{\Omega})_{n}=\frac{1}{\sqrt{N_{\rm
runs}}}\sqrt{\langle(N_{\Omega})_{n}^{2}\rangle-\langle(N_{\Omega})_{n}\rangle^{2}}$
(7)

##  Appendix C Model for the atom number statistics at $\vec{k}=\vec{0}$ 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
${\bm{k}}={\bf 0}$. We assume that the statistical properties of the BEC mode
are those of a coherent state, $g^{(n)}_{\rm BEC}(0)=1$. In contrast, the
statistical properties of the depletion are assumed to be those of a thermal
state, $g^{(n)}_{\rm dep}(0)=n!$. This latter assumption holds because we
consider a small voxel of size $dk=0.025k_{d}\ll l_{c}$ 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
$\hat{N}_{\Omega}=a_{\Omega}^{\dagger}a_{\Omega}$ falling in the considered
voxel, where $a_{\Omega}=a_{\rm BEC}+a_{\rm dep}$ is the sum of the operators
associated with the BEC and the depletion. Since we assume the operators
$a_{\rm BEC}$ and $a_{\rm dep}$ are uncorrelated, one simply obtains

$\langle(a_{\Omega}^{\dagger})^{n}a_{\Omega}^{n}\rangle=\sum_{p=1}^{n}\binom{n}{p}^{2}\langle(a_{\rm
BEC}^{\dagger})^{p}a_{\rm BEC}^{p}(a_{\rm dep}^{\dagger})^{n-p}a_{\rm
dep}^{n-p}\rangle$ (8)

where by definition $a_{\Omega}^{n}=\sum_{p=1}^{n}\binom{n}{p}a_{\rm
BEC}^{p}a_{\rm dep}^{n-p}$. This leads to the formula given in the main text,

$\displaystyle g^{(n)}(0)-1$ $\displaystyle=$
$\displaystyle\frac{\langle(a_{\Omega}^{\dagger})^{n}a_{\Omega}^{n}\rangle}{\langle
a_{\Omega}^{\dagger}a_{\Omega}\rangle^{n}}-1$ $\displaystyle=$
$\displaystyle\sum_{p=1}^{n-1}\left[(n-p)!\binom{n}{p}^{2}-\binom{n}{p}\right]f_{\rm
coh}^{p}(1-f_{\rm coh})^{n-p}$

##  Appendix D Relation $f_{\rm coh}$ versus $f_{c}$.

The coherent fraction is approximatively given by

$f_{\rm coh}\simeq\frac{\rho_{\rm BEC}(0)}{\rho_{\rm BEC}(0)+\rho_{\rm
dep}(0)}$ (9)

where $\rho_{\rm BEC}(0)$ (resp. $\rho_{\rm dep}(0)$) is the momentum density
of the BEC (resp. of the depletion) at the center of the Brillouin zone,
$\vec{k}=\vec{0}$. The condensate fraction is approximately given by

$f_{c}\simeq\frac{\rho_{\rm BEC}(0)\mathcal{V}_{c}}{\rho_{\rm
BEC}(0)\mathcal{V}_{c}+\rho_{\rm dep}(0)\mathcal{V}_{d}}$ (10)

where $\mathcal{V}_{c}$ (resp. $\mathcal{V}_{d}$) is the total volume of the
momentum-space occupied by the BEC (resp. by the depletion) normalized to the
density at $\vec{k}=\vec{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 $\mathcal{V}_{c}$ and $\mathcal{V}_{d}$. We describe the
BEC density profile with an isotropic function [42]

$n_{\rm BEC}({\bf k})=\rho_{\rm BEC}(0)e^{-k^{2}/2\sigma_{\rm bec}^{2}}$ (11)

and the density profile of the depletion with a lorentzian function along each
(uncoupled) lattice axis,

$n_{\rm dep}({\bf k})=\rho_{\rm dep}(0)\ \Pi_{j=x,y,z}\left(\frac{\sigma_{\rm
dep}^{2}/4}{k_{j}^{2}+(\sigma_{\rm dep}/2)^{2}}\right).$ (12)

We obtain $\mathcal{V}_{c}=(\sqrt{2\pi}\sigma_{\rm bec})^{3}$ and
$\mathcal{V}_{d}=\left(\sigma_{\rm dep}\arctan\left[\frac{k_{d}}{\sigma_{\rm
dep}}\right]\right)^{3}$. From fitting the measured density profiles with the
above dependency, we obtain $\mathcal{V}_{c}/\mathcal{V}_{d}=0.03(1)$, a value
which we use to plot the theoretical prediction in Fig. 4.