The contact in the unitary Fermi gas across the superfluid phase transition

TL;DR

This paper uses advanced quantum Monte Carlo simulations to accurately determine how the contact parameter in a unitary Fermi gas varies with temperature, especially across the superfluid transition, aligning well with experimental data.

Contribution

It provides the first ab initio calculation of the temperature-dependent contact in the homogeneous unitary Fermi gas across the superfluid transition.

Findings

Contact decreases sharply near the critical temperature from below.

Gradual decrease of contact in the normal phase with increasing temperature.

Results agree with recent ultracold atomic gas experiments.

Abstract

A quantity known as the contact plays a fundamental role in quantum many-body systems with short-range interactions. The determination of the temperature dependence of the contact for the unitary Fermi gas of infinite scattering length has been a major challenge, with different calculations yielding qualitatively different results. Here we use finite-temperature auxiliary-field quantum Monte Carlo (AFMC) methods on the lattice within the canonical ensemble to calculate the temperature dependence of the contact for the homogeneous spin-balanced unitary Fermi gas. We extrapolate to the continuum limit for 40, 66, and 114 particles. We observe a dramatic decrease in the contact as the superfluid critical temperature is approached from below, followed by a gradual weak decrease as the temperature increases in the normal phase. Our results are in excellent agreement with the most recent precision ultracold atomic gas experiments. We also present results for the energy of the unitary gas as a function of temperature in the continuum limit.

Tables

Table 1. Table 1: Various estimates of the Bertsch parameter ξ 𝜉 \xi

Method	$\xi$	error
Fixed-node diffusion Monte Carlo Astrakharchik et al. (2004)	0.42	0.01
Duke experiment Luo and Thomas (2009)	0.39	0.02
ENS experiment Nascimbéne et al. (2010); Navon et al. (2010)	0.41	0.01
Ground-state fixed-node Monte Carlo Forbes et al. (2011)	$\le 0.383$	0.001
Ground-state AFMC Carlson et al. (2011)	0.372	0.005
MIT experiment Ku et al. (2012)	0.376	0.005
Lattice quantum Monte Carlo Endres et al. (2013)	0.366	${}_{-0.011}{}^{}{}_{}{}^{+0.016}$
AFMC (this work)	0.367	0.007

Equations

$$\int d^{3}R{g}_{\uparrow ,\downarrow }^{(2)}(\mathbb{R}+r/2,\mathbb{R}-r/2)\underset{r\to 0}{\sim }\frac{C}{(4\pi r{)}^2},$$

$$C=\frac{4\pi m}{{\hslash }^2}{\frac{\partial E}{\partial (-1/a)}|}_S.$$

$$\hat{H}=\sum _{\mathbf{k},\sigma }{ϵ}_k{\hat{a}}_{\mathbf{k},\sigma }^{†}{\hat{a}}_{\mathbf{k},\sigma }+g\sum _{\mathbf{x}}{\hat{n}}_{\mathbf{x},\uparrow }{\hat{n}}_{\mathbf{x},\downarrow },$$

$$\frac{1}{V_0}=\frac{m}{4\pi {\hslash }^{2}a}-{\int }_B\frac{d^{3}k}{(2\pi {)}^{3}2{ϵ}_k}$$

$$C=\frac{4\pi m}{{\hslash }^2}{\frac{\partial F}{\partial (-1/a)}|}_T,$$

$$C=\frac{m^2{V}_0⟨\hat{V}⟩}{{\hslash }^4},$$

$$e^{-\beta \hat{H}}=\int D[\sigma ]G_{\sigma }{\hat{U}}_{\sigma }+O((\Delta \beta {)}^2),$$

$$⟨\hat{O}⟩=\frac{\text{Tr}(\hat{O}{e}^{-\beta \hat{H}})}{\text{Tr}(e^{-\beta \hat{H}})}=\frac{\int D[\sigma ]⟨\hat{O}{⟩}_{\sigma }{W}_{\sigma }{\Phi }_{\sigma }}{\int D[\sigma ]W_{\sigma }{\Phi }_{\sigma }},$$

$${ϵ}_k^{(2)}=\frac{{\hslash }^2{k}^2}{2m},$$

$${ϵ}_{\mathbb{k}}^{(h)}=\frac{{\hslash }^2}{m\delta x^2}[3-\sum _i\text{cos}(k_i\delta x)],$$

$${ϵ}_k^{(3)}=\frac{{\hslash }^2{k}^2}{2m}\left[1-\alpha {\left(\frac{k\delta x}{\pi }\right)}^2\right],$$

Full text

Introduction.— The unitary Fermi gas (UFG) describes a system of spin-1/2 particles with a zero-range interaction and a diverging s-wave scattering length $a$ which saturates the upper bound on the modulus of the scattering amplitude imposed by the unitarity condition. This system is of interest for understanding the properties of other systems such as high-$T_{c}$ superconductors Randeria (2010); Mueller (2017) and nuclear matter Carlson et al. (2012); Gandolfi et al. (2015), and has been realized experimentally with 6Li and 40K ultracold atomic Fermi gases Regal et al. (2005); Ketterle and Zwierlein (2008); Mukherjee et al. (2017). Its quantitative understanding presents a challenge to theorists and experimentalists.

A quantity called the contact $C$ describes the short-range correlations of particles of opposite spin and is defined by

[TABLE]

where $g_{\uparrow,\downarrow}^{(2)}(\mathbb{r}_{\uparrow},\mathbb{r}_{\downarrow})=\langle\hat{n}_{\uparrow}({\bf r}_{\uparrow})\hat{n}_{\downarrow}({\bf r}_{\downarrow})\rangle$ is the density-density correlation function, with $\hat{n}_{\sigma}({\bf r})$ the density of particles at position ${\bf r}$ and spin $\sigma$. Several exact relations involving the contact, known as Tan’s relations, were derived in Refs. Tan (2008a, b, c). In particular, the contact characterizes the high-momentum tail of the normalized momentum distribution $n_{\sigma}(\mathbb{k})$ through the relation $n_{\sigma}(\mathbb{k})\underset{k\rightarrow\infty}{\sim}C/k^{4}$, where $k$ is the wavenumber and the distribution is normalized with $N_{\sigma}=\int\frac{d^{3}k}{(2\pi)^{3}}n_{\sigma}(\mathbb{k})$ ($N_{\sigma}$ being the total number of particles with spin $\sigma$) Tan (2008a). The contact also characterizes the high-frequency tail of the shear viscosity spectral function Taylor and Randeria (2010); Enss et al. (2011). It can be expressed in terms of the adiabatic derivative (at constant entropy $S$) of the thermal energy $E$ with respect to the inverse scattering length Tan (2008b)

[TABLE]

Other relations involving the contact were introduced in Refs. Combescot et al. (2006); Punk and Zwerger (2007); Baym et al. (2007); Braaten and Platter (2008); Braaten et al. (2008); Zhang and Leggett (2008); Werner et al. (2009); Zhang and Leggett (2009); Pieri et al. (2009); Combescot et al. (2009); Schneider et al. (2009); Hu et al. (2010); Braaten et al. (2010); Son and Thompson (2010); Nishida (2012); Hofmann and Zwerger (2017); see Ref. Braaten (2012) for a review.

Tan’s relations were verified experimentally in the ultracold atomic gas experiments of Refs. Stewart et al. (2010); Kuhnle et al. (2010). Soon after, the temperature dependence of the contact for the UFG was measured in a trap Kuhnle et al. (2011), followed by the measurement for the homogeneous system Sagi et al. (2012). Ref. Sagi et al. (2012) observed a sharp decrease in the contact as the temperature was lowered below the superfluid critical temperature. Recently, two independent precision experiments Carcy et al. (2019); Mukherjee et al. (2019) were able to address quantitatively the temperature dependence of the contact across the superfluid phase transition. Both experiments agree well with each other and show a dramatic increase in the contact as the temperature is lowered below the superfluidity transition temperature.

Calculating the temperature dependence of the contact for the UFG has proven challenging, and published results differ widely Palestini et al. (2010); Enss et al. (2011); Drut et al. (2011); Hu et al. (2011); Enss et al. (2011); Goulko and Wingate (2016); Rossi et al. (2018). This is not surprising given that many of the theoretical results were derived using uncontrolled approximations. However, two recent works are based on methods that have, in principle, controlled errors. Ref. Goulko and Wingate (2016) used a diagrammatic Monte Carlo approach on a lattice Burovski et al. (2006a) both in the superfluid and in the normal phases. Ref. Rossi et al. (2018) used the bold diagrammatic Monte Carlo method of Ref. Van Houcke et al. (2019), and was limited to the normal phase.

Here we use canonical-ensemble auxiliary-field quantum Monte Carlo (AFMC) methods Alhassid (2017); Jensen et al. (2018) on a spatial lattice to calculate the temperature dependence of the contact across the superfluid transition for $N=40,66$, and $114$ particles. For each of these particle numbers, we extrapolate to the continuum limit with no remaining systematic errors due to a finite filling factor (or equivalently finite effective range $r_{e}$ Werner and Castin (2012)).

Our continuum limit results differ substantially from the grand-canonical AFMC results of Ref. Drut et al. (2011), which were carried out at a finite filling factor. The temperature dependence we find is qualitatively similar to that found in the diagrammatic Monte Carlo approach Goulko and Wingate (2016) at temperatures below the critical temperature $T\lesssim T_{c}\simeq 0.15\,T_{F}$ (where $T_{F}$ is the Fermi temperature), but exhibits a different behavior above $T_{c}$. Our results for the contact show a similar qualitative behavior to the results of the bold diagrammatic Monte Carlo method Rossi et al. (2018) at temperatures $T>T_{c}$, but are systematically below them. Our calculations of the contact are in remarkable agreement with the recent precision experiments of Refs. Carcy et al. (2019); Mukherjee et al. (2019) both below and above $T_{c}$. Among available theoretical results for the contact, our calculations provide the best quantitative agreement with these experiments across the superfluid phase transition.

We also calculate the temperature dependence of the thermal energy in the continuum limit for $N=40$ and $66$ particles, and compare it with the experimental results of Ref. Ku et al. (2012). Taking the zero-temperature limit of the thermal energy, we estimate the Bertsch parameter to be $\xi=0.367(7)$, in agreement with the experimental value $\xi=0.376(5)$ of Ref. Ku et al. (2012).

Lattice formulation.— We discretize space with a cubic lattice of linear size $L=N_{L}\delta x$, where $\delta x$ is the lattice spacing. We use periodic boundary conditions and take a zero-range interaction of strength $V_{0}$, i.e., $V=V_{0}\delta(\mathbb{r}-\mathbb{r^{\prime}})$. The corresponding lattice Hamiltonian is given by

[TABLE]

where $g=V_{0}/(\delta x)^{3}$ is the coupling constant determined by the condition

[TABLE]

so as to produce the given scattering length $a$ on the lattice ($a\to\infty$ for the UFG). The integral over the wavevector $\bf{k}$ is restricted to the first Brillouin zone $B$ of the reciprocal lattice in momentum space of a spatial cubic lattice $\mathbb{x}=(n_{x},n_{y},n_{z})\delta x$, $n_{i}\in\{-M,-M+1,...,M\}$ where $M=(N_{L}-1)/2$ (we use odd $N_{L}$). The operators $\hat{a}^{\dagger}_{\bf{k},\sigma}$ and $\hat{a}_{\bf{k},\sigma}$ are, respectively, the creation and annihilation operators of a particle with wavevector $\bf{k}$ and spin $\sigma=\pm 1/2$ obeying fermionic anti-commutation relations $\{\hat{a}^{\dagger}_{\bf{k},\sigma},\hat{a}_{\bf{k}^{\prime},\sigma^{\prime}}\}=\delta_{\bf{k},\bf{k}^{\prime}}\delta_{\sigma,\sigma^{\prime}}$. The operator $\hat{n}_{\bf{x},\sigma}=\hat{\psi}^{\dagger}_{\bf{x},\sigma}\hat{\psi}_{\bf{x},\sigma}$ is the number operator of particles at lattice site $\bf{x}$ with spin $\sigma$, where $\hat{\psi}^{\dagger}_{\bf{x},\sigma}$ and $\hat{\psi}_{\bf{x},\sigma}$ are the creation and annihilation operators satisfying $\{\hat{\psi}^{\dagger}_{\bf{x},\sigma},\hat{\psi}_{\bf{x^{\prime}},\sigma^{\prime}}\}=\delta_{\bf{x},\bf{x}^{\prime}}\delta_{\sigma,\sigma^{\prime}}$. Here we use a quadratic single-particle dispersion relation $\epsilon_{k}={\hbar^{2}\bf{k}}^{2}/2m$. In the supplemental material we show that dispersion relations used in other works Burovski et al. (2006b, a); Goulko and Wingate (2010); Carlson et al. (2011); Goulko and Wingate (2016) lead to similar results after extrapolation to the continuum limit.

For a given lattice size $N_{L}^{3}$ and particle number $N$, there is a systematic error that arises from the finite lattice filling factor $\nu=N/N_{L}^{3}$, and an extrapolation $\nu\rightarrow 0$ is necessary to obtain the continuum limit for the given particle number. In the limit of low filling factor, the many-body energies scale as $\nu^{1/3}$ Burovski et al. (2006a); Pricoupenko and Castin (2007); Werner and Castin (2012). We therefore use a linear fit in $\nu^{1/3}$ for our low-filling-factor simulations to extract the continuum results.

Results.— We performed AFMC simulations in the canonical ensemble as described in Ref. Jensen et al. (2018). The simulations are carried out for $N=40,66$, and $114$ particles, on lattices of size $N_{L}^{3}=5^{3},7^{3},9^{3},11^{3},13^{3}$ and $15^{3}$. We divide the inverse temperature $\beta=1/T$ into discrete time slices of length $\Delta\beta$ (using the Trotter product for the propagator $e^{-\beta\hat{H}}$) and perform the simulations for several values of $\Delta\beta$. We then extrapolate to the limit $\Delta\beta\rightarrow 0$ using a quadratic $\Delta\beta$ dependence that characterizes the symmetric Trotter decomposition, thus removing the systematic error introduced by the finite $\Delta\beta$. Results for multiple lattice sizes $N_{L}^{3}$ for a given particle number $N$ are used to extrapolate to the continuum limit $\nu\rightarrow 0$ (see the supplemental material for detailed extrapolation results). In the following we discuss results for two measurable thermal observables: the contact and the thermal energy.

(i) Contact: The expression (2) for the contact can also be written as

[TABLE]

where $F$ is the free energy and the derivative is evaluated at constant temperature $T$. In the lattice formulation the contact can then be calculated from

[TABLE]

where $\langle\hat{V}\rangle$ is the thermal expectation value of the potential energy $\hat{V}=g\sum_{\bf{x}}\hat{n}_{\bf{x},\uparrow}\hat{n}_{\bf{x},\downarrow}$. In Fig. 1 we show our AFMC results for the temperature dependence of the contact $C$ calculated from (6) in the continuum limit in units of $Nk_{F}$ for $N=40$ (solid blue squares). The temperature $T$ is expressed in units of the Fermi temperature $T_{F}=\varepsilon_{F}/k_{B}$, where $k_{B}$ is the Boltzmann constant and $\varepsilon_{F}=(\hbar^{2}/2m)(3\pi^{2}\rho)^{2/3}$ is the Fermi energy for a free gas of density $\rho=\nu/(\delta x)^{3}$. Our results are in excellent agreement with the recent experimental results of the Swinburne group Carcy et al. (2019) (solid purple diamonds) and of the MIT group Mukherjee et al. (2019) (solid red up triangles), both above and below the critical temperature $T_{c}$.

We also compare our results with the theoretical calculations of Refs. Palestini et al. (2010); Gandolfi et al. (2011); Enss et al. (2011); Goulko and Wingate (2016); Rossi et al. (2018); Hu et al. (2011); Leyronas (2011); Liu (2013); Sun and Leyronas (2015) and the low-temperature experimental result of Ref. Hoinka et al. (2013).

Our results for the contact show similar qualitative behavior to those of the lattice diagrammatic Monte Carlo method of Ref. Goulko and Wingate (2016) (open gray diamonds) in the low-temperature regime, but have markedly different qualitative behavior for $T>T_{c}$. Our results above $T_{c}$ are more consistent with the bold diagrammatic Monte Carlo results of Ref. Rossi et al. (2018) (open black circles), but they are systematically lower.

In Fig. 1, we also compare our AFMC results for the contact with those of Ref. Enss et al. (2011) (solid pink line), where good overall qualitative agreement is seen for the entire temperature range. This is somewhat surprising since the work of Ref. Enss et al. (2011) used the Luttinger-Ward approach with uncontrolled systematic errors. However, this method has been shown to produce reliable results for other observables of the UFG Zwerger (2016); Jensen et al. (2018). Quantitatively, our results are above those of Ref. Enss et al. (2011) at low temperatures, and significantly below them for $T>T_{c}$.

Ref. Drut et al. (2011) used an AFMC approach similar to the current work but in the grand-canonical ensemble, and extracted the contact above $T_{c}$ from the tail of the momentum distribution at a finite filling factor. The calculated temperature dependence of the contact in Ref. Drut et al. (2011) is substantially different from our results. As can be seen in Fig. 2 of the supplemental material, the contact is very sensitive to the filling factor, particularly at temperatures $T>T_{c}$, and the continuum extrapolation leads to qualitatively different results.

We tested our continuum extrapolations by comparing the results of different dispersion relations for the single-particle energy. For a finite filling factor $\nu$, the contact depends on the dispersion relation but similar results should be obtained in the limit $\nu\to 0$. In Fig. 4 of the supplemental material, we show the contact for multiple dispersion relations for $N=40$ particles at $T/T_{F}\simeq 0.24$ and demonstrate that they extrapolate to similar values (within statistical errors) in the continuum limit. In the comparison we use a quadratic dispersion (the one implemented in our calculations), the hopping dispersion $\epsilon^{(h)}_{\mathbb{k}}=\frac{\hbar^{2}}{m\delta x^{2}}[3-\sum_{i}\textrm{cos}(k_{i}\delta x)]$ (used in Ref. Goulko and Wingate (2016)), and the dispersion $\epsilon^{(3)}_{k}=\frac{\hbar^{2}k^{2}}{2m}[1-\alpha(\frac{k\delta x}{\pi})^{2}]$ with $\alpha=0.257022$ Carlson et al. (2011).

The inset of Fig. 1 shows the continuum contact results for $N=40,66$, and $114$ particles. The results for $N=66$ and $114$ particles show little systematic difference from the $N=40$ particle results, although the results for the latter have smaller statistical errors. This suggests that our results for the contact are close to the thermodynamic limit.

Our calculations are limited to $T\lesssim 0.45\,T_{F}$. Large lattice simulations with lower filling factors are necessary to determine the contact at higher temperatures up to $T/T_{F}\approx 1$, where a meaningful comparison with the virial expansion results can be made.

(ii) Thermal energy: We also calculated the thermal energy $E=\langle\hat{H}\rangle$ of the UFG (in units of the non-interacting Fermi gas energy at zero temperature $E_{\rm FG}=\frac{3}{5}N\varepsilon_{\rm F}$) as a function of temperature $T$ (measured in units of the Fermi temperature $T_{F}$). In Fig. 2 we show our AFMC results for $E/E_{FG}$ as a function of $T/T_{F}$ in the continuum limit for $N=40$ (solid squares) and $N=66$ (solid circles) particles. We compare our results with the experimental results of Ref. Ku et al. (2012) (open circles), the AFMC results of Ref. Drut et al. (2012) (open squares) and the zero-temperature quantum Monte Carlo result of Ref. Carlson et al. (2011) (open triangle).

In the high-temperature regime we find good quantitative agreement between our results and those of Refs. Drut et al. (2012) and Ku et al. (2012). Below the critical temperature $T_{c}\simeq 0.15\;T_{F}$, the AFMC results of Ref. Drut et al. (2012) are systematically above our results. This is anticipated since the results of Ref. Drut et al. (2012) were calculated at a finite filling factor of $\nu\simeq 0.040-0.045$ (corresponding to a non-negligible effective range parameter $k_{F}r_{e}\simeq 0.36$ for the quadratic dispersion relation), while in the current work we use a continuum extrapolation to remove the systematic error associated with a finite filling factor. When comparing to the experimental results of Ref. Ku et al. (2012), our results for $N=40$ and $66$ particles are systematically lower in the superfluid regime.

We can use our low-temperature results to extract the Bertsch parameter $\xi$ defined by $E(T=0)=\xi E_{FG}$. Taking an average of its values for our lowest two temperatures and for both $N=40$ and $N=66$ particles, we find $\xi=0.367(7)$. In Table 1 we compare values of the Bertsch parameter determined from recent experimental and theoretical works. Our results are in agreement with the value $\xi=0.372(5)$ found in the $N=66$ ground-state quantum Monte Carlo calculation of Ref. Carlson et al. (2011), and with the lattice quantum Monte Carlo result $\xi=0.366^{+0.016}_{-0.011}$ of Ref. Endres et al. (2013). Our value for $\xi$ also agrees with the experimental value $\xi=0.376(5)$ of Ref. Ku et al. (2012).

Conclusions.— We carried out canonical-ensemble AFMC simulations for the UFG on a lattice using a quadratic single-particle dispersion relation for $N=40,66$ and $114$ particles. Our results for each particle number include extrapolations to the continuum limit of zero filling factor $\nu\rightarrow 0$. In particular, we have calculated the temperature dependence of the contact across the superfluid phase transition, and find excellent agreement with the recent experimental results of Refs. Carcy et al. (2019); Mukherjee et al. (2019). Among various existing calculations of the temperature dependence of the contact, our AFMC results provide the best quantitative agreement with these recent experiments. We also calculated the thermal energy as a function of temperature and estimated a value of $\xi=0.367(7)$ for the Bertsch parameter, in agreement with the experimental value and with zero-temperature quantum Monte Carlo calculations.

Acknowledgments.— We thank K. Van Houcke, N. Navon, and F. Werner for useful discussions. We also thank J. Carlson, J.E. Drut, T. Enss, O. Goulko, M.J.H. Ku, B. Mukherjee, P. Pieri, K.E. Schmidt, C. J. Vale, and M.W. Zwierlein for providing the data shown in Fig. 1 and Fig. 2.

This work was supported in part by the U.S. DOE grants Nos. DE-FG02-91ER40608, DE-SC0019521, and DE-FG02-00ER41132. The research presented here used resources of the National Energy Research Scientific Computing Center, which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231. We also thank the Yale Center for Research Computing for guidance and use of the research computing infrastructure.

I Supplemental material: The contact in the unitary Fermi gas across the superfluid phase transition

II Finite-temperature AFMC

We use auxiliary-field quantum Monte Carlo (AFMC) methods SAlhassid2017 ; SJensen2018 ; SJensen2018-2 on a spatial lattice to calculate thermal expectation values of observables in the canonical ensemble. The method is based on a Hubbard-Stratonovich representation of $e^{-\beta\hat{H}}$, where $\beta=1/k_{B}T$ is the inverse temperature (with $k_{B}$ the Boltzmann constant). Dividing the imaginary time $\beta$ into $N_{\tau}$ time slices of length $\Delta\beta$, we use a symmetric Trotter decomposition of $e^{-\beta\hat{H}}$ and a Gaussian Hubbard-Stratonovich transformation for each lattice site ${\bf x}$ and discretized imaginary time $\tau_{n}=n\Delta\beta$ ($n=1,2,...,N_{\tau})$. This results in a path integral over auxiliary fields $\sigma_{\bf{x}}(\tau_{n})$:

[TABLE]

where $G_{\sigma}$ is a Gaussian weight and $\hat{U}_{\sigma}$ is a propagator of non-interacting particles moving in external auxiliary fields $\sigma_{\bf x}(\tau)$. The thermal expectation value of an observable $\hat{O}$ is then given by

[TABLE]

where $\Phi_{\sigma}=\textrm{Tr}(\hat{U}_{\sigma})/|\textrm{Tr}(\hat{U}_{\sigma})|$ is the Monte Carlo sign, $W_{\sigma}=G_{\sigma}|\textrm{Tr}(\hat{U}_{\sigma})|$ is a positive-definite weight, and $\langle\hat{O}\rangle_{\sigma}=\textrm{Tr}(\hat{O}\hat{U}_{\sigma})/\textrm{Tr}(\hat{U}_{\sigma})$ is the thermal expectation value of the observable $\hat{O}$ for the auxiliary-field configuration $\sigma$. Here we use the canonical ensemble, so the traces are evaluated for fixed particle numbers $N_{\sigma}$  SAlhassid2017 ; SGilbreth2013 ; SJensen2018 using the method of Ref. SGilbreth2015 .

III Data Analysis

The symmetric Trotter decomposition we use produces an error $O((\Delta\beta)^{2})$ for small imaginary time step $\Delta\beta$. In Fig. 1, we show extrapolations in $\Delta\beta$ for the contact with $N=40$ particles and lattice size $9^{3}$ at temperatures (a) $T/T_{F}=0.353$, (b) $T/T_{F}=0.202$, and (c) $T/T_{F}=0.149$, where a linear fit has been carried out in $(\varepsilon_{F}\Delta\beta)^{2}$ for small $\Delta\beta$ ($\varepsilon_{F}$ is the Fermi energy of the free gas).

A significant systematic error is due to the finite filling factor $\nu$ of the simulations. In panels (d)-(f) of Fig. 1 we show the continuum extrapolations $\nu\to 0$ of the contact at several temperature (after the $\Delta\beta\rightarrow 0$ extrapolation), where a linear fit in $\nu^{1/3}$ is carried out for low values of the filling factor $\nu$. In Fig. 2 we show the contact as a function of temperature for several values of the filling factor $\nu$ at constant number of particles $N=40$ (panel (a)) and $N=66$ (panel (b)). We observe that the contact is particularly sensitive to finite filling factor effects. The extrapolated values for $\nu\to 0$ are also shown by the solid squares in panel (a) and solid circles in panel (b).

IV Comparison of various dispersion relations

Several dispersion relations for the dependence of the single-particle energy on momentum were used in the literature SBulgac2006 ; SBurovski2006 ; SBurovski2006-2 ; SBulgac2008 ; SGoulko2010 ; SCarlson2011 ; SGoulko2016 ; SJensen2018 for the UFG. The results shown in the main text use a quadratic dispersion relation as in Refs. SBulgac2006 ; SBulgac2008 ; SJensen2018 . In Figs. 3 and 4 we compare results obtained for different dispersion relations to further test our continuum limit extrapolations. We consider the following dispersion relations

[TABLE]

where $\epsilon^{(2)}_{k}$ is the quadratic dispersion, $\epsilon^{(h)}$ is the standard hopping relation used in Refs. SBurovski2006 ; SBurovski2006 ; SBurovski2006-2 ; SGoulko2010 ; SCarlson2011 ; SGoulko2016 ($\delta x$ is the lattice spacing), and $\epsilon^{(3)}_{k}$ is a quartic dispersion introduced in Ref. SCarlson2011 with $\alpha=0.257022$. Each dispersion relation has a different dependence on the filling factor with different effective range parameters $r_{e}$ and $R_{e}$ SWerner2012 .

Using the method of Ref. SPricoupenko2007 , we calculated the two-particle energies with center-of-mass wavevector $\mathbb{K}=0$ for lattices of size up to $41^{3}$. In Fig. 3, we show the lowest two such energies as a function of $\nu^{1/3}$ for the dispersion relations in Eqs. (9). We see that various dispersion relations exhibit a different dependence on $\nu^{1/3}$ but they all extrapolate to the same energies in the continuum limit.

In Fig. 4 we show continuum extrapolations of the contact for $N=40$ particles and $T/T_{F}\simeq 0.235$ using the dispersions $\epsilon^{(2)}_{k}$, $\epsilon^{(h)}_{\mathbb{k}}$ and $\epsilon^{(3)}_{k}$ in Eqs. (9). Carrying out AFMC calculations on lattices of size $7^{3},9^{3},11^{3},13^{3}$, and $15^{3}$, and performing a linear extrapolation in $\nu^{1/3}$ for values of $\nu^{1/3}$ below $\sim 0.4$, we find that the extrapolated values for the different dispersions agree within their statistical errors.

V Momentum distribution

The momentum distribution $n_{k}=\langle\hat{a}^{\dagger}_{k}\hat{a}_{k}\rangle$ (we suppress the spin index $\sigma$ as the distribution is independent of spin for the spin-balanced case) is shown in Fig. 5(a) for $N=40$ particles and temperature of $T/T_{F}=0.235$ for lattice sizes $7^{3},11^{3}$ and $15^{3}$ (open symbols). The momentum distribution is broadened by both the interaction and temperature.

In Fig. 5(b) we show the scaled momentum distributions $3\pi^{2}(k/k_{F})^{4}n_{k}$ of Fig. 5(a). For reference we also show the values of the contact $C/(Nk_{F})$ for lattice sizes of $7^{3},11^{3}$ and $15^{3}$, calculated from the expectation value of the potential energy $\langle\hat{V}\rangle$ using Eq. (6) (horizontal lines). We observe that for the smaller lattice size of $7^{3}$ there is a substantial difference between the scaled tail of the momentum distribution and the value of the contact extracted from the potential energy, while this difference becomes much smaller for larger lattice sizes. This stronger lattice size dependence of the tail makes reliable extraction of the contact from the tail of the momentum distribution challenging. In this work, we therefore extracted the contact from used the average potential energy.