The imbalanced Fermi gas at unitarity
Abstract
Lattice field theory is a useful tool for studying strongly interacting theories in condensed matter physics. A prominent example is the unitary Fermi gas: a twocomponent system of fermions interacting with divergent scattering length. With Monte Carlo methods this system can be studied from first principles. In the presence of an imbalance (unequal number of particles in the two components) a sign problem arises, which makes conventional algorithms inapplicable. We will show how to apply reweighting techniques to generalise the recently developed worm algorithm to the imbalanced case, and present results for the critical temperature, the energy per particle, the chemical potential and the contact density for equal, as well as unequal number of fermions in the two spin components.
The imbalanced Fermi gas at unitarity
M. Wingate
DAMTP, University of Cambridge, Wilberforce Road, Cambridge CB3 0WA, UK
Email: M.Wingate@damtp.cam.ac.uk
\abstract@cs
1 Introduction
With the discovery of the renormalization group came the birth of both lattice gauge theory and statistical field theory as we know them today. Nowhere is this kinship more obvious than in the study of critical phenomena. For example, the numerical studies of deconfinement in gauge theories are identical in most respects to those of demagnetization in spin systems. Now a challenge in both fields is to develop more efficient methods for systems with strongly interacting fermions. Then we could better understand the phases of gauge theories with matter and of superfluids and superconductors.
Given that the common problem lies in dealing with fermionic degrees of freedom, we are motivated to consider the unitary Fermi gas, which is beautiful in its simplicity. This is a system of 2component nonrelativistic fermions (typically ultracold Li or K atoms in experiments) interacting via a shortrange potential. The gas is so dilute that only swave scattering is relevant, thus it is sufficient to treat the interaction as a local 4fermion coupling. On a spatial lattice, the simplest Hamiltonian is that of the attractive Hubbard model.
The lattice theory describing unitary Fermi gases avoids some of the complications of lattice gauge theory. With nonrelativistic fermions there is no chiral symmetry, so there is no fermion doubling problem. Without nonabelian gauge fields, there are no topological sectors to worry about sampling with the correct distribution. This is an ideal system for concentrating on Monte Carlo methods for fermions. In addition to both being strongly coupled theories, QCD and unitary Fermi gas both exhibit spontaneous symmetry breaking at low temperatures, restored at some critical temperature. Like finite density QCD, the fermion matrix of the unitary gas can have a nonpositive determinant.
2 Setup
We start from the FermiHubbard model, which in the grand canonical ensemble reads,
(2.0) 
where is the discrete dispersion relation, the chemical potential and () the timedependent fermionic creation (annihilation) operator. We use . The index labels the fermionic species. The coupling constant corresponding to attractive interaction can be tuned so that the scattering length becomes infinite. The corresponding value is . We work on a 3D periodic lattice with sites. The continuum limit can be taken by extrapolation to vanishing filling factor .
According to [1] the partition function for this model can be written as a series of products of two matrix determinants built of free finitetemperature Green’s functions. If (the balanced case) these determinants are equal so that all terms in the series are positive and it can be used a probability distribution for Monte Carlo sampling. The order parameter for the phase transition is given by a twopoint correlation function of the operator . To obtain from the numerical data, previous work [2, 3] used a procedure involving an approximation which introduced a systematic error. We have improved the data analysis method so that this approximation is no longer necessary [4]. Our final result will be the dimensionless quantity , where the Fermi energy is the only energy scale of the system. The continuum limit is taken by obtaining for different values of and extrapolating to vanishing filling factor [2, 5].
A detailed description of our numerical setup is given in [4] and [6]. We use the DDMC algorithm as introduced in [2] with a few modifications which increase the efficiency by reducing autocorrelation effects that are present in the original setup. Also we generalise the algorithm to the spinimbalanced case. If , a sign problem arises, since the distribution function is no longer positive for all configurations. To deal with this problem we make use of the “sign quenched method”, which is based on the “phase quenched method” known from lattice QCD [7]. The idea is to write the thermal distribution as a product of its modulus and its sign, and to use the positive function given by the modulus as the new probability distribution. This implies a reweighting of each MC estimator, in this case simply a multiplication with the relative sign of the two matrix determinants. Additionally, each expectation value must be divided by the expectation value of the sign. If the latter is is close to zero, numerical errors will be very large, as it happens in QCD. However, for the unitary Fermi gas the sign remains very close to unity for small imbalances, as shown in Fig. 1, so that sign quenching is applicable for imbalances up to approximately . Our method can provide a useful tool to examine the trend of the critical temperature for small deviations from the balanced limit.
3 Results
We obtained data at different values , of which were at . The lattice sizes varied between for the highest filling factor and for the lowest, so that the volume range in physical units was approximately constant. As discussed in [2] for the balanced case, the dimensionless physical observables scale linearly with for sufficiently small . With our data this behaviour is seen for . This condition was fulfilled for out of the points and in particular for out of the balanced points.
Since the chemical potential difference is less prone to numerical errors, we use instead of the relative density difference to quantify imbalance. For the values of imbalance considered in our study these two quantities are proportional to each other, with , see [4]. Every physical observable is a function of filling factor and imbalance . We fit our data to a three dimensional surface, where the following assumptions are made for the form of the fitted function:

At fixed imbalance, is a linear function of with slope : . This is a generalisation of the relation valid in the balanced case.

and viewed as functions of the imbalance can be Taylor expanded.

Due to symmetry in all odd powers in the expansions of and have to vanish.
If we expand and to leading order in the fitted function becomes
(3.0) 
In the following we will present our fit results for several physical quantities. The results from sections 3.1, 3.2 and 3.3 were presented in more detail in [4].
3.1 Critical Temperature
Fitting a line through the points with results in with d.o.f , as shown in Fig. 3. For comparison we also fit a quadratic through all data points, resulting in a continuum value of , which is in excellent agreement with the linear extrapolation. This confirms that subleading corrections proportional to can indeed be neglected for sufficiently small .
Now we also include data with . The best fit according to (3) yields , , and in units of , with d.o.f.. Note that the value corresponding to the minimal is positive, which is forbidden by physical arguments – the critical temperature can only decrease with increasing imbalance. Since the function is very flat along the direction, forcing results in d.o.f.. From the error on we derive the lower bound . The best fit values for and are in excellent agreement with the ones obtained from the fit of the balanced data only.
By simplifying the fit model setting , the lower bound on can be tightened to . The other parameters and agree with the results from the previous fit. This fit has d.o.f. and is still consistent with . For this reason we also perform a fit to constant and , and obtain with d.o.f.. Again the result agrees with the previous fits. We also performed fits using the jackknife method and several robust fits and obtained consistent results. A three dimensional plot of the data together with the constant surface fit is presented in Fig. 4 (left).
3.2 Energy per particle
The total energy is composed of the kinetic energy and the interaction energy . For the explicit MC estimators see [4]. Our results are obtained at , but the temperature dependence of the energy per particle was found to be weak. Using only balanced data we obtain the continuum value , or , in units of the ground state energy of the free gas . The goodness of fit is d.o.f. . Since with increasing imbalance interactions become suppressed, the absolute value of the interaction energy must decrease. This in turn means an increase of the total energy, since the interaction energy is negative. As we did for the critical temperature we fit the energy in units of to the function (3) and obtain the best fit parameters , , and , with d.o.f.. These results are consistent with the balanced fit. Forcing yields a best fit result of with d.o.f., which agrees with the previous result. For a plot of the data see Fig. 4 (right).
3.3 Chemical potential
For the chemical potential at we obtain the continuum value with d.o.f. using only balanced data. A similar analysis can be performed for the average chemical potential in presence of an imbalance. Since this quantity is not expected to depend on the imbalance we fit our data to a constant function and obtain in units of with d.o.f.. This is in very good agreement with our balanced result. A plot of the data and the fit are shown in Fig. 3.
3.4 Contact density
Another important quantity is the contact density, which can be interpreted as a measure of the local pair density [8]. The contact plays an important role for several universal relations derived by Tan [9]. The definition of the contact is , where is the physical coupling constant [8, 10], and it is related to the contact density via , or for homogeneous systems simply . The dimensionless quantity can be expressed as
(3.0) 
Using only balanced data the best fit is with d.o.f. . In the presence of an imbalance we expect the interaction energy and hence the contact density to decrease. A three dimensional fit of the data to (3) yields , , and , with d.o.f.. This is consistent with the balanced fit. The parameter is negative as required. Forcing yields , and , with d.o.f.. Fitting to a constant function yields with d.o.f.. Figure 5 summarises our results.
3.5 Comparison with the literature
Our final result for using both balanced and imbalanced data, , is significantly higher than the previous result from [2], where . In [11] the result of [2] was found to agree with a continuous spacetime DDMC method. The authors of [3] found an upper bound of . They used an auxiliary field Monte Carlo approach and extracted using the same approximation as [2] and [11], which might explain the discrepancy between our results. Through extrapolating Monte Carlo results of lowdensity neutron matter, the authors of [12] found , which agrees with our result. There are also results obtained with the Restricted Path Integral Monte Carlo method [13], , and an upper bound of obtained with a hybrid Monte Carlo method [14]. Results from an expansion are also available [15]. For comparison, the critical temperature in the BEC limit is .
Our result for the energy per particle shows excellent agreement with the value at quoted in [3]. The value quoted in [2] is , which roughly corresponds to . Our result for the chemical potential differs from quoted in [2], but is consistent with the value quoted in [3]. Some theoretical predictions for the contact density are also available, but to our knowledge only at temperatures much lower than (see [8] and references therein).
There are several recent experimental studies of the homogeneous unitary Fermi gas. The direct measurement presented in [16] is , which agrees well with our result, and at , which differs from our value. The values from [17], and at , show excellent agreement with our results. Their result for the energy per particle at is higher than our value. In another experimental work [18] an estimate for at zero imbalance is extrapolated from data at higher values of imbalance. An experimental value for the contact density of the homogeneous unitary Fermi gas at zero temperature, , is presented in [19].
4 Outlook and Acknowledgements
We have presented a Monte Carlo calculation of several thermodynamic observables of the unitary Fermi gas with equal and unequal chemical potentials in the two spin components. The improved DDMC algorithm with sign quenching also offers the intriguing possibility to explore the case of unequal masses of the two species. If the dispersion relations are different for the two components and the mass ratio enters as a new parameter. As in the spinimbalanced case the two matrix determinants no longer need to be identical, so that sign quenching is required. The mass ratio is expected to influence the behaviour of the system significantly, so that exploring the phase diagram promises many new interesting insights.
This work used resources provided by the Cambridge High Performance Computing Facility. OG is supported by the German Academic Exchange Service (DAAD), the EPSRC and the CET.
References
 [1] A. N. Rubtsov, V. V. Savkin, A. I. Lichtenstein, Phys. Rev. B 72, 035122 (2005).
 [2] E. Burovski, N. Prokof’ev, B. Svistunov, M. Troyer, New J. Phys. 8, 153 (2006).
 [3] A. Bulgac, J. E. Drut, P. Magierski, Phys. Rev. A 78, 023625 (2008).
 [4] O. Goulko, M. Wingate, [arXiv:1008.3348], to appear in Phys. Rev. A.
 [5] J.W. Chen, D. B. Kaplan, Phys. Rev. Lett. 92, 257002 (2004).
 [6] O. Goulko, M. Wingate, PoS(LAT2009)062, [arXiv:0910.3909].
 [7] D. Toussaint, Nucl. Phys. B Proc. Suppl. 17, 248 (1990).
 [8] E. Braaten, [arXiv:1008.2922].
 [9] S. Tan, Ann. Phys. 323, 2952 (2008); ibid 2971; ibid 2987.
 [10] F. Werner, Y. Castin, [arXiv:1001.0774].
 [11] E. Burovski, E. Kozik, N. Prokof’ev, B. Svistunov, M. Troyer, Phys. Rev. Lett. 101, 090402 (2008).
 [12] T. Abe, R. Seki, Phys. Rev. C 79, 054003 (2009).
 [13] V. K. Akkineni, D. M. Ceperley, N. Trivedi, Phys. Rev. B 76, 165116 (2007).
 [14] D. Lee, T. Schäfer, Phys. Rev. C 73, 015202 (2006).
 [15] Y. Nishida, Phys. Rev. A 75, 063618 (2007).
 [16] S. Nascimbène, N. Navon, K. J. Jiang, F. Chevy, C. Salomon, Nature 463, 1057 (2010).
 [17] M. Horikoshi, S. Nakajima, M. Ueda, T. Mukaiyama, Science 327, 442 (2010).
 [18] Y.I. Shin, C. H. Schunck, A. Schirotzek, W. Ketterle, Nature 451, 451 (2008).
 [19] N. Navon, S. Nascimbène, F. Chevy, C. Salomon, Science 328, 729 (2010).