Effective-range dependence of resonant Fermi gases

TL;DR

This study investigates how the effective range of interactions influences the properties of a resonant Fermi gas, providing new computational methods and identifying stability limits.

Contribution

We developed a pseudopotential formalism for smooth potentials with variable effective range and applied Monte Carlo methods to analyze the ground state of a unitary Fermi gas.

Findings

Universal constants: ξ = 0.388(1), ζ = 0.087(1)

Computed condensate fraction, momentum distribution, and pair correlations

Identified thermodynamic instability for k_F R_eff ≥ 1.9

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.

Equations

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

$$\hat{H}=-\frac{1}{2}\sum _{i=1}^N{\nabla }_i^2+\sum _{i\ne j}^{N}V(r_{ij}).$$

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

$${\psi }_{lm}(r,\theta ,\varphi )=Y_{lm}(\theta ,\varphi )R_l(r),$$

$$R_l(r)={\mathcal{A}}_l(k)j_l(kr)+{\mathcal{B}}_l(k)n_l(kr),$$

$$R_l(r)={\mathcal{A}}_l^{\prime }(k)h_l^{(1)}(kr)+{\mathcal{B}}_l^{\prime }(k)h_l^{(2)}(kr).$$

$$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)],$$

$$\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{]},$$

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

$$f_l(k)=\frac{1}{\cot ({\delta }_l(k))-\mathrm{i}k}.$$

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

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

$$\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}$$

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

$$V(r)=\{\begin{array}{ll}-U_0,& r\le R_0,\\ U_1,& R_0<r\le R_1,\\ 0,& r>R_1,\end{array}$$

$$V^{\mathrm{UTP}}(r)=\{\begin{array}{ll}{\left(1-\frac{r}{r_{\mathrm{c}}}\right)}^2\left[u_1\left(1+\frac{2r}{r_{\mathrm{c}}}\right)+\sum _{i=2}^{N_{\mathrm{u}}}{u}_i{\left(\frac{r}{r_{\mathrm{c}}}\right)}^i\right],& r\le r_{\mathrm{c}},\\ 0,& r>r_{\mathrm{c}},\end{array}$$

$$⟨{\left|{\delta }_l^{\mathrm{UTP}}(k)-{\delta }_l(k)\right|}^2⟩$$

$$\sum _l{\int }_0^{k_{\mathrm{F}}}$$

$$\varphi ({\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},$$

$$\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},$$

$$\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},$$

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

$${\rho }_{\downarrow \uparrow }^{(2)}({\mathbf{r}}_1^{\prime },{\mathbf{r}}_2^{\prime };{\mathbf{r}}_1,{\mathbf{r}}_2)=⟨c_{\uparrow }^{†}({\mathbf{r}}_1^{\prime })c_{\downarrow }^{†}({\mathbf{r}}_2^{\prime })c_{\downarrow }({\mathbf{r}}_2)c_{\uparrow }({\mathbf{r}}_1)⟩,$$

$$\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}|),$$

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

$$E_{\mathrm{HF}}$$

$$\displaystyle=\frac{3E_{\mathrm{F}}}{5}\bigg{(}1-\frac{10\pi}{108}k_{\mathrm{F}}R_{\mathrm{eff}}\bigg{)},$$

$$|{\Psi }_{\mathrm{BCS}}⟩=\prod _{\mathbf{k}}(u_{\mathbf{k}}^{\ast }+v_{\mathbf{k}}^{\ast }{c}_{\uparrow \mathbf{k}}^{†}{c}_{\downarrow -\mathbf{k}}^{†})|0⟩,$$

$$F(\mathbf{r})=\frac{1}{\Omega }\sum _{\mathbf{k}}{\mathrm{e}}^{\mathrm{i}\mathbf{k}\cdot \mathbf{r}}{u}_{\mathbf{k}}{v}_{\mathbf{k}}^{\ast },$$

$$c$$

$$=\frac{2}{N}\sum _{\mathbf{k}}|u_{\mathbf{k}}{|}^2|v_{\mathbf{k}}{|}^2.$$

Full text

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

[TABLE]

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,

[TABLE]

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

[TABLE]

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

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

[TABLE]

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

[TABLE]

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)$,

[TABLE]

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,

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

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),

[TABLE]

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 [math], 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

[TABLE]

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,

[TABLE]

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

[TABLE]

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,

[TABLE]

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}}$,

[TABLE]

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,

[TABLE]

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

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.

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

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

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)

[TABLE]

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

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

[TABLE]

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)

[TABLE]

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

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

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)

[TABLE]

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

[TABLE]

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).

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.

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.

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)

[TABLE]

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

[TABLE]

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)

[TABLE]

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)

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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.