Exact Thermodynamics For Weakly Interacting Normal-Phase Quantum Gases: Equations of State For All Partial Waves

While the thermodynamics for bosonic systems with weak $s$-wave interactions has been known for decades, a general and systematic extension to higher partial waves has not yet been reported. We provide closed-form expressions for the equations of state for weakly interacting systems with arbitrary partial waves in the normal phase. Thermodynamics, including contact, loss rate, and compressibility, are derived over the entire temperature regime. Our results offer an improved thermometer for ultracold atoms and molecules with weak high-partial wave interactions.

This article is devoted to the EOS of single-component Bose and Fermi gases with weak interactions in the normal phase.The EOS is well understood and available in an analytical form for single-species bosons with weak s-wave interactions [34].In contrast, for singlecomponent Fermi gases with weak p-wave interactions, the contacts and EOS have only been studied in the lowand high-temperature regimes [32,33,35], even though higher partial-wave physics has attracted increased attention recently [36][37][38][39].Analytical expressions for the EOS of single-component gases beyond the s-wave case (i.e., for p-wave Fermi gases, d-wave Bose gases, f -wave Fermi gases, etc.)-applicable over the entire temperature regime-do not exist.
Within the Hartree-Fock framework, we derive analytical closed-form expressions for the EOS, applicable to all temperatures, of single-component atomic or molecular quantum gases with weak l-wave interactions in the * yqyan@cuhk.edu.hknormal phase.We calculate the contact, which determines the chemical reaction rate of ultracold gases.Using the virial expansion, we find that, while the contact of weakly-interacting s-wave Bose gases in the normal phase is a pure two-body quantity, that of weakly-interacting p-wave Fermi gases displays pronounced three-body effects even at temperatures as high as the degeneracy temperature.This effect is shown to arise from many-body dressing, i.e., the emergence of quasi-particles at leading order in the interaction strength.We also discuss the relation between the resulting reaction rate and that obtained through a simple thermal average over the inelastic cross-section.Applying the local-density approximation (LDA), we calculate the contacts of the harmonically trapped systems.Our results show that the trapped system needs to be cooled to rather low temperatures to probe the "low-temperature" portion of the EOS of the homogeneous system.
The article is arranged as follows: Section II introduces the l-wave low-energy two-body interaction potential employed in Sec.III to derive the normal-phase l-wave EOS in the weak-interaction limit.Section IV applies the EOS of the homogeneous system to deduce explicit, yet general, expressions for the two-body contact and two-body loss rate coefficient, which are interpreted using the virial expansion (see Sec. V).Section VI focuses on the homogeneous l = 0 and l = 1 systems.The loss rate coefficients derived in this work are compared with heuristic thermal averages in Sec.VII.Section VIII applies the homogeneous EOS to harmonically trapped systems using the LDA.Finally, Sec.IX discusses the applicability regime of the theory results derived in this work, while Sec.X concludes.Technical details are relegated to several appendices.

II. INTERACTION MODEL
The low-energy two-body potential for arbitrary partial-wave channel l reads U l (q, q ′ ) = 4πg l q l (q ′ ) l l m=−l Y lm (q)Y * lm (q ′ ), (1) where Y lm (q) is the spherical harmonic, and q and q ′ are the incoming and outgoing relative momenta [40].
The two-body phase shifts δ l are given by k 2l+1 cot(δ l ) = −1/a l + O(k 2 ), where a l is the scattering length in the lth partial-wave channel.To describe the binding energy of shallow two-body bound states, the leading-order effective-range correction needs to be included, and δ l needs to be expanded up to order k 2 [41].However, effective-range corrections can be excluded since we work in a weakly interacting regime where bound states do not contribute.A standard renormalization procedure gives (see Appendix A) Noticing that since we are using first-order perturbation theory where ultraviolet divergencies are absent, renormalization is not required, implying that the bare coupling g l and scattering length a l are related by g l = 4πℏ 2 a l /M , where M denotes the mass of the gas constituents (atoms or molecules).For l = 0, the interaction U 0 is, as expected, equal to g 0 [42].The next section uses the interaction potential U l (q, q ′ ) to derive perturbative results for the normal-phase EOS.

III. EQUATIONS OF STATE IN NORMAL PHASE
To include the two-body interactions in the EOS, we account for the mean-field corrections to the Bose-Einstein distribution function (l even) and Fermi-Dirac distribution function (l odd) in momentum space [42], where ϵ 2M denotes the single-particle kinetic energy, ℏΣ l (k) the self-energy, and ℏk the momentum.In Eq. ( 3) and in what follows, the upper sign is for even l (single-component bosons) and the lower sign for odd l (single-component fermions).The self-energy reads from which we can obtain the normal-phase grand potential Ω by the "generalized Hellman-Feynman theorem" [42] Ω − where n k (λ) and Σ(k, λ) are defined through Eqs. ( 3) and ( 4) with the two-body potential U l (q, q ′ ) scaled by λ 1 .Here, is the non-interacting grand potential, where is the thermal wavelength, and z is the fugacity.At the leading order in the scattering length, we find where the indices i, j, and n start from 0, and Here, Li s and Γ are the polylogarithm and gamma functions, respectively.To construct the full EOS, the mean particle density n needs to be expressed in terms of z.
We achieve this by treating the self-energy as a small parameter and integrating Eq. (3) in momentum space: One can check that Eqs. ( 7) and ( 9) fulfill the thermo- , where µ denotes the chemical potential.Equations ( 7) and ( 9) are the first main result of this article.From Eqs. ( 7) and ( 9), one can-at least formally-calculate all thermodynamic quantities.Fully analytical expressions for the isothermal compressibility, entropy, and isochoric heat capacity are given in Appendix B.

IV. CONTACT AND TWO-BODY LOSS
In addition to the observables considered in Appendix B, we consider the contact C l , which is conjugate to the inverse scattering length.The contact has been discussed extensively for the two-component Fermi gas at unitary [15][16][17][18][19][20][21].Working in the grand canonical ensemble, where the fugacity z is a thermodynamic variable, C l is defined in terms of the grand potential Ω by the adiabatic relation This definition of C l generalizes the definition of the pwave contact C 1 [24].For s-wave interacting Bose gases, the most commonly employed definition of the contact C 0 differs from Eq. ( 10) by a factor of 2π [12][13][14].The description of higher-partial wave systems typically requires a second contact, namely the conjugate to the effective range [23][24][25][26][27]. Since we find that it affects the thermodynamics of weakly interacting systems at sub-leading order, we exclude it from our discussion.The contact of weakly-interacting systems is a fascinating thermodynamic quantity since it determines the loss rate due to chemical reactions between two particles.Examples of reactions in molecular NaRb and KRb gases are: When the two incoming reactants are "scattered" into final products, the (typically large) binding energy is converted to the kinetic energy of the products.Consequently, the products have so much energy that they are not held in place by the comparatively shallow trapping potential.Since the reaction time is short compared to the typical time scale of experimental observations, a non-hermitian Hamiltonian with a complex interaction potential can effectively describe the process.For the single-component p-wave gas, it was shown that the change of the number N of constituents is related to the imaginary part of the scattering length [33], where H is the effective Hamiltonian with complex interaction and ⟨•⟩ denotes the thermal average.Since the derivation in Ref. [33] was done in real space, without making any assumptions about the form of the interaction, the result can be straightforwardly generalized to arbitrary partial-wave channels: From Eq. (11) and the definition of the contact, Eq. ( 10), one obtains where n denotes the particle density; the loss-rate coefficient β l can be measured experimentally [43][44][45].The loss-rate coefficient characterizes-due to the n 2 termlosses that arise from two-body processes.In general, though, the loss-rate coefficient may be n-dependent, implying that dn/dt may effectively scale with n 3 or n to some other power.
Combining the EOS and the definition of the contact, we find the contact in the canonical ensemble: (±z (0) ), (14) where the fugacity z (0) of the non-interacting system is implicitly determined by Li 3/2 (±z (0) ) = ±nλ 3  T .Equation (14) and its interpretation and implications (see below) are the second main result of this paper.

V. VIRIAL EXPANSION ANALYSIS
To unravel how the many-body thermodynamics emerges from the two-body scattering length and fewbody correlations, we employ the virial expansion, which provides a systematic expansion in terms of one-, two-, three-, and higher-body clusters [1].Formally, we expand Ω in terms of the fugacity z, where Z 1 = V /λ 3 T is the canonical partition function for a single constituent in a box with volume V .The determination of the virial coefficient b j requires information up to the canonical partition function Z j for j constituents [46], i.e., b j contains one-, two-, • • • , j-body physics; Z j with j > 1 accounts for interactions as well as exchange statistics.Since we have an analytical expression for Ω, the virial coefficients b j can be calculated analytically up to arbitrarily large j by Taylor expanding Eq. ( 7) around z = 0. We provide expressions for ∆b j = b j − (±1) j−1 j 5/2 up to j = 4: The expressions for ∆b j will be interpreted below.The following section applies the l-wave result for C l to two commonly investigated systems, namely s-wave Bose (in this case, our virial coefficients agree with the literature [47]) and p-wave Fermi gases.

VI. HOMOGENEOUS SYSTEMS A. Single-component s-wave Bose gas
The contact C 0 for the weakly-interacting singlecomponent Bose gas, applicable to any temperature T above the transition temperature T C , is directly proportional to n 2 : Since the quantity n 2 can be interpreted as the semiclassical pair density, the thermodynamic variable C 0 is a two-body quantity in the weak-interaction limit; in other words, many-body dressing is absent.As a consequence, the corresponding loss-rate coefficient β 0 is independent of n, Even above degeneracy, two-body chemical reactions of the weakly-interacting single-component s-wave gas do not exhibit three-or higher-body correlations.This behavior can be traced back to how the self-energy modifies the momentum distribution Eq. ( 3).At the mean-field level, the s-wave interactions lead to a self-energy Σ 0 that is independent of k, i.e., Σ 0 (k) = Σ 0 (see Appendix B 1).According to Eq. ( 3), the interactions can thus be interpreted as modifying the chemical potential without modifying the character of the constituents, i.e., the constituents remain free particles, and each two-body collision involves exactly two "physical" constituents.The virial expansion formalism can further confirm the interpretation.By self-consistently calculating the contact with truncated virial expansion at j, we find that j = 2 is enough to produce the exact results and that choosing a higher j does not introduce new terms (see Appendix C).

B. Single-component p-wave Fermi gas
Setting l = 1 in Eq. ( 14), we find Since the polylogarithm on the right-hand side of Eq. ( 19) has the index 5/2 as opposed to 3/2, the polylogarithm cannot, contrary to the s-wave case, be directly converted to n.Consequently, C 1 features a non-trivial dependence on n and T .At high temperatures (z (0) → 0), Eq. ( 19) becomes In this regime, the contact C 1 has-similar to the contact C 0 -a two-body nature.However, unlike in the s-wave 1. Contact (or two-body loss-rate coefficient), both in scaled dimensionless form, for single-component p-wave gas as a function of scaled temperature.The solid line shows Eq. ( 19); dotted, dashed, and dash-dotted lines show the second-, third-, and forth-virial expansions.Inset (i): Extension to larger T /TF , illustrating that the second-order virial expansion converges to the exact result at relatively high temperatures.Inset (ii): Contact for-from bottom to top at T = T C/F -p-wave (black), d-wave (cyan), s-wave (yellow), and f -wave (magenta).
case, the high-temperature p-wave contact has an explicit temperature dependence.Since λ −2 T is directly proportional to T , C 1 increases linearly with temperature, i.e., reactions become slower as the gas is getting colder.The corresponding β 1 at high temperatures is independent of n and linearly dependent on T , In the zero temperature limit, Eq. ( 19) approaches Appendix D discusses how to evaluate the zerotemperature limits of some of the functions that enter into the l = 1 EOS.Since the n-dependence where denotes the Fermi temperature.The black solid line in Fig. 1 shows Eq. ( 19) as a function of T /T F .The T 1 -and T 0 -scalings in the highand low-temperature regimes fully agree with previous works [32,33].
To interpret the change of the dependence of C 1 from being proportional to n 2/3 at low temperatures to being proportional to n 2 at high temperatures, we first note that p-wave interacting gases may exist in the normal phase approximately all the way down to zero temperature since the superfluid transition temperature is exponentially small [40].It is then natural to assume that many-body effects will modify the reaction rate in the low-temperature limit as the incoming and outgoing momenta are expected to be constrained due to the fermionic exchange statistics, i.e., intuitively, one expects some dressing of the constituents due to many-body effects.A careful analysis of the self-energy confirms this picture.Substituting Σ 1 (k) = A 1 + B 1 k 2 (see Appendix B 1) into Eq.( 3), the constant A 1 can be shown to modify, just as in the s-wave case, the chemical potential.The B 1 k 2 term, in contrast, modifies the single-particle energies ϵ (0) k , thereby effectively renormalizing the mass of the physical constituents.When a chemical reaction happens at low temperatures, two quasi-particles with effective mass interact instead of two physical constituents.Since the renormalization of the mass is due to manybody dressing, the chemical reaction involves more than two physical constituents.
The above analysis is complemented by the virial expansion up to the fourth order in z. Figure 1 compares the contact C 1 , calculated up to second, third, and fourth order, with the exact result, Eq. ( 19). Figure 1 shows that the second-order expansion agrees with the exact expression at T ≫ T F [see inset (i)].Importantly, the second-order virial expansion deviates notably from the exact result for temperatures as high as T /T F = 2.The third-order virial expansion provides an excellent description down to T /T F ≈ 0.25.Interestingly, the fourthorder virial expansion does not yield much improvement over the third-order expansion, indicating that threebody processes are essential in chemical reactions of weakly interacting p-wave gases for T /T F ≈ 0.25 − 2. At higher temperatures, three-body processes contribute very little.At lower temperatures, the chemical reactions acquire many-body characteristics.
Extending the analysis to higher partial waves, we find that the C l for l > 1 also have non-negligible three-body contributions in the vicinity of the degeneracy temperature (T C for even l and T F for odd l).In the high-T limit, C l is directly proportional to T l ; this scaling is consistent with the two-particle Bethe-Wigner threshold law [50][51][52].The inset (ii) of Fig. 1 plots our analytical expressions for C l for l = 0 − 3 as a function of temperature.

VII. STATISTICAL AVERAGE OF INELASTIC CROSS SECTION
In the literature, the loss-rate coefficients β l have been calculated by thermally averaging the two-body inelastic cross sections σ in,l (E).In what follows, we review the steps taken within this approach to derive β l [53][54][55].According to the definition of the scattering length a l for the l-th partial-wave channel, namely k 2l+1 cot(δ l ) = −1/a l , the scattering matrix element S l in the low-energy threshold limit reads The inelastic partial-wave cross-section σ in,l is related to the scattering matrix element through [66][67][68][69] where the factor 2l + 1 originates from the fact that the lth partial-wave channel has a (2l + 1)-fold degeneracy.
Assuming |Im(a l )|k 2l+1 ≪ 1 and |Re(a l )|k 2l+1 ≪ 1, we find Utilizing the definition of the scattering energy E = ℏ 2 k 2 /2µ, where µ = M/2 is the two-body reduced mass, one obtains The loss-rate coefficient β l is then found by thermally averaging the two-body inelastic cross section σ in,l (E) over the Boltzmann distribution function [53][54][55]: where the factor 2 reflects that one inelastic collision process eliminates two particles.Since Eq. ( 28) employs the Boltzmann distribution function, it is instructive to compare it with the high-temperature limit of the expression for β l derived in this work within the thermodynamic formalism.Our exact result and its high-temperature limit read It can be checked that Eq. ( 29) agrees with Eq. ( 28) for each partial wave channel.This can be understood because the two particles' center-of-mass and relative momenta obey the Boltzmann distribution separately.
At lower temperatures, however, the thermal average needs to be generally performed over the product of two Bose-Einstein or two Fermi-Dirac distribution functions (the three-body analog is discussed in Ref. [56]).Since the product of two such distribution functions does not, heuristic thermal average [Eq.(30)] thermodynamics [Eq.( 29)] Comparison between results from Eq. ( 30) and Eq. ( 29) up to l = 3.
unlike in the case of the Boltzmann distribution function, separate in relative and center-of-mass coordinates, the thermal-average approach does not straightforwardly extend to the low-temperature regime.By naively replacing the classical Boltzmann distribution with the quantum version (Bose-Einstein or Fermi-Dirac distribution), Eq. ( 28) becomes (30) The question mark over the equal sign indicates that the expression is not rigorous but instead, deduced heuristically.The integral in Eq. ( 30) can be evaluated analytically, and the results for l = 0 − 3 are reported in the second column of Table I.Curiously, a comparison of the thermal-average approach and our exact results (third column of Table I) shows that the heuristic thermal-average approach does yield the same expressions for l = 0 and l = 1 as the rigorous thermodynamic framework developed in this work.For the higher partial wave channels (l = 2 and l = 3), however, the heuristic approach yields a different temperature dependence.The correction of the heuristic expressions of the loss rate coefficient constitutes the third main result of this paper.

VIII. HARMONICALLY TRAPPED SYSTEMS
We now apply our results to harmonically trapped Nparticle systems, which are being studied extensively experimentally.To account for the trap-induced inhomogeneity of the density, we convert our homogeneous EOS, namely Eqs. ( 7) and ( 9), to those for the trapped system via the LDA [57].Within this framework, the local density at position r determines the EOS of the homogeneous system to be used at that point: Ω trap = d 3 rΩ[n(r)]/V.To obtain the contact C trap l of the trapped system, Eq. ( 10) is evaluated numerically.The black solid lines in Fig. 2 show the result for a spherically symmetric harmonic trap with angular frequency ω.To gain physical insight, the EOS of the inhomogeneous system can be described through the virial expansion.In an isotropic harmonically trapped system, the relation between the virial coefficients b trap j of the trapped system and those of the homogeneous system is b trap j = b j /j 3/2 [46].The explicit description of the thermodynamics of harmonically trapped single-component gases with weak interactions and the interpretation thereof (see below) constitute the fourth main result of this paper.
Figures 2(a) and 2(b) compare the contacts of the harmonically trapped s-wave Bose and p-wave Fermi gas, obtained from the virial expansion up to fourth order, with the full numerical results.We make two observations: (i) For the s-wave Bose gas, the contact of the trapped system does not coincide with the second-order virial expansion, indicating that the contact of the trapped system is not, unlike that of the homogeneous system, a twobody quantity.This is because each position in the trap has a distinct local self-energy.Correspondingly, the fugacity of the trapped system cannot be interpreted as a globally shifted fugacity of the non-interacting gas.(ii) The second-order viral expansion for the p-wave system works quite well at T ∼ T trap F ; apparently, three-body corrections play a rather small role near the degeneracy temperature.The two-body (high-temperature) approximation works better for the inhomogeneous system than the homogeneous system since the former is much hotter than the latter for the same scaled temperature (e.g., T /T trap F = 1 and T /T F = 1 correspond to z ≈ 0.17 and z = 0.98, respectively).This validates previous works on loss processes of harmonically trapped p-wave systems [32,33].Finally, we note that while the adiabatic relation, Eq. ( 10), holds for the trapped system, the relation between the loss rate and the contact differs from that for the homogeneous system since the loss-rate coefficient of the trapped system is not only governed by particles being lost from the trap but also by a so-called deformation effect [33].

IX. VALIDITY REGIME OF THE THEORY BASED ON b2
The results presented in this paper employ the meanfield framework and expansions applicable to the weakly interacting regime.It is thus natural to ask what the validity regime of the theory is and whether the theory covers the operating regime of typical state-of-the-art experiments.
The applicability regime can be estimated using the virial EOS of the homogeneous system that accounts for the second-order virial coefficient.At this level, the virial EOS can be analytically tackled for an arbitrary interaction strength.The virial EOS up to b 2 allows us to obtain exact reference results that can be used to assess the accuracy of Eq. ( 7).To facilitate the comparison, we compare the second-order virial coefficient b 2 , obtained within the mean-field framework, directly with its exact counterpart.
We consider the s-wave Bose and p-wave Fermi gases as examples.Their exact b 2 are where ãs is equal to a s k C and T is equal to T /T C , and where ṽp is equal to v p k 3 F and T is equal to T /T F .erf, erfc, and p F q denote the error, complementary error, and generalized hypergeometric functions, respectively.Equations ( 31) and (32) are obtained from the Beth-Uhlenbeck formula [46].The expressions are expected to apply to relatively large |ã s | or |ṽ p |.In writing Eq. ( 32), we ignored the contribution from the shallow bound state, i.e., we assumed that the effective range is infinitely large.Using Eq. ( 16), the second-order meanfield level virial coefficients, extracted from the EOS for l = 0 and l = 1 derived in this work, read Figure 3  2 |/|b 2 |between the mean-field coefficient and the exact virial coefficient is small, the mean-field treatment is expected to provide a faithful description.The gray-shaded area in Fig. 3 demarcates the parameter combinations for which the normalized difference |b 2 − b mf 2 |/|b 2 | between the mean-field virial coefficient and the exact virial coefficient is smaller than 0.05.A difference below 5 % is interpreted as indicating that the mean-field description provides an accurate description of the system.At T ≃ 2, which is typical for experiments, the normalized difference is smaller than 5 % for |ã s | ≲ 0.15 and |ṽ p | ≲ 0.1 for the s-wave Bose gas and p-wave Fermi gas, respectively.For lower values of T , the validity regime of the mean-field treatment extends, according to our criteria, to larger reduced interaction strengths.We caution that-even though the validity regime, as determined by |b 2 −b mf 2 |/|b 2 |-is quite large for temperatures that are much lower than the degeneracy temperature, Fig. 3 may not accurately reflect the validity regime of the full EOS.The reason is that the virial expansion fails at these low temperatures for which the fugacity is large (in this regime, virial coefficients other than b 2 come into play).We expect the "true" validity range for the mean-field EOS to be comparable to that at temperatures notably above the transition temperature.
Ultracold molecular gas experiments typically work with samples that are characterized by extremely weak interactions, which are thus expected to be well described by the theory developed in this work.Table II shows three examples: one for a s-wave Bose gas and two for p-wave Fermi gases.The molecular gas is loaded into an approximately harmonic trap in the experiments.We use the typical densities reported in the experimental works [11,45,58] to calculate the reduced interaction strength and subsequently calculate the transition temperature T C and Fermi temperature T F assuming that the systems are homogeneous.Table II shows that the reduced interaction strength is well within the parameter regime where the mean-field-based theory developed in this work is applicable.II.Species used in molecular gas experiments, their statistics, their interaction strength, and their reduced interaction strength.

X. CONCLUSION
In summary, we theoretically derived the EOS for single-component normal-phase quantum gases, which can be used to obtain all thermodynamic quantities.We focused on the behavior of the contact of two commonly produced systems-s-wave Bose and p-wave Fermi gases.We showed that the contact is purely a two-body quantity in the former system and exhibits many-body characteristics in the latter.We analyzed the behavior of the p-wave contact in the near-degenerate regime and found that the three-body contribution plays a vital role.The discussion was extended to harmonically trapped systems, where we analyzed the contacts under the LDA.
Our study provides critically needed guidance for recent ultracold molecular gas experiments, i.e., for weakly interacting molecules in the deeply degenerate regime where the virial expansion fails.Specifically, our results can be used to calibrate loss rate and temperature measurements.Moreover, our results also apply to single-component Fermi gases such as 6 Li and 40 K [60,61], and provide a reference for studying crossover from weakly to strongly interacting systems.To obtain the relation between the bare couplings g l and the scattering lengths a l for each partial wave channel l, we follow the standard renormalization procedure, i.e., we compare the T -matrix element T l (k, k ′ ) and the partial wave scattering amplitude f l [57]: where In Eq. (A1), M , V , and k denote the mass of the constituent, volume, and relative wave vector, respectively.P l is the Legendre polynomial of degree l.The T -matrix elements are defined by the Schwinger-Dyson equation [62] T where ϵ has an infinitesimally small positive real value that ensures retarded propagation.Making use of the form of interaction in the main text, Eq. (A3) can be worked out explicitly: We find that the infinite sum inside the square brackets forms a geometric sequence of the form 1 where c is equal to the second term in square brackets in the last line of Eq. (A4).Hence, we can write Because of energy conservation, the magnitudes of k 1 and k 2 should be the same.Setting |k 1 | = |k 2 | = k, denoting the angle between k 1 and k 2 by θ, and using the addition theorem of spherical harmonics, we obtain Combining Eq. (A1) and the definition of the partialwave phase shifts δ l (k), we get (A8) The integral on the left-hand side of Eq. (A8) diverges.It can be reexpressed using the well-known low-energy relation [62]: where P denotes the Cauchy principal value.Using Eq. (A9) in Eq. (A8), we finally obtain the renormalization condition Eq. ( 2).It is known that the derived renormalization condition cannot eliminate all diverging terms that may arise in the many-body treatment of the singlespecies bosonic system, especially when the interaction is strong [63].For example, in single-component bosonic systems, Efimov physics requires one to introduce an additional parameter for renormalization [64].However, since the lowest-order mean-field interactions dominate the many-body physics for the weakly-interacting systems of interest in this work, the renormalization condition is not needed to eliminate divergencies, i.e., we can use the non-integral part of the bare coupling constant.
Appendix B: Details of Calculation on Thermodynamics 1. Self-Energy and Grand Potential Substituting Eq. (1) into Eq.( 4), the self-energy is found explicitly to be where θ denotes the angle between k and k ′ .To proceed, we work in the weakly-interacting limit and assume that n k is equal to the non-interacting distribution functions n Inserting Eq. (B2) into Eq.(B1), we note that the θ dependence only appears in the term |k−k ′ | 2l .To evaluate the integral over θ, the usual binomial theorem for scalars cannot be applied to |k − k ′ | 2l .Instead, we first take the square and then construct a series expansion using the trinomial theorem [65]: where the indices i, j, and n each go from 0 to l. Letting, as before, the angle between k and k ′ be θ and using Using the integral expression of the polylogarithm function, the self-energy becomes Li 2j+n+3 2 (±z).
(B7) Equation (B7) reveals the structure of the self-energy for the lth partial-wave channel clearly: Σ l (k) = Σ l (k) is a polynomial of even degree in k (see the factor of k 2i+n in the summand) since 2i is always even and the summand is zero when n is odd.For example, Σ 0 = A 0 for the s-wave channel, Σ for the d-wave channel, etc., where A l , B l , and C l are constants that depend on the temperature T .Since Σ l (k, λ) is directly proportional to λ and a l [see Eq. (B7)] and since n k (λ) has no dependence on λ at leading order [n k (λ) ≈ n (0) k ], the leading-order modification of the grand potential based on Eq. ( 5) is given by Using Eqs.(B5) and (B6), we arrive at Eq. ( 7) in the main text.

Isothermal Compressibility
This section considers the isothermal compressibility, which is defined through In terms of the grand potential Ω, we find Using Eq. ( 7), we find the isothermal compressibility in terms of the fugacity: This expression is not yet in the "standard form" of the compressibility, as it depends on the fugacity (a thermodynamic variable in the grand canonical potential) rather than the density (a thermodynamic variable in the canonical ensemble).To convert z to n, we use the Gibbs-Duhem relation We know that n is equal to ±Li 3/2 (±z (0) )/λ 3 T at leading order.If we write z in terms of z (0) , z = z (0) + δz, then we find to first order, Substituting z = z (0) + δz into Eq.(B15), we have The explicit expressions for l = 0 and l = 1 read The divergence of the compressibility of the bosonic s-wave system indicates the Bose-Einstein condensate (BEC) transition.The first equation of Eqs.(B19) diverges at z (0) = 1 or T = T C , i.e., the BEC transition temperature of the non-interacting system.The argument extends to Bose gases with higher partial-wave interactions because Li s (x) diverges at x = 1 for all s, s ≤ 1.It is worthwhile pointing out that T C does, in fact, have a correction of order n 1/3 a 0 , which arises from higher-order fluctuations that are not captured by the mean-field approach considered here [34].

Entropy and Isochoric Heat Capacity
This section determines the isochoric heat capacity.We start with the entropy, which is defined through To obtain the entropy in the canonical ensemble in terms of n and T , we change z to z (0) and keep terms up to first order in a l /λ 2l+1 T : It can be noted that the entropy for l = 0 is quite special since it is independent of a 0 .After simplification, we find that the entropy of the weakly-interacting s-wave Bose gas is identical to that of the non-interacting Bose gas: ) ) .
(B22) The independence of a 0 is a consequence of the fact that the s-wave mean-field interaction simply shifts the chemical potential, making the interacting system resemble a non-interacting gas.The p-wave interactions, in contrast, modify the entropy of the non-interacting system: We can now calculate the isochoric heat capacity directly from the entropy, Explicitly, the expression is Consistent with the discussion above, one can check that the weakly interacting s-wave system has the same C V as the corresponding non-interacting system.In contrast, the isochoric heat capacity of the weakly-interacting pwave gas contains a correction due to the interactions.

Appendix C: Calculation of Contacts via Truncated Virial Series
This work uses the virial expansion, truncated at order j, to analyze the j-body contribution to the contact and other thermodynamic quantities.A crucial point of the derivations is to consistently account for terms that contribute to the j-body physics.In particular, this implies that higher-order terms that go beyond j-body physics need to be excluded.
We start from the general virial expansion in the grand canonical ensemble, truncated at order j, Ω = −k B T Z 1 j i=1 b i z i .Here, Z 1 denotes the one-body partition function for an as-of-yet unspecified system (could be the homogeneous or inhomogeneous system).Using the standard thermodynamics relation From the expression for N and the definition of the contact, it is straightforward to calculate the contact in the grand canonical ensemble at order j: In what follows, we are interested in the contact in the canonical ensemble.To convert Eq. (C2) from the grand canonical to the canonical ensemble, we need to express z in terms of N .This is achieved by applying the Lagrange inversion theorem (see Sec. 3.6.6 of Ref. [70]) to Eq. (C1).
To obtain the forth-order results presented in the main text, we need to include terms up to the fourth power of N , When substituting Eq. (C3) into Eq.(C2) and working at order j, one needs to truncate all terms at order N j .Even though the z 2 term, e.g., generates terms that scale as N 2 , • • • , N 8 , only the terms proportional to N 2 , N 3 , and N 4 are kept to obtain consistent forth-order results.
The results reported below are obtained by additionally taking, consistent with the weak interaction regime assumption, the |a l |/λ 2l+1 T → 0 limit.The abovementioned strategy is critical for determining that the homogeneous system's s-wave contact has a pure two-body character.If we-incorrectly so-kept all terms in N when going from the grand canonical to the canonical ensemble, we would obtain + O(N 5 ) (truncated at j = 2) On the other hand, when following the correct approach that consistently keeps terms up to order N j , we find: (C7) The analysis above leads to the following important conclusion: When one includes terms up to j = 2, the s-wave result is exact; contributions from larger clusters do not lead to any corrections.The main text characterizes weakly interacting singlecomponent systems in the normal phase.This appendix provides limiting expressions for several quantities of fermionic systems in the T /T F → 0 limit.In the T /T F → 0 limit, it is computationally inefficient to evaluate expressions that explicitly or implicitly contain polylogarithm functions numerically.The reason is that the direct use of T = 0 yields expressions like Li s (−∞) cannot be evaluated numerically.Instead, it is generally more convenient to use asymptotic analytic expressions.The large-argument asymptote of fermionic-type polylogarithm functions is [71] −Li s (−x)

(D3)
One also needs to change the integration limits accordingly for harmonically trapped systems.To illustrate this, we consider the p-wave system.The exact integral for the contact is , which shows that the integrand vanishes at r = 1.This is consistent with the LDA, where the density goes to zero at r = R F .Since the density of the cloud is zero for r > R F , the upper limit of the integration can be changed from r = ∞ to r = 1, (D5) We note that our work does not rule out the existence of other phases in the extremely low-temperature regime.If, e.g., a BCS phase existed, the BCS transition temperature would be exponentially small in the weakly-interacting regime (|a l |k 2l+1 F ≪ 1) considered in this work.Evaluating the low-temperature normal-phase behavior would still be justified by considering the limiting T /T F → 0 expressions discussed above.

FIG. 2 .
FIG. 2.Normal-phase contacts, in scaled dimensionless units, for harmonically trapped (a) s-wave Bose and (b) pwave Fermi gas.The gray-shaded region in (a) denotes the BEC phase where our calculation does not apply.The exact results (solid lines) are compared with the virial expansions up to fourth order (see legend).Here, T trap C = (N/ζ(3)) 1/3 ℏω/kB [48] [T trap F = (6N ) 1/3 ℏω/kB [49]] is the transition [Fermi] temperature of the non-interacting trapped Bose [Fermi] gas; ζ(s) denotes the Riemann Zeta function, and k trap C and k trap F denote the momentum scales of the corresponding energy scales.

shows the virial coefficients b bose 2 (left panel) and b fermi 2 (
Figure 3 shows the virial coefficients b bose 2 (left panel) and b fermi 2 (right panel) as functions of the reduced interaction strength and reduced temperature.If the normalized difference |b 2 − b mf2 |/|b 2 |between the mean-field coefficient and the exact virial coefficient is small, the mean-field treatment is expected to provide a faithful description.The gray-shaded area in Fig.3demarcates the parameter combinations for which the normalized difference |b 2 − b mf 2 |/|b 2 | between the mean-field virial coefficient and the exact virial coefficient is smaller than 0.05.A difference below 5 % is interpreted as indicating that the mean-field description provides an accurate description of the system.At T ≃ 2, which is typical for experiments, the normalized difference is smaller than 5 % for |ã s | ≲ 0.15 and |ṽ p | ≲ 0.1 for the s-wave Bose gas and p-wave Fermi gas, respectively.For lower values of T , the validity regime of the mean-field treatment extends, according to our criteria, to larger reduced interaction strengths.We caution that-even though the validity regime, as determined by |b 2 −b mf 2 |/|b 2 |-is quite large for temperatures that are much lower than the degeneracy temperature, Fig.3may not accurately reflect the validity regime of the full EOS.The reason is that the virial expansion fails at these low temperatures for which the fugacity is large (in this regime, virial coefficients other than b 2 come into play).We expect the "true" validity range for the mean-field EOS to be comparable to that at temperatures notably above the transition temperature.

6 FIG. 3 .
FIG. 3. Validity regime determined by comparing the mean-field b mf 2 and the exact b2.Left panel: s-wave Bose gas.Right panel: p-wave Fermi gas.The heatmaps show the value of the exact b2 [Eqs.(31) and (32)]; the color scheme is defined in the bars that are shown to the right of the main panels.The gray-shaded areas denote the parameter regime for which |b2 − b mf 2 |/|b2| ≤ 0.05.A normalized difference below 5 % is interpreted as an indicator that the theory framework developed in this work is applicable.

truncated at j = 2 )
units of k C and T C , we find:C 0 N [Re(a 0 )] 2 k 3 Appendix D: Zero-Temperature Limit in Fermionic Systems

1 N
[Re(a 1 )] 2 (k trap F ) .05.A normalized difference below 5 % is interpreted as an indicator that the theory framework developed in this work is applicable.