Evidences for interaction-induced Haldane fractional exclusion statistics in one and higher dimensions

Haldane fractional exclusion statistics (FES) has a long history of intense studies, but its realization in physical systems is rare. Here we study repulsively interacting Bose gases at and near a quantum critical point, and find evidences that such strongly correlated gases obey simple non-mutual FES over a wide range of interaction strengths in both one and two dimensions. Based on exact solutions in one dimension, quantum Monte Carlo simulations and experiments in both dimensions, we show that the thermodynamic properties of these interacting gases, including entropy per particle, density and pressure, are essentially equivalent to those of non-interacting particles with FES. Accordingly, we establish a simple interaction-to-FES mapping that reveals the statistical nature of particle-hole symmetry breaking induced by interaction in such quantum many-body systems. Whereas strongly interacting Bose gases reach full fermionization in one dimension, they exhibit incomplete fermionization in two dimensions. Our results open a route to understanding correlated interacting systems via non-interacting particles with FES in arbitrary dimensions.

chanics. However, they are not the only possible forms of quantum statistics [1]. In two dimensions (2D), anyonic excitations can carry fractional charges and obey fractional statistics [2][3][4][5][6][7]. To extend the concept of fractional statistics, Haldane generalized the Pauli exclusion principle and formulated a theory of fractional exclusion statistics (FES) that continuously interpolates between the Bose and Fermi statistics in arbitrary spatial dimensions [8]. This theory breaks particle-hole symmetry [9] and defines a FES parameter g αβ by counting how much the dimensionality d α of the Hilbert space of available single-particle states, namely the "number of holes" N h,α of species α decreases as particles of various species β are added to a system of fixed size and boundary conditions: where α and β are "labels of species" consisting of a certain set of quantum numbers (such as the quasi-momentum), N P,β is the particle number in species β, and g αβ is independent of the particle numbers. Fig. 1(a) illustrates a simplified example of non-mutual FES with g αβ = gδ αβ and D α = max{d α }, where Bose and Fermi statistics correspond to g = 0 and 1, respectively.
The statistical distribution of particles in an ideal gas with FES can be derived via the standard methods in statistical mechanics [10,11].
Haldane's FES approach reveals the statistical nature of a physical system with respect to its energy spectrum rather than the exchange statistics of wave functions. Therefore, it applies to generic quantum matters regardless of whether the constituent particles are interacting or not. As a result, FES provides a powerful approach for studying interacting quantum many-body systems, and has found realizations in a few physical systems. In one dimension (1D), the Calogero-Sutherland model of particles interacting through a 1/r 2 potential [12][13][14][15], Lieb-Liniger Bose gases [15,16], and anyonic gases with delta-function interaction [17,18] have been exactly mapped onto ideal gases with FES. In three dimensions, FES was assumed to be valid and used to analyze the equation of states of unitary Fermi gases [19]. On the other hand, it remains challenging to find evidences for FES in generic interacting quantum systems with varied spatial dimensionality and interaction strength.
In this letter, we consider repulsively interacting Bose gases at and near a quantum critical point in 1D and 2D. Under zero temperature, a system undergoes a quantum phase transition from a vacuum to a quantum liquid when the chemical potential µ exceeds a critical value µ c ( Fig. 1(b)). Here "quantum liquid" denotes Tomonaga-Luttinger liquid (TLL) [20] in 1D or su- perfluid in 2D [21]. We study non-mutual Haldane FES in such systems and show the relation between the interaction strength and FES parameter g ( Fig. 1(c)).
We report evidences for interaction-induced FES in the thermodynamic properties of these Bose gases based on exact solutions in 1D and high-precision quantum Monte Carlo (QMC) simulations in 2D. Our numerical data are confirmed by existing experiments [20][21][22][23][24][25]. We establish a one-toone mapping between a transformed interaction parameter C tr and the FES parameter g over a wide range of interaction strengths. Under this mapping, we observe agreements on the entropy per particle, density, and pressure between interacting Bose and non-interacting FES systems at and near the quantum critical point. In 1D, interaction drives the system into a full fermionization with g max,1D = 1, whereas in 2D, g max,2D = 0.432 (14) reveals an incomplete fermionization.
Here we formulate non-interacting particles with non-mutual FES parameter g. The occupation number f in a state with energy is given by [10] f ( ) = 1 w(ζ) + g Refs. [20,22]). (b) Power-law scaling of S c /N with respect to C tr . Dotted line denotes the fermionization limit where T is the temperature and µ the chemical potential. The number density and energy density are given by n = G( )f ( )d and e = G( )f ( ) d , where the density of states per volume is given by G( ) = 1 2π √ in 1D and 1 4π in 2D for non-relativistic particles ( = k 2 ; k is the momentum). In this work, we set 2m = k B = = 1, where m is the particle mass, k B is the Boltzmann constant, and is the reduced Planck constant. We aim to search for such simple non-mutual FES in strongly correlated matters of ultracold atoms.
The 1D repulsively interacting Bose gases (with no inelastic losses) are described by the Hamil- where c is the interaction strength and N is the number of particles. In its dilution limit, the discrete 1D Bose-Hubbard model used in QMC simulations relates to Eq. 3 via where U and t are the onsite interaction and tunneling parameters, respectively.
We show S c /N increases with the growth of a scaled interaction strengthc = c/ √ T ( Fig. 2(a)).
Atc → ∞, it reaches A ∞,1D ≈ 1.89738 (dotted line) [26] that exactly matches the S c /N of noninteracting fermions [27] (corresponding to g = 1). Our TBA solutions agree with data extracted from experiments performed by the Kaiserslautern group [22] and the USTC group [20], and agree with our 1D QMC simulations [26].
We further observe that S c /N obeys an empirical power-law scaling ( Fig. 2 with respect to a transformed interaction parameter C tr : where β 1D = 0.298(1) andc 1 = 0.772(5) in 1D are determined by a two-parameter fit. The fit agrees with the TBA data within 1% over a large range of interaction strengths (0.002 <c < ∞).

Equation 6
is inspired by a Ginzburg-Landau theory [28] for 2D superfluid [24]. The dependence of S c /N on interaction, observed in a previous experiment [21], is accurately described here as a power-law scaling (Eq. 5) that accordingly signals interaction-induced FES, as explained below.
We find a simple and explicit interaction-to-FES mapping by comparing Eq. 5 with the behavior of non-interacting particles with non-mutual FES. The TBA equation is a consequence of breaking particle-hole symmetry in excitations determined by the Bethe Ansatz equations [29]. Such particlehole symmetry breaking can in general be quantified by momentum-dependent mutual FES [15], and can be described by non-mutual FES in strongly interacting systems [29]. At and near a quantum critical point where the correlation length is large, the underlying FES physics can be greatly simplified into non-mutual FES. To demonstrate this point, we compute the "critical" entropy per particle (at µ c = 0), S c,FES /N , of a non-interacting gas with a non-mutual FES parameter g (Fig. 2(c), blue curve), and find that S c,FES /N exhibits an approximate power-law scaling with respect to g: S c,FES /N = A ∞,1D g β FES,1D , with β FES,1D = 0.298(2) fitted for 0.05 < g ≤ 1. This second power-law scaling and Eq. 5 agree very well, and the corresponding numerical data agree within 4% (Fig. 2(c)), which strongly suggests a one-to-one mapping between an interacting Bose gas with C tr and a non-interacting FES gas with g: with g max = 1 in 1D.
The agreements are within 15% forñ c and 8% forp c . The overall good agreements on S c /N , n c , andp c provide a comprehensive set of evidences for interaction-induced FES in a 1D Bose gas at a quantum critical point. Here, g max = 1 corresponds to a full fermionization of a 1D Bose gas atc = ∞, which was predicted and observed for quantum gases in the Tonks-Girardeau regime [16,[30][31][32][33].
Having established Eq. 7 in 1D, we now investigate whether interaction-induced FES [35][36][37] exists in 2D Bose gases at and near a quantum critical point. Based on QMC simulations [38,39], we study a 2D Bose-Hubbard lattice gas that has a vacuum-to-superfluid quantum phase transition at µ c = −4t [21]. The Bose-Hubbard model with no loss is described by a Hamiltonian: whereb † i andb i are the creation and annihilation operators at site i,n i =b † ib i , and < i, j > runs over all nearest neighboring sites. For this 2D lattice gas, we identify a scaled interaction strength c 2D = U/(2t) [26] that is the lattice-gas equivalence [21,24] of the interaction parameter √ 8πa/l z for a weakly interacting 2D Bose gas without lattices, where a is the scattering length and l z is an oscillator length [25].
To obtain physical properties that are insensitive to the lattice structure, we perform QMC simulations for eachc 2D at a series of temperatures down to T = 0.1t. We extract scaled thermodynamic quantities S c /N ,ñ c,2D = n c /T , andp c,2D = p c /T 2 at each T , and then perform extrapolation towards T = 0 for each quantity [26]. We test our extrapolation protocol by studying a 1D Bose-Hubbard system [26] and find excellent agreements with solutions to the TBA equation (Fig. 2). and experiments [21,23]. (b) Power-law scaling of S c /N with respect to C tr . (c), (d), (e) Evidences for equivalence between a 2D interacting Bose gas and 2D non-interacting particles with non-mutual FES, under the mapping g = g max,2D C tr with g max,2D = 0.432 (14), regarding three thermodynamic observables: (c) critical entropy per particle S c /N ; (d) scaled critical densityñ c ; (e) scaled critical pressurep c , with the dotted line denoting the non-interacting boson limitp c0 = π 24 [34]. Our results agree with existing experiments [21,23,24]. Horizontal and vertical gray bands mark A ∞,2D and g max,2D , respectively.
To provide a statistical interpretation for Eqs. 6 and 7 based on particle-hole symmetry breaking analysis [9], we rewrite these two equations as follows: On the left hand side, the numerator equals the dimensionality of Hilbert space occupied by one single particle; the denominator equals the dimensionality of Hilbert space that is "unoccupied under interaction strengthc" but still "occupiable under infinite interaction" by one single particle.
Our work empirically validates Eqs. 6 and 7 for both 1D and 2D systems at and near a quantum critical point, and hence reveals a phenomenological relation that the dimensionality ratio of the above occupied / unoccupied-but-occupiable Hilbert space is proportional to the scaled interaction strengthc.
Finally, we further explore the scope of application of our non-mutual FES mapping formula, Eq. 7, by comparing the scaled equations of stateñ(μ) =ñ µ−µc T of interacting Bose gases with those of non-interacting particles with FES in a finite range ofμ besides the quantum critical pointμ c . With no additional adjustable parameters, a strongly interacting 1D Bose gas with c 1D = 100 shows excellent equivalence to 1D non-interacting particles with g = 0.992 ( Fig. 4(a)).
As interaction weakens, the equivalence at µ ≤ µ c is still good, whereas deviations become more significant asμ exceeds 0. Here we presentñ because under the samec andμ, the relative deviations forp and S/N are smaller than that forñ [26]. Fig. 4(b) shows similar condition of equivalence between 2D interacting Bose gases and 2D ideal gases with FES. We attribute the deviations at positive finiteμ (Fig. 4), as well as the residual small deviations forμ ≤ 0 (Figs. 2 ∼ 4), primarily to the need of including more complex mutual FES effects [15,26], which is subject to future research.
To conclude, we find strong evidences for interaction-induced non-mutual Haldane fractional exclusion statistics in both 1D and 2D Bose gases at and near the vacuum-to-quantum-liquid transition. Our unified, non-perturbative mapping approach can be generalized by including mutual where ε(k) is called "dressed energy", k is quasi-momentum, µ is chemical potential, T is temperature and c = −2/a 1D with a 1D being 1D scattering length. Thus the grand thermodynamic potential of unit length, namely, pressure p can be obtained by Other thermodynamic properties, such as particle density n, entropy density s are given by the derivatives of the pressure, namely, For our convenience in analysis of critical phenomenon, we define the following dimensionless and the dimensionless thermodynamic properties Consequently, the dimensionless Yang-Yang TBAE is given bỹ This serves as an equation of states for a whole temperature regime. This integral equation can be numerically solved by iteration method. In order to make a discretization in the variable spacek, we need to find a proper cutoffk c . The cutoffk c is determined by choosing ln(1 + e −ε(kc) ) < 10 −9 because this term decreases quickly with increasingk, as it is shown in the Fig. 5. The number of discretization fromk = 0 tok =k c is N k = 1000. Once we getε(k) the dimensionless pressurep can be obtained byp Furthermore, we compare the pressure by taking different N k and find the difference can be negligible, as shown in the Fig.6.
Based on the above numerical method, the densityñ and entropys are obtained by the derivatives of pressure [45]ñ where the derivatives of ε µ ≡ ∂ε(k) ∂µ , ε T ≡ ∂ε(k) ∂T are given by These two equations may be numerically solved by iteration.
As an example, whenc approaches +∞, the last integral term in the dimensionless TBA, Eq. A6, can be ignored and the gas behaves like free fermions [27]. In particular, at the critical pointμ = 0, the entropy per particle can be computed analytically: where Li s is a polylogarithmic function of order s.

Appendix B: Fractional Exclusion Statistics
In 1991, Haldane [8] formulated a description of the fractional exclusion statistics (FES) based on a generalized Pauli exclusion principle. This FES was further formulated by Wu [10,15] and others [9,11]. It has been proved that the 1D δ-function interacting Bose gas can be mapped onto ideal particles with FES, see [10,15,29]. In this sense, the dynamical and statistical interactions are transmutable. This allows one to deal with the thermodynamic properties through Haldane's FES.
In general, the relation between the interaction strength and FES parameter is very complicated.
However, under a strong interaction strength, the system may be equivalent to an ideal gas with a non-mutual FES, i.e. the FES parameter does not depends on the momenta of particles. For such a non-mutual FES with a parameter g, the occupation number f in a state with energy = 2 k 2 2m is given by where w obeys The thermodynamic properties in D dimension, such as pressure p, energy density E, and particle density n are given by and the entropy density s can be obtained from thermodynamic relation E = −p + µn + sT .
The left hand side of Eq. B2 is monotonically increasing with w. For a given g,μ andk, Eq.(B2) can be numerically solved by bisection method. The relative error of w in our numerical calculation is less than 10 −6 . The cutoffk c and discrete number N k here are the same with the one for solving Yang-Yang equation.
where σ =ñ,p, S/N and the subscripts "Y" and "F" denote the property obtained by Yang-Yang equation and FES, respectively. Our numerical results are shown in the Fig. 7. The contour plots of the deviations of density and pressure are respectively shown near the critical point, where 15%, 8% and 3% derivations are marked.
The power of simple, non-mutual FES under strong interactions: For a strong interaction, i.e.c 1, these thermodynamic properties obtained from the interacting system and from the ideal particles with FES are in excellent agreement, even for largeμ.
The bottom part of each panel in Fig. 7 illustrates this point.

The need of mutual FES under weak interactions:
On the other hand, under weak interactions, even when S/N shows fairly good agreement  Fig. 8 here, respectively. Thus these discrepancies are not caused by inaccuracy of our mapping formulae (Eqs. 7 and 6 in the main text) and cannot be alleviated by choosing a different "effective FES parameter g". Rather, such discrepancies inñ andp are primarily due to the need of including more complex mutual FES effects [15].

Appendix D: Dimension analysis and dimensionless quantities
In this section and the following two sections, we will explicitly explain how we derive the dimensionless observables (particle density, pressure and entropy per particle, especially at the quantum critical point) in Bose gas in the continuous space based on the simulations for the Bose-Hubbard model on the lattice.
For D-dimensional ultracold Bose gas, the Hamiltonian could be written as [46,47] This Hamiltonian has an equivalent field theory form whereψ(r) is the wave function operator at the position r.  By doing dimension analysis to Eq. D2, we obtain the dimension of the wavefunction operator and the parameters, that is where E and L represent the dimension of energy and length, respectively.
If we choose E to be the unit of energy and λ to be the unit of length, we have the following relation equations between the quantities and their corresponding dimensionless quantities Here, we add a tilde to each of the symbols to denote their dimensionless quantities; for convenience, the dimensionless quantity forψ is denoted asψ with the hat dropped.
From Eq. D6, the Jacobian determinant is So, the integrals have the following relation equation Substituting these relation equations into Eq. D2, we can finally obtain the dimensionless form of the Hamiltoniañ with the dimensionless parametersc = cλ −D E −1 andμ = µE −1 .
Furthermore, if we take E = T and λ = λ dB /(2 √ π) = 1/ √ T , where λ dB = 2 π/T is the thermal de Broglie wavelength, we can get Eλ 2 = 1, and the dimensionless Hamiltonian becomes withc = cT D/2−1 andμ = µT −1 . The corresponding dimensionless forms of those observables we are interested in, the particle number density, pressure and entropy per particle, becomẽ As for the Bose-Hubbard model simulated in our numerical part of work, its Hamiltonian reads where x denotes the lattice site vector,b x (b † x ) is the annihilation (creation) operator for bosons on the site x,n x =b † xb x is the particle number operator, and < x, x > indicates the summation runs over all the nearest neighbor sites. The parameter t is the tunneling parameter, U is the onsite interaction strength, and µ BH is the chemical potential. Following the same approach, we can also obtain the dimensionless Hamiltonian for this model, whereb x is the dimensionless quantity forb x , which is exactlyb x itself sinceb x is already dimensionless, andt = t/T ,Ũ = U/T andμ BH = µ BH /T are the corresponding dimensionless quantities for each parameter, and the dimensionless forms of the observables are as follows Appendix E: Mapping between Bose gases and the discrete Bose-Hubbard model Since d Drψ † (r)∇ 2ψ (r) = − d Dr ∇ψ † (r) · ∇ψ(r), the dimensionless form of the Hamiltonian for Bose gas, Eq. D15, can also be written as In order to investigate the mapping relation between Bose gas and Bose-Hubbard model, we discretize the space into N D small cells with the side length∆ =L/N , whereL is the size of the system, making the integral in Eq. E1 into sums: whereψ x ≡ψ( x∆), and x = (x 1 , x 2 , ..., x D ) is the index of the cells with the integer components x i ranging from 1 to N .∇ψ x are numerical differences ofψ x defined as Thus, the discretized Hamiltonian becomes If we make the replacementψ we can getH By comparing Eqs. E6 and D20, we obtain the mapping relations between Bose gas in the continuous space and Bose-Hubbard model in the lattice as follows or reversely,∆ From the last equation, we can know that the critical point in Bose gas µ = 0 corresponds to µ BH = −2Dt, that is, the critical point for the phase transition from a vacuum phase to a quantum liquid phase in Bose-Hubbard model.
With some more analysis, we can also obtain the mapping relationships between the observables in both models as followsñ What we need to notice is that the corresponding dimensionless observables in Bose-Hubbard model forñ andp are not simplyñ BH andp BH .

The extrapolation protocol and its application to 1D Bose-Hubbard systems
According to the discretization approximation above, we would expect that when the dimensionless parametersc,μ in Bose gas model and ratios U/t, µ BH /t, T /t in Bose-Hubbard model are associated by Eq. E11 and Eq. E12, the corresponding dimensionless observables in two systems are equivalent to each other except a correction brought by the finite dimensionless spacing∆, which could be expressed by the following equatioñ whereÕ andÕ BH are two corresponding dimensionless observables in interacting Bose gas and Bose-Hubbard model respectively, andf is the dimensionless correction function decaying to zero when∆ is approaching to zero, that is From Eq. E10, we know that the correction function could be rewritten as a function of the temperature, that isf (g,μ, T /t), and∆ → 0 is equivalent to T /t → 0. Thus Eq. F1 and Eq. F2 becomesÕ (g,μ) =Õ BH (U/t, µ BH /t, T /t) +f (g,μ, T /t), Concluding from the analysis above, in order to obtain a dimensionless observable in Bose gas atc andμ, we can first compute the corresponding dimensionless observables in Bose-Hubbard model at different temperatures with the parameter ratios U/t and µ BH /t determined by Eq. E11 and Eq. E12, and then extrapolate the results towards zero temperature, which corresponds to the results in a Bose gas with no lattices.
As a typical example shown in Fig. 9, we measure and compute the dimensionless observables

A discussion on scale invariance and the extrapolation towards zero temperature
As shown by Eq. A6, the 1D interacting boson system studied in this work satisfies "scale invariance", namely, numerical or experimental data taken at different temperatures can "collapse" onto universal functions (S/N versusμ,ñ versusμ,p versusμ) when the parameters (quasimomentum k, interaction strength c) and thermodynamic quantities are scaled properly according to the thermal de Broglie wavelength. On the other hand, a lattice gas system does not strictly satisfy scale invariance. To approach these universal functions using 1D Bose-Hubbard model, we need to reach a parameter regime where the dimensionless lattice spacing∆ is much smaller than all other relevant dimensionless length scales -and thus decoupled from the physical properties of the system. This is equivalent to requiring that the thermal de Broglie wavelength be much larger than all other relevant length scales. We further convert this requirement into a final criterion that the temperature be much lower than all other relevant energy scales -and thus decoupled from the physical properties of the system. Mathematically this is satisfied when T /t becomes sufficiently small. This is the physical motivation of our protocol of "extrapolation towards zero temperature".
In a practical simulation, all our numerical data are obtained under finite temperatures. So our extrapolation results represent the physical properties when the system temperature is much lower than all other relevant energy scales. As long as the temperature does not play a role in influencing the physical properties of the system, we consider the purpose of the extrapolation protocol to be fulfilled.
In 1D, we know beforehand that the system satisfies scale invariance, such that the extrapolation towards zero temperature must have a well-defined limit. We indeed observe convergence in extrapolation (see Fig. 9), and observe that the QMC results (after extrapolation) agree excellently with the solutions to TBA equations (main text, Fig. 2). The statistical uncertainties of the extrapolated results reflect the accuracy of our QMC simulations and our confidence in using the QMC data to reveal the known physical and scaling properties of 1D Bose gases in continuous space. This computation serves as a calibration of our extrapolation protocol. When we apply this extrapolation protocol to 2D Bose-Hubbard lattice gases, where the corresponding continuousspace 2D Bose gases don't have an explicit equation like the TBA equation, we do not presume a prerequisite that the continuous-space system satisfies scale invariance. Instead, we rely on the statistical uncertainties of the extrapolated results to provide understanding on the scaling behaviors and physical properties of the continuous-space Bose gases under sufficiently low temperatures.
In an example shown in Fig. 10, In our work, we mainly focus on the following observables: particle density n, pressure p and entropy per particle S/N . In this section, we will show how we derive these observables in Bose-Hubbard model.
For both one-dimensional and two-dimensional systems, we apply worm algorithm in path- method [38,39]. During the simulations, we can directly measure the particle number density n BH and the grand energy density BH ≡ 1 V H BH , where V is the volume of the system, and in lattice model is just the total number of the sites. In order to ensure the data we obtained are in the thermodynamic limit, we keep the size of the system we simulated to be at least L = 40t/T , except when T /t = 0.1 we choose L = 20t/T for some of them. A typical example in two dimension is shown in Fig. 11, and it presents that the system sizes we choose are large enough to allow us to neglect the errors brought by finite system sizes.
Since the particle number density for Bose-Hubbard model n BH could be measured directly in our simulations, we will principally introduce how we derive the pressure p BH and the entropy S BH for Bose-Hubbard model.
According to the Gibbs-Duhem equation, dp = ndµ + sdT, (G1) if the temperature is fixed, we have Considering that when µ → −∞, the system will become vacuum with n(µ → −∞) = 0 and p(µ → −∞) = 0, we set µ 0 = −∞ and Eq. G2 becomes that is  11. The particle number density and the grand energy density in Bose-Hubbard model, n BH and BH , vary with LT /t. Both curves could be well fitted by the formula f (x) = a + be −cx , and according to the fitting results, even when LT /t = 20, the relative errors brought by the finite size of the system are of the order 10 −6 , which is pretty small compared with the statistical relative errors (typically 10 −4 ∼ 10 −3 ). the region below µ min BH , that is (−∞, µ min BH ). Therefore, we first measure n BH at different chemical potentials ranging from µ min BH to µ BH , with other parameters U , t and T keeping fixed, and then calculate the integral in Eq. G4 numerically by trapezoidal rule in the region [µ min BH , µ BH ]. And as for the region below µ min BH , we apply the distribution function for ideal Bosons to fit the tail of our data for n BH with the only fitting parameter ε, and based on the fitting result, we estimate the integral in the region (−∞, µ min BH ]. An example is shown in Fig. 12. To calculate the pressure at the critical point p c BH , we typically set µ min BH = µ c BH − 10T and the interval for the numerical integral ∆µ BH = 0.04T , where µ c BH = −2Dt is the location of the critical point for the phase transition from a vacuum to a quantum liquid in Bose-Hubbard model. This corresponds toμ ranging from −10 to 0 with the dimensionless interval ∆μ = 0.04. As the example shown in Fig. 13, the errors coming from the finite interval during the numerical integrating can be ignored compared with the statistical errors.
In order to calculate the entropy of the system, we apply the following quasi-static process: keep In order to calculate the pressure at, for example, µ BH /t = −4, we measure n BH in the region µ BH /t ∈ [−9, −4] with the interval ∆µ BH /t = 0.02, and calculate the integral in Eq. G4 numerically by trapezoidal rule and get 0.02686 (6). We also fit the tail (here, we select µ BH /t ∈ [−9, −6.5]) of the data by the formula Eq. G5 with the fitting result ε = −2.4280 (5), and according to this, we obtain the estimate for the integral in the region (−∞, −9], which is around 1 × 10 −6 . Thus, in fact, the relative deviation induced by the truncation is of the order 10 −5 which is pretty small compared with the relative statistical errors(∼ 10 −3 ). We have taken account of the integral from −∞ to µ min BH in our results, but since µ min BH we choose are small enough, we could safely ignore the effect brought by the truncation in principle.
all the parameters (including t, U and T ) fixed except the chemical potential µ varying from µ 0 to µ extremely slowly, so that we could regard the system as always in equilibrium states. Following this process, we have On the other hand, the entropy of a system is defined by the following formula where ρ is the density matrix with the explicit expression Here, Z = T r e −βH is the partition function, H is the Hamiltonian in grand ensemble and typically could be written as H = H 0 − µN . According to Eq. G8, we could rewrite Eq. G7 in a more explicit form, Therefore, Substitute this into Eq. G6 and we can get For Bose-Hubbard model, it becomes where BH (µ BH ) could be directly measured in our QMC simulations and p BH (µ BH ) could be obtained via the formula Eq. G4 by numerical integral techniques. For the 1D gas experiment by the Kaiserslautern group, we directly obtain the S c /N andñ c data from Ref. [22]. For the 1D gas experiment by the USTC group [20], we obtain the original data of equation of state (n(µ)) that are taken under multiple temperatures, and re-analyze them to obtain the critical pressure p c and the critical entropy density s c based on s = ∂P ∂T µ . For the 2D gas experiment by the ENS, we obtain the S c /N data directly from Ref. [23]. For the 2D gas experiments by the Chicago group [21,25] where the 2D lattice gases / continuous-space gases are experimentally observed to satisfy scale invariance, we either obtain data from these two references or process the original data to extract the needed thermodynamic observables.
For another experiment by the Chicago group [24], we obtain theñ c andp c data from Ref. [24], but do not have sufficient data to extract S c /N because the authors neither test scale invariance in the strong coupling regime nor perform groups of measurements for the samec 2D under multiple temperatures.
Finally, while our theoretical and numerical data agree quite well with existing experiments [20][21][22][23][24][25] in Figs. 2, 3, and 4 of the main text, we note that there are experimental factors that in principle can lead to the break-down of scale invariance to some extent in experimental systemsin particular, in the strongly interacting experimental systems. For example, three-body loss effects scale as the fourth power of the atomic scattering length. In addition, when one prepare strongly interacting experimental systems by approaching the unitarity limit (Feshbach resonance) or using deep optical lattices, the scattering amplitude at finite temperatures can become momentumdependent. These and other experimental factors could in principle cause the break-down of scale invariance to some extent in experimental systems, but the quantitative characterization of such break-down still remains sparse. The comparison of our theoretical and numerical data with existing experiments provide new insights regarding this issue. The overall good agreements suggest that in existing experiments, a description based on scale invariance is fairly consistent with the underlying physics within the current level of experimental uncertainties. The residual small discrepancies between our results and existing experiments may come from practical factors including inelastic losses and finite-temperature effects in experiments.