Observation of universal non-Gaussian statistics of the order parameter across a continuous phase transition

TL;DR

This study experimentally observes universal non-Gaussian fluctuations of the order parameter in a continuous phase transition of an interacting Bose gas, revealing critical scaling and challenging traditional homogeneous models.

Contribution

It provides the first direct measurement of non-Gaussian order parameter statistics and demonstrates their universality and critical scaling across a phase transition.

Findings

Non-Gaussian order parameter fluctuations are observed near the transition.

Fluctuations exhibit universal oscillations in high-order cumulants.

Effective potential analysis reveals a non-trivial minimum in the superfluid phase.

Abstract

Second-order phase transitions are characterised by critical scaling and universality. The singular behaviour of thermodynamic quantities at the transition, in particular, is determined by critical exponents of the universality class of the transition. However, critical properties are also characterised by the probability distribution of order parameter across the transition, where non-Gaussian statistics are expected, but remain largely unexplored. Here, making use of single-atom-resolved detection in momentum space, we measure the full probability distribution of the amplitude of the order parameter across a continuous phase transition in an interacting lattice Bose gas. We find that fluctuations are captured by an effective potential -- reconstructed from the measured probability distribution by analogy with Landau theory -- displaying a non-trivial minimum in the superfluid (ordered) phase, which vanishes at the transition point. Additionally, we observe non-Gaussian statistics of the order parameter near the transition, distinguished by non-zero and oscillating high-order cumulants. We provide direct experimental evidence that these oscillations are universal, and show numerically that they exhibit critical scaling. Our experiments are conducted in inhomogeneous systems, challenging the conventional understanding of criticality, which is primarily based on homogeneous models. Our results underscore the crucial role of order parameter statistics in probing critical phenomena and universality.

Equations

$$\hat{H}=-J\sum _{⟨\mathbf{i},\mathbf{j}⟩}[{\hat{a}}_{\mathbf{i}}^{†}{\hat{a}}_{\mathbf{j}}+h.c.]+\frac{U}{2}\sum _{\mathbf{i}}{\hat{n}}_{\mathbf{i}}({\hat{n}}_{\mathbf{i}}-1)+\sum _{\mathbf{i}}{V}_{\mathbf{i}}{\hat{n}}_{\mathbf{i}}.$$

$$\mathcal{F}(|\psi |)\simeq a_0+a_2|\psi {|}^2+a_4|\psi {|}^4,$$

$$P[\varphi (\mathbf{r})]\propto {\mathrm{e}}^{-\beta \mathcal{H}[\varphi (\mathbf{r})]},$$

$$P_{\mathrm{GL}}[\psi (\mathbf{r})]\propto {\mathrm{e}}^{-{\mathcal{F}}_{\mathrm{GL}}[\psi (\mathbf{r})]},$$

$${\mathcal{F}}_{\mathrm{GL}}[\psi (\mathbf{r})]=\int \mathrm{d}\mathbf{r}\left[|\nabla \psi (\mathbf{r}){|}^2+{\tilde{a}}_2|\psi (\mathbf{r}){|}^2+{\tilde{a}}_4|\psi (\mathbf{r}){|}^4+\dots \right].$$

$$P(\psi )\simeq P_{\mathrm{L}}(\psi )=e^{-V{\mathcal{F}}_{\mathrm{L}}(\psi )}.$$

$$P(\psi )=L^{2\beta /\nu }\tilde{P}(L^{\beta /\nu }\psi ,L^{1/\nu }\delta u).$$

$$\sum _i{\mathrm{e}}^{-(a_0+a_2|\psi {|}_i^2+a_4|\psi {|}_i^4)}=1,$$

$$\sum _{i}f({\psi }_i)=1,$$

$${\hat{a}}_i^{†}{\hat{a}}_j+h.c.\approx 2n_{\mathrm{QR}}\cos ({\varphi }_i-{\varphi }_j)$$

$$H_{\mathrm{QR}}=-2J{n}_{\mathrm{QR}}\sum _{⟨ij⟩}\cos ({\varphi }_i-{\varphi }_j)-\frac{U}{2}\sum _i\frac{{\partial }^2}{\partial {\varphi }_i^2}.$$

$${\kappa }_n(u,L)=L^{-n\beta /\nu }{f}_n(\delta ),$$

Full text

Continuous phase transitions driven by spontaneous symmetry breaking [14, 15] are fundamental phenomena covering all branches of physics [16]. The concept of a fluctuating order parameter $\psi$ plays an essential role in their description. Its expectation value $\braket{\psi}$ was shown by Landau [14] to distinguish the ordered phase, where it assumes a non-zero value in the thermodynamic limit, from the disordered phase where it vanishes. Landau’s mean-field theory, however, neglects fluctuations in the order parameter and fails to capture the singular behaviour at the transition, i.e. the universal critical regime.

These phenomena are in fact driven by the large fluctuations associated with the divergence of the correlation length $\xi$ [8, 16]. As this length grows beyond microscopic scales (related to the interparticle distance), the system enters the critical regime (Fig. 1a) characterised by power laws with universal critical exponents. These features are well understood in the framework of Wilson’s renormalisation group [17, 18]. In finite-size systems, though, a peculiar regime appears within the critical one, when the correlation length grows larger than the linear system size $L$. Under this condition, the central limit theorem breaks down, leading to significant changes in the order parameter’s statistics, which become non-Gaussian [19, 20].

The probability distribution function (PDF) of the order parameter, which fully captures its statistical properties, is a fundamental quantity in the description of phase transitions. At the transition point, its non-Gaussian shape is universal, i.e. it only depends on the universality class of the transition. It has received extensive theoretical and numerical attention over the past few decades [2, 21, 22, 23, 24, 25, 26, 27], focusing in the thermodynamic limit where $\xi$ and $L$ are sent to infinity at fixed $\xi/L$. From an experimental perspective, however, measuring the non-Gaussian PDF of an order parameter at criticality poses major challenges. The non-Gaussian regime is typically a narrow region of the control parameter near the transition point (Fig. 1a) that shrinks with increasing system size. In a macroscopic system with Avogadro number of particles, it reduces to the critical point, requiring therefore an unrealistic precision in the tuning of the control parameter. Furthermore, characterising critical non-Gaussian statistics is inherently difficult [28], and, to our knowledge, has only been reported in a classical, mesoscopic system of liquid crystals [9, 29]. More generally, non-Gaussian statistics are a diagnostic of complexity and criticality, not just in physics but across many real-world phenomena [30], and their characterisation is fundamental to the understanding of complex systems.

Modern experimental platforms implement quantum many-body models in which phase transitions with substantial non-Gaussian critical regimes are expected, thanks to small system sizes of tens to hundreds of particles [31, 32, 33, 34, 35, 36, 37, 38, 39, 40]. In these platforms, single-particle-resolved detection methods [41, 42] provide access to the full counting statistics of observables [43, 44, 45, 36, 37, 46, 47]. These advances enable the exploration of the universal order-parameter statistics, which, as argued above, is complementary to the standard measurements of critical exponents [48, 49, 50, 51].

Here, we leverage these techniques to measure the full probability distribution of the condensate order parameter across the interaction-driven superfluid transition, using an ultracold Bose gas of $\sim 4000$ atoms in optical lattices (Fig. 1c). Our main results are three-fold. First, we observe non-Gaussian order-parameter statistics in the critical regime near the transition point, thereby providing an experimental confirmation of renormalisation group predictions. Second, we show that the measured PDFs are well approximated by a phenomenological approach similar in spirit to Landau’s, which constitutes the basic framework to the traditional approach of phase transitions [14]. We additionally explain why this simple picture effectively describes our measurements. Third, we characterise the non-Gaussian nature of the statistics through high-order cumulants of the order parameter, revealing their large oscillations across the transition. We provide experimental evidence that these features are universal, and show in numerical simulations of homogeneous models that these oscillations obey critical scaling. The observation of cumulant oscillations stems from the unique level of diagnostics offered by quantum simulators; it represents a new perspective on the fundamental nature of critical fluctuations [6, 52, 53, 54].

Our experiment investigates the superfluid-to-normal fluid transition in a lattice gas of bosons. More specifically, we use an ultracold gas of metastable helium (${}^{4}\mathrm{He}^{*}$) atoms adiabatically loaded in the lowest band of a cubic optical lattice [13], realising the Bose-Hubbard model (BHM) [55]

[TABLE]

Here, $J$ and $U$ designate the tunnelling and on-site interaction energies respectively, $\hat{a}_{\bm{i}}^{{\dagger}}$ ($\hat{a}_{\bm{i}}$) is the bosonic creation (annihilation) operator on site $\bm{i}$, $\hat{n}_{\bm{i}}=\hat{a}_{\bm{i}}^{\dagger}\hat{a}_{\bm{i}}$ is the number operator and $\braket{\bm{i},\bm{j}}$ indicates nearest-neighbour lattice sites. The effect of the harmonic trapping potential is encoded in the spatial dependence of $V_{\bm{i}}$ (Methods). At zero temperature, the homogeneous BHM ($V_{\bm{i}}=0$) features a quantum phase transition between a superfluid state for weak interactions and a Mott insulating state for large interactions [56], occurring for a critical value $u_{c}$ of the ratio of the interaction and tunnelling energies $u=U/J$. At finite temperature, this transition acquires a classical character, i.e. it is governed by thermal fluctuations in a narrow region around a temperature-dependent critical interaction strength $u_{\rm c}(T)$ (dashed line in Fig. 1c). This transition is associated with a spontaneous breaking of the $U(1)$ symmetry and belongs to the universality class of the 3D classical XY model. The presence of a harmonic potential $V_{\bm{i}}$ implies that two different phases may co-exist, and leads to a modification of the critical value (red curve in Fig. 1c). However, this does not prevent the observation of critical scaling [49, 50, 57, 58, 59]. In our experiment, the total atom number is $N_{\mathrm{tot}}=$4.0(5)\text{\times}{10}^{3} (Methods), yielding $\braket{\hat{n}_{\bm{i}}}\approx$1 at the centre of the system; the control parameter $u$ is tuned between $u\approx$5 to $u\approx$70 by adjusting the optical lattice depth.

The order parameter of the superfluid-to-normal fluid transition is the condensate complex amplitude $\psi=\sqrt{N_{0}}\mathrm{e}^{\mathrm{i}\phi}$ [60, 61]. The ordered phase is associated with long-range order which manifests itself through a macroscopic population $N_{0}\gg 1$ of the zero quasi-momentum mode $\bm{q}=\bm{0}$, corresponding to the mode of the condensate. In the following we focus on the amplitude of the order parameter $|\psi|=\sqrt{N_{0}}$. Our detection method is ideally suited to monitor the full statistics of $|\psi|$ as it reconstructs the 3D momentum distribution atom-by-atom. It relies on a long time-of-flight expansion towards a pair of micro-channel plates, allowing to measure the momentum $\hbar\bm{k}$ of individual atoms with a detection efficiency of about $50\text{\,}\mathrm{\char 37\relax}$ (Fig. 2a,b). More details on the experimental apparatus can be found in our previous works [10, 62, 63] and in the Methods. We show in Fig. 2b exemplary single-shot experimental realisations. For $u<u_{\rm c}$, we measure the typical 3D diffraction pattern associated with a superfluid with long-range phase coherence. Deep in the Mott insulator regime ($u>u_{\rm c}$), where phase coherence is lost, the 3D momentum distribution is featureless. In the following, for each experimental shot we denote $N$ the total detected atom number in the the first Brillouin zone (blue cube in Fig. 2b corresponding to all momentum states $k_{i}<k_{\rm d}/2$, $i=x,y,z$) and $N_{0}=|\psi|^{2}$ the detected atom number in a sphere of radius $|\delta\bm{k}|=$0.03$k_{\rm d}$ centred on $\bm{k}=\bm{0}$ (Methods). Here, $k_{\rm d}=2\pi/d$ is the lattice momentum, with $d$ the lattice constant.

Our experiment is naturally subject to fluctuations of the atom number $N$, associated with the preparation of ultracold lattice gases. The probability distribution of such fluctuations, here obtained from about 1200 experimental repetitions, is represented in Fig. 2c in the vicinity of $u_{\rm c}$. A key asset of our detection method is that it offers the full counting statistics (FCS) of both $N$ and $N_{0}$ [46]. While the statistics of $N$ are well described by a bell-shaped (Gaussian) distribution, the statistics associated with $N_{0}=|\psi|^{2}$ (Fig. 2d) shows a completely different behaviour, with large relative fluctuations. Such a feature is characteristic of the phase transition between the superfluid state with sizeable $\braket{|\psi|}/\sqrt{N}$ and the normal state with $\braket{|\psi|}/\sqrt{N}\to 0$. For any value of $u$, the measured statistics of $N_{0}$ is slightly influenced by the fluctuations of the total number $N$ (see SI, section .1). Post-selecting on $N$ (see Fig. 2c) further reduces this effect.

We reconstruct the probability distribution function $P(|\psi|)$ of the amplitude of the order parameter $|\psi|$ from the measured FCS of $N_{0}$. The PDFs shown in Fig. 2e exhibit drastic differences depending on the value of $u$. In the ordered (superfluid) regime, the distribution has a maximum around $|\psi_{0}|\approx 10$, capturing the macroscopic population of the lowest energy (condensate) mode. Its shape is well matched by a Gaussian distribution, as expected from the Poissonian statistics for $N_{0}=|\psi|^{2}$ that result from sampling randomly a small fraction of the condensate atoms [46]. In the disordered (Mott) phase, the PDF is different, taking a maximum value for $|\psi_{0}|=0$. Once more, it is well adjusted by a Gaussian distribution, as expected from the expected thermal statistics [64, 46]. The PDF near the phase transition $u\approx 36$ strongly contrasts with both the ordered and disordered regimes. It features a relatively large width, with significant values in the range $0\leq|\psi|\lesssim 7$, capturing the large fluctuations of the order parameter near the phase transition. The overall shape of the distribution is peculiar and is distinct from Gaussian.

The behaviour of the order-parameter PDF near the transition is captured by introducing a function $\mathcal{F}(|\psi|)=-\ln P(|\psi|)$ playing the role of an effective thermodynamic potential. Close enough to the transition, it is well approximated by a low-order polynomial

[TABLE]

similar in spirit to Landau’s approach to phase transitions. As detailed in the Supplemental Information (SI), this thermodynamic potential is expected to efficiently describe the statistics of order-parameter amplitude regardless of the universality class. It has a non-trivial minimum for $a_{2}<0$ (ordered phase) and a minimum at $\psi=0$ for $a_{2}>0$ (disordered phase). Near the transition, $P(|\psi|)$ becomes highly non-Gaussian when $a_{2}^{2}/a_{4}\lesssim 1$.

In Fig. 3a, we plot the PDF $P(|\psi|)$ measured near the transition in the range $33\leq u\leq 40$, while Fig. 3b shows the corresponding potentials. We observe a continuous change from a distribution similar to that of the ordered phase deep in the superfluid regime ($u=25$), to a distribution closer to the disordered phase deep in the Mott regime ($u=55$), with flatter distributions in between. Equation (2) fits well our experimental data (solid lines in Fig. 3a,b). The fitted coefficient $a_{2}$ is shown in Fig. 3c. As anticipated, $a_{2}$ increases from negative to positive values as the system evolves from the ordered to the disordered phase. Empirically, we identify a critical value $u_{\rm c,0}=$37.1(2)$$ where $a_{2}$ is equal to zero. We stress that, away from the thermodynamic limit, $u_{\rm c,0}$ is different from $u_{\rm c}$, the critical point associated with universal scaling that we discuss below. This result is well-established in the phenomenology of finite-size systems [25] (see SI). We also show in Fig. 3c the expected asymptotic behaviours of $a_{2}$, for $u\ll u_{\rm c,0}$ and $u\gg u_{\rm c,0}$ (dashed lines). In the former case, the Poissonian statistics of $|\psi|^{2}$ yield $a_{2}=-1$, while in the latter case, where statistics are thermal, one finds $a_{2}=-\ln[\braket{N_{0}}/(1+\braket{N_{0}})]$.

Even though the PDF is well approximated by the exponential of a Landau-like free energy which is analytic in $|\psi|^{2}$, we emphasise here that we observe beyond-Landau — i.e. beyond mean-field — criticality. Our measurement of the average value $\braket{|\psi|}$ as a function of $u-u_{\rm c}$ confirms this picture (Fig. 3d). First, we find a critical exponent $\beta=0.35(1)$ compatible with that of the 3D XY universality class, $\beta_{\rm 3DXY}=0.3485$ [65], and distinct from the mean-field Landau prediction $\beta_{\mathrm{MF}}=0.5$. Second, a critical scaling for $\braket{|\psi|}$ is found even for small values of $u$, suggesting that the correlation length $\xi$ decreases slowly with $u_{c}-u$ and remains larger than the Ginzburg length $\ell_{\rm G}$ (see Fig. 1) in the ordered (superfluid) phase. An algebraic scaling is thus observed in the regime where the order-parameter PDFs are well captured by Gaussian statistics. This finding supports the idea that the range of the control parameter ($u$ here) for non-Gaussian statistics is narrower than that of the universal regime, defined by scaling laws (see Fig. 1a). Finally, as previously mentioned (see Fig. 1c), the harmonic trap leads to sizeable corrections to $u_{\rm c}$, as compared to predictions for a homogeneous system in the thermodynamic limit. The extracted value of $u_{\rm c}=$38.1(2)$$ agrees well with predictions based on Quantum Monte-Carlo calculations in a harmonic trap [66]. The fact that $u_{c}$ exceeds the critical value ($u\approx 29.3$ [67]) for the superfluid-Mott-insulator transition in a bulk 3D lattice with filling $n=1$ highlights the fact that the superfluid-to-normal transition we observe is occurring in the outer region of the trap, where the local filling is $n<1$ (see Methods and Fig. 5).

In view of the above observations, it might seem surprising that our experimental results align with the Landau-like theory. However, at fixed $u$, beyond-mean-field effects associated with the 3D XY universality class manifest themselves only in the tails of the order-parameter PDF, and correspond to rare events [68] (see SI, section .5). Experimentally, the typical number ($\sim 700$) of shots used to construct a PDF, although large, remains insufficient to resolve these rare events. Consequently, a single measured distribution does not allow to distinguish between the Landau-like potential introduced in Eq. (2) and the expected distribution obtained from 3D XY numerical simulations (SI, section .7).

Equivalently to the PDF, full statistical information of the order parameter is found in its moments, or cumulants, at all orders. To quantify explicitly the non-Gaussian nature of the PDFs, we focus on the evolution of its high-order cumulants across the transition. This methodology offers a dual advantage. First, cumulants are a standard tool in statistics for capturing non-Gaussian behaviour; the presence of any cumulant $\kappa_{n\geq 3}\neq 0$ is necessary and sufficient for non-Gaussian statistics. Second, as we shall see, the dependence of cumulants on the control parameter proves to be extremely rich, providing new and surprising insight into the critical behaviour.

The measured cumulants $\kappa_{n}(|\psi|)$ of the amplitude of the order parameter $|\psi|$ are shown in Fig. 4. Cumulants $\kappa_{1}=\braket{|\psi|}$ and $\kappa_{2}=\braket{|\psi|{}^{2}}-\braket{|\psi|}^{2}$ behave as expected from the phenomenology of a continuous phase transition; $\kappa_{1}$ is finite in the ordered phase ($u\ll u_{\rm c}$) and vanishes in the disordered phase ($u\gg u_{\rm c}$), while $\kappa_{2}$, which quantifies the fluctuations of the order parameter and can be associated with a susceptibility [16], exhibits a peak in the critical regime ($u\approx u_{\rm c}$).

The evolution of higher-order cumulants with the control parameter $u$ is instead more instructive. First, we find that cumulants of order $n\geq 3$ vanish in both the ordered ($u\ll u_{c}$) and disordered ($u\gg u_{c}$) phases, while they take large values near the phase transition $u\approx u_{\rm c}$. These variations unambiguously reveal the crossover of the order-parameter statistics from Gaussian to non-Gaussian across a continuous phase transition. A similar phenomenology has been theoretically proposed as a key signature to identify the phase transition from hadron gas to quark-gluon plasma in quantum chromodynamics(QCD) [69, 70, 71]. Second, we observe fast oscillations of the high-order cumulants as a function of $u$ near the phase transition. We discuss below the origin of this unexpected observation.

In Fig. 4a, we compare the measured cumulants with those obtained numerically (insets). In the numerical simulation for various models belonging to the 3D XY universality class, we find similar oscillations of the cumulants (SI). Furthermore, a finite-size scaling analysis performed on these models, which are cast on homogeneous lattices with periodic boundary conditions, demonstrates an excellent collapse of each cumulant onto a universal curve (see Fig. S7 in the SI). This proves that the oscillations of the order-parameter cumulants are a signature of universality across the phase transition. Whether such features can survive in a harmonic trap is not straightforward. Remarkably, our results suggest that the coexistence in the trap of different phases, which implies that the entire system is not uniformly critical, does not hinder the manifestation of universal behaviour through high-order cumulants. A possible interpretation for this observation is the fact that only the non-Gaussian critical region contributes to these high-order cumulants. Interestingly, while a qualitative agreement is found between the experiment and the simulations on homogeneous models, we observe quantitative differences for cumulants of order $n\geq 4$. These deviations can be interpreted as originating from the harmonic trap. Theoretically predicting the detailed behaviour of these cumulants in non-homogeneous models remains a challenge that this work motivates to address.

To find direct experimental evidence of universality in the behaviour of the high-order cumulants, we measure them while crossing the superfluid transition in a different configuration. More specifically, we decrease the total atom number (from $4.0(5)\text{\times}{10}^{3}$ to $2.6(3)\text{\times}{10}^{3}$) and slightly increase the entropy per atom. This leads to a substantial shift of the position of the transition, now observed at $u_{\rm c,0}=26.8(1)$, and compatible with critical behaviour developing at the trap center with nearly unit filling $n\approx 1$ (see Methods and Fig. 5). This shift highlights the different microscopic configuration of the quantum gas as it crosses the transition, whose universality class remains however unchanged. Crucially, the results in Fig. 4b show identical oscillations of the high-order cumulants of $|\psi|$ as those displayed in Fig. 4a. This observation strongly suggests that the measured oscillations are universal features of criticality in a harmonically trapped $U(1)$-symmetric system.

In conclusion, we measured the statistics of the amplitude of the order parameter across a continuous phase transition using ultracold lattice bosons. Our experiment is performed in a harmonic trap and realises an inhomogeneous system of mesoscopic size. Interestingly, we observe many of the emblematic predictions for homogeneous systems in the thermodynamic limit, including the change in the effective thermodynamic potential — similar in spirit to Landau theory — and the non-Gaussian critical regime.

In addition, we experimentally unveiled universal oscillations of the order-parameter cumulants across the transition. These oscillations are general features of criticality, which are not specific to the 3D XY universality class explored here. They have been observed numerically in the $O(N)$ universality classes [72] and stem from the shape of the critical PDF across the transition [27]. Furthermore, they are considered as a promising method to identify phase transitions in the QCD phase diagram [69, 73] and negative cumulants were recently reported in this context [74, 75]. The oscillations of the cumulants of order $n\geq 3$ are confined to a smaller region around the critical point as the system size increases. Hence they offer a tool to resolve the critical point with high resolution, which can be leveraged for metrological applications, and possibly extended to complex systems in general to identify critical transitions [76].

The standard approach of probing critical scaling is experimentally achievable across large dynamic ranges only in macroscopic systems where finite-size effects are negligible. Conversely, our measurement of the order-parameter statistics demonstrates a fully complementary approach that is well suited to mesoscopic systems. This approach, probing the universal properties of continuous phase transitions through their non-Gaussian critical regime, is very promising for several quantum-simulation platforms [35, 33, 34, 36, 37, 38]. In the future, probing systems with homogeneous densities and varying sizes would provide experimental validation of the universal finite-size scaling of high-order cumulants. Progress in theoretical predictions for inhomogeneous systems would furthermore provide crucial insight in our understanding of universality.

Methods

Production of low-entropy and low-number BECs. Our experiment uses standard laser cooling and evaporative cooling steps to produce a spin-polarized Bose-Einstein condensate (BEC) of $\sim$4\text{\times}{10}^{5}$$ metastable helium-4 (${}^{4}\mathrm{He}^{*}$) in the sub-level $m_{J}=1$ of the state $2^{3}$S1 [77]. The atom number is then reduced to $4.0(5)\text{\times}{10}^{3}$ using state-dependent inelastic losses: we use two-photon Raman transitions to transfer a controlled fraction of the atoms to the sub-level $m_{J}=0$ where two-body Penning ionization leads to losses. The small BEC is then loaded adiabatically in the cubic optical lattice. In the 3D lattice, the trapping frequency of the harmonic potential $V_{\bm{i}}$ is $\omega/2\pi=140(5)\sqrt{V_{0}}~$Hz, with $V_{0}$ the lattice amplitude in unit of the recoil energy $E_{\rm R}=h^{2}/8md^{2}$, and $m$ the mass of an atom. A thermometry based on the comparison with ab initio Monte-Carlo calculations [13] yields a temperature in the range $k_{\rm B}T/J\sim 1.5-2$. This is compatible with an entropy per atom of $S/N\sim$0.4(1)$k_{\rm B}$ extracted from the measured condensate fraction before loading the BEC in the optical lattices [13].

Volume used to compute the statistics of $|\psi|$. In translationally invariant lattice systems where the quasi-momentum ${\bm{q}}$ is a good quantum number, the number of condensate atoms $N_{0}$ is equal to the atom number in the mode ${\bm{q}}={\bm{0}}$. In a finite-size and harmonically-trapped gas, $N_{0}$ remains well approximated by the atom number at ${\bm{q}}={\bm{0}}$. In this work, we probe the gas in the momentum basis where the modes located at $\bm{n}k_{\rm d}$, with $\bm{n}\in\mathbb{Z}^{3}$, are copies of the mode ${\bm{q}}={\bm{0}}$. We use a sphere of radius $0.03k_{\rm d}$ centred on the central peak (${\bm{n}}={\bm{0}}$) or first-order diffraction peaks ($\|\bm{n}\|=1$) to measure $N_{0}$ in each shot. This radius is unambiguously determined from the width of the bosonic bunching signal which is set by the size of one mode [64, 78]. As shown in the SI, the exact value of the radius we use in the analysis does not affect our observations.

Bootstrap procedure. The errorbars on experimental cumulants, probability distribution functions and fit parameters are $68\%$ confidence intervals obtained from a bootstrap method. 500 pseudosamples are generated from sampling the experimental shots with replacement. The observables are computed for each pseudosample. The resulting histograms take a typical bell-shape curve, from which the mean, the $16^{\mathrm{th}}$ and $84^{\mathrm{th}}$ percentiles are extracted. The percentiles provide an estimate of the statistical uncertainty on the mean.

Trapped atomic densities at different total atom number. In Fig. 4, we compare the order-parameter cumulants measured in lattice gases with different total atom numbers and slightly different entropy per atom.

As discussed in the main text, the difference in the measured critical values $u_{\rm c,0}$ highlights the distinct configurations of in-trap densities close to $u_{\rm c,0}$: at small atom number (Fig. 4b) the condensate develops at the trap center; in contrast the condensate appears in the outer region of the trap at large atom number (Fig. 4a). These differences are confirmed by numerical calculations of the harmonically trapped Bose-Hubbard Hamiltonian using quantum Monte Carlo (QMC) [13], see Fig. 5. There we contrast the local density $n_{i}=\braket{\hat{a}_{i}^{\dagger}\hat{a}_{i}}$ with the condensate density $n^{(0)}_{i}=\frac{1}{V}\sum_{j}\braket{\hat{a}_{i}^{\dagger}\hat{a}_{j}}$, where $V$ is the volume of the simulation box. All the parameters used in these numerical calculations are calibrated from the experiment, except for the temperature, which is estimated by matching the momentum distribution calculated numerically via QMC with the one measured experimentally, as described in [13].

Acknowledgements — We thank I. Balog, A. Browaeys, Z. Hadzibabic 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 QuanTiP, 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 numbers ANR-22-PETQ-0004). A.R. has benefited from the financial support of the Grant No. ANR-24-CE30-6695 FUSIoN. Numerical simulations were conducted on the CPBsmn cluster at the ENS of Lyon.

.1 Effect of the total atom number on the order-parameter statistics

The statistics of the order parameter could in principle be affected by technical noise, i.e., by the shot-to-shot fluctuations of the total atom number. To mitigate this contribution we post-select experimental shots on the total atom number $N=N_{\rm FBZ}$ in the first Brillouin zone. Our data-taking procedure consists of measuring at different values of $u$ in consecutive shots, so that the atom number statistics is guaranteed to be independent of $u$. For instance, we observe in the experiment the same variation of $N$ with $u$ as obtained in Quantum Monte Carlo simulations of the Bose-Hubbard Hamiltonian in the presence of an harmonic trap. These variations simply reflect the change in the shape of Wannier functions with $u$, which in turn affects the number of atoms measured in the first Brillouin zone. Note that $N$ is slightly under-estimated in the superfluid regime as the central peak saturates the detector due to its high atomic density. The scaling analysis of Fig. 3d, in particular, is performed on the diffracted peaks at $\bm{k}=\bm{k_{\rm d}}$ to avoid the saturation effects. We discuss this in more detail in section .3.

Fig. S1 presents the cumulants of the order parameter obtained when varying the amplitude $\Delta N$ of the fluctuations around a fixed mean value $\braket{N}$ in the post-selection on the total atom number. This illustrates that our results are only marginally affected by the fluctuation amplitude $\Delta N$.

.2 Effect of the measurement volume on the order-parameter statistics

For each experimental shot, we count the atom number $N_{0}$ in the condensed mode (from which we deduce $|\psi|$), in a spherical volume centred on $\bm{k}=\bm{0}$. We choose the radius of this sphere to match the size of a mode in momentum space, as obtained from the width of the bunching peak [62]. We show in Fig. S2 that small variations of the radius around this value does not modify the measured cumulants significantly. The small variations are instead expected: the increase of $\kappa_{1}$ with volume reflects the increase in the atom number; the contrast of higher-order cumulants decreases with volume as other modes are mixed with the condensate mode.

.3 Momentum vs quasi-momentum

The condensate order parameter $|\psi|\sim\sqrt{N_{\bm{q}=\bm{0}}}$ is associated with the mode with quasi-momentum $\bm{q}=\bm{0}$. Time-of-flight measurements as performed in this work give instead access to populations in momentum modes, the in-trap quasi-momenta being projected onto the plane wave basis as the lattice is abruptly switched off. The quasi-momentum mode $\bm{q}=\bm{0}$ is replicated at the reciprocal lattice vectors $k_{\rm d}\bm{n}$, with $k_{\rm d}=2\pi/d$ ($d=$775\text{,}\mathrm{nm}$$ is the lattice constant) and $\bm{n}\in\mathbb{Z}^{3}$. This mechanism is analogous to a beam-splitter with multiple outputs and leads to a multinomial noise of each of the output modes $\bm{k}=\bm{n}k_{\rm d}$.

In addition, the average populations of the momentum peaks $\bm{k}=\bm{n}k_{\rm d}$ are not all proportional to that of the quasi-momentum mode $\bm{q}=\bm{0}$. This effect results from the variation of the Wannier functions with $u$, whose envelop in momentum space modifies the amplitudes of the diffracted peaks. Numerically, we find that the average population of the first-order peaks $|\bm{n}|=1$ is proportional to that of the quasi-momentum mode $\bm{q}=\bm{0}$, independently of $u$. In other words, the ratio $\braket{N_{\bm{k}=\pm\bm{k}_{\mathrm{d}}}}/\braket{N_{\bm{q}=\bm{0}}}$ is constant. In contrast, the ratio $\braket{N_{\bm{k}=\bm{0}}}/\braket{N_{\bm{q}=\bm{0}}}$ varies with $u$ similarly to $\braket{N}$. As a result, one must consider the normalised variable $N_{\bm{k}=\bm{0}}/\braket{N}$ to obtain a variation with $u$ that faithfully reflects that expected for the variable $N_{\bm{q}=\bm{0}}$. To verify this assertion, we compare the cumulants of $\sqrt{N_{\bm{k}=\bm{0}}}/\sqrt{\braket{N}}$ with those of $\sqrt{N_{\bm{k}=\pm\bm{k}_{\mathrm{d}}}}$ in Fig. S3.

In Fig. S3, we plot the cumulants normalised to the value of their average $\kappa_{1,\rm c}$ measured at $u=u_{\rm c,0}$. The high-order cumulants ($n\geq 3$) are almost identical. The first and second order cumulants slightly differ instead, with differences which are well understood. The lower values of $\kappa_{1}/\kappa_{1,\rm c}$ of the central peak $\bm{k}=\bm{0}$ are due to a saturation of the Micro-Channel Plate detector occurring at small values of $u$ where peak densities are the largest. The second order cumulants $\kappa_{2}/\kappa_{1,\rm c}^{2}$ differ mainly due to stronger effect of the multinomial process for the diffracted peaks at $\pm\bm{k}_{\rm d}$, with respect to that occurring for the peak $\bm{k}=\bm{0}$.

The analysis of the cumulants of the central and first-order peaks shows how reliable the observation of oscillating high-order cumulants near the transition point is. Furthermore, it validates the study of the cumulants of the variable $\sqrt{N_{\bm{k}=\bm{0}}}/\sqrt{\braket{N_{\rm FBZ}}}$ in order to quantitatively reflect the properties of the quasi-momentum mode $\bm{q}=\bm{0}$. All the cumulants shown in the main text and in the SI for the central peak $\bm{k}=\bm{0}$ refer to this variable.

.4 Critical behaviour of $\braket{|\psi|}$

Here we complement the analysis of the critical scaling of the order-parameter amplitude, shown in Fig. 3d, by fitting the data over different ranges of $u$ values. In particular, we fit the data over windows $u_{\rm min}\leq u\leq 37$ with a varying minimum value $u_{\rm min}$. The results are illustrated in Fig. S4.

For $u_{\rm min}=5$, corresponding to using all the available data with $u\leq 37$, we observe a small discrepancy between the experimental value of $\beta$ and the expected value in the 3D XY universality class. For $u_{\rm min}=15,20$ or $25$, the measured and 3D XY critical exponents are compatible, suggesting that the critical regime, described by scaling and universality, extends at least down to $u=15$. For a larger value $u_{\rm min}=30$, the measured exponent $\beta$ disagrees significantly with the 3D XY universality class, although with much larger errorbars. This discrepancy occurs in the range of $u$ where high-order cumulants become non-zero, signalling non-Gaussian statistics associated with finite-size effect ($\xi$ is of the order of $L$). It could therefore originate from finite-size effects that are known to affect critical scalings.

.5 Critical probability distribution functions

In this section, we discuss the Probability Distribution Functions (PDFs) of critical systems and contrast them with those obtained from the Landau theory of critical phenomena. Even though the system of interest is quantum, the transition we are studying here occurs at finite temperature. Therefore, it can be described by an effective classical theory when focusing on large-scale properties [79]. We therefore focus on the PDF of classical critical systems.

Assume that the microscopic system is described in terms of a complex field $\phi(\bm{r})$, which can be defined on a lattice or in the continuum. Its thermodynamics is described by a Hamiltonian $\mathcal{H}[\phi(\bm{r})]$ and the Boltzmann weight

[TABLE]

which gives the probability of a configuration of the microscopic field. We assume that a control parameter $u$ drives the transition at some critical value $u_{\rm c}$.

The physics of the system is usually more conveniently expressed in terms of a coarse-grained order-parameter field $\psi(\bm{r})$ (which can, for instance, correspond to the spatial average of $\phi(\bm{r})$ over a mesoscopic region). The probability distribution of a configuration $\psi(\bm{r})$ is then given by a Ginzburg-Landau (GL) functional

[TABLE]

with

[TABLE]

Working in the vicinity of the phase transition, we can assume that $\psi(\bm{r})$ has a small amplitude and varies slowly in space, so that it is legitimate to expand the functional in powers of the field and its gradients. In the case of systems with $U(1)$ symmetry, as in the case of the Bose-Hubbard model, the order-parameter field $\psi$ is complex. Since the Ginzburg-Landau functional is an analytic function of $\psi(\bm{r})$ and $\psi^{*}(\bm{r})$ and invariant under $U(1)$ transformations of the field, only even powers of $|\psi(\bm{r})|$ and $|\nabla\psi(\bm{r})|$ are allowed in its expansion.

Focusing on the spatially-averaged/zero-momentum field $\psi=V^{-1}\int\mathrm{d}\bm{r}\,\psi(\bm{r})$ (with $V$ the volume of the system), its statistical average $\braket{\psi}$, taken using $P_{\rm GL}[\psi(\bm{r})]$, corresponds to the order parameter that discriminates between different phases. The statistical properties of $\psi$ are given in terms of its PDF $P(\psi)$, which is the marginal distribution of $P_{\rm GL}$ after integrating out the modes with non-zero momenta. We will use the common abuse of language and refer to $\psi$ as the order parameter.

Since all momentum modes interact with each other through the non-linear coupling (approximated by a quartic term in Eq. (S3)), the computation of $P(\psi)$ from $P_{\rm GL}[\psi(\bm{r})]$ is a priori highly non-trivial. Landau’s theory of phase transitions amounts to performing a mean-field calculation and approximating $P(\psi)$ by

[TABLE]

with $\mathcal{F}_{\rm L}(\psi)=\tilde{a}_{0}+\tilde{a}_{2}|\psi|^{2}+\tilde{a}_{4}|\psi|^{4}$ the Landau free energy ($\tilde{a}_{0}$ is related to the normalization of $P_{\rm L}(\psi)$). In the thermodynamic limit, $V\to\infty$, one finds that $|\braket{\psi}_{\rm L}|^{2}\propto-\tilde{a}_{2}/\tilde{a}_{4}$ if $\tilde{a}_{2}<0$ and $|\braket{\psi}_{\rm L}|^{2}=0$ otherwise, recovering the standard Landau phenomenology. Note that $\tilde{a}_{2}$ depends on $u$, but will generically vanish (linearly) at a value $u=u_{\rm L}\neq u_{\rm c}$ where $u_{\rm c}$ is the exact critical point. In this approximation and for a finite system, $\psi$ will have non-Gaussian statistics for a sufficiently small $\tilde{a}_{2}$ (such that $\tilde{a}^{2}_{2}V/\tilde{a}_{4}\lesssim 1$).

However, it is well-known that in dimensions lower than four, the above mean-field picture breaks down. The coupling between small and large scales (or, alternatively, large- and small-momentum modes of $\psi(\bm{r})$) strongly renormalises the PDF, which, near the critical point, takes the scaling form

[TABLE]

Here $L$ is the linear system size ($V=L^{3}$), $\beta$ and $\nu$ are the order-parameter and correlation length critical exponents, respectively, and $\delta u=u-u_{\rm c}$ is the distance to the critical point.

The scaling function $\tilde{P}(x,y)$ is universal, in the sense that it does not depend on microscopic details of the system. These details enter only in so-called non-universal amplitudes, which amount to fixing the units of the variables $x$ and $y$ (for instance, the amplitude of the field $\psi$, and the link between $\delta u$ and the correlation length). The scaling function $\tilde{P}$ depends on the symmetry of the order parameter (here $U(1)$), the dimensionality of space (here $d=3$), and, importantly, the boundary conditions [21]. The critical exponents also depend on the former two, but not on the latter.

Independently of the universality class, the shape of $\tilde{P}(x,y)$ is qualitatively of the following form [25, 27]: i) for large $x$, it takes the form of a compressed exponential, $e^{-c(y)x^{\delta+1}}$, with $\delta$ the isothermal critical exponent (that depends on $\beta$ and $\nu$, with $\delta\simeq 4.78$ in the case of the 3D XY universality class); ii) at small $x$, it can be approximated by $e^{-a_{2}x^{2}-a_{4}x^{4}-\ldots}$, where $a_{2}$ depends on $y$ and thus on $\delta u$. Importantly, $a_{2}$ changes sign at some $y_{0}\neq 0$ that depends on the boundary conditions. At fixed system size, this change of sign corresponds to $u=u_{\rm c,0}\equiv u_{\rm c}+y_{0}L^{-1/\nu}$, implying that $u_{\rm c,0}$ is generically different from $u_{\rm c}$ (but $u_{\rm c,0}\to u_{\rm c}$ as $L\to\infty$). On the other hand, we expect $a_{4}$ to weakly depend on $\delta u$. The region where $\psi$ assumes small values is therefore well captured by a naive Landau form Eq. (S4) with $\tilde{a}_{i}\to a_{i}$ ($i=2,4$), even though the parameters and critical scaling are very different. Only the behaviour in the tail, through the exponent $\delta$, discriminates between the different universality classes.

.6 Fitting of the probability distribution functions

In the main text, the measured probability distribution functions are fitted with the ansatz $\mathrm{e}^{-(a_{0}+a_{2}|\psi|^{2}+a_{4}|\psi|^{4})}$. We show in Fig. S5b the fitted values for the coefficient $a_{4}$. The form used to fit the experimental data is expected to be valid close to the transition point and whether it remains valid away the transition is not granted.

In the specific case of the Mott insulator phase whose statistics in momentum space is well known [64], we expect a distribution whose quartic term is equal to zero ($a_{4}=0$). As a matter of fact, when $u\geq 40$ the shape of the curves in the disordered phase does not allow to easily distinguish between a quadratic and a quartic form, leading to large correlations between the fitted $a_{2}$ and $a_{4}$ as well as large errorbars displayed in Fig. S5a (blue disks, see also Fig. 3c).

Figure S5a (grey squares) shows the fit results obtained by enforcing $a_{4}=0$ for $u\geq 40$. Under this assumption, we note that the fitted values for $a_{2}$ deep in the normal (Mott) phase ($u\gg u_{\rm c}$) better match the expected value from the Mott statistics (dashed line).

In all cases, the coefficient $a_{0}$ is fixed by normalisation of the probability by

[TABLE]

where the $|\psi|_{i}$ are the values sampled by $\sqrt{N_{0}}$ experimentally, i.e. $|\psi|_{i}^{2}$ takes integer values.

.7 Comparison of experimental distributions with 3D XY simulations and Landau predictions

In this section, we compare the experimental PDF of the order parameter with the PDF predicted by a Landau-like approximation, and by Monte Carlo (MC) simulations of the homogeneous classical 3D XY model with either free or periodic boundary conditions. The Landau-like approximation is $P(|\psi|=x)\propto e^{-(a_{2}x^{2}+a_{4}x^{4})}$. The MC simulations are performed using Wolff’s algorithm, sampling more than a million order-parameter configurations. We have taken a sufficiently large system to ensure convergence of the PDFs with size. They obey the scaling form and behaviour discussed in Sec. .5. Note in particular that even though their shapes look qualitatively similar, they are in fact different, especially when compared at the same value of $L^{1/\nu}\delta u$ (with $\delta u$ being related to the inter-spin coupling in this context).

Fig. S6 shows the experimental data across the transition, compared to the best fit to a Landau-like distribution, and to the PDF of the 3D XY model at the distance to the critical point that gives the best overlap with the experimental PDF. We observe that the three PDFs are quite similar, and that the differences between them cannot be distinguished experimentally with our signal-to-noise ratio (using $\sim 1000$ shots). Given the finite statistics of the experiments, in particular the lack of data for rare events, we find that the Landau-like approximation of the PDF works very well. The compressed exponential behaviour expected for rare events in the exact distribution only appears in the tail of the distribution itself, and its observation would require in principle a much larger number of experimental shots than that used in this work. The finite size of our sample introduces another fundamental limit to the reconstruction of the tail of the distribution.

The three theoretical predictions shown Fig. S6 are continuous PDFs, while we measure discrete PDFs associated with discrete atom numbers. This leads to an issue regarding the normalization of the PDF when the bin size (=1) is not small compared to the variations of the PDF. To solve this, we normalize the continuous PDFs $f$ such that,

[TABLE]

where $\psi_{i}$ are the bins sampled by $\sqrt{N_{0}}$ experimentally.

.8 Finite-size scaling of 3D quantum rotors.

To reconstruct the non-Gaussian statistics of critical fluctuations in a quantum-mechanical model exhibiting the superfluid-to-normal-gas transition, we have performed simulations of the 3D quantum-rotor model. The quantum-rotor model is the limit of the Bose-Hubbard model for large, integer filling $n_{\rm QR}\gg 1$ [80]. Decomposing the bosonic operator in terms of amplitude and phase, $\hat{a}_{i}=e^{i\hat{\phi}_{i}}\sqrt{\hat{n}_{i}}$, with $[\hat{n}_{i},\hat{\phi}_{i}]=i$, and rewriting the density as $\hat{n}_{i}=n_{\rm QR}+\delta n_{i}$, if $\sqrt{\langle(\delta n_{i})^{2}\rangle}\ll n_{\rm QR}$ then one can write

[TABLE]

upon neglecting density fluctuations. Density-density interactions introduce a quantum kinetic term in the Hamiltonian, since $n_{i}^{2}=-\frac{\partial^{2}}{\partial\phi_{i}^{2}}$. At integer filling, all linear terms in the density drop from the Hamiltonian, leading to the quantum rotor Hamiltonian

[TABLE]

As done in Ref. [63], the momentum distribution of the quantum rotor model can be sampled via a coherent-state path-integral approach [80, 81, 82], leading to the possibility of reconstructing the FCS of the condensate population $N_{0}=\hat{a}_{0}^{\dagger}\hat{a}_{0}$.

At finite temperature, the phase transition is in the 3D XY universality class. From Eq. (S5), we expect the cumulants of the order parameter to obey the finite-size scaling law

[TABLE]

with $u=U/(2Jn_{\rm QR})$ the coupling constant, $\delta\propto(u/u_{\rm c}-1)\times L^{1/\nu}$, $u_{\rm c}$ the (temperature-dependent) critical coupling, and the critical exponents $\beta$ and $\nu$ are taken as those of the 3D XY model, $\nu\simeq 0.6718$ and $\beta\simeq 0.3485$. For a given temperature, the critical value $u_{\rm c}$ is found by finding the point where all the different-$L$ curves $L^{\beta/\nu}\kappa_{1}(u,L)$ cross, since $L^{\beta/\nu}\kappa_{1}(u_{\rm c},L)=f_{1}(0)$ for all $L$. Then, for a given temperature and $n$, the curves $L^{n\beta/\nu}\kappa_{n}(u,L)$ plotted as a function of $\delta$ are expected to collapse for all $L$ on a single curve.

However, to obtain the collapse of curves taken at different temperatures, we need to fix two non-universal amplitudes, associated with the amplitude of the order parameter and to a length scale. Note that $\kappa_{n}$ has the dimension of the order parameter to the power $n$, which implies that $f_{n}$ has the dimension $\textrm{(order parameter)}^{n}\textrm{(length)}^{n\beta/\nu}$, and that $\delta$ has the dimension of a length to the power $1/\nu$. These two amplitudes could be obtained, for instance, from $f_{1}(0)$ and $f_{1}^{\prime}(0)$, and using them to rescale all cumulants and $\delta$ appropriately.

In practice, it is more convenient to define $\tilde{\kappa}_{n}=\kappa_{n}(u,L)/\kappa_{1}(u_{\rm c},L)^{n}$ and find $\alpha(T)$ such that $\tilde{\kappa}_{1}$ plotted as a function of $\delta=\alpha(T)(u/u_{\rm c}-1)\times L^{1/\nu}$ collapse on a single curve for all $T$. Then, all $\tilde{\kappa}_{n}$ plotted as a function of their corresponding $\delta$ also collapse on a single curve (which depends on $n$).

The curves shown in the inset of Fig. 4 correspond to the universal curves shown in Fig. S7. The lines correspond to the average of the universal curves over the different sizes and temperatures. The shaded areas represent the dispersion of the different curves obtained with this method.