# Effective-range dependence of resonant Fermi gases

L.M. Schonenberg  Cavendish Laboratory, J.J. Thomson Avenue, Cambridge, CB3
0HE, United Kingdom  G.J. Conduit  Cavendish Laboratory, J.J. Thomson Avenue,
Cambridge, CB3 0HE, United Kingdom

###### Abstract

A Fermi gas of cold atoms allows precise control over the dimensionless
effective range, $k_{\mathrm{F}}R_{\mathrm{eff}}$, of the Feshbach resonance.
Our pseudopotential formalism allows us to create smooth potentials with
effective range, $-2\leq k_{\mathrm{F}}R_{\mathrm{eff}}\leq 2$, which we use
for a variational and diffusion Monte Carlo study of the ground state of a
unitary Fermi gas. We report values for the universal constants of
$\xi=0.388(1)$ and $\zeta=0.087(1)$, and compute the condensate fraction,
momentum distribution, and pair correlations functions. Finally, we show that
a gas with $k_{\mathrm{F}}R_{\mathrm{eff}}\gtrsim 1.9$ is thermodynamically
unstable.

##  I Introduction

Cold atom gases have delivered a series of surprises and insights, including
polaron physics Schirotzek  _et al._ (2009), the realization of the Bose-
Hubbard model Greiner  _et al._ (2002), and the BEC-BCS crossover Leggett
(1980); Bourdel  _et al._ (2004); Giorgini  _et al._ (2008). The development
of uniform trapping potentials has enabled the experimental realization of
particles in a box Gaunt  _et al._ (2013), while the Feshbach resonance offers
a unique level of control of the inter-particle interactions Chin  _et al._
(2010). Fermi gases interacting via zero-range contact interactions offer
scale invariant physics in the unitary limit of diverging scattering length,
captured by the Bertsch parameter Baker (1999). However, despite their
universal physics, contact interactions do not represent finite range
interactions seen in nature, e.g. screened Coulomb forces, neutron-neutron
interactions, and narrow Feshbach resonances. In this paper we present a study
of the consequences of finite ranged interactions in a unitary Fermi gas.

The scattering of two particles at low energies is described by the scattering
phase shift Landau and Lifshitz (1981), which up to first order in the wave
vector $k$ is given by

$\cot(\delta(k))=-\frac{1}{ka}+\frac{1}{2}kR_{\mathrm{eff}},$

where $a$ is the scattering length and $R_{\mathrm{eff}}$ the effective range.
In the limit of zero interaction range $R_{\mathrm{eff}}=0$, a vanishing
scattering length $a=0$ corresponds to a noninteracting gas, while the unitary
limit of infinite scattering length $a^{-1}=0$ results in scale-invariance. We
use this scale-invariance as a solid basis to investigate the effects of the
length scale introduced by the effective range term $R_{\mathrm{eff}}$.
Typical values for the effective range are
$k_{\mathrm{F}}R_{\mathrm{eff}}\approx 3$ Miller  _et al._ (1990); Forbes  _et
al._ (2012) for neutron matter, and $k_{\mathrm{F}}R_{\mathrm{eff}}\gtrsim-4$
for the $543.25\mathrm{G}$ narrow Feshbach resonance of 6Li Gurarie and
Radzihovsky (2007); Wang  _et al._ (2011). There is a wide variety of Feshbach
resonances available Chin  _et al._ (2010) and several of those exhibit large
negative effective ranges, summarized in Ref. Wang  _et al._ (2011).

So far most quantum Monte Carlo (QMC) studies of finite-range interactions
have used the Pöschl-Teller interaction potential for
$0<k_{\mathrm{F}}R_{\mathrm{eff}}<0.4$, and then extrapolate effective range
effects to zero to study the ground state of the unitary Fermi gas Li  _et
al._ (2011); Gandolfi  _et al._ (2011); Forbes  _et al._ (2011). Forbes  _et
al._ (2012) purposefully consider the effect of small positive effective
ranges up to $k_{\mathrm{F}}R_{\mathrm{eff}}=0.35$ in the context of neutron
matter. Negative effective ranges have been studied in an Eagles-Leggett mean-
field theory using a well-barrier interaction potential at zero temperature De
Palo  _et al._ (2004); Jensen  _et al._ (2006), at finite-temperature De Palo
_et al._ (2005), and also using the two-channel model of the Feshbach resonant
interaction at both zero and finite temperatures Timmermans  _et al._ (1999);
Gurarie and Radzihovsky (2007).

Here we study a gas at unitarity across a broad spread of effective ranges
$-2\leq k_{\mathrm{F}}R_{\mathrm{eff}}\leq 2$. Many-body physics arises from
repeated two-body scattering events, so a Hamiltonian where the opposite spin
fermions interact via a pseudopotential that exactly reproduces the scattering
phase shift with $a^{-1}=0$ and $R_{\mathrm{eff}}$ is the ideal platform for
an accurate many-body simulation. To smoothly connect positive and negative
effective ranges we develop a new pseudopotential following Refs. Bugnion  _et
al._ (2014); Whitehead  _et al._ (2016a). The pseudopotential is smooth and
extended in space, making it easy to sample with the variational and diffusion
Monte Carlo methods that we use to calculate ground state properties Foulkes
_et al._ (2001).

In Section II we use two-body scattering theory to understand the properties
of our potential. In Section III we evaluate four possible choices for the
interaction potential, including the newly proposed Ultra-Transferable
Pseudopotential (UTP) and select the UTP as the potential of choice for
numerical studies. In Section IV we discuss the quantum Monte Carlo formalism
and present results for the ground state energy including values for the
universal constants, condensate fraction, momentum distribution and Tan’s
contact, and pair correlation functions. Finally, we consider the
thermodynamic stability of the system in Section V and find that gases with
$k_{\mathrm{F}}R_{\mathrm{eff}}\gtrsim 1.9$ are unstable.

##  II Formulation of the problem

We study the Hamiltonian for spin $1/2$ fermions in three dimensions with
resonant interactions between opposite spins,

$\hat{H}=-\frac{1}{2}\sum_{i=1}^{N}\nabla_{i}^{2}+\sum_{i\neq
j}^{N}V(r_{ij}).$

Atomic units ($\hbar=m=1$) are used throughout. $\nabla_{i}^{2}$ is the
Laplacian with respect to the coordinates of particle $i$, $N$ is the total
number of particles and we study equal numbers of up and down spin particles.
$r_{ij}$ is the distance between particles $i$ and $j$, and $V$ is an
interaction potential that acts between particles with opposite spins,
characterized by the idealized scattering phase shift
$\cot(\delta(k))=kR_{\mathrm{eff}}/2$. To understand the form of this
interaction potential we first summarize some important results from
scattering theory and in particular consider the possible emergence of bound
states, before we discuss the explicit forms of the potentials used.

###  II.1 Scattering theory

We consider two identical distinguishable fermions in a vacuum. In their
center-of-mass frame, the Schrödinger equation for particles interacting via a
radially symmetric potential $V(r)$ is given in spherical coordinates by

$[-\nabla^{2}+V(r)]\psi(r,\theta,\phi)=E\psi(r,\theta,\phi),$

where $E$ is the energy of the relative motion.

The analytic solution for noninteracting particles, $V(r)=0$, takes the form

$\psi_{lm}(r,\theta,\phi)=Y_{lm}(\theta,\phi)R_{l}(r),$

with $l$ the angular momentum and $m$ the component of the angular momentum
along the quantization axes. $Y_{lm}$ are the spherical harmonics, and the
radial function $R_{l}$ is given by

$R_{l}(r)=\mathcal{A}_{l}(k)j_{l}(kr)+\mathcal{B}_{l}(k)n_{l}(kr),$ (1)

where $k=\sqrt{E}$ is the wave vector in the center-of-mass frame, and the
coefficients $\mathcal{A}_{l}(k)$ and $\mathcal{B}_{l}(k)$ are set by the
boundary conditions. $j_{l}(kr)$ and $n_{l}(kr)$ are the spherical Bessel and
spherical Neumann functions respectively. To connect to scattering waves we
rewrite the radial function in terms of spherical Hankel functions
$h_{l}^{(1,2)}(kr)=j_{l}(kr)\pm\mathrm{i}n_{l}(kr)$,

$R_{l}(r)=\mathcal{A}_{l}^{\prime}(k)h_{l}^{(1)}(kr)+\mathcal{B}_{l}^{\prime}(k)h_{l}^{(2)}(kr).$

The Hankel functions $h_{l}^{(1,2)}(kr)$ behave as spherical waves at large
radii $\sim\exp[\pm\mathrm{i}(kr-l\pi/2)]/r$.

The effect of a spherically symmetric interaction potential $V(r)$ on the wave
function is limited by angular momentum conservation and causality to the
introduction of a phase shift $\delta_{l}(k)$ in the outgoing wave
$h_{l}^{(1)}(kr)$ of the radial wave function,

$R_{l}^{\mathrm{int}}(r)=\mathcal{C}_{l}(k)[\mathrm{e}^{\mathrm{i}2\delta_{l}(k)}h_{l}^{(1)}(kr)+h_{l}^{(2)}(kr)],$
(2)

with $\mathcal{C}_{l}(k)$ a normalization constant. At large radii
$R_{l}^{\mathrm{int}}(r)\sim\sin(kr+\delta_{l}(k)-l\pi/2)/r$, verifying the
interpretation of $\delta_{l}(k)$ as a phase shift. $\delta_{l}(k)$ is related
to the coefficients in Equation (1) as
$\delta_{l}(k)=\arctan[-\mathcal{B}_{l}(k)/\mathcal{A}_{l}(k)]$. The ratio
$-\mathcal{B}_{l}(k)/\mathcal{A}_{l}(k)$ can be expressed in terms of the
logarithmic derivative of the interacting radial wave function by matching
$R_{l}^{\mathrm{int}}(r)$ and $R_{l}(r)$ at the cutoff radius $r_{\mathrm{c}}$
beyond which the interaction potential vanishes. Combining both results, the
phase shift can be expressed as

$\delta_{l}(k)=\arctan\bigg{[}\frac{kj_{l}^{\prime}(kr_{\mathrm{c}})-\gamma_{l}j_{l}(kr_{\mathrm{c}})}{kn_{l}^{\prime}(kr_{\mathrm{c}})-\gamma_{l}n_{l}(kr_{\mathrm{c}})}\bigg{]},$

where
$\gamma_{l}=(R_{l}^{\mathrm{int}})^{\prime}(r_{\mathrm{c}})/R_{l}^{\mathrm{int}}(r_{\mathrm{c}})$.

At large radii, the interacting wave function can also be written as the sum
of an incoming plane wave and a spherical outgoing scattered wave

$\lim_{\mathbf{r}\to\infty}\psi^{\mathrm{int}}(\mathbf{r})=\mathrm{e}^{\mathrm{i}\mathbf{k}\cdot\mathbf{r}}+\frac{f(k,\theta)}{r}\mathrm{e}^{\mathrm{i}kr},$

with $f$ the scattering amplitude and $\theta$ the scattering angle. By
equating the radial component of this expression in angular momentum channel
$l$ with Equation (2), the scattering amplitude can be related to the phase
shift

$f_{l}(k)=\frac{1}{\cot(\delta_{l}(k))-\mathrm{i}k}.$

This expression reveals bound states of the interaction potential because they
introduce poles into the scattering amplitude $f_{l}(k)$ Landau and Lifshitz
(1981). From now on we focus on the $l=0$ channel that dominates interactions
between opposite spin fermions, starting with an examination of possible bound
states in the next section.

###  II.2 Bound states

Figure 1: (Color online) Scattering phase shift
$\delta(k)=\operatorname{arccot}(kR_{\mathrm{eff}}/2)$ for the three different
cases of $R_{\mathrm{eff}}$. For $R_{\mathrm{eff}}<0$ the phase shift of a
realistic potential with the same low energy scattering properties is
indicated by the red dotted line. The noninteracting phase shift is shown by
the gray dashed line.

Each time the phase shift accumulates a factor of $\pi$, a node is introduced
in the wave function of the scattered wave $\sim\sin(kr+\delta(k))/r$. Since
each node introduced into the wave function by the potential corresponds to an
additional bound state, this establishes the link between the scattering phase
shift and the number of bound states $n\geq 0$ for any well-behaved potential,
which is formalized in Levinson’s theorem Chadan and Sabatier (1989),

$\displaystyle\delta(0)-\delta(\infty)=\begin{cases}n\pi,&a^{-1}\neq 0,\\\
(n+\frac{1}{2})\pi,&a^{-1}=0.\end{cases}$

We are interested in the latter case $a^{-1}=0$. As is evident from Fig. 1,
for $R_{\mathrm{eff}}>0$ the phase shift decreases from $\pi/2$ to $0$, and
there is no bound state. For $R_{\mathrm{eff}}<0$,
$\delta(0)-\delta(\infty)=-\pi/2$, which gives $n=-1$. Because the number of
bound states cannot be negative, this phase shift does not correspond to a
physical potential. However, potentials with the same low-energy scattering
properties may be obtained from a phase shift with additional contributions at
higher order in $k$. Provided these contributions occur at momenta beyond the
largest momentum scale in the system, i.e. the Fermi momentum $k_{\mathrm{F}}$
for a fermionic many-body system, they do not affect the physics of the system
as the interacting particles cannot probe these high momentum features. As
seen in the figure, the effect of the higher order term is to introduce a
phase winding of $\pi$ so that $\delta(0)-\delta(\infty)=\pi/2$, which
corresponds to a physical potential with no bound state. We conclude that in
both cases there is no bound state and the potential is therefore completely
characterized in terms of its scattering phase shift Chadan and Sabatier
(1989).

Despite the absence of a true bound state with negative energy, virtual bound
states may exist. The scattering amplitude for our idealized phase shift reads

$f=\frac{1}{\frac{1}{2}k^{2}R_{\mathrm{eff}}-\mathrm{i}k},$

which has a pole at zero energy $k=0$. In the zero-range limit,
$R_{\mathrm{eff}}=0$, this pole corresponds to a zero energy virtual bound
state, which is the $a\to\infty$ limit of the familiar bound state with energy
$E=-1/(2a^{2})$ Giorgini  _et al._ (2008); Gurarie and Radzihovsky (2007).
Because the pole in the scattering amplitude extends to finite effective range
$R_{\mathrm{eff}}$, so does the virtual bound state, which will be important
for our discussion of the many-body system in Section IV.1.

##  III Pseudopotentials

Having defined the interaction potential in terms of scattering properties, we
evaluate four possible real space interaction potentials for use in our many-
body simulations. For positive effective range we consider the potential well
and Pöschl-Teller interactions; for negative effective range we consider the
well-barrier potential. Furthermore, we propose the Ultra-Transferable
Pseudopotential (UTP) Bugnion  _et al._ (2014); Lloyd-Williams  _et al._
(2015); Whitehead and Conduit (2016); Whitehead  _et al._ (2016a), which is
equally applicable for both positive and negative effective ranges. After
comparing all four potentials, we select the UTP for our numerical study. The
software used to generate the UTP is available online Schonenberg and Conduit
(2017).

###  III.1 Positive effective range

Positive effective ranges $R_{\mathrm{eff}}>0$ for attractive interactions are
usually obtained from uniformly attractive potentials, $V(r)\leq 0$ for all
$r$. In this case, the effective range is approximately equal to the physical
interaction range Landau and Lifshitz (1981), while the depth of the potential
can be used to tune the scattering length.

####  III.1.1 Potential well

A spherical potential well interaction was used in Refs. Astrakharchik  _et
al._ (2004, 2005); Jensen  _et al._ (2006) as a model for contact
interactions,

$\displaystyle
V(r)=\begin{cases}-\bigg{(}\dfrac{\pi}{2R_{\mathrm{eff}}}\bigg{)}^{2},&r\leq
R_{\mathrm{eff}},\\\ 0,&r>R_{\mathrm{eff}},\end{cases}$

tuned to have scattering length $a^{-1}=0$, and effective range
$R_{\mathrm{eff}}$. The scattering phase shift of this potential is correct at
low incident energies, but is incorrect at intermediate energies where higher
order terms start to contribute.

####  III.1.2 Pöschl-Teller

The Pöschl-Teller interaction gives the exact phase shift with scattering
length $a$ and effective range $R_{\mathrm{eff}}$ in the lowest angular
momentum channel Chadan and Sabatier (1989), and has been used in several
studies Morris  _et al._ (2010); Carlson  _et al._ (2011); Forbes  _et al._
(2012); Li  _et al._ (2011). At unitarity, $a^{-1}=0$, the potential can be
written in terms of its effective range $R_{\mathrm{eff}}$ as

$V(r)=-\frac{8R_{\mathrm{eff}}^{-2}}{\cosh^{2}(\frac{2r}{R_{\mathrm{eff}}})}.$

###  III.2 Negative effective range

Scattering phase shifts with negative effective range result from potentials
with an attractive well hosting a (virtual) bound state at short radii, and a
potential barrier at intermediate radii. Quantum tunneling through the
potential barrier couples the (virtual) bound state with the continuum of
scattering states at large radii. When a rising barrier suppresses quantum
tunneling, the (virtual) bound and scattering states become uncoupled.

These potentials are called Shape resonances and exhibit the same physics as
Feshbach resonances. In the Feshbach resonance model the (virtual) bound state
in the well is represented by the closed channel, and the tunneling through
the potential barrier is described by a hybridization term that mixes the
closed channel with the open channel that describes the continuum of
scattering states Wang  _et al._ (2011).

####  III.2.1 Well-barrier potential

Following Refs. De Palo  _et al._ (2004, 2005); Jensen  _et al._ (2006) we
consider a well-barrier potential,

$\displaystyle V(r)=\begin{cases}-U_{0},&r\leq R_{0},\\\ U_{1},&R_{0}<r\leq
R_{1},\\\ 0,&r>R_{1},\end{cases}$

with $U_{0},U_{1}>0$ and $R_{1}>R_{0}>0$. This potential reduces to the
potential well for $U_{1}=0$. A potential with scattering length $a$ and
effective range $R_{\mathrm{eff}}$ for given radii $\\{R_{0},R_{1}\\}$ can be
obtained by suitably tuning the well depth and barrier height
$\\{U_{0},U_{1}\\}$ as described in Ref. Jensen  _et al._ (2006). We discuss
our choice for $\\{R_{0},R_{1}\\}$ in Section III.4.

As discussed in Section II.2, the scattering phase shift of physical
potentials with negative effective range include a phase winding by $\pi$ at
some high momentum $k$. Dimensional analysis confirms that this momentum may
be pushed to arbitrarily high momentum where it does not affect the scattering
of low-energy particles by reducing $\\{R_{0},R_{1}\\}$, at the expense of
diverging $\\{U_{0},U_{1}\\}$.

###  III.3 UTP

We now propose a pseudopotential that describes both positive and negative
effective ranges. It is also smooth and extended in space, easing the
application of numerical methods. Following Bugnion  _et al._ (2014);
Whitehead  _et al._ (2016a) we propose a UTP that takes a polynomial form
within a cutoff radius $r_{\mathrm{c}}$,

$\displaystyle
V^{\mathrm{UTP\\!}}(r)\\!=\\!\begin{cases}\\!\left(\\!1\\!-\\!\frac{r}{r_{\mathrm{c}}}\\!\right)^{2}\\!\left[u_{1}\\!\left(1+\frac{2r}{r_{\mathrm{c}}}\right)\\!+\\!\displaystyle\sum_{i=2}^{N_{\mathrm{u}}}u_{i}\left(\\!\frac{r}{r_{\mathrm{c}}}\\!\right)^{\\!i}\right]\\!,&\\!\\!\\!\\!r\leq
r_{\mathrm{c}},\\\ 0,&\\!\\!\\!\\!r>r_{\mathrm{c}},\end{cases}$

where the $u_{i}$ are the $N_{\mathrm{u}}=5$ optimizable coefficients. The
term $(1-r/r_{\mathrm{c}})^{2}$ ensures that the pseudopotential goes smoothly
to zero at $r=r_{\mathrm{c}}$, and the component $u_{1}(1+2r/r_{\mathrm{c}})$
constrains the pseudopotential to have zero gradient at particle coalescence,
to ensure that the wave function in the potential is smooth.

We calibrate the potential to deliver the correct scattering phase shift for
particles with momenta up to the characteristic momentum scale of our many-
body system, the Fermi momentum $k_{\mathrm{F}}$. To determine the
coefficients $\\{u_{i}\\}$ we numerically solve the scattering problem,
extract the scattering phase shift $\delta_{l}^{\mathrm{UTP}}(k)$, and then
minimize the total squared error in the phase shift over angular momentum
channels $l$ and relative scattering wave vectors $k$ of particles in the
Fermi sea,

$\displaystyle\langle\left|\delta_{l}^{\mathrm{UTP}}(k)-\delta_{l}(k)\right|^{2}\rangle$
$\displaystyle=$ $\displaystyle\sum_{l}\int_{0}^{k_{\mathrm{F}}}$
$\displaystyle\left|\delta_{l}^{\mathrm{UTP}}(k)-\delta_{l}(k)\right|^{2}g(k/k_{\mathrm{F}})\mathrm{d}k,$

where the weighting is given by the density of scattering wave vectors in the
center of mass frame $g(x)=12x^{2}(2-3x+x^{3})$. The virtue of a large cutoff
radius $r_{\mathrm{c}}$ is that it leads to potentials that are more extended
in space. On the other hand, $r_{\mathrm{c}}$ should be smaller than the
inter-particle spacing so that three-body scattering events are rare. We
therefore choose $r_{c}=1/k_{\mathrm{F}}$, except for large positive effective
ranges where we need a cutoff radius of the order of $R_{\mathrm{eff}}$, so we
adopt $r_{c}=\max(1/k_{\mathrm{F}},2R_{\mathrm{eff}})$.

###  III.4 Comparison of potentials

Figure 2: (Color online) Plots of different potentials normalized by the
reciprocal Fermi energy $E_{\mathrm{F}}$ as a function of radius for effective
ranges $k_{\mathrm{F}}R_{\mathrm{eff}}=\\{-1,0.1,1\\}$. For the case
$k_{\mathrm{F}}R_{\mathrm{eff}}=0.1$, the potentials are shown using solid
lines, whereas the red and blue dashed lines show the Pöschl-Teller and
potential well interactions with $k_{\mathrm{F}}R_{\mathrm{eff}}=0.2$, and the
shaded purple region shows the variation of the UTP for effective range $0\leq
k_{\mathrm{F}}R_{\mathrm{eff}}\leq 0.2$. The parameters of the well-barrier
potential potential were chosen so that its depth and height correspond to
those of the UTP.

Having introduced four possible interaction potentials, we now compare how
accurately they recover the correct scattering phase shift and their numerical
efficiency. The latter is a combination of two factors: numerical convergence
is aided by smooth potentials as they produce smooth wave functions, and also
by potentials with a wide spatial extent as they occupy a larger volume of
configuration space so are more rapidly sampled.

To visualize the smoothness and extent of the interaction potentials we plot
potentials with effective ranges
$k_{\mathrm{F}}R_{\mathrm{eff}}=\\{-1,0.1,1\\}$ in Fig. 2. For positive
effective range the UTP is similar to the Pöschl-Teller interaction. The
potential well is of similar spatial extent, but shows a discontinuity. The
diverging depth of the Pöschl-Teller and potential well interactions in the
zero-range limit, $k_{\mathrm{F}}R_{\mathrm{eff}}=0$, is illustrated by the
deepening of those potentials as the effective range decreases from
$k_{\mathrm{F}}R_{\mathrm{eff}}=0.2$, indicated by the dotted line, to
$k_{\mathrm{F}}R_{\mathrm{eff}}=0.1$, indicated by the solid line. In this
limit the interaction becomes momentum independent, which for the potential
well and Pöschl-Teller interactions implies that they also become short-
ranged. This is not the case for the UTP, as we calibrate the potential only
for wave vectors up to an intermediate momentum scale $k_{\mathrm{F}}$. The
shaded region shows the variation of the UTP with effective range $0\leq
k_{\mathrm{F}}R_{\mathrm{eff}}\leq 0.2$, demonstrating its shape remains
similar even in the zero-range limit. The numerical advantage of the UTP
becomes clear: it remains smooth and extended in space. For negative effective
range the UTP displays a barrier at intermediate distances, like the well-
barrier potential. For the well-barrier potential we set
$\\{k_{\mathrm{F}}R_{0}=0.2,k_{\mathrm{F}}R_{1}=0.4\\}$, so that its depth and
height are similar to that of the UTP. Many-body simulations will benefit from
the smoothness of the UTP compared to the two discontinuities for the well-
barrier. We conclude that the UTP is the only potential that is of finite
depth at all effective ranges, is smooth and extended in space, and is
therefore well-suited for use in a QMC simulation.

Figure 3: (Color online) Root Mean Square (RMS) phase shift error for the
potential well, well-barrier and UTP for $k$ ranging from $0$ to
$k_{\mathrm{F}}$. On average, the RMS phase shift error for the UTP is about
two orders of magnitude less than that for the potential well and well-barrier
potentials.

Having examined the numerical advantages of the UTP compared to the other
potentials, we now evaluate the accuracy of their scattering phase shifts. In
Fig. 3 we plot the root mean square (RMS) phase shift error of the potentials
summed over angular momentum channels and integrated over scattering wave
vectors up to the Fermi wave vector $k_{\mathrm{F}}$. For positive effective
range the UTP is over two orders of magnitude more accurate than the potential
well. Its maximum error of less than $10^{-3}$ also renders it equivalent to
the exact Pöschl-Teller interaction for all practical purposes. For negative
effective ranges we note that even though the depth and height of the well-
barrier potential have been chosen to mimic the UTP, its phase shift
properties are almost two orders of magnitude worse. As the UTP is the easiest
to work with in a many-body simulation, accurate, and applies at all effective
ranges, we select it for our QMC many-body study.

##  IV Quantum Monte Carlo

To calculate the ground state properties of the Fermi gas we use a quantum
Monte Carlo (QMC) method that is a tandem of the variational Monte Carlo (VMC)
and fixed-node diffusion Monte Carlo (DMC) techniques Ceperley and Alder
(1980); Umrigar  _et al._ (1993); Foulkes  _et al._ (2001). We use the casino
implementation Needs  _et al._ (2010) with a Slater–Jastrow trial wave
function $\Psi=\mathrm{e}^{J}D$, where $D$ is a Slater determinant of $N/2$
pairing orbitals $\phi(\mathbf{r}_{ij})$, each holding an up and down spin
particle, and $J$ a Jastrow factor that we optimize first using VMC, before
using DMC to further relax the wave function to its ground state. DMC is an
accurate Green’s function projector method for determining ground state
energies and other expectation values, and is well-suited to investigating
homogeneous gaseous phases.

The pairing orbitals Carlson  _et al._ (2003); Morris  _et al._ (2010) are
formed of a linear combination of plane-waves and polynomials

$\displaystyle\phi(\mathbf{r}_{ij})=$
$\displaystyle\sum_{n=0}^{N_{\mathrm{PW}}}a_{n}\sum_{\mathbf{G}\in
S_{n}}\exp(\mathrm{i}\mathbf{G}\cdot\mathbf{r}_{ij})$
$\displaystyle+\Theta(L_{\mathrm{P}}-r_{ij})\bigg{(}1-\frac{r_{ij}}{L_{\mathrm{P}}}\bigg{)}^{3}\sum_{n=0}^{N_{\mathrm{P}}}b_{n}r_{ij}^{n},$

with $\mathbf{G}$ an element of the set of symmetry related reciprocal lattice
vectors $S_{n}$, $N_{\mathrm{PW}}$ is the number of sets to include, and
$N_{\mathrm{P}}$ is the order of the polynomial.
$\mathbf{r}_{ij}=\mathbf{r}_{i}-\mathbf{r}_{j}$ is the separation between two
particles with opposite spins at positions $\mathbf{r}_{i}$ and
$\mathbf{r}_{j}$, and $r_{ij}=|\mathbf{r}_{ij}|$ its magnitude. The term
$\Theta(L_{\mathrm{P}}-r_{ij})(1-r_{ij}/L_{\mathrm{P}})^{3}$, where $\Theta$
is the Heaviside step function and $L_{\mathrm{P}}$ is an optimizable cutoff
length, ensures that the polynomial orbital smoothly approaches zero before
the edge of the cell. The coefficients $a_{n}$ and $b_{n}$ are optimizable,
with the exception of $a_{0}$, which we set to 1, and $b_{1}$, which is fixed
by requiring that the orbital is cuspless at the origin. The Slater
determinant of these orbitals contains both the noninteracting limit where the
particles fill the $N_{\mathrm{PW}}$ shells of plane-waves, and a
superconducting state of Cooper pairs captured by the polynomial series.

As superconductivity is a collective phenomenon, it is important to to capture
many-body correlations in the pairing orbitals. We therefore use a backflow
transformation López Ríos  _et al._ (2006) that replaces the particle
coordinates $\mathbf{r}_{i}$ by collective coordinates
$\mathbf{x}_{i}(\mathbf{R})=\mathbf{r}_{i}+\mathbf{\zeta_{i}}(\mathbf{R})$
with

$\mathbf{\zeta_{i}}(\mathbf{R})=\sum_{i\neq
j}\mathbf{r}_{ij}\Theta(L_{\mathrm{B}}-r_{ij})\bigg{(}1-\frac{r_{ij}}{L_{\mathrm{B}}}\bigg{)}^{3}\sum_{n=0}^{\mathrm{N_{B}}}\eta_{s_{i}s_{j},n}r_{ij}^{n},$

where $s_{i},s_{j}$ are magnetic quantum numbers of particles $i$ and $j$,
$N_{\mathrm{B}}=5$ the order of the polynomial, and $L_{\mathrm{B}}$ is an
optimizable cutoff length. The optimizable coefficients $\eta_{\alpha\beta,i}$
have to obey the symmetry requirements
$\eta_{\uparrow\uparrow,i}=\eta_{\downarrow\downarrow,i}$ and
$\eta_{\uparrow\downarrow,i}=\eta_{\downarrow\uparrow,i}$. We find backflow
corrections between particles of equal spin to be insignificant and therefore
set $\eta_{\uparrow\uparrow,i}=\eta_{\downarrow\downarrow,i}=0$.

The Slater determinant is multiplied by a Jastrow factor $\mathrm{e}^{J}$, to
capture the short-distance behavior of the pairwise interaction potential. We
use

$\displaystyle J=\sum_{i\neq
j}\Theta(L_{\mathrm{J}}-r_{ij})\bigg{(}1-\frac{r_{ij}}{L_{\mathrm{J}}}\bigg{)}^{3}\sum_{n=0}^{N_{\mathrm{J}}}u_{s_{i}s_{j},n}r_{ij}^{n},$

where $N_{\mathrm{J}}=8$ is the order of the polynomial and $L_{\mathrm{J}}$
is an optimizable cutoff length that we choose in accordance with the cutoff
radius of the pseudopotential Drummond  _et al._ (2004). The optimizable
coefficients $u_{\alpha\beta,i}$ have to obey the symmetry requirements
$u_{\uparrow\uparrow,i}=u_{\downarrow\downarrow,i}$ and
$u_{\uparrow\downarrow,i}=u_{\downarrow\uparrow,i}$, and $u_{\alpha\beta,1}$
is fixed by requiring zero gradient at the origin. Similar results are
obtained with a Jastrow factor optimized for periodic systems Whitehead  _et
al._ (2016b).

In the zero-range limit $k_{\mathrm{F}}R_{\mathrm{eff}}=0$ the Slater-Jastrow
trial wave function captures 93% of the correlation energy, defined as the
difference in ground state energy between the Hartree-Fock and DMC results.
The backflow transformation captures another 3.5%, raising the total to 96.5%.
Backflow transformations are especially important for negative effective
range, where the amount of correlation energy captured without backflow
transformations is only 85% at $k_{\mathrm{F}}R_{\mathrm{eff}}=-2$, while a
trial wave function with backflow transformations captures 92% of the
correlation energy.

Observables other than the ground state energy are computed using the
extrapolated estimator
$\langle\hat{A}\rangle=2\langle\hat{A}\rangle_{\mathrm{DMC}}-\langle\hat{A}\rangle_{\mathrm{VMC}}$,
which combines the DMC and VMC expectation values of the operator $\hat{A}$ to
reduce the bias from linear to quadratic in the difference between the VMC and
DMC wave functions Ceperley and Kalos (1986). In agreement with Refs. Morris
_et al._ (2010); Gandolfi  _et al._ (2011); Li  _et al._ (2011) we find the
results of this extrapolation to be within the statistical error bar of the
DMC estimate and therefore expect residual errors to be small. We extrapolate
to zero DMC time step and infinite number of walkers to obtain accurate ground
state energies following the procedure detailed in Section A.1. We expect that
the use of a quadratic DMC algorithm would give similar results Mella  _et
al._ (2000); Sarsa  _et al._ (2002). We calculate the ground state wave
function in the thermodynamic limit using datapoints for systems with {66,
114, 162, 186, 294} particles; technical details of the extrapolation to
infinite system size are provided in Section A.2. For the smallest system, we
use $N_{\mathrm{PW}}=5$ plane-waves to accommodate the $2\times 33$ spin-up
and down particles, while for the largest system we use $N_{\mathrm{PW}}=10$.
We set $N_{\mathrm{P}}=6$, allowing us to accurately describe particles in the
virtual bound state for negative effective range. In total our trial wave
function includes 34-39 parameters that we optimize using VMC, before using
the trial wave function as a starting point for our DMC calculations.

###  IV.1 Ground state energy

Figure 4: (Color online) Ground state energy per particle of the unitary Fermi
gas as a fraction of that of a noninteracting gas with effective range
$k_{\mathrm{F}}R_{\mathrm{eff}}$. Results obtained using the UTP are shown in
purple, and results obtained using the Pöschl-Teller interaction for positive
effective range in red. For large negative effective range the mean-field (MF)
theory of Ref. Gurarie and Radzihovsky (2007) is indicated by the dotted green
line.

Having developed a pseudopotential that smoothly connects positive and
negative effective ranges and outlined our trial wave function, we are well
positioned to study the ground state properties of resonant Fermi gases with
effective ranges $-2\leq k_{\mathrm{F}}R_{\mathrm{eff}}\leq 2$. We first study
the ground state energy per particle $E$, plotted as a fraction of that of a
noninteracting gas $E_{0}=\frac{3}{5}E_{\mathrm{F}}$, with $E_{\mathrm{F}}$
the Fermi energy, in Fig. 4. Starting from the zero-range case
$k_{\mathrm{F}}R_{\mathrm{eff}}=0$, we discuss the large negative and positive
effective range limits. The zero-range case itself will be discussed in the
next section.

We observe a decreasing energy as the effective range tends to
$k_{\mathrm{F}}R_{\mathrm{eff}}=-2$. The potential barrier we saw in Fig. 2
rises and decouples the virtual bound sate inside the barrier from the Fermi
sea. This reduces the energy of a pair of opposite spin particles in the bound
state towards the zero energy of the bare virtual bound state, causing more
particles at the Fermi surface to pair and the energy to approach zero. This
behavior is qualitatively the same as that from the BCS mean-field calculation
of Ref. Gurarie and Radzihovsky (2007), while quantitatively our DMC energy
approaches their mean-field energy, which is exact only in the limit
$k_{\mathrm{F}}R_{\mathrm{eff}}\to-\infty$.

For positive effective ranges, we observe a maximum value in the ground state
energy $E=0.432(1)\,E_{0}$ at $k_{\mathrm{F}}R_{\mathrm{eff}}=0.8$. For larger
effective ranges, the physical range of the interaction increases so that one
particle can now interact simultaneously with several opposite spin particles,
causing the energy to fall. The rapid decrease in energy in the
$k_{\mathrm{F}}R_{\mathrm{eff}}\to\infty$ limit gives rises to a thermodynamic
instability that will be discussed in Section V.

We also calculate the ground state energy using the alternative Pöschl-Teller
interaction available for positive effective ranges. The results for the UTP
and Pöschl-Teller interaction coincide, demonstrating the universality of the
many-body ground state energy for potentials with equivalent scattering
properties in the Fermi sea.

###  IV.2 Zero-range limit

Figure 5: (Color online) Ground state energy per particle of the unitary Fermi
gas as a fraction of that of a noninteracting gas with effective range
$k_{\mathrm{F}}R_{\mathrm{eff}}$. The UTP shown in purple and Pöschl-Teller in
red are results from this work. Also shown are other numerical results from
Refs. Astrakharchik  _et al._ (2004); Morris  _et al._ (2010); Carlson  _et
al._ (2011); Forbes  _et al._ (2012); Li  _et al._ (2011); Pessoa  _et al._
(2015), where for Ref. Carlson  _et al._ (2011) we extrapolated their results
to infinite system size as shown by the solid gray line, with the uncertainty
in our extrapolation indicated by the dashed lines. Experimental results at
$k_{\mathrm{F}}R_{\mathrm{eff}}=0$ from Refs. Luo and Thomas (2009); Navon
_et al._ (2010); Zürn  _et al._ (2013); Ku  _et al._ (2012) are shown in the
box on the left, slightly offset from unitarity for improved readability.

Having studied the variation of the ground state energy over the full extent
of effective ranges we now focus on the behavior near the zero range limit
$k_{\mathrm{F}}R_{\mathrm{eff}}=0$ in Fig. 5. For small effective range the
ground state energy per particle can be parameterized as Werner and Castin
(2012)

$\frac{E}{\frac{3}{5}E_{\mathrm{F}}}=\xi+\zeta
k_{\mathrm{F}}R_{\mathrm{eff}}+O((k_{\mathrm{F}}R_{\mathrm{eff}})^{2}),$

where the Bertsch parameter $\xi$ and $\zeta$ are universal constants for
Galilean invariant continuous space models. From Ref. von Keyserlingk and
Conduit (2013) we note that the effective coupling is stronger for more
negative effective range and we therefore expect $\zeta>0$. We report
$\xi=0.388(1)$, which agrees with the experimental result of Luo and Thomas
(2009). Our result is two times the experimental standard error lower than the
result of Navon  _et al._ (2010), while it is approximately two standard
errors higher than the results of Ku  _et al._ (2012) and Zürn  _et al._
(2013). Our statistical error estimates are negligible in comparison with the
experiments, but the fixed node constraint on the variational wave function
introduces a systematic error, which could explain why our value is higher
than the experimental measurements of Refs. Ku  _et al._ (2012); Zürn  _et
al._ (2013). Our result agrees with $\xi=0.390(2)$ from a DMC calculation by
Pessoa et al. Pessoa  _et al._ (2015). For the slope we find $\zeta=0.087(1)$,
in agreement with the auxiliary field result $\zeta=0.11(3)$ from Carlson  _et
al._ (2011), but in disagreement with the DMC result $\zeta=0.127(4)$ of
Forbes  _et al._ (2012). Before discussing how this deviating value may be
understood as the result of the computational method employed, we first
illuminate how the choice of pseudopotential influences the results.

To compare the UTP with the Pöschl-Teller interaction used in Refs. Morris
_et al._ (2010); Forbes  _et al._ (2012); Li  _et al._ (2011) we calculate the
ground state energy for both potentials using the same trial wave function. As
we have seen before the equivalent phase shift of the UTP and Pöschl-Teller
interaction guarantees the same ground state energy for both potentials, but
as the effective range is reduced the DMC energy calculated using the Pöschl-
Teller interaction overestimates that of the UTP and the error bars for the
Pöschl-Teller interaction become larger. As QMC is a variational method, it is
important to use an accurate trial wave function. This is especially true for
attractive interactions where the BCS instability requires that the nodal
surface is optimized during a VMC calculation, before fixing the nodes and
further reducing the energy using DMC Carlson  _et al._ (2003). The quality of
the trial wave function is described by the variance of the local energy
$E_{\mathrm{L}}=\Psi^{-1}\hat{H}\Psi$, which is zero for the true ground
state. Because the depth of the of the Pöschl-Teller interaction diverges, the
local energy variance for $k_{\mathrm{F}}R_{\mathrm{eff}}\leq 0.2$ calculated
with the Pöschl-Teller interaction is more than four times that calculated
with the UTP, explaining the overestimate of the ground state energy and
larger error bars, and confirming the numerical advantage of our wide and
smooth UTP.

Having understood how the smooth UTP results in a lower variational estimate
of the energy, we now compare our results with other DMC studies Astrakharchik
_et al._ (2004); Morris  _et al._ (2010); Forbes  _et al._ (2012); Li  _et
al._ (2011); Pessoa  _et al._ (2015). We find a lower variational energy
because, besides the smooth UTP, our study employs a trial wave function that
includes more variational freedom over previous studies. In particular, we
have combined the flexible Jastrow factor, polynomial pairing orbitals and
backflow transformation of Ref. Morris  _et al._ (2010) with the plane-wave
orbitals used in Refs. Carlson  _et al._ (2003); Forbes  _et al._ (2012). This
could explain why our reported value for the slope $\zeta=0.087(1)$ is lower
than $\zeta=0.127(4)$ from Forbes  _et al._ (2012).

We have also compared our DMC with the auxiliary field QMC study that is free
from the sign problem for a spin balanced system with attractive interactions
and therefore does not require the fixed node approximation Carlson  _et al._
(2011). The results in their Fig. 2 display finite size effects, leading to
uncertainty in our extrapolation of their results to infinite system size.
Nevertheless, the extrapolated ground state energy for effective range
$0.15<k_{\mathrm{F}}R_{\mathrm{eff}}<0.3$ agrees with our result within
$0.01\,E_{0}$.

###  IV.3 Condensate fraction

Figure 6: (Color online) Condensate fraction as a function of effective range
calculated using the UTP in purple and Pöschl-Teller interaction in red for
comparison. We also plot the results obtained near the zero-range limit by
Refs. Astrakharchik  _et al._ (2005); Morris  _et al._ (2010); Li  _et al._
(2011).

Having studied the variation of the ground state energy we now examine other
expectation values starting with the condensate fraction. A defining feature
of a superconductor is the existence of a condensate that introduces
correlations between Cooper pairs of opposite spin particles irrespective of
their separation. Correlations between pairs of opposite spins are naturally
captured by the off-diagonal two-body density matrix

$\rho_{\downarrow\uparrow}^{(2)}(\mathbf{r}_{1}^{\prime},\mathbf{r}_{2}^{\prime};\mathbf{r}_{1},\mathbf{r}_{2})=\langle
c_{\uparrow}^{\dagger}(\mathbf{r}_{1}^{\prime})c_{\downarrow}^{\dagger}(\mathbf{r}_{2}^{\prime})c_{\downarrow}(\mathbf{r}_{2})c_{\uparrow}(\mathbf{r}_{1})\rangle,$

where $c_{\alpha}^{\dagger}(\mathbf{r})$ is the fermionic creation and
$c_{\alpha}(\mathbf{r})$ the annihilation operator for a particle with spin
$\alpha$ at position $\mathbf{r}$. It is convenient to work in coordinates
that make the separation between two pairs,
$\mathbf{R}=\frac{1}{2}(\mathbf{r}_{1}^{\prime}+\mathbf{r}_{2}^{\prime})-\frac{1}{2}(\mathbf{r}_{1}+\mathbf{r}_{2})$,
and the size of the pairs $\mathbf{r}=\mathbf{r}_{1}-\mathbf{r}_{2}$, and
$\mathbf{r}^{\prime}=\mathbf{r}^{\prime}_{1}-\mathbf{r}^{\prime}_{2}$
explicit. In the limit $R=|\mathbf{R}|\to\infty$ the two-body density matrix
is proportional to the condensate fraction $c$ Leggett (2006); Astrakharchik
_et al._ (2005)

$\rho_{\downarrow\uparrow}^{(2)}\bigg{(}\mathbf{R}+\frac{\mathbf{r}^{\prime}}{2},\mathbf{R}-\frac{\mathbf{r}^{\prime}}{2};\frac{\mathbf{r}}{2},-\frac{\mathbf{r}}{2}\bigg{)}\to
c\frac{N}{2}\phi^{*}(|\mathbf{r}^{\prime}|)\phi(|\mathbf{r}|),$ (3)

where $\phi(r)$ is the complex pair wave function normalized to reciprocal
volume $1/\Omega$, and $N$ is the number of particles. In the normal phase
where correlations between pairs vanish as $R\to\infty$, $c=0$, while for a
superconducting phase the pairs remain correlated however far apart they are
and the fraction of particles in the condensate is $0<c\leq 1$.

The numerical computation of the condensate fraction using the relation above
is complicated as it requires an extrapolation to the $R\to\infty$ limit
Astrakharchik  _et al._ (2005); Morris  _et al._ (2010); Li  _et al._ (2011).
This motivates us to use a Fourier transform to capture the long-distance
behavior as a zero-momentum mode, in order to accurately compute the
condensate fraction using all accumulated samples of the two-body density
matrix across the entire simulation cell, following the procedure outlined in
Appendix B.

The condensate fractions calculated with the UTP and Pöschl-Teller
interactions agree as seen in Fig. 6, confirming the accuracy of the UTP. We
also show data from references Astrakharchik  _et al._ (2005); Morris  _et
al._ (2010); Li  _et al._ (2011) for comparison and, after taking into account
the effective ranges used, observe good agreement between results. In the
zero-range limit $k_{\mathrm{F}}R_{\mathrm{eff}}=0$, we report $c=0.56(2)$ and
the negative slope for the condensate fraction is consistent with the positive
slope for the energy encountered earlier, as the breaking of Cooper pairs
increases the energy. There is a maximum in the condensate fraction at
$k_{\mathrm{F}}R_{\mathrm{eff}}\approx-1.8$ of $0.83(1)$. For more negative
effective range, the virtual bound states decouple from the Fermi sea, and so,
although the particles remain paired, they no longer interact with each other
and correlations between pairs vanish, causing the condensate fraction to
decrease Gurarie and Radzihovsky (2007).

###  IV.4 Momentum distribution

Figure 7: (Color online) (Top) Momentum distribution $n(k)$ for
$k_{\mathrm{F}}R_{\mathrm{eff}}=\\{-1,0,1\\}$. The noninteracting distribution
is indicated by the gray dashed line. (Bottom) Tail of the momentum
distribution on logarithmic axes. The lines indicate a weighted least squares
fit of $n(k)=C/k^{\alpha}$ to the momentum tail $k>2.2k_{\mathrm{F}}$.

Having surveyed how the ground state energy and condensate fraction vary with
effective range, we now select three characteristic effective ranges
$k_{\mathrm{F}}R_{\mathrm{eff}}=\\{-1,0,1\\}$ to study one- and two-body
correlation functions. In this section we study the momentum distribution
shown in Fig. 7. In the limit $R_{\mathrm{eff}}\to\infty$ the physical range
of the potential diverges, approaching a constant background potential and so
the momentum distribution approaches that of a noninteracting system. When the
effective range is decreased from positive to negative, the sharp cutoff at
the Fermi momentum disappears as weight is moved from low momenta to the high
momentum tail characteristic of a state of paired particles.

The tail of the momentum distribution at unitarity in the zero-range limit is
$n(k)\to C/k^{4}$, where $C$ is Tan’s contact Tan (2008a); *Tan2008a;
*Tan2008b. As shown by Ref. Braaten  _et al._ (2008) this result extends to
$R_{\mathrm{eff}}<0$ and the contact becomes a function of effective range
$C(k_{\mathrm{F}}R_{\mathrm{eff}})$. In the zero-range limit we report
$C(0)/k_{\mathrm{F}}^{4}=0.119(1)$, which is in reasonable agreement with
0.1147(3) Gandolfi  _et al._ (2011) and 0.0961(1) Pessoa  _et al._ (2015)
computed using different trial wave functions. Our results for
$k_{\mathrm{F}}R_{\mathrm{eff}}=-1$ are consistent with a $\sim 1/k^{4}$ tail,
and we find an increased value for the contact
$C(-1)/k_{\mathrm{F}}^{4}=0.157(3)$ as expected when more particles pair. For
positive effective ranges we observe a more rapidly decaying tail.

###  IV.5 Pair-correlation function

Figure 8: (Color online) Pair-correlation function for opposite (top) and
equal (bottom) spins for $k_{\mathrm{F}}R_{\mathrm{eff}}=\\{-1,0,1\\}$. The
noninteracting correlation function is indicated by the gray dashed line.

To better understand the two-body interactions that cause the deformation of
the Fermi surface, we show the pair-correlation function for opposite and
equal spins in Fig. 8. For opposite spins, we correct the pair-correlation
function for short-range effects due to the particular form of the
pseudopotential Lloyd-Williams  _et al._ (2015)

$g_{\uparrow\downarrow}(r)=\frac{g_{\uparrow\downarrow}^{\mathrm{2-body,exact}}(r)}{g_{\uparrow\downarrow}^{\mathrm{2-body,pseudo}}(r)}g_{\uparrow\downarrow}(r),$

where $g_{\uparrow\downarrow}^{\mathrm{2-body,\\{exact,pseudo\\}}}(r)$ are the
pair-correlation functions for the two-body problem computed using the exact
and pseudopotential wave functions respectively. Since our UTP is norm-
conserving Bugnion  _et al._ (2014); Whitehead  _et al._ (2016a) no correction
is necessary outside of the interaction region.

For opposite spins, the pair-correlation function naturally shows that due to
the attractive interaction the particles are more likely to be found in close
proximity compared to the noninteracting case. In the zero-range limit the
pair correlation diverges at short inter-particle distances as $\sim 1/r^{2}$
Tan (2008a); *Tan2008a; *Tan2008b. This divergence becomes stronger for
negative effective range where particles of opposite spins are more likely to
be found in pairs, whereas for positive effective range the particles are
further apart compared to the zero-range case. The dominant contribution to
the correlations between equal spins at positive and zero effective range is
the exchange-correlation hole due to the Pauli exclusion principle. The volume
of the exchange-correlation hole diminishes as the effective range becomes
negative, because the fermions are more likely to be paired in the virtual
bound state and behave as composite bosons.

##  V Thermodynamic stability

We saw in Fig. 4 that the ground state energy of a Fermi gas at unitarity
falls rapidly with increasing effective range. This raises the possibility of
a thermodynamic instability towards phase separation into a high density phase
with a high value of $k_{\mathrm{F}}R_{\mathrm{eff}}$, and so a large negative
energy, and a low density phase. To analyze this possibility we first assess
the behavior of the ground state energy using a mean-field approximation
before investigating the thermodynamic instability.

###  V.1 Hartree-Fock theory

The proposed collapse into the dense phase means the dimensionless physical
interaction range $k_{\mathrm{F}}R_{\mathrm{eff}}$ diverges and the
interaction potential approaches a constant background potential. In this
limit the wave function approaches that of a noninteracting system so we can
use the Hartree-Fock approximation to estimate the ground state energy per
particle as

$\displaystyle E_{\mathrm{HF}}$
$\displaystyle=\frac{3E_{\mathrm{F}}}{5}+\frac{1}{4}\int\mathrm{d}^{3}r\,V(\mathbf{r})n(\mathbf{r})$
$\displaystyle=\frac{3E_{\mathrm{F}}}{5}\bigg{(}1-\frac{10\pi}{108}k_{\mathrm{F}}R_{\mathrm{eff}}\bigg{)},$
(4)

with $n(\mathbf{r})$ the density, which in our case is uniform so $n(r)=n$,
and the factor of $1/4$ accounts for the fact that interactions act only
between particles of opposite spin. The dependence on the explicit form of the
interaction potential only enters via the integral, and the result is
independent of the choice for the potential-well, Pöschl-Teller, or UTP.

###  V.2 Stability

To assess the thermodynamic stability we consider the Helmholtz free energy
density $F/\Omega$, with $\Omega$ the volume. The Helmholtz free energy is
$F=E-TS$, with temperature $T=0$ and $S$ the entropy. For a thermodynamically
(meta)stable phase the free energy density is required to be a convex function
of density, so $d^{2}F/dn^{2}>0$, whereas if $d^{2}F/dn^{2}<0$ the system
phase separates Landau and Lifshitz (2001).

Figure 9: (Color online) (Top) Energy density as a function of the Fermi gas
density on the bottom axes and Fermi wave vector on the top axes. HF is the
Hartree-Fock theory Equation (V.1) shown by the green dashed line. We also
show our DMC results obtained using the UTP in purple, fitted with cubic
splines. (Bottom) Curvature for the Hartree-Fock theory and UTP. The spinodal
point estimated from UTP data at $nR_{\mathrm{eff}}^{3}\approx 0.25$
($k_{\mathrm{F}}R_{\mathrm{eff}}\approx 1.9$) is near the spinodal point
calculated using the Hartree-Fock theory at $nR_{\mathrm{eff}}^{3}=0.235$
($k_{\mathrm{F}}R_{\mathrm{eff}}=1.91$).

In Fig. 9 we show the free energy density derived from Equation (V.1) as a
function of the Fermi gas density, as well as the curvature derived from this
free energy density. A spinodal point where the free energy density turns from
convex to concave occurs at $nR_{\mathrm{eff}}^{3}=0.235$. This signals the
onset of phase separation into an infinitely dense phase and a low density
phase for effective range $nR_{\mathrm{eff}}^{3}>0.235$, that is
$k_{\mathrm{F}}R_{\mathrm{eff}}>1.91$.

In the same figure we also show our DMC data obtained using the UTP. Our DMC
simulations are insensitive to phase separation as the trial wave function has
insufficient overlap with that of a phase separated state, so we can address
the full extent of effective ranges. We interpolate our DMC results with cubic
splines to determine $d^{2}F/dn^{2}$, which show that the spinodal point is at
$nR_{\mathrm{eff}}^{3}\approx 0.25$, $k_{\mathrm{F}}R_{\mathrm{eff}}\approx
1.9$, consistent with the Hartree-Fock result.

We establish that resonant gases with effective range
$k_{\mathrm{F}}R_{\mathrm{eff}}\gtrsim 1.9$ are thermodynamically unstable to
particle collapse. The instability of a Fermi gas with finite range
interactions at unitarity is reminiscent of the instability for neutron matter
at finite scattering lengths, where a three-body repulsive force is necessary
to ensure thermodynamic stability Friedman and Pandharipande (1981). The
concerns of Ref. Forbes  _et al._ (2012) that the Pöschl-Teller interaction
might harbor a many-body bound state is a precursor to the instability
presented here.

##  VI Conclusion

We have proposed the UTP as an interaction potential to model resonant
scattering of fermions with varying effective interaction range. Unlike other
potentials, the UTP smoothly connects the positive and negative effective
range regimes. Moreover, at the midpoint between those regimes where the
effective range is zero, the UTP remains smooth, extended in space, and of
finite depth. This allows us to perform an accurate calculation of the ground
state properties as we can directly simulate the zero-range limit, with no
need for extrapolations.

Exploiting the numerical advantages of the UTP and an improved estimator for
the condensate fraction, we have performed DMC ground state calculations for
resonant gases as a function of effective interaction range. In the zero-range
limit, we report values for the universal constants of $\xi=0.388(1)$ and
$\zeta=0.087(1)$, for the contact $C/k_{\mathrm{F}}^{4}=0.119(1)$, and for the
condensate fraction $c=0.56(2)$. Furthermore, by studying the momentum
distribution and pair correlation functions, we have demonstrated how the
system evolves from a state of independent pairs of opposite spin particles
for negative effective range, to the strongly interacting state in the zero-
range limit $k_{\mathrm{F}}R_{\mathrm{eff}}=0$, and finally to the weakly
interacting BCS superconductor for positive effective range. We find resonant
gases with effective range $k_{\mathrm{F}}R_{\mathrm{eff}}\gtrsim 1.9$ are
unstable to phase separation into an infinitely dense phase and a vacuum phase
containing no particles.

Having covered the complete gamut of effective interaction ranges, we expect
our results will be relevant for cold atom gases with both broad and narrow
Feshbach resonances. On the positive effective range side we also expect our
results to be relevant for neutron matter. Furthermore, the UTP formalism,
extended here to include the effective range term, will be useful for future
studies of both contact and finite ranged interactions.

###### Acknowledgements.

The authors thank Pablo López Ríos, Thomas Whitehead, Neil Drummond, Richard
Needs, Stefano Giorgini, and Jordi Boronat for useful discussions. The authors
acknowledge the financial support of the EPSRC [EP/J017639/1]; LMS
acknowledges financial support from the Cambridge European Trust, VSB Fonds,
and the Prins Bernhard Cultuurfonds; and GJC acknowledges the financial
support of the Royal Society and Gonville & Caius College. Computational
facilities were provided by the University of Cambridge High Performance
Computing Service. There is Open Access to this paper and data Schonenberg and
Conduit (2017).

##  Appendix A Details of QMC extrapolations

In this section we provide technical details of the extrapolations employed to
acquire accurate QMC data. To accurately extract the ground state energy it is
important to extrapolate to zero time step and infinite walker population,
discussed in the next section, and to the thermodynamic limit, discussed in
the second section.

###  A.1 Time step and walker population extrapolation

In DMC the imaginary time evolution operator $\mathrm{e}^{-\hat{H}\Delta\tau}$
is applied at each time step $\Delta\tau$ to a trial wave function,
represented by a finite number of walkers, to project out the ground state
Foulkes  _et al._ (2001). The gradient of the energy as a function of time
step is expected to be proportional to the local energy variance Lloyd-
Williams  _et al._ (2015); Whitehead and Conduit (2016) and the true ground
state is recovered by extrapolating to zero time step and infinite walker
population Lee  _et al._ (2011).

In practice we perform these extrapolations simultaneously by reducing the
time step by a factor of two, while increasing the walker population by the
same factor, as in Fig. 10. For optimal efficiency the computational effort
should be increased by $2\sqrt{2}$ for each division of the time step by 2 Lee
_et al._ (2011). For effective range $R_{\mathrm{eff}}\geq 0$ our trial wave
function optimized using VMC results in small local energy variances. The
variation with energy and walker population is therefore small, less than
$10^{-3}E_{0}$ in the linear regime, which for
$k_{\mathrm{F}}R_{\mathrm{eff}}=1$ extends to time steps $\Delta\tau
E_{\mathrm{F}}\leq 10^{-2}$ and for $k_{\mathrm{F}}R_{\mathrm{eff}}=0$ to time
steps $\Delta\tau E_{\mathrm{F}}\leq 2.5\times 10^{-3}$. However, for negative
effective range $R_{\mathrm{eff}}<0$ the local energy variance is larger and
we observe a significant variation with time step and walker in the linear
regime extending up to time steps $\Delta\tau E_{\mathrm{F}}\leq 2.5\times
10^{-3}$. We have performed additional tests to show that the variation of the
energy originates from the reduction in time step and not from the increase in
walker population. Extrapolating to zero time step is essential as even the
smallest time step used $\Delta\tau E_{\mathrm{F}}=0.625\times 10^{-3}$
introduces a systematic error to the ground state energy of $3\times
10^{-3}E_{0}$, and smaller time steps would require a large number of steps to
exceed the auto-correlation time of the random Monte Carlo walk. We conclude
that extrapolating to zero time step and infinite walker population is
essential to accurately calculate ground state energies.

Figure 10: (Color online) Ground state energy for
$k_{\mathrm{F}}R_{\mathrm{eff}}=\\{-1,0,1\\}$ as a function of time step
(bottom axes) and walker population (top axes). The straight lines show a
weighted least squares fit to the data points in the linear regime.

###  A.2 System size extrapolation

Having extrapolated to zero time step and infinite walker population, we now
extrapolate to infinite system size. As the length scale associated with the
features of our interaction potential is less than the average inter-particle
separation, we expect system size effects to be dominated by the kinetic
energy term in the Hamiltonian. The finite-size error in the kinetic term
originates from the discretization of the plane-wave wave vectors, which in
three dimensions is proportional to the reciprocal number of particles $1/N$
Drummond  _et al._ (2008).

Exploiting our smooth pseudopotential to create low-variance trial wave
functions, we can study systems with up to 294 particles. Results are shown in
Fig. 11, where we observe the expected linear regime for systems with more
than 162 particles. For $R_{\mathrm{eff}}\leq 0$, finite size effects are
$<2\times 10^{-3}E_{0}$ for systems with more than 162 particles, where the
particles are bound in pairs described by the polynomial term in the pairing
orbitals. In contrast, for $k_{\mathrm{F}}R_{\mathrm{eff}}=1$ the plane-wave
term dominates and the trial wave function is closer to that of a
noninteracting system, thus displaying larger finite-size effects with
variations in energy up to $5\times 10^{-3}E_{0}$ as the system size is
decreased from an infinite number of particles to 162 particles. The residual
errors in the extrapolation are $<10^{-3}E_{0}$ and we conclude that
extrapolating to infinite system size using systems with at least 162
particles is essential to obtain accurate predictions for the ground state
energy.

Figure 11: (Color online) Ground state energy for
$k_{\mathrm{F}}R_{\mathrm{eff}}=\\{-1,0,1\\}$ as a function of the number of
particles. The straight lines show a weighted least squares fit to the data
points with $N\geq 162$.

##  Appendix B Evaluating the condensate fraction

A central property of a superconductor is the existence of a condensate of
pairs of particles. The condensate manifests itself as a macroscopic
eigenvalue of the two-body density matrix for opposite spins irrespective of
the distance between the pairs of opposite spin particles, i.e. off-diagonal
long range order Yang (1962). In practice the limiting behavior of the two-
body density matrix is often used to compute the condensate fraction
Astrakharchik  _et al._ (2005); Morris  _et al._ (2010); Li  _et al._ (2011),
thereby ignoring available knowledge of the two-body density matrix at short
distances and in the corners of the simulation cell. Here, we propose a
Fourier transform to exploit knowledge of a modified two-body density matrix
over the entire simulation cell to accurately estimate the condensate
fraction. We show that this improved estimator gives direct access to the
condensate fraction.

We consider the BCS wave function Annett (2004)

$|\Psi_{\mathrm{BCS}}\rangle=\prod_{\mathbf{k}}(u_{\mathbf{k}}^{*}+v_{\mathbf{k}}^{*}c^{\dagger}_{\uparrow\mathbf{k}}c^{\dagger}_{\downarrow-\mathbf{k}})|0\rangle,$

with $u_{\mathbf{k}},v_{\mathbf{k}}$ the usual complex coherence factors to
evaluate the expectation values in this section. As demonstrated by the
Eagles-Leggett mean-field theory of the BEC-BCS crossover this wave function
is qualitatively correct even in the strong coupling limit Eagles (1969);
Leggett (1980). We introduce the order parameter $F(\mathbf{r})=\langle
c_{\downarrow}(\mathbf{r})c_{\uparrow}(\mathbf{r})\rangle$, related to the
pair wave function introduced earlier as $F(r)=\sqrt{cN/2}\phi(r)$. Expressed
in terms of the coherence factors $F(\mathbf{r})$ is

$F(\mathbf{r})=\frac{1}{\Omega}\sum_{\mathbf{k}}\mathrm{e}^{\mathrm{i}\mathbf{k}\cdot\mathbf{r}}u_{\mathbf{k}}v_{\mathbf{k}}^{*},$

with $\Omega$ the volume. $F(\mathbf{r})$ is the eigenfunction of the off-
diagonal two-body density matrix at large inter-pair separation $R$ and the
condensate fraction $c$ is defined in terms of its macroscopic eigenvalue
Leggett (2006)

$\displaystyle c$
$\displaystyle=\frac{2\Omega}{N}\int_{\Omega}\mathrm{d}^{3}r|F(\mathbf{r})|^{2}$
$\displaystyle=\frac{2}{N}\sum_{\mathbf{k}}|u_{\mathbf{k}}|^{2}|v_{\mathbf{k}}|^{2}.$

To compute the condensate fraction numerically we introduce the spatially
averaged one- and two-body density matrices. We use the projected two-body
density matrix obtained by setting the separation between particles in each
pair equal to each other, $\mathbf{r}=\mathbf{r}^{\prime}$, thereby
eliminating one coordinate to integrate over De Palo  _et al._ (2002); Ortiz
and Dukelsky (2005)

$\displaystyle\bar{\rho}^{(1)}_{\alpha}(\mathbf{R})=$
$\displaystyle\frac{1}{\Omega}\int_{\Omega}\mathrm{d}^{3}\bar{r}\rho^{(1)}_{\alpha}(\bar{\mathbf{r}}+\mathbf{R},\bar{\mathbf{r}}),$
$\displaystyle\bar{\rho}_{\downarrow\uparrow}^{(\mathrm{2})}(\mathbf{R})=$
$\displaystyle\frac{1}{\Omega^{2}}\int_{\Omega}\mathrm{d}^{3}\bar{r}\,\mathrm{d}^{3}r$
$\displaystyle\rho_{\downarrow\uparrow}^{(2)}\bigg{(}\bar{\mathbf{r}}+\mathbf{R}+\frac{\mathbf{r}}{2},\bar{\mathbf{r}}+\mathbf{R}-\frac{\mathbf{r}}{2};\bar{\mathbf{r}}+\frac{\mathbf{r}}{2},\bar{\mathbf{r}}-\frac{\mathbf{r}}{2}\bigg{)},$

where $\rho^{(1)}_{\alpha}(\mathbf{r}^{\prime},\mathbf{r})=\langle
c_{\alpha}^{\dagger}(\mathbf{r}^{\prime})c_{\alpha}(\mathbf{r})\rangle$ is the
one-body density matrix for spin $\alpha$.

To remove known short-ranged one-body contributions from the two-body density
matrix we follow Ref. Astrakharchik  _et al._ (2005) and introduce an
estimator for the condensate fraction

$c_{\uparrow\downarrow}(\mathbf{R})=\frac{2\Omega^{3/2}}{N}\big{(}\bar{\rho}^{(2)}_{\uparrow\downarrow}(\mathbf{R})-\bar{\rho}^{(1)}_{\uparrow}(\mathbf{R})\bar{\rho}^{(1)}_{\downarrow}(\mathbf{R})\big{)}.$

Using Equation (3) we find at large radius $R$

$\lim_{R\to\infty}c_{\downarrow\uparrow}(R)=c/\sqrt{\Omega}.$

The extrapolation to the large $R$ limit is problematic in numerical studies
where information is available for finite $R$ values only, and would neglect
information available in the simulation cell at smaller distances and further
out in the corners of the simulation cell. Instead, we propose a Fourier
transform to capture the long distance behavior of the two-body density matrix
as a discontinuity at small momentum. Defining the Fourier transform pair as

$\displaystyle f(\mathbf{r})$
$\displaystyle=\frac{1}{\sqrt{\Omega}}\sum_{\mathbf{k}}\mathrm{e}^{-\mathrm{i}\mathbf{k}\cdot\mathbf{r}}f_{\mathbf{k}},$
$\displaystyle f_{\mathbf{k}}$
$\displaystyle=\frac{1}{\sqrt{\Omega}}\int_{\Omega}\mathrm{d}^{3}r\,\mathrm{e}^{\mathrm{i}\mathbf{k}\cdot\mathbf{r}}f(\mathbf{r}),$

we compute the Fourier transform of the modified two-body density matrix to
give us direct access to the condensate fraction $c$

$\displaystyle c_{\uparrow\downarrow\mathbf{q}}$
$\displaystyle=\frac{2\Omega^{3/2}}{N}\big{(}\bar{\rho}^{(2)}_{\uparrow\downarrow\mathbf{q}}-\sum_{\mathbf{k}}\bar{\rho}^{(1)}_{\uparrow\mathbf{k}}\bar{\rho}^{(1)}_{\downarrow\mathbf{q}-\mathbf{k}}\big{)}$
$\displaystyle=\frac{2}{N}\sum_{\mathbf{k}}\big{[}\langle
c^{\dagger}_{\uparrow\mathbf{k}}c^{\dagger}_{\downarrow\mathbf{q}-\mathbf{k}}c_{\downarrow\mathbf{q}-\mathbf{k}}c_{\uparrow\mathbf{k}}\rangle$
$\displaystyle\hskip 50.0pt-\langle
c^{\dagger}_{\uparrow\mathbf{k}}c_{\uparrow\mathbf{k}}\rangle\langle
c^{\dagger}_{\downarrow\mathbf{q}-\mathbf{k}}c_{\downarrow\mathbf{q}-\mathbf{k}}\rangle\big{]}$
$\displaystyle=\frac{2\delta_{\mathbf{q}\mathbf{0}}}{N}\sum_{\mathbf{k}}|v_{\mathbf{k}}|^{2}|u_{\mathbf{k}}|^{2}$
$\displaystyle=\delta_{\mathbf{q}\mathbf{0}}c,$

with $\delta_{ab}$ the Kronecker delta function. The condensate fraction
exists as a discontinuous peak at zero momentum as expected, and due to the
subtraction of the one-body density matrix the condensate fraction estimator
contains no other contributions. We exploit this relation to accurately
compute the condensate fraction using all accumulated samples of the modified
two-body density matrix in the simulation cell.

## References

  * Schirotzek  _et al._ (2009) A. Schirotzek, C.-H. Wu, A. Sommer, and M. W. Zwierlein, Phys. Rev. Lett. 102, 230402 (2009).
  * Greiner  _et al._ (2002) M. Greiner, O. Mandel, T. Esslinger, T. W. Hänsch, and I. Bloch, Nature (London) 415, 39 (2002).
  * Leggett (1980) A. J. Leggett, in  _Modern Trends in the Theory of Condensed Matter_ , edited by A. Pekalski and J. Przystawa (Springer, Berlin, 1980) pp. 13–27.
  * Bourdel  _et al._ (2004) T. Bourdel, L. Khaykovich, J. Cubizolles, J. Zhang, F. Chevy, M. Teichmann, L. Tarruell, S. J. J. M. F. Kokkelmans, and C. Salomon, Phys. Rev. Lett. 93, 050401 (2004).
  * Giorgini  _et al._ (2008) S. Giorgini, L. P. Pitaevskii, and S. Stringari, Rev. Mod. Phys. 80, 1215 (2008).
  * Gaunt  _et al._ (2013) A. L. Gaunt, T. F. Schmidutz, I. Gotlibovych, R. P. Smith, and Z. Hadzibabic, Phys. Rev. Lett. 110, 200406 (2013).
  * Chin  _et al._ (2010) C. Chin, R. Grimm, P. Julienne, and E. Tiesinga, Rev. Mod. Phys. 82, 1225 (2010).
  * Baker (1999) G. A. Baker, Phys. Rev. C 60, 054311 (1999).
  * Landau and Lifshitz (1981) L. D. Landau and E. M. Lifshitz,  _Quantum Mechanics: Non-Relativistic Theory_ (Butterworth-Heinemann, Oxford, UK, 1981).
  * Miller  _et al._ (1990) G. Miller, B. Nefkens, and I. Šlaus, Phys. Rep. 194, 1 (1990).
  * Forbes  _et al._ (2012) M. M. Forbes, S. Gandolfi, and A. Gezerlis, Phys. Rev. A 86, 053603 (2012).
  * Gurarie and Radzihovsky (2007) V. Gurarie and L. Radzihovsky, Ann. Phys. (NY) 322, 2 (2007).
  * Wang  _et al._ (2011) Y. Wang, J. P. D’Incao, and B. D. Esry, Phys. Rev. A 83, 042710 (2011).
  * Li  _et al._ (2011) X. Li, J. Kolorenč, and L. Mitas, Phys. Rev. A 84, 023615 (2011).
  * Gandolfi  _et al._ (2011) S. Gandolfi, K. E. Schmidt, and J. Carlson, Phys. Rev. A 83, 041601 (2011).
  * Forbes  _et al._ (2011) M. M. Forbes, S. Gandolfi, and A. Gezerlis, Phys. Rev. Lett. 106, 235303 (2011).
  * De Palo  _et al._ (2004) S. De Palo, M. Chiofalo, M. Holland, and S. Kokkelmans, Phys. Lett. A 327, 490 (2004).
  * Jensen  _et al._ (2006) L. M. Jensen, H. M. Nilsen, and G. Watanabe, Phys. Rev. A 74, 043608 (2006).
  * De Palo  _et al._ (2005) S. De Palo, M. L. Chiofalo, M. J. Holland, and S. Kokkelmans, Laser Phys. 15, 376 (2005).
  * Timmermans  _et al._ (1999) E. Timmermans, T. Tomassini, M. Hussein, and A. Kerman, Phys. Rep. 315, 199 (1999).
  * Bugnion  _et al._ (2014) P. O. Bugnion, P. López Ríos, R. J. Needs, and G. J. Conduit, Phys. Rev. A 90, 033626 (2014).
  * Whitehead  _et al._ (2016a) T. M. Whitehead, L. M. Schonenberg, N. Kongsuwan, R. J. Needs, and G. J. Conduit, Phys. Rev. A 93, 042702 (2016a).
  * Foulkes  _et al._ (2001) W. M. C. Foulkes, L. Mitas, R. J. Needs, and G. Rajagopal, Rev. Mod. Phys. 73, 33 (2001).
  * Chadan and Sabatier (1989) K. Chadan and P. C. Sabatier,  _Inverse Problems in Quantum Scattering Theory_ (Springer, New York, 1989).
  * Lloyd-Williams  _et al._ (2015) J. H. Lloyd-Williams, R. J. Needs, and G. J. Conduit, Phys. Rev. B 92, 075106 (2015).
  * Whitehead and Conduit (2016) T. M. Whitehead and G. J. Conduit, Phys. Rev. A 93, 022706 (2016).
  * Schonenberg and Conduit (2017) L. M. Schonenberg and G. J. Conduit, Cambridge Univ. Apollo Repos. https://doi.org/10.17863/CAM.7090 (2017).
  * Astrakharchik  _et al._ (2004) G. E. Astrakharchik, J. Boronat, J. Casulleras, and S. Giorgini, Phys. Rev. Lett. 93, 200404 (2004).
  * Astrakharchik  _et al._ (2005) G. E. Astrakharchik, J. Boronat, J. Casulleras, and S. Giorgini, Phys. Rev. Lett. 95, 230405 (2005).
  * Morris  _et al._ (2010) A. J. Morris, P. López Ríos, and R. J. Needs, Phys. Rev. A 81, 033619 (2010).
  * Carlson  _et al._ (2011) J. Carlson, S. Gandolfi, K. E. Schmidt, and S. Zhang, Phys. Rev. A 84, 061602 (2011).
  * Ceperley and Alder (1980) D. M. Ceperley and B. J. Alder, Phys. Rev. Lett. 45, 566 (1980).
  * Umrigar  _et al._ (1993) C. J. Umrigar, M. P. Nightingale, and K. J. Runge, J. Chem. Phys. 99, 2865 (1993).
  * Needs  _et al._ (2010) R. J. Needs, M. D. Towler, N. D. Drummond, and P. López Ríos, J. Phys. Condens. Matter 22, 023201 (2010).
  * Carlson  _et al._ (2003) J. Carlson, S.-Y. Chang, V. R. Pandharipande, and K. E. Schmidt, Phys. Rev. Lett. 91, 050401 (2003).
  * López Ríos  _et al._ (2006) P. López Ríos, A. Ma, N. D. Drummond, M. D. Towler, and R. J. Needs, Phys. Rev. E 74, 066701 (2006).
  * Drummond  _et al._ (2004) N. D. Drummond, M. D. Towler, and R. J. Needs, Phys. Rev. B 70, 235119 (2004).
  * Whitehead  _et al._ (2016b) T. M. Whitehead, M. H. Michael, and G. J. Conduit, Phys. Rev. B 94, 035157 (2016b).
  * Ceperley and Kalos (1986) D. M. Ceperley and M. H. Kalos, in  _Monte Carlo Methods in Statistical Physics_, edited by K. Binder (Springer, Berlin; Heidelberg, 1986) 2nd ed., Chap. 4, pp. 145–194.
  * Mella  _et al._ (2000) M. Mella, G. Morosi, and D. Bressanini, Phys. Rev. E 61, 2050 (2000).
  * Sarsa  _et al._ (2002) A. Sarsa, J. Boronat, and J. Casulleras, J. Chem. Phys. 116, 5956 (2002).
  * Pessoa  _et al._ (2015) R. Pessoa, S. Gandolfi, S. A. Vitiello, and K. E. Schmidt, Phys. Rev. A 92, 063625 (2015).
  * Luo and Thomas (2009) L. Luo and J. E. Thomas, J. Low Temp. Phys. 154, 1 (2009).
  * Navon  _et al._ (2010) N. Navon, S. Nascimbene, F. Chevy, and C. Salomon, Science 328, 729 (2010).
  * Zürn  _et al._ (2013) G. Zürn, T. Lompe, A. N. Wenz, S. Jochim, P. S. Julienne, and J. M. Hutson, Phys. Rev. Lett. 110, 135301 (2013).
  * Ku  _et al._ (2012) M. J. H. Ku, A. T. Sommer, L. W. Cheuk, and M. W. Zwierlein, Science 335, 563 (2012).
  * Werner and Castin (2012) F. Werner and Y. Castin, Phys. Rev. A 86, 053633 (2012).
  * von Keyserlingk and Conduit (2013) C. W. von Keyserlingk and G. J. Conduit, Phys. Rev. B 87, 184424 (2013).
  * Leggett (2006) A. J. Leggett,  _Quantum Liquids: Bose Condensation and Cooper Pairing in Condensed-Matter Systems_ (Oxford University Press, Oxford, UK, 2006).
  * Tan (2008a) S. Tan, Ann. Phys. (NY) 323, 2971 (2008a).
  * Tan (2008b) S. Tan, Ann. Phys. (NY) 323, 2987 (2008b).
  * Tan (2008c) S. Tan, Ann. Phys. (NY) 323, 2952 (2008c).
  * Braaten  _et al._ (2008) E. Braaten, D. Kang, and L. Platter, Phys. Rev. A 78, 053606 (2008).
  * Landau and Lifshitz (2001) L. Landau and E. M. Lifshitz,  _Statistical Physics. Part 1_ (Butterworth-Heinemann, Oxford, UK, 2001).
  * Friedman and Pandharipande (1981) B. Friedman and V. Pandharipande, Nucl. Phys. A 361, 502 (1981).
  * Lee  _et al._ (2011) R. M. Lee, G. J. Conduit, N. Nemec, P. López Ríos, and N. D. Drummond, Phys. Rev. E 83, 066706 (2011).
  * Drummond  _et al._ (2008) N. D. Drummond, R. J. Needs, A. Sorouri, and W. M. C. Foulkes, Phys. Rev. B 78, 125106 (2008).
  * Yang (1962) C. N. Yang, Rev. Mod. Phys. 34, 694 (1962).
  * Annett (2004) J. F. Annett,  _Superconductivity, Superfluids and Condensates_ (Oxford University Press, Oxford, UK, 2004).
  * Eagles (1969) D. M. Eagles, Phys. Rev. 186, 456 (1969).
  * De Palo  _et al._ (2002) S. De Palo, F. Rapisarda, and G. Senatore, Phys. Rev. Lett. 88, 206401 (2002).
  * Ortiz and Dukelsky (2005) G. Ortiz and J. Dukelsky, Phys. Rev. A 72, 043611 (2005).