The Fermi–Hubbard model at unitarity

We simulate the dilute attractive Fermi–Hubbard model in the unitarity regime using a diagrammatic determinant Monte Carlo (MC) algorithm with worm-type updates. We obtain the dependence of the critical temperature on the filling factor ν and, by extrapolating to ν → 0, determine the universal critical temperature of the continuum unitary Fermi gas in units of Fermi energy: Tc/εF = 0.152(7). We also determine the thermodynamic functions and show how the MC results can be used for accurate thermometry of a trapped unitary gas.


Introduction
In recent years, ultracold atomic systems have served as a controlled and tunable toolbox for studying many-body quantum phenomena. The continuous tunability of the interaction by Feshbach resonances makes these systems ideal candidates to study the crossover from momentum space pairing in the Bardeen-Cooper-Schrieffer (BCS) theory to Bose-Einstein condensate (BEC) of fermions bound into molecules. This BCS-BEC crossover has been one of the most studied problems in recent experiments in both magnetic or optical traps [1]- [6] and optical lattices [7].
Tuning across the Feshbach resonance, one traverses the whole range of the gas parameter k F a s , where k F is the Fermi momentum and a s is the s-wave scattering length. The regime of k F |a s | 1 (negative a s ) is described by BCS theory. At k F a s 1 (positive a s ) fermions pair into bosonic molecules and form a BEC. At the microscopic scale, the BCS and BEC regimes are radically different; however, the macroscopic and, in particular, critical behaviour is supposed to be qualitatively the same for the whole range of k F a s : the system undergoes the superfluid (SF) phase transition at a certain critical temperature. Separating the BCS and BEC extremes is the socalled unitarity point (k F a s ) −1 → 0. It is worth noting that the unitarity regime is approximately realized in the inner crust of the neutron stars, where the neutron-neutron scattering length is nearly an order of magnitude larger than the mean interparticle separation [8].
The Fermi gas at unitarity is a peculiar case of a strongly interacting system with no interaction-related energy scale: the divergent scattering length and any related energy scale drop out completely. This gives rise to universality of the dilute gas properties, in the sense that the only relevant energy scale left in the system is given by the density, n. Because of this 3 Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT universality one obtains a unified description of such diverse systems as cold atoms in magnetic or optical traps, Fermi-Hubbard model in optical lattices and inner crusts of neutron stars.
The theoretical description of the Fermi gas in the BCS-BEC crossover regime is a major challenge, since the system features no small parameter on which one could build a theory in a rigorous way. The original analytical treatments were confined to zero temperature and were based on the extension of the BCS-type many-body wavefunction [9]. Most of the subsequent elaborations are also of mean-field type (with or without remedies for the effects of fluctuations) [10]- [16]. The accuracy and reliability of such approximations is questionable since they inevitably involve an uncontrollable approximation.
Numerical investigations of the unitary Fermi gases are hampered by the sign problem, inherent to any Monte Carlo (MC) simulations of fermion systems [17,18]. One way of avoiding the sign problem at the expense of a systematic error is the fixed-node MC framework, which has been used to study the ground state [19]. The systematic error of the fixed-node MC depends on the quality of the variational ansatz for the nodal structure of the many-body wavefunction and is not known precisely. Only in a few exceptional cases can the sign problem be avoided without incurring systematic errors. One of such cases is given by fermions with attractive contact interaction, for which a number of sign-problem-free schemes has been introduced [20]- [23]. Fortunately, this system can be tuned to the unitarity regime. Still, despite a number of calculations at finite temperature [24]- [26], an accurate description of the finite-temperature properties of the unitary Fermi gas is missing.
In the present paper, we simulate the Fermi-Hubbard model in the unitary regime by means of a determinant diagrammatic MC method. By studying the dilute limit of the model, we extract properties of the homogeneous continuum Fermi gas. A brief summary of the main results has been given in [27]. Here, we provide a detailed description of the MC scheme and methods of data analysis. We also report new results relevant to experimental realizations of the Fermi-Hubbard model in optical lattices and trapped Fermi gases.
The Fermi-Hubbard model is defined by the Hamiltonian Here, c † kσ is a fermion creation operator, n xσ = c † xσ c xσ , σ =↑, ↓ is the spin index, x enumerates L 3 sites of the three-dimensional (3D) simple cubic lattice with periodic boundary conditions, the quasimomentum k spans the corresponding Brillouin zone, k = 2t 3 α=1 (1 − cos k α a) is the single-particle tight-binding spectrum, a and t are the lattice spacing and the hopping amplitude, respectively, µ stands for the chemical potential and U < 0 is the on-site attraction. Without loss of generality we henceforth set a and t equal to unity; the effective mass at the bottom of the band is then m = 1/2.
In section 2, we will study the two-body problem at zero temperature, show how the Hubbard model can be used to study the continuum unitary gas and investigate the functional structure of lattice corrections to the continuum behaviour. In section 3, we discuss the finitetemperature diagrammatic expansion for the Hubbard model (subsection 3.1), and give a 4 Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT Figure 1. Diagrammatic series for the vertex insertion (ξ, p) (heavy dot). Small dots represent the bare Hubbard interaction U, and lines are the single-particle propagators.
qualitative description of the MC procedure to sum the diagrammatic series (subsection 3.2), with details of the updating procedures given in appendix A. In order to extract the thermodynamic limit properties from MC data, we use finite-size scaling analysis described in section 4. Section 5 gives an overview of the scaling functions describing thermodynamics of the unitary gas, and results are presented and discussed in section 6.

Two-body problem
Consider a quantum mechanical problem of two fermions at zero temperature described by the Hamiltonian (1.1a)-(1.1c) with µ = 0. The most straightforward way to tackle this problem is within the diagrammatic technique in the momentum-frequency representation [28,29], which, in the present case, is built on four-point vertices, U, with two incoming (spin-↑ and spin-↓) and two outgoing (spin-↑ and spin-↓) ends, connected by single-particle propagators. The scattering of two particles is then described by a series of ladder diagrams [16,29] shown in figure 1. Ladder diagrams can be summed by introducing the vertex insertion (ξ, p), which depends on frequency ξ and momentum p. Since (0, 0) is proportional to the scattering amplitude, the unitarity limit corresponds to (ξ → 0, p → 0) → ∞. The summation depicted in figure 1 leads to −1 (ξ, p) = U −1 + (ξ, p), (2.1) where (ξ, p) is the polarization operator (the integration is over the Brillouin zone) It immediately follows from equations (2.1) and (2.2) that the unitary limit corresponds to U = U * , where A straightforward numeric integration yields U * ≈ −7.915t.
In the limit of vanishing filling factor ν → 0 for the many-body problem of (1.1a)-(1.1c), the typical values of ξ and p are related to the Fermi energy ξ ∼ F ∼ ν 2/3 and Fermi 5 Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT momentum p ∼ k F ∼ ν 1/3 and are small compared to the bandwidth and reciprocal lattice vector, respectively. In zeroth order with respect to ν, the lattice system is identical to the continuum one. Indeed, by combining (2.1), (2.2) and (2.3), we get 4) and observe that for small ξ and p, one can replace ε(k) with k 2 and extend integration over dk to the whole momentum space with the result With this form of cont , and particle propagators based on the parabolic dispersion relation one recovers the continuum limit behaviour. Now we are in a position to treat the lattice corrections. The first correction should come from , not from propagators, since only in large momenta play a special role due to resonance in the two-particle channel. In the lowest non-vanishing order in ξ and p, we have (the summation over repeating subscripts is implied) with the difference between the lattice and continuous model given by (2.9) In the limit of ξ → 0 and p → 0, we have −1 ≈ −1 cont ∼ k F ∼ ν 1/3 , and −1 − −1 cont ∼ k 2 F ∼ ν 2/3 . Hence, the leading lattice correction is of the form Incidentally, equations (2.8) and (2.9) hint into an intriguing possibility of completely suppressing the leading-order lattice correction by tuning the single-particle spectrum k , so that A = B = 0. We did not explore this possibility in the present study.

Determinant diagrammatic MC
The diagrammatic technique employed in the previous section is not particularly convenient for numerical studies. In this section, we briefly review the Matsubara technique and then present a MC scheme of summing the resultant diagrammatic series.

Rubtsov's representation
To construct a diagrammatic expansion for the model (1.1a)-(1.1c), we follow [22,30] and consider the statistical operator in the coordinate-imaginary time representation. In the interaction picture, we get where β is an inverse temperature, H 1 (τ) = e τH 0 H 1 e −τH 0 , and T τ stands for the imaginary time ordering. Expanding equation (3.1) in powers of H 1 , one obtains for the partition function Expansion (3.2) generates the standard set of Feynman diagrams. Graphically, the diagrams are similar to those of section 2, and consist of the four-point vertexes, U, connected by the singleparticle propagators The pth order of the perturbation theory is then graphically given by a set of (p!) 2 possible interconnections of vertices by propagators, see figure 2.
The diagrammatic expansion (3.2) is unsuitable for the direct MC simulation since it has a sign problem built in: different terms in the series have different signs-a closed fermion loop brings in an extra minus sign [28]. The trick is to consider all diagrams of a given order p with the fixed vertex configuration and A σ (S p ) are the p × p matrices built on the single-particle propagators For equal number of spin-up and spin-down particles det A ↑ det A ↓ = |detA| 2 , and the sign problem is absent. 6  The following two-point pair correlation function will prove useful where P(x, τ) and P † (x , τ ) are the pair annihilation and creation operators in the Heisenberg picture, respectively: P(x, τ) = c x↑ (τ)c x↓ (τ). The nonzero asymptotic value of G 2 (xτ; x τ ) as |x − x | → ∞ is proportional to the condensate density. Feynman diagrams for g 2 (xτ; x τ ) are similar to those for Z, but contain two extra elements: a pair of two-point vertices with two incoming (outgoing) ends which represent P(x, τ) (P † (x , τ )), see figure 3. The vertex configurations for the correlation function (3.7) slightly differ from those for the partition function The partially summed diagrammatic expansion for g 2 (xτ; x τ ) is similar to equation (3.4) whereÃ σ (S p ) is a (p + 1) × (p + 1) matrix which differs from equation (3.6) only by that it has an extra row i 0 and an extra column j 0 such thatÃ Below we deal only with equal number of spin-up and spin-down fermions and (for the sake of brevity) suppress the spin indexes wherever possible. We also peruse the generic notation (with superscripts) D(S p ) for pth order terms of the diagrammatic expansions similar to (3.4), e.g., D (Z) (S p ) = (−U ) p |det A(S p )| 2 . To simplify the notation, we also omit superscripts if this does not lead to ambiguity.

Diagrammatic MC and Worm algorithm
Equations (3.4) and (3.10) have similar general structure of a series of integrals and sums with ever increasing number of integration variables and summations. In [31,32], it has been shown how to arrange a numerical procedure that sums such convergent series. To this end one considers the space of all possible vertex configurations S (for the series (3.4) S ≡ S (Z) , while for the series (3.10) S ≡ S (G) ), with the 'weight' D(S p ) associated with each element of the space. One then uses the Metropolis principle [33] to arrange a stochastic Markov process which sequentially generates vertex configurations S p according to their weights D(S p ), thus sampling the space S . In the course of sampling, one also collects statistics for observables in the form of MC estimators, see subsection 3.3.
The stochastic process consists of elementary MC updates performed on vertex configurations S p . The set of updates is problem-specific being restricted only by the requirements of (i) the ergodicity, i.e., given a particular diagram S p , it takes a finite number of steps to convert it into any other diagram S q , and (ii) detailed balance, i.e., relative contributions of diagrams S p and S q to the statistics is given by the ratio of their weights, D(S p )/D(S q ). The set of updates satisfying these requirements is not unique, the freedom is used to maximize the efficiency of simulations, as explained in detail in appendix A.
In view of close similarity between the diagrammatic expansions (3.4) and (3.10), it is advantageous to construct a MC process which samples these two series in a single simulation. This way one has access to both diagonal, e.g., energy, and off-diagonal properties, e.g., the superfluid response. An efficient way of performing such concurrent simulation is provided by the worm algorithm, which was originally devised for the worldline MC simulations [34]. In the context of the diagrammatic determinant MC, the generic worm algorithm principles imply the following. Firstly, one works in the joint configuration space S (Z) ∪ S (G) , accommodating diagrams (3.3) and (3.9). Secondly, all the updates are performed exclusively in terms of the two-point vertices P(x, τ) and P † (x, τ)-through their creation/annihilation, motion, and 'interactions' with adjacent vertexes. The worm-type updating procedures are further detailed in appendix A. Within the worm-algorithm framework, the configuration spaces S (Z) and S (G) are disjoint subsets of one extended configuration space. In what follows, we will refer to them as Z-(or 'diagonal') and G-(or 'off-diagonal') sectors of the configuration space.

9
Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT Formally, the extended configuration space corresponds to the generalized partition function (3.11) where the value of ζ is arbitrary-it controls the relative statistics of Z-and G-sectors and the efficiency of simulation.

MC estimators
Suppose we have an observable X(α), which depends on a set of variables α, e.g. temperature and chemical potential. A MC estimator for the observable X is an expression which, upon averaging over the sequence of MC configurations converges to the expectation value of X(α).
In accordance with equation (3.11), the simplest worm-algorithm MC estimators are and Their MC averages are and where . . . MC denotes averaging over the set of stochastically generated configurations. In particular, for our purposes, it will be quite useful that The general rules for constructing an estimator for a quantity X(α) specified by the diagrammatic expansion (3.17) are standard. We adopt a convenient convention: if the actual summation in (3.17) involves only a subset S 0 of vertex configurations-a typical example is an expansion defined within Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT the Z-sector only-then we extend summation over the entire configuration space by simply defining D (X) (S p ∈ S 0 ) ≡ 0. If vertex configurations S p are sampled from the probability density D (Z W ) (S p ) coming from the expansion for the generalized partition function: as In what follows, by the estimator for a quantity x(α) = X(α)/Z(α), we basically mean corresponding function Q (X) (α; S p ).

Estimators for number density and kinetic energy.
We start with the estimator for the number density. The expectation value of the number density reads Here, (x, τ) is an arbitrary space-time point (the system is space/time translational invariant) and σ is one of the two spin projections; the factor of 2 comes from the spin summation. The diagrammatic expansion of the numerator is similar to that for Z, with the diagram weight given by is a similar (p + 1) × (p + 1) matrix with an extra row and a column, corresponding to the extra creation and annihilation operators in the numerator of (3.21), respectively. This immediately leads to the following estimator (3.20) for the number density We utilize the freedom of choosing (σ, x, τ) to suppress autocorrelations in measurements. The density measurement starts with randomly generated (σ, x, τ). The estimator for kinetic energy is derived similarly. One employs the coordinate-space expression for the kinetic energy in terms of hopping operators and deals with a slightly generalized version of equation (3.21) The rest is identical to the previous discussion up to replacement B σ p+1 (S p , x, τ) → B σ p+1 (S p , x, g, τ), since now the spatial position of the creation operator is shifted from that of the annihilation operator by the vector g. In our case, only the nearest-neighbour correlator (3.24) has to be computed.

11
Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT

Estimator for the interaction energy. The estimator for the interaction energy
is readily constructed using a generic trick of finding the expectation value of operator in terms of which the perturbative expansion is performed. Consider the Hamiltonian H(λ) = H 0 + λH 1 and observe that Tr The differentiation of equation (3.2) for Z = Z(λ) and letting λ = 1 afterwards is straightforward since the diagram of order p is proportional to λ p . Hence

Estimator for the integrated correlation function.
Following the general treatment of [32], one can construct an estimator for the correlation function (3.7). In this work, we just need to sum and integrate this correlator over all its variables (see section 4) (3.28) and the estimator for K(L, T ) is particularly simple

Extrapolation towards macroscopic continuum system
The MC setup discussed in section 3 works in the grand canonical ensemble with external parameters {L, T, µ}. In order to extract the critical temperature of a continuum gas from MC data, one has to perform the two-step extrapolation. (i) Upon taking the limit of L → ∞, one obtains T c (µ), the critical temperature of a lattice system at a given chemical potential, and translates it into T c (ν) by extrapolating the measured filling factor to the infinite system size: The extrapolation towards the continuum limit is then done by taking the limit of ν → 0. The latter procedure is based on results presented in section 2.
The finite-size extrapolation is performed by considering a series of system sizes L 1 < L 2 < L 3 · · · . At the critical point the correlation function (3.7) decays at large distances as a power-law: where η is the anomalous dimension [35]. Since we expect the transition to belong to U(1) universality class, we take η ≈ 0.038. If one re-scales the summed correlator (3.28) according to the corresponding quantity is supposed to become size-independent at the critical point, i.e., the crossing of R(L i , T ) and R(L j , T ) curves can be used to obtain an estimate T L i ,L j for the critical temperature T c (µ) [36]. Indeed, for temperatures in the vicinity of the critical point the correlation length diverges as ξ corr ∝ |t| −ν ξ , where t = (T c (µ) − T )/T c (µ), and ν ξ ≈ 0.671 for the U(1) universality class. In the renormalization group (RG) framework [35], the finite-size scaling of the re-scaled correlator R obeys the relation is the universal scaling function analytic at x = 0, c is a non-universal constant, ω ≈ 0.8 is the critical exponent of the leading irrelevant field, see e.g. [37], and dots represent higher-order corrections to scaling. If the irrelevant field corrections were not present, all R(L i , T ) curves would intersect at a unique point, T c (µ). Expanding equation (4.2) up to terms linear in t one obtains for the crossing T L i ,L j To employ equation (4.3), one performs a linear fit of the sequence of T L i ,L j against the right-hand side of equation (4.3) for several pairs of system sizes. The intercept of the best-fit line yields the thermodynamic limit critical temperature T c (µ). Note that if the universality class is not U(1) and the values of η, ν ξ , ω are different, or system sizes are too small to justify the scaling limit, the whole procedure fails. Hence, the adopted scheme of pinpointing T c features a built-in consistency check.

Thermodynamic scaling functions of a unitary Fermi gas
As has been noted in section 1, the only relevant microscopic energy scale in the continuum unitary Fermi gas is given by the Fermi energy ε F = κh 2 n 2/3 /m, where κ = (3π 2 ) 2/3 /2 for a two-component Fermi gas. Therefore, as was first noticed in [38], all thermodynamic potentials feature self-similarity properties and can be written in terms of dimensionless scaling functions of the dimensionless temperature x = T/ε F . All scaling functions are mutually related; it is sufficient to know just one of them to unambiguously restore the rest. Apart from the shape of scaling functions, the self-similarity at the unitary point is identical to that of a non-interacting Fermi gas [39], including functional relations between different thermodynamic potentials. A derivation of scaling functions and relations between them can be found in [38]. In this section, we render the scaling analysis in a form convenient for our MC study. In terms of the dimensionless chemical potential y = µ/ε F , the dimensionless equation of state reads y = f µ (x). The f µ function can be calculated numerically. Another quantity which is also available in our simulation is the dimensionless energy per particle E/(Nε F ) = f E (x). The scaling relations for other thermodynamic quantities are defined likewise. For instance, the entropy and pressure read S/N = f S (x), and P/(nε F ) = f P (x). Though f S and f P are not directly calculated in our simulation, and we can relate them to f E . It is also important to relate f µ to f E , since this yields a consistency check for the numerical results.
To establish desired relations, we start with the scaling relation for the Helmholtz free energy F/(Nε F ) = f F (x) which in canonical variables reads where γ = κh 2 /m. Then, for the entropy and pressure, we have (the prime stands for the derivative) The expression for energy in terms of f F is We thus see that On may also consider equation (5.4) as a differential equation to be integrated with respect to f F from x = ∞ down to finite x, taking advantage of known asymptotic behaviour of f F and f E for the weakly interacting two-component gas. Now, we note that from the general thermodynamic relation E = −PV + TS + µN, it immediately follows that which, in turn, leads to the following relations that allow one to express f S and f F functions through numerically available functions f E and f µ . Another useful relation is which allows us to extract f F directly from f E and f µ , and thus provides a simple check for the data consistency: the result (5.9) for the f F curve should be consistent with the derivative deduced from (5.10). By integrating equation (5.6), we get Here, we took into account the asymptotic ideal-gas behaviour of f E and introduced the corresponding term into the integral to render the latter convergent. The free constant of integration, C 0 , can be restored from (5.9)-(5.12) combined with the asymptotic ideal-gas behaviour of f µ , The result is Note that if higher-order terms in the asymptotic (x → ∞) behaviour of f E (x) are also known, then (5.11) can be used to establish the corresponding corrections for f F (x) and other scaling functions. For instance, it has been found in [40] that, as In accordance with (5.11), this implies (x → ∞) Finally, the scaling functions W 0 and G 0 defined in [38] are related to f E and f µ as follows

Trapped Fermi gas
So far we have considered the uniform Fermi gas, while in experimental realizations [1]- [6], one has to deal with the parabolic trapping potential. The standard procedure, especially in systems with short 'healing length' is to use the local density approximation (LDA), i.e., to replace the chemical potential with its coordinate-dependent counterpart µ(r) = µ − V(r). This procedure can be easily combined with MC results as follows. We introduce the dimensionless variable u = µ/T = f µ (x)/x and define the scaling function for the number density as w n = x −3/2 . This is equivalent to the parametric {u(x), w n (x)} dependence of n on u. The scaling functions for other thermodynamic quantities are defined in a similar manner, e.g. w E (u) ≡ f E (x(u)) for the energy, and w S ≡ f S (x(u)) for the entropy. Within LDA, u acquires the coordinate dependence u(r) = (µ − V(r))/T which translates into the density profile n(r; µ, T ) = w n (u(r))(mT/κh 2 ) 3/2 .
Likewise, other thermodynamic functions are to be understood as local, coordinate-dependent quantities. Consider the case of N particles in a cigar-shaped parabolic trap, characterized by the axial and radial frequencies ω and ω ⊥ . The characteristic energy in this case is E F = (3N ) 1/3h3 (ω 2 ω ⊥ ) 1/3 , which would coincide with the Fermi energy for a non-interacting gas in the trap. Note, that we denote the Fermi energy in the trap by an upright capital E F in order to avoid confusion with the uniform system Fermi energy ε F . By integrating over the radial coordinates, one obtains the axial density profile n a (z) (z is an axial coordinate) in the form where L = λ 1/3 (3N) 1/6 l , the aspect ratio λ = ω ⊥ /ω , the oscillator length l 2 =h/mω , and Obtaining the temperature dependence of thermodynamic functions for a non-uniform system within LDA is also straightforward. For the total energy of a cloud E tot , we obtain and likewise for the entropy (5.25)

Results and discussion
We performed simulations outlined in previous sections for filling factors ranging from 0.95 down to 0.06 with up to about 300 fermions on lattices with up to 16 3 sites. The typical rank of determinants involved in computations of acceptance ratios (appendix A) and estimators Since we only need ratios of determinants, we use fastupdate formulae [22] to reduce the computational complexity of updates from M 3 down to M 2 . We validate our numerical procedure by comparing results against the exact diagonalization data for a 4 × 4 cluster [41], and other simulations of the critical temperature at quarter filling ν = 0.5 [42,43] and ν = 0.25 [43]. In all cases, we find agreement within statistical errors of a few per cent. Figure 4 shows a typical example of the finite-size analysis outlined in section 4. Despite the fact that numerically accessible system sizes are quite small, figure 4 (and similar analysis for the whole range of filling factors) supports expectations that the universality class for the SF phase transition is U(1). The finite-size analysis allows us to pinpoint the phase transition temperature to within a few per cent.

Critical temperature
Shown in figure 5 is the dependence of the critical temperature on the lattice filling factor. The critical temperature is measured in units of the Fermi energy, as is natural for the unitarity limit. We define the Fermi momentum for a lattice system with filling factor ν as k F = (3π 2 ν) 1/3 and the Fermi energy ε F = k 2 F , as those of a continuum gas with the same effective mass and number density n = ν.
It is clearly seen that the presence of the lattice suppresses the critical temperature considerably, nearly by a factor of 4, depending on the filling factor. Strong dependence of T c on ν is in apparent contradiction with [24], which claims weak or no ν-dependence. This disagreement might be due to the difference in the single-particle spectra ε k used: reference [24] employs the parabolic spectrum with spherically symmetric cutoff, while we use the tightbinding dispersion law. Indeed, equations (2.8)-(2.9) indicate that a particular choice of ε k does influence lattice corrections to T c , which may even have different signs for different ε k . It is also clear from figure 5 that close to half-filling T c is essentially constant, as expected (see, e.g. [44]). The predicted ∼ ν 1/3 scaling (2.10) sets in at about ν ≈ 0.5. We thus use a linear fit T c (ν)/ε F (ν) = T c /ε F − const · ν 1/3 to eliminate lattice corrections in the final result. Such fitting procedure results in the best-fit line given by 0.152(7) − 0.13(2)ν 1/3 . We further analyse the fit residues in order to estimate the effect of the sub-leading lattice corrections which are expected to be proportional to ν 2/3 . As shown in figure 6, such corrections, if any, are smaller than the uncertainty of the ν 1/3 fit. The finite-size scaling of the condensate fraction data from [24]. Raw data points are re-scaled similar to equation (4.1) by the L 1+η factor. Shaded vertical strips represent results for T c /ε F of this work and [24] respectively, solid lines are drawn to guide the eye.
This analysis yields the final result T c /ε F = 0.152 (7) for the continuum uniform gas, which is noticeably below the transition temperature in the BEC limit T BEC = 0.218ε F . Various approximate analytical treatments led in the past to T c either above [10,12,13,15], or below [11,14,16] T BEC .
It is instructive to compare our results for T c to other numerical calculations available from the literature. The simulations of [25] yield T c = 0.05ε F , but at the value of the scattering length which has not been determined precisely. This result most probably corresponds to a deep BCS regime, where the transition temperature is exponentially suppressed. Lee and Schäfer [26] report an upper limit T c < 0.14ε F , based on a study of the caloric curve of a unitary Fermi gas down to T/ε F = 0.14 for filling factors down to ν = 0.5. The caloric curve of [26] shows no signs of divergent heat capacity which would signal the phase transition. This upper limit is consistent with T c (ν = 0.5)/ε F ≈ 0.054, see figure 5.
The Seattle group has performed simulations of the caloric curve and condensate fraction, n 0 , of the unitary gas, [24]. Using 'visual inspection' of the caloric curve shape the critical temperature was estimated in [24] to be T c = 0.22(3)ε F . Unfortunately, the authors did not perform the finite-size analysis and ν → 0 extrapolation. The overall shape of the caloric curve seem to be little affected by the finite volume of the system. This is hardly surprising since even in the thermodynamic limit E(T ) and its derivative dE/dT are continuous at the transition point. These properties also make it hard to use non-quantitative measures for reliable estimates of critical parameters from the E(T ) curve. On the other hand, the condensate fraction which has singular properties at T c does show sizable finite-size corrections, see figure 1 of [24]. At this point, we note that scaling of the condensate fraction is identical to that for K(L, T ). In figure 7, we plot the data of the Seattle group as n 0 L 1+η versus temperature. The intersection of scaled curves turns out to be inconsistent with the estimate for T c derived from the caloric curve inspection.

Thermodynamic functions
The filling factor dependence of thermodynamic quantities is similar to that of T c figure 8 displays the behaviour of energy and chemical potential along the critical line T = T c (ν). The extrapolation towards ν → 0 yields for the continuum gas The numerical values for other thermodynamic functions at criticality can be easily restored using the formulae of section 5. In order to elucidate the thermodynamic behaviour of the unitary gas, we performed simulations for a range of temperatures T > T c . Shown in figure 9 are the simulation results for energy and chemical potential for the continuum gas as functions of temperature. Each point was obtained using data analysis similar to that depicted in figure 8. In the high-temperature region, we simulated up to 80 fermions on lattices with up to 32 3 sites. In this region, the condition ν 1 is necessary but not sufficient for extrapolation to the continuum limit, for it is crucial to keep temperature much smaller than the bandwidth: T 6t. As can be seen from figure 9, our results for both energy and chemical potential approach the virial expansion [40] as T/ε F → ∞. For T/ε F 0.5, our data are not far from the curve of [24]. Though we do not have data points for T < T c , there is still a reasonable agreement even at T c with the T → 0 fixed-node MC values [19]. In this region, our results are consistent with a very weak dependence of energy and chemical potential on temperature, and the numeric values of both are consistent with the experimental results [3]- [5].
Using equation (5.9) and data from figure 9, we deduce the dependence of free energy on temperature, see figure 10. We also use equation (5.10) to make sure that our MC data for energy and chemical potential are consistent with each other.  [40] (plus the first virial Fermi correction), black triangles are the MC results of [24], and the purple stars denote the ground-state fixed-node MC results [19].

Trapped gas
As discussed in subsection 5.1, the thermodynamic functions for the uniform case can be used for analysis of experimental system within the LDA. In this section, we report preliminary results of our ongoing study of trapped gas experiments. It starts with the interpolation procedure which produces continuous functional behaviour for thermodynamic functions, consistent with the discrete set of simulated points. We use a piecewise-cubic ansatz with a smooth crossover to the virial expansion, equation (5.16), for the free energy. Temperature dependence of both energy and entropy are then deduced using numerical integration of equations (5.24) and (5.25) respectively.
As the trapped gas is cooled down, the superfluidity first sets in at the centre of the trap, where the density is the highest. Equation (5.22) can be used to pinpoint this onset temperature: at T = T c , µ/T = (µ/ε 0 F )/(T/ε 0 F ), where ε 0 F is the Fermi energy of the uniform gas with the density equal to the density at the trap centre. Using T c /ε 0 F = 0.152(7) and equation (6.2), one obtains T c /E F = 0.20 (2). We quote here a conservative estimate for the uncertainty, which incorporates both the uncertainty of the critical temperature itself, and a systematic uncertainty which stems  [24], purple star denotes the ground-state fixed-node MC result [19], black dotted line shows the Boltzmann gas curve, and the blue dashed is the asymptotic prediction of [40].
from restoring the continuous functional dependence of the chemical potential out of the finite set of the MC calculated points with finite errorbars. Experimentally, the temperature of the strongly interacting Fermi gas is not easily accessible. On the contrary, thermometry of the non-interacting Fermi gases is well established. In the adiabatic ramp experiments, one starts from the non-interacting gas at some temperature (in units of Fermi energy) (T/T F ) 0 , and slowly ramps magnetic field towards the Feshbach resonance [6], thus adiabatically connecting the system at unitarity to a non-interacting one. Assuming the entropy conservation during the magnetic field ramp, equation (5.25) can be employed for the thermometry of the interacting gas: by matching the entropy of a non-interacting gas at the temperature (T/T F ) 0 with the entropy calculated via equation (5.25), one relates the initial temperature (before the magnetic field ramp) to the final temperature (after the ramp). We find that the onset of the superfluidity corresponds to (T/T F ) 0 = 0.12 ± 0.015 (again, we quote here the most conservative estimate for the errorbar). This value seems to be somewhat lower than the value suggested by the experimental results [6]. Nevertheless, given the level of the noise in figure 4 of [6], the consistency is reasonable.
An alternative thermometry can be build on recent advances in the experimental technique [3,4] which made it possible to directly image the in situ density profiles of the interacting system. Such density profiles can be directly fit to equation (5.20), which gives the shape of the cloud depending on µ/T and T/E F . By relating the chemical potential to the temperature using equation (5.22), one is left with only one fitting parameter, T/E F (apart from a trivial fitting parameter z 0 which accounts for the overall shift of the cloud image off the trap centre).
As an illustrative example of such procedure, we have analysed the experimental density profiles measured by the Rice group [4], as depicted in figure 11. From this analysis, we deduce the upper bound for the temperature in the experiments [4] T < 0.1E F , which is consistent with the results of the measurements of the condensate fraction [2]. Since in the experiments [2,4] the gas is very degenerate, one is able to put an upper limit on temperature only.
Note that if the temperature is known from, e.g. the adiabatic ramp experiments, equations (5.20)-(5.22) must reproduce the cloud shape without free parameters (apart from z 0 ).

Conclusions
We have developed a worm-type scheme within the systematic-error-free determinant diagrammatic MC approach for lattice fermions. We applied it to the Hubbard model with attractive interaction and equal number of spin-up and spin-down particles. At finite densities, the model describes ultracold atoms in optical lattice. In the limit of vanishing filling factor, ν → 0, and fine-tuned (to the resonance in the two-particle s-wave channel) on-site attraction, a universal regime sets in, which is identical to the BCS-BEC crossover in the continuous space. In the present work, we confined ourselves to a special value of the on-site interaction, U = U * ≈ −7.915t, corresponding to the divergent s-scattering length. At U = U * and ν → 0, the system reproduces the unitary point of the BCS-BEC crossover. The unitary regime is scaleinvariant: all thermodynamic potentials are expressed in terms of dimensionless scaling functions of the dimensionless ratio T/ε F (temperature in units of Fermi energy). We obtained these scaling functions by extrapolating results for the Hubbard model to ν → 0. For the critical temperature of the superfluid-normal transition in the uniform case, we found T c /ε F = 0.152 (7). Our results form a basis for an unbiased thermometry of trapped fermionic gases in the unitary regime: in particular, we found (within the LDA) that for the parabolic confinement, the critical temperature in units of the characteristic trap energy E F is T/E F = 0.20 (2). For the experimentally relevant case of an isentropic conversion of a gas from the non-interacting regime to the unitary regime, we find that the onset of the superfluidity corresponds to the initial temperature (before the magnetic field ramp) (T/T F ) 0 = 0.12 ± 0.015, which is reasonably consistent with the experimental result [6], to within the experimental noise.

Acknowledgments
We appreciate the generosity of A Bulgac, P Magierski, and J Drut, who kindly provided us with their numeric data. We are also indebted to W Li and R Hulet for providing us with their unpublished experimental data. This research was enabled by computational resources of the Center for Computational Sciences and in part supported by the Laboratory Research and Development program at Oak Ridge National Laboratory. Part of the simulations were performed on the 'Hreidar' cluster of ETH Zürich. We also acknowledge partial support by NSF grants numbers PHY-0426881 and PHY-0456261.

Appendix A. Updating procedures
The generic diagrammatic MC rules for constructing updates are as follows [45]. Let an update B transform a diagram D(S p ) into diagram D(S q ). Here, configurations S p and S q may have different order and may or may not feature the 'worm' two-point vertices. The update B involves two steps. First, a modification of the configuration is proposed, with some probability density W(S p → S q ) for new continuous/discrete variables. There are no strict requirements fixing the form of W(S p → S q ), it should rather be chosen on physical grounds to maximize the efficiency of the algorithm, and be simple enough to allow analytical evaluation of the normalization integral. Then, the update is either accepted, with probability P S p →S q , or rejected. The complementary update C, which transforms D(S q ) into D(S q ) proposes the modification with the probability density W(S p ← S q ), which, in principle, may differ from W(S p → S q ), and is accepted with the probability P S p ←S q . For a pair of updates B and C to be balanced, the acceptance probabilities must obey the Metropolis relations [45] where R is a solution of the detailed-balance equation If the probabilities of selecting updates B and C are not equal, they must also be included into the definition of W(S p → S q ) and W(S p ← S q ) respectively. a boson is zero. Thus, multi-ladder diagrams still dominate, but the xand τ-span of individual ladders becomes finite and related to the inter-particle distance, see figure A.1(B). Resonant ladders also become more 'dilute' since the average distance between vertices diverges, but they still contain many vertices since the resonance formation involves distances much smaller than ν −1/3 .
Proposing diagrams which disregard the nature of the resonant interaction results in small acceptance ratios. To significantly reduce this type of slowing down, we resort to a wormtype updating scheme which allows one to, loosely speaking, directly 'draw' the multi-ladder diagrams.
A.2.1. Worm creation/annihilation updates. These updates are the necessary ingredient of the worm-type scheme, since they switch between the diagonal (3.3) and off-diagonal (3.9) diagrams. These updates transform diagrams as follows To create a pair of two-point vertices P(x, τ) and P † (x , τ ), one first selects x and τ uniformly over the spatial lattice and (0, β) interval, respectively. One then selects x among the lattice sites in the spatial cube with the edge length l centred around x with equal probabilities, and selects τ uniformly on the interval (τ − τ/2, τ + τ/2). The values of l and τ are arbitrary. The overall probability density is thus given by The reverse update attempts to remove P(x, τ) and P(x , τ ) provided |x α − x α | < l/2 (α = 1, 2, 3) and |τ − τ | < τ/2, as prescribed by the balance requirement. The solution to the detailed balance equation (A.3) then reads In order to maximize the efficiency of the numerical scheme, it is advantageous to have the acceptance probabilities (A.1) and (A.2) to be of the order unity. On the other hand, the explicit macroscopically large factors in the right-hand side of equation (A.9) render the acceptance probabilities macroscopic. The use of the extra weighting factor ζ is now clear: by choosing 1/ζ = βL d τl d /ζ, one removes undesired factors from the right-hand side of equation (A.9). The remaining freedom of choosingζ has to be employed to further fine-tune the acceptance probabilities. The required values ofζ are typically of the order of unity.
A.2.2. Adding/removing vertices within the worm framework. These updates are central for the method since they change the diagram order. The main idea of these updates is: given a pair of 'worm' two-point vertices P and P † , use, say, P † (x, τ) as a tip of a magic pen to add or remove vertices from the diagram (note that it suffices to use either P or P † as a dynamical variable, the choice being only a matter of taste) Low-density version. The strategy just described works well close to half-filling. In the lowdensity regime, the multi-ladder diagrams dominate and the above scheme is not optimal for a number of reasons: (i) it does not respect the diagram structure which should feature (at least locally) well-defined vertex chains; (ii) the probability density according to which the value of τ new is proposed should favour shifting P † towards smaller τ's. To see this, consider two particles in vacuum: the matrix elements of P † (x, τ 1 )P(x, τ 2 ) simply equal zero for τ 2 < τ 1 .]; (iii) both large and small values of δτ = τ − τ new have to be accounted for. Consider two fermions in vacuum again. The local probability density for a shift of length δτ and x new = x is proportional to the free diffusion propagator: w(δτ) ∝ δτ −3/2 , which is singular as δτ → 0. However, the mean value of δτ for this distribution, w(δτ)δτ d(δτ), diverges at the upper limit. These properties of resonant ladders are ignored when shifts are proposed uniformly over some ± τ interval. The requirements (ii) and (iii) are best met if, when adding a vertex, a new position of P † is proposed according to the probability density (cf equation (A.11)) Since the free single-particle propagator is known only numerically, we use equation (A.13) on a uniform mesh with some step σ: δτ j = σj. That is, we tabulate the values w j = w(y; δτ j ) prior to the start of computations; we then select y and δτ j according to the (normalized) distribution w(y; δτ j ). The proposed new position of P † is then chosen as x new = x + y, and τ new = τ + δτ j + , where is uniform over the interval [−σ/2, σ/2]. We find that σ ∼ 1/5U produces good results. In order to meet the requirement (i), the notion of the closest-in terms of a certain distance-neighbouring vertex is introduced: the reverse update S p+1 → S p deterministically deletes the closest neighbour of P † . The detailed balance than requires that the forward update where w j is given by equation (A.13).
The scheme (A.14) in the dilute regime typically yields an order of magnitude gain in efficiency, as compared to the scheme (A.12).

A.2.3. Shifting the worm two-point vertices.
Clearly, the scheme (A.14) should be supplemented by an update which allows for an easy change of the closest neighbour. The following update performs the task We select x among the nearest neighbours of the lattice site x with equal probabilities and τ uniformly in some interval around τ.
This update is self-balanced, and the acceptance ratio function R is simply (A. 16)