Fermions in Two Dimensions: Scattering and Many-Body Properties

TL;DR

This paper investigates the properties of two-dimensional ultracold Fermi gases across different interaction regimes using ab initio Monte Carlo methods, focusing on scattering theory and many-body ground-state characteristics.

Contribution

It provides a detailed analysis of 2D scattering parameters and computes ground-state properties of strongly interacting Fermi gases with variational wave functions.

Findings

Calculated scattering length and effective range in 2D.

Determined ground-state energy and equation of state.

Analyzed chemical potential in strongly interacting regime.

Abstract

Ultracold atomic Fermi gases in two-dimensions (2D) are an increasingly popular topic of research. The interaction strength between spin-up and spin-down particles in two-component Fermi gases can be tuned in experiments, allowing for a strongly interacting regime where the gas properties are yet to be fully understood. We have probed this regime for 2D Fermi gases by performing T=0 ab initio diffusion Monte Carlo calculations. The many-body dynamics are largely dependent on the two-body interactions, therefore we start with an in-depth look at scattering theory in 2D. We show the partial-wave expansion and its relation to the scattering length and effective range. Then we discuss our numerical methods for determining these scattering parameters. We close out this discussion by illustrating the details of bound states in 2D. Transitioning to the many-body system, we use variationally optimized wave functions to calculate ground-state properties of the gas over a range of interaction strengths. We show results for the energy per particle and parametrize an equation of state. We then proceed to determine the chemical potential for the strongly interacting gas.

Figures

Tables

Table 1. Table 1: The final parameters for our EOS, where we fit to ( E TL / N − ϵ b / 2 ) / E FG subscript 𝐸 TL 𝑁 subscript italic-ϵ 𝑏 2 subscript 𝐸 FG (E_{\rm TL}/N-\epsilon_{b}/2)/E_{\rm FG} . A 7th-order polynomial f ​ ( η ) 𝑓 𝜂 f(\eta) is used in the crossover regime and we fit to f BEC ​ ( η ) superscript 𝑓 BEC 𝜂 f^{\rm BEC}(\eta) and f BCS ​ ( η ) superscript 𝑓 BCS 𝜂 f^{\rm BCS}(\eta) on either side. Their functional forms, Eq. ( 36 ) and Eq. ( 37 ), are based on the limiting behaviour on each side of the crossover and they each contain 3 parameters that are determined by continuity restrictions. In an effort to include as much of the intermediate polynomial as possible and minimize the overall variance of our fit, we select the matching points as η = − 0.25 𝜂 0.25 \eta=-0.25 and η = 2.5 𝜂 2.5 \eta=2.5 .

$c_i$	$f^{\mathrm{BEC}}\left(\eta \right)$	$f\left(\eta \right)$	$f^{\mathrm{BCS}}\left(\eta \right)$
$c_0$	28.545	0.18181
$c_1$	-42.648	0.13334
$c_2$	15.555	0.076788	-0.50515
$c_3$		-0.0099012	2.9215
$c_4$		-0.017582	-2.8151
$c_5$		0.0072822
$c_6$		-0.0010898
$c_7$		0.000057991

Equations

$$-\frac{{\hslash }^2}{2m_r}{\nabla }^2\psi (r,\theta )+V(r)\psi (r,\theta )=E\psi (r,\theta ),$$

$$\frac{1}{r}\frac{\partial }{\partial r}(r\frac{\partial \psi (r,\theta )}{\partial r})+\frac{1}{r^2}\frac{{\partial }^2\psi (r,\theta )}{\partial {\theta }^2}+k^2\psi (r,\theta )-\frac{2m_r}{{\hslash }^2}V(r)\psi (r,\theta )=0,$$

$$\psi (r,\theta )=\sum _{l=0}^{\infty }{a}_l{R}_l(r)T_l(\theta ).$$

$$-\bigg{[}k^{2}-\frac{2m_{r}}{\hbar^{2}}V(r)\bigg{]}=\frac{1}{r}\frac{\partial}{\partial r}\bigg{(}r\frac{\partial R_{l}(r)}{\partial r}\bigg{)}\frac{1}{R_{l}(r)}+\frac{1}{r^{2}}\bigg{(}\frac{1}{T_{l}(\theta)}\bigg{)}\frac{\partial^{2}T_{l}(\theta)}{\partial\theta^{2}}\,,$$

$$f(r)+g(\theta )=0,$$

$$\frac{{\partial }^2{T}_l(\theta )}{\partial {\theta }^2}=-c_l{T}_l(\theta ),$$

$$T_l(\theta )=a_l\sin (\sqrt{c_l}\theta )+b_l\cos (\sqrt{c_l}\theta ).$$

$$\frac{{\partial }^2{T}_l(\theta )}{\partial {\theta }^2}=-l^2{T}_l(\theta ),$$

$$T_l(\theta )=\frac{1}{\sqrt{\pi }}\cos (l\theta ).$$

$$\begin{split}-\,\bigg{[}k^{2}-\frac{2m_{r}}{\hbar^{2}}V(r)-\frac{l^{2}}{r^{2}}\bigg{]}&=\bigg{(}\frac{1}{u_{l}(r)\sqrt{r}}\bigg{)}\frac{\partial}{\partial r}\bigg{[}r~{}\frac{\partial}{\partial r}\Big{(}\frac{u_{l}(r)}{\sqrt{r}}\Big{)}\bigg{]}\,\\ &=\bigg{(}\frac{1}{u_{l}(r)\sqrt{r}}\bigg{)}\frac{\partial}{\partial r}\bigg{[}\sqrt{r}~{}\Big{(}\frac{\partial u_{l}(r)}{\partial r}-\frac{u_{l}(r)}{2r}\Big{)}\bigg{]}\,\\ &=\frac{1}{u_{l}(r)\sqrt{r}}\bigg{[}\frac{\partial^{2}u_{l}(r)}{\partial r^{2}}\sqrt{r}+\frac{u_{l}(r)}{4r^{3/2}}\bigg{]}\,\\ &=\frac{1}{u_{l}(r)}\frac{\partial^{2}u_{l}(r)}{\partial r^{2}}+\frac{1}{4r^{2}}\,.\\ \end{split}$$

$$-\,\frac{\partial^{2}u_{l}(r)}{\partial r^{2}}=u_{l}(r)\bigg{[}k^{2}-\frac{2m_{r}}{\hbar^{2}}V(r)-\frac{l^{2}-1/4}{r^{2}}\bigg{]}\,,$$

$$-\frac{\partial^{2}u_{0}(r)}{\partial r^{2}}=u_{0}(r)\bigg{[}k^{2}-\frac{2m_{r}}{\hbar^{2}}V(r)+\frac{1}{4r^{2}}\bigg{]}\,.$$

$$\frac{\partial^{2}u_{0}(r)}{\partial r^{2}}+u_{0}(r)\bigg{[}k^{2}+\frac{1}{4r^{2}}\bigg{]}=0\,.$$

$$u_{0}(r)\propto\sqrt{r}\big{[}J_{0}(kr)\cos\delta_{0}-N_{0}(kr)\sin\delta_{0}\big{]}\,.$$

$$\cot {\delta }_0=\frac{\tilde{K}{N}_0(k{r}_1)-N_0(k{r}_2)}{\tilde{K}{J}_0(k{r}_1)-J_0(k{r}_2)},$$

$$\frac{{\partial }^2{y}_0(r)}{\partial r^2}=-\frac{1}{4r^2}{y}_0(r),$$

$$y_0(r)=-\sqrt{r}\log (r/a_{2D}).$$

$$r_{e}^{2}=4\int_{0}^{\infty}\big{[}y_{0}^{2}(r)-u_{0}^{2}(r)\big{]}_{k\rightarrow 0}\,dr\,,$$

$$\cot\delta_{0}\big{|}_{k\rightarrow 0}\,\approx\frac{2}{\pi}\bigg{[}\gamma+\log\Big{(}\frac{ka_{\rm 2D}}{2}\Big{)}\bigg{]}+\frac{k^{2}r_{e}^{2}}{4}\,,$$

$$\begin{array}{r}{u}_0(r_1)=\xi \sqrt{r_1}\log (r_1/a_{2D})\mathrm{and}{u}_0(r_2)=\xi \sqrt{r_2}\log (r_2/a_{2D}),\end{array}$$

$$\displaystyle\begin{split}a_{\rm 2D}=r_{1}\exp\Bigg{[}\frac{-u_{0}(r_{1})}{\xi\sqrt{r_{1}}}\Bigg{]}\,,~{}~{}~{}{\rm where}~{}~{}\xi=\frac{\sqrt{r_{1}}\,u_{0}(r_{2})-\sqrt{r_{2}}\,u_{0}(r_{1})}{\sqrt{r_{1}r_{2}}\,\log(r_{2}/r_{1})}\,.\end{split}$$

$${ϵ}_b=-\frac{4{\hslash }^2}{m{a}_{2D}^2{e}^{2\gamma }},$$

$$V(r)=-v_0\frac{{\hslash }^2}{m_r}\frac{{\mu }^2}{{\cosh }^2(\mu r)},$$

$$k_F=\sqrt{2\pi n}.$$

$$\eta =\log (k_F{a}_{2D}),$$

$$\begin{array}{r}{\Phi }_{\mathrm{BCS}}(\mathbf{R})=\mathcal{A}[\varphi ({\mathbf{r}}_{11^{\prime }})\varphi ({\mathbf{r}}_{22^{\prime }})...\varphi ({\mathbf{r}}_{N_{\uparrow }{N}_{\downarrow }^{\prime }})],\\ {\Psi }_T(\mathbf{R})=\prod _{i{j}^{\prime }}{f}_J(r_{i{j}^{\prime }}){\Phi }_{\mathrm{BCS}}(\mathbf{R}),\end{array}$$

$$\varphi (\mathbf{r})=\sum _n{\alpha }_n{e}^{i{\mathbf{k}}_{\mathbf{n}}\cdot \mathbf{r}}+\tilde{\beta }(r),$$

$$\tilde{\beta }(r)=\beta (r)+\beta (L-r)-2\beta (L/2)\text{mbox}forr\le L/2,$$

$$=0\text{mbox}forr>L/2,$$

$$\beta (r)=[1+cbr][1-e^{-dbr}]\frac{e^{-br}}{dbr},$$

$$\begin{array}{rl}\frac{\partial \tilde{\beta }(r)}{\partial {\alpha }_i}& =\frac{\partial \tilde{\beta }(r)}{\partial r}\frac{\partial r}{\partial {\alpha }_i}\\ & =\frac{\partial \tilde{\beta }(r)}{\partial r}\frac{1}{2}\frac{2\Delta {\alpha }_{ij}}{\sqrt{\Delta x_{ij}^2+\Delta y_{ij}^2}}\frac{\partial \Delta {\alpha }_{ij}}{\partial {\alpha }_i}\\ & =\frac{\partial \tilde{\beta }(r)}{\partial r}\frac{\Delta {\alpha }_{ij}}{r}\frac{\partial \Delta {\alpha }_{ij}}{\partial {\alpha }_i}\\ & =\pm \frac{\partial \tilde{\beta }(r)}{\partial r}\frac{\Delta {\alpha }_{ij}}{r},\end{array}$$

Full text

††journal: Journal of Low Temperature Physics

1 Introduction

In recent years cold atomic gas experiments have seen novel developmentsBloch:2008 ; Giorgini:2008 ; Levinsen:2015 . With the advent of Feshbach resonances it is now possible to probe the interactions of these systems in regimes ranging from tightly bound Bose-Einstein Condensate (BEC) dimers to weakly interacting Bardeen-Cooper-Schrieffer (BCS) pairs. This has allowed for the possibility of experimental verification of ground state properties, providing a strong motivation for further investigation into these cold dilute gas systems. The BEC-BCS crossover of cold atomic gases is of particular interest due to the existence of a scale independent unitary regime. In application to the crossover, mean-field theory provides a quantitatively inadequate description, which is, however, useful as a guidepost. As a result, the unitary regime has been the target of several first-principles Quantum Monte Carlo attempts to determine the ground state properties of these dilute cold atomic gas systems Carlson:2003 ; Chang:2004 ; Astrakharchik:2004 ; Forbes:2011 ; Gandolfi:2011 ; Forbes:2012 . Intriguingly, BEC-BCS crossover of cold atomic gases is closely related to the physics of neutron matter in compact stars, which is found on the BCS side of the crossover. Gezerlis:2008 ; Carlson:2012 ; Stein:2012 ; Gandolfi:2015 ; Buraczynski:2016 ; Lacroix:2017 .

A rich area of study subject to ongoing investigation is that of low dimensionality Gunter:2005 ; Liu:2010 ; Martiyanov:2010 ; Valiente:2011 ; Frohlich:2011 ; Feld:2011 ; Orel:2011 ; Makhalov:2014 ; Bauer:2014 ; Mulkerin:2015 ; He:2015 ; Klawunn:2016 ; Anderson:2015 ; He:2016 ; Wong:2015 ; Murthy:2015 ; Ries:2015 ; Fenech:2016 ; Boettcher:2016 ; Rammelmuller:2015 ; Martiyanov:2016 ; Cheng:2016 ; Luciuk:2017 . These dilute cold atomic gases have been trapped using anisotropic potentials, resulting in a quasi-2D pancake-shape gas cloud. Very recent experiments have used box potentials to directly probe the physics of homogeneous two-dimensional Fermi gases Hueck:2017 . These systems of reduced dimensionality display properties that are distinct from the analogous in 3D phenomena. As will be discussed in detail in the following section, the main difference arises due to the logarithmic dependence on the coupling that appears in 2D. Mean-field theory BCS has also been applied to 2D Miyake:1983 ; Randeria:1990 . As in 3D, the 2D regime is not well described by the two-dimensional BCS theory in the crossover region. As a result, the determination of ground-state properties of strongly interacting Fermi gases has been attempted with Quantum Monte Carlo methods starting with a pioneering calculation using DMC Bertaina:2011 , which was later updated using Auxiliary-Field Quantum Monte Carlo (AFQMC) Shi:2015 , as well as DMC using a more sophisticated wave function that included several variational parameters Galea:2016 . These previous works focused on the determination of the contact parameter and ground state energies throughout the crossover.

In this paper we start with the two-body interaction from which our effective range and scattering length are determined. This is of particular importance to low-energy scattering phenomena, and is directly applied to the many body 2D s-wave problem. In section 2 we first provide a self-contained discussion of scattering in 2D using the partial-wave expansion. Then we define the effective range and scattering length parameters and go over their determination. Finally we discuss the formation of bound states in 2D, noting the substantial differences between 2D scattering and the well-known 3D scattering phenomena. In section 3 we discuss the strongly interacting 2D Fermi gas in the BEC-BCS crossover. First, we give a detailed discussion of the pairing function used in the many-body wave function. We explicitly determine the variational parameters used in the pairing function. A brief overview of the DMC method is given, then explicit variationally optimized DMC results for a range of interactions are shown. The equation of state is fit to these DMC results and then the corresponding chemical potential is calculated.

2 Scattering

The two-body problem is the starting point for our many-body study. Since three-dimensional (3D) scattering is a more familiar topic, we make sure to highlight crucial differences between this theory and the 2D one. We start with a partial-wave expansion of the time-independent Schrödinger equation and then define relevant scattering parameters, namely the scattering length $a_{2D}$ and effective range $r_{e}$. Next we show our techniques for determining $a_{2D}$ and $r_{e}$. Finally we describe and illustrate bound states in 2D.

2.1 Partial Wave Expansion

Describing the interaction of two particles in their center of mass reference frame, the two-body problem is reduced to that of one body with mass $m_{r}$ in the presence of a potential that represents the interaction. This $m_{r}$ is the reduced mass of the two-body system, and simplifies to $m/2$ if the particles have equal mass. The problem is further simplified by taking the interaction to be spherically symmetric. In this case the potential $V({\bf r})\rightarrow V(r)$, where $r=|{\bf r}|=\sqrt{r^{2}_{x}+r^{2}_{y}}\,$ is the separation distance between particles. The time-dependent wave functions for such a system can be written as $\psi({\bf r},t)=\psi({\bf r})e^{-iEt/\hbar}$, where $\psi({\bf r})$ is an eigenfunction of the time-independent Schrödinger equation:

[TABLE]

and the scattering energy $E$ is the eigenvalue. The second angle (usually denoted with $\phi$ in physics) which ranges from 0 to $\pi$ in 3D has no meaning in 2D. Inserting the Laplacian for 2D polar coordinates, we can rewrite Eq. (1) as

[TABLE]

where $k^{2}=2m_{r}E/\hbar^{2}$. Now we perform a partial-wave expansion, separating the wave function into a sum of the product of radial and angular terms ($l$ will turn out to be the orbital angular momentum):

[TABLE]

After this substitution and some simple rearranging we find the following equation, which must be satisfied for all $l$:

[TABLE]

where we have removed the summation over partial waves. This can be rewritten as

[TABLE]

where $f$ depends only on $r$ and $g$ includes all the angular dependence. If we imagine keeping $r$ fixed, then g(\theta)=\big{(}1/T_{l}(\theta)\big{)}\partial^{2}T_{l}(\theta)/\partial\theta^{2} must be the same for all values of $\theta$ if Eq. (5) is to be consistent. The angular functions $T_{l}(\theta)$, therefore, must satisfy the wave equation:

[TABLE]

where $c_{l}$ is some constant (which can be different for each value of $l$). The general solution can be written as

[TABLE]

Considering, for example, an incident beam at $\theta=0$, we assume symmetric scattering and can thereby eliminate the odd $\sin$ function dependence. Then, with the simple periodic condition $T_{l}(0)=T_{l}(2\pi)$, we can determine that $\sqrt{c_{l}}$ must be an integer $l$ and Eq. (6) can be rewritten as Adhikari:1986a :

[TABLE]

with the normalized solution

[TABLE]

Like the Legendre polynomials involved in the analogous 3D partial wave expansion, the $T_{l}(\theta)$ functions are linearly independent and therefore we were justified in removing the summation for Eq. (4). Using Eq. (8) to simplify the angular term and making the substitution $R_{l}(r)\rightarrow u_{l}(r)/\sqrt{r}$ to simplify the radial term, Eq. (4) becomes

[TABLE]

Rearranging, we find the Schrödinger equation for the 2D reduced radial wave function $u_{l}(r)$:

[TABLE]

Differing from the 3D case, the solution to the 2D equation is related to the radial wave function by $u_{l}(r)=\sqrt{r}R_{l}(r)$ instead of $u_{l}(r)=r\,R_{l}(r)$.

For purely s-wave scattering we consider only the $l=0$ partial wave:

[TABLE]

The wave function $u_{0}(r)$ can be solved for numerically, although analytic solutions exist for simple potentials such as the square well. The singularity at $r=0$ in Eq. (12) means that boundary conditions must be carefully selected. In the asymptotic region, when the radial separation is larger than the range of the potential, Eq. (12) simplifies to

[TABLE]

We stress at this point the significance of the $1/r^{2}$ term in this equation: note that it is present even for purely $s$-wave scattering. This is to be compared with the 3D case, where a centifugal barrier is only present for beyond-$s$-wave partial waves. As a result of this $1/r^{2}$, together with the different relationships between $u_{0}(r)$ and $R_{0}(r)$, the 2D problem always supports a bound state for pairwise attractive interactions (regardless of the strength of the interaction), in direct contradistinction to what goes on in the 3D problem.

A solution for $u_{0}(r)$ can be written in terms of Bessel functions of the first ($J_{n}$) and second ($N_{n}$) kind Khuri:2009 :

[TABLE]

By calculating $u_{0}(r)$ at two radial separations $r_{1}$ and $r_{2}$ beyond the range of the potential (always for a small $k$), the s-wave phase shifts $\delta_{0}$ can be calculated with the relation

[TABLE]

where $\tilde{K}=u_{0}(r_{2})\sqrt{r_{1}}/u_{0}(r_{1})\sqrt{r_{2}}$.

2.2 Scattering length and effective range

To define the $s$-wave scattering length $a_{2D}$ and the effective range $r_{e}$ we first look at the zero-energy Schrödinger equation in the asymptotic region. We take Eq. (12) outside the range of the potential and set $k=0$ to find:

[TABLE]

where $y_{0}(r)$ represents the asymptotic form of $u_{0}(r)$. In general, we can write the solution $y_{0}(r)=\sqrt{r}(\alpha+\beta\log(r))$ Khuri:2009 . In this and all other instances we take $\log$ to represent the natural logarithm. The 2D scattering length $a_{\rm 2D}$, in analogy to the 3D interpretation, is defined as the $r$-intercept of $y_{0}(r)$. For sufficiently short range potentials, such as those used in Section 3, $a_{\rm 2D}$ can also be given by the $r$-intercept of $u_{0}(r)$ for even strongly bound states. Using the condition $y_{0}(a_{\rm 2D})=0$, where $a_{\rm 2D}$ is the 2D scattering length, we set $\alpha=-\beta\log(a_{\rm 2D})$ and the solution becomes $y_{0}(r)=\beta\sqrt{r}\log(r/a_{\rm 2D})$. The choice for $\beta$ will influence the effective range; we set $\beta=-1$ and use the solution

[TABLE]

An analogous parameter to the above $\beta$ is encountered for the 3D case, where it can be determined by setting $y_{0}(0)=1$. Ideally, we would use this condition for 2D scattering as well, however $y_{0}(r)$ in 2D does not extrapolate to $r<0$ and always approaches 0 at the origin. This is related to the fact that $a_{\rm 2D}$, unlike the 3D scattering length, can never be negative. In other words, a 2-particle bound state exists for even arbitrarily weak attraction. Our choice of $\beta=-1$ for the 2D case is consistent with work done by Adhikari et al.Adhikari:1986 .

The effective range is related to the area between $u_{0}(r)$ and $y_{0}(r)$. In 2D this is defined by the integral Adhikari:1986 :

[TABLE]

and is the second-order term in the effective-range expansion relating low-energy phase shifts $\delta_{0}$ to the scattering parameters $a_{\rm 2D}$ and $r_{e}$. In 2D for small values of $k$: Khuri:2009

[TABLE]

where $\gamma\approx 0.577215$ is Euler’s constant. Although differing from the 3D expansion, this relationship has the same implication; low-energy scattering is independent of the details of the potential. The logarithmic scattering-length dependence in Eq. (19) is characteristic of 2D interactions and also appears in the asymptotic solution $y_{0}(r)$, Eq. (17), which describes the kinetic energy alone.

To determine $a_{2D}$ and $r_{e}$ for an arbitrary potential we solve for $u_{0}(r)$ numerically and determine the asymptotic form $y_{0}(r)$ by extrapolating from two points, $r_{1}$ and $r_{2}$, outside the range of the potential. We consider two generalizations of Eq. (17), and use the fact that $u_{0}(r)=y_{0}(r)$ for both $r_{1}$ and $r_{2}$ to write

[TABLE]

where $\xi$ has been introduced (in place of $\beta$) for the purposes of determining $a_{\rm 2D}$ by extrapolation. This is done using the following equations, which we find by working from Eq. (20):

[TABLE]

The scattering length is then used to determine $y_{0}(r)$ as given in Eq. (17), which is then in turn used to scale $u_{0}(r)$ such that $u_{0}(r)=y_{0}(r)$ at $r_{1}$ and $r_{2}$. The effective range $r_{e}$ can then be determined by solving Eq. (18), which depends on the difference between $y^{2}_{0}(r)$ and $u^{2}_{0}(r)$.

In order to find the correct asymptotic form $y_{0}(r)$, it’s important to solve $u_{0}(r)$ up to sufficiently large $r$, outside the range of the potential. Another technical detail has to do with the scattering energy $E=\hbar^{2}k^{2}/2m_{r}$. The parameters $a_{\rm 2D}$ and $r_{e}$ are defined in the limit of $E\rightarrow 0$. To make sure that finite-energy scattering effects are not influencing our determination of the parameters, we reduce $E$ until the results for $a_{\rm 2D}$ and $r_{e}$ have converged.

2.3 Bound States in 2D

By plotting the scattering length as a function of the potential depth, we can visualize the formation of bound states. In this respect, a difference exists between the 2D and 3D scattering theories. As the depth of $V(r)$ is increased, the scattering length $a_{\rm 2D}$ approaches [math] then diverges to $+\infty$ when a new bound state is created; whereas $a_{\rm 3D}$ changes from $-\infty$ to $+\infty$ when a new bound state is formed. In each case a scattering length of $+\infty$ corresponds to a weakly bound state that becomes tighter as scattering length decreases. The binding energy of the particle pair in 2D is given by

[TABLE]

where there would also be small correction terms for a finite effective range, but we use this zero-range expression for concreteness. This compares to the 3D case where, for equal mass particles, the binding energy of the two particle state is $\epsilon_{b,{\rm 3D}}=-\hbar^{2}/(ma^{2}_{\rm 3D})$. We note that an alternate definition of $a_{\rm 2D}$ is sometimes used in other work: $a^{\prime}_{\rm 2D}=a_{\rm 2D}e^{\gamma}/2$, such that $\epsilon_{b}=-\hbar^{2}/(ma^{\prime 2}_{\rm 2D})$.

For the examples in this section and the DMC results presented in Section 3, we use the modified Pöschl-Teller potential:

[TABLE]

where $r$ is the interparticle spacing. This potential is purely attractive and continuous. The parameters $v_{0}$ and $\mu$ roughly correspond to the depth and inverse width respectively and are tuned such that $V(r)$ reproduces the desired scattering parameters $a_{\rm 2D}$ and $r_{e}$.

For a given $\mu$, the scattering length and effective range exhibit a repetitive pattern of spiking up and then decaying as a function of $v_{0}$. This is illustrated in Fig. 1, where we have set $\mu/k_{F}=100$. Here we have introduced the Fermi wave vector $k_{F}$ that has units of inverse length (as does $\mu$) and is related to the 2D number density $n$ by:

[TABLE]

When $a_{\rm 2D}$ diverges to $+\infty$, a new (and initially, arbitrarily weak) bound state is created. In other words, at these values of $v_{0}$, the potential becomes deep enough to support an additional bound state. Near these locations, in this example, we see the effective range become very large. The $v_{0}$ values where $a_{\rm 2D}\rightarrow\infty$ depend on the specific $\mu$ value selected. In the case of an attractive square well potential, the effective range integral in Eq. (18) will give a negative number for very strongly bound states, causing the effective range to be imaginary. If plotting scattering parameters for the square well in the same style as Fig. 1, we would see a very similar plot. A major distinction is that $r_{e}$ would continue decreasing to zero and then become imaginary as $v_{0}$ is increased. Instead of gradually increasing as we approach a new bound state, $r_{e}^{2}$ becomes increasingly large and negative before diverging to a large positive value after the bound state threshold is surpassed.

Starting in the region where only one bound state can exist (plotted in the inset of Fig. 1), we will look at the wave function evolution as the depth of $V(r)$ is increased. In the top panel of Fig. 2, we show the reduced radial wave function $u_{0}(r)$ (solid line) and the asymptotic form $y_{0}(r)$ (dashed line), as defined by Eq. (13) in the limit as $k\rightarrow 0$ and Eq. (16) respectively. These wave functions are plotted for $v_{0}=0.1$, 0.2, and 0.4, with the corresponding potentials shown in the bottom panel. We have expressed $V(r)$ in units of the Fermi energy $\epsilon_{F}=\hbar^{2}k^{2}_{F}/2m$ and set $\mu/k_{F}=100$. The 3 states in this figure are weakly bound and have large scattering lengths which are identified as the points where $y_{0}(r)$ would become zero. These values are determined by extrapolating $u_{0}(r)$ as described in Section 2.2. In this example, the asymptotic zone is reached far before the wave functions cross the $r$-axis and therefore $u_{0}(a_{\rm 2D})=y_{0}(a_{\rm 2D})=0$.

As $v_{0}$ is increased further, the non-zero node of $y_{0}(r)$ becomes increasingly central (i.e., $a_{\rm 2D}\rightarrow 0$). As depicted in Fig. 3, we find that $y_{0}(r)$ and $u_{0}(r)$ cross the $r$-axis at dramatically different locations for very strongly bound states and look far more distinct than at smaller $v_{0}$. In the bottom panel, we plot the effective range integrand in Eq. (18). Here the dotted line simply marks the $r$-axis. Curves correspond to wave functions plotted in the top panel, which can be distinguished by line thickness. For $v_{0}=1.2$, the integrand is almost completely positive and peaks where the difference between $u_{0}(r)$ and $v_{0}(r)$ is maximum. As $v_{0}$ is increased, the features become more pronounced and we find large negative contributions to the effective range integral. For this example we find the positive contributions are dominant for any $v_{0}$. As discussed above, this is not generally true (e.g., for the square well potential where $r_{e}$ becomes imaginary).

When a new bound state is formed, the scattering length diverges discontinuously to $+\infty$ and an extra node exists. In Fig. 4, we show the wave function behaviour past this threshold, where the potential supports two bound states. We plot $V(r)/\epsilon_{F}$ in the bottom panel as was done in Fig. 2. In this figure, however, the scale has increased by an order of magnitude. The scattering length of each state is roughly the same as $a_{\rm 2D}$ for the equivalent state in Fig. 2.

3 Strongly Interacting 2D Fermi gases in the BEC-BCS Crossover

Now we shift to the many-body context of dilute Fermi gases with tunable interactions. We study interaction strengths where the gas is in between a BEC state of tightly bound pairs (dimers) and a weakly paired BCS superfluid. In this regime, where the coupling of opposite-spin particles is intermediate, the gases are said to be strongly interacting and their properties are not fully understood. As expected, we find that the mean-field BCS calculation, which gives the correct energy on each side of the crossover, is unreliable in between.

In this section we describe our many-body system including the interaction parametrizations and many-body wave function. We briefly introduce our DMC method before showing ground-state energy results for a range of interaction strengths. We first calculate the energy per particle and then parametrize an equation of state (EOS) in order to determine the chemical potential.

3.1 The BEC-BCS Crossover in 2D

Due to the omnipresent bound state, and therefore a positive scattering length for all interaction strengths, identifying the exact region of BEC-BCS crossover point is not as obvious in 2D as in 3D. Here the crossover interaction strength is chosen to be the value at which the chemical potential switches signs; this is a reasonably intuitive choice. For $k_{F}a_{2D}\gg 1$ we encounter the BCS limit and for $k_{F}a_{2D}\ll 1$ we have the corresponding BEC limit.

Calculations are done for a range of interaction strengths, defined as

[TABLE]

in order to determine the gas properties for a large fraction of the crossover. The number density of the many-body system is fixed such that the Fermi wave vector $k_{F}$, as defined in Eq. (24), is constant. To change the interaction strength we vary $a_{\rm 2D}$. Given that $n=N/A$ is the number density of the system (where $N$ is the number of particles and $A$ the area of the periodic box) and $r_{0}=1/\sqrt{\pi n}$ is the mean interparticle spacing, the diluteness requirement is satisfied by taking $r_{e}\ll r_{0}$. We maintain a constant effective range of $k_{F}r_{e}=0.006$ by adjusting $\mu$ as $v_{0}$ is varied.

3.2 Many-Body Wave Function

To describe the strongly interacting Fermi gas for any attraction strength, we use the Jastrow-BCS many-body trial wave function Forbes:2011 ; Gandolfi:2011 ; Chang:2004 ; Astrakharchik:2004 ; Forbes:2012 :

[TABLE]

where the anti-symmetry requirement of $\Psi_{T}({\bf R})$ for the Fermi gas is enforced by the operator $\cal A$. Correlations between interacting particles are accounted for through the Jastrow terms $f_{J}(r_{ij^{\prime}})$. The pairing functions $\phi({\bf r})$ are expressed as

[TABLE]

which contains variational parameters $\alpha_{n}$ for each momentum state up to some level $n_{max}$ and the $\beta(r)$ function to account for higher-momentum contributions. This two-body function encodes details of the many-body system which vary with the interaction strength.

The spherically symmetric short-range function is given by:

[TABLE]

which contains variational parameters $b$, $c$ and $d$. This form of the beta function has been used for 3D calculations, and we have explicitly checked its behaviour in 2D. Specifically, when calculating the local energy given by $\Psi^{-1}_{T}({\bf R})\hat{H}\Psi_{T}({\bf R})$, we need to evaluate terms of the form $\partial\tilde{\beta}(r)/\partial\alpha_{i}$, where $\alpha_{i}$ is the coordinate of a specific particle (e.g., $x_{2},\,y_{3}\ldots$) and $r=\sqrt{\Delta x^{2}_{ij}+\Delta y^{2}_{ij}}$ is the radial separation between two particles. Defining $\Delta\alpha_{ij}$ as the projection of $r$ along a coordinate (e.g., $\Delta x_{ij}$ or $\Delta y_{ij}$) we can write:

[TABLE]

where $\partial\Delta\alpha_{ij}/\partial\alpha_{i}$ can be positive or negative 1. For example, $\partial(x_{2}-x_{5})/\partial x_{2}=1$ and $\partial(x_{2}-x_{5})/\partial x_{5}=-1$.

The result in Eq. (29) has a singularity at $r=0$ due to the $1/r$ term. This can cause large fluctuations in the local energy for small $r$, therefore we define the variational parameter $c$ in Eq. (28) such that

[TABLE]

Making use of L’Hôpital’s rule, it is mostly straightforward to show that

[TABLE]

In this work, we set $b=0.5k_{F}$ and $d=5$, as done by Gandolfi et al. Gandolfi:2011 for the 3D unitary Fermi gas. With these values of $b$ and $d$, we find $c\simeq 3.5$.

3.3 DMC

To determine ground-state properties of Fermi gases we use DMC to project the ground state $\Phi_{0}$ from the trial wave function $\Psi_{T}({\bf R})$. This is done by propagating in imaginary time $\tau=it$:

[TABLE]

where the trial energy $E_{T}$ is a constant offset applied to the Hamiltonian.

DMC expectation values are determined by averaging over a set of equilibrated configurations. In this work we use the mixed estimate to calculate the energy:

[TABLE]

In the second line we wrote the explicit form of the imaginary-time evolved ket. The propagator is then used to act on the trial wave function bra straightforwardly due to the fact that it commutes with the Hamiltonian. Finally taking $\tau$ $\rightarrow$ $\infty$ we see that this is the energy of the ground state.

3.4 Equation of State

A mean-field calculation Randeria:1990 gives a ground-state energy per particle of

[TABLE]

for Fermi gases in the BEC-BCS crossover, where the binding energy $\epsilon_{b}$ is given by Eq. (22). This is expected to be accurate for weakly paired systems in the BCS limit in which $E/N\rightarrow E_{\rm FG}$. The BEC limit of tightly bound pairs is also expected to be reasonably well described by mean field as the energy scale grows rapidly by many orders of magnitude due to large binding energies. The QMC results vary dramatically from the mean-field description in the crossover but become increasingly similar to mean-field predictions in each limit.

Our DMC calculations for a range of interactions strengths are shown in Fig.. 5. The errors represent statistical uncertainty, which becomes larger as the energy scale increases. We use the Jastrow-BCS wave function $\Phi_{\rm BCS}(\mathbf{R})$, Eq. (26), which contains parameters that are optimized for each $\eta$ independently. We have compared our results to previous ab-initio work, finding significantly lower energies than prior ground-state DMC results Bertaina:2011 in the crossover regime and excellent agreement with AFQMC Shi:2015 . More recently another QMC study has emerged Madeira:2017 that finds DMC results in agreement with ours. Note that we here provide more results in the BCS side than where available previously Galea:2016 .

In order to calculate other ground-state properties of strongly interacting Fermi gases, we calculate an EOS for our thermodynamic-limit energy results $E_{\rm TL}$. This quantity is a finite-size corrected version of the results in Fig. 5. The correction ranges from zero on the BEC side of the crossover to $\sim 0.041~{}E_{\rm FG}$ in the BCS regime. Using similar methods as previous ab initio studies Bertaina:2011 ; Shi:2015 ; Galea:2016 , we parametrize the EOS using three functions. In the crossover regime we fit to a 7th-order polynomial:

[TABLE]

This is joined by a dimer form in the BEC regime Levinsen:2015 and an expansion in $1/\eta$ in the BCS regime Engelbrecht:1992 . The dimer form is given by:

[TABLE]

where $x=\log[4\pi/(k_{F}a_{d})^{2}]\approx 3.703-2\eta$ (for the dimer scattering length $a_{d}\approx 0.557a_{\rm 2D}$ Petrov:2003 ) and $d=\log\pi+2\gamma+0.5$. The BCS form is given by:

[TABLE]

Values of $c_{i}$ in Eq. (36) and Eq. (37) are determined using continuity conditions for $f$, $\partial f/\partial\eta$, and $\partial^{2}f/\partial\eta^{2}$ at the matching points. Our EOS parameters are included in Tab. 1. The matching points were selected as $\eta=-0.25$ and $\eta=2.5$. We found these values result in the most optimal overall fit while including as much of the crossover polynomial function as possible. Also, we ensure that our matching point for Eq. (36) is selected on the BEC side of the crossover.

3.5 Chemical Potential

We have determined the chemical potential $\mu$ using our EOS. In order to derive a relationship, we define

[TABLE]

where $E_{\rm TL}(\eta)$ is the total ground-state energy. This quantity is related to our parametrized EOS, where we fit to $f(\eta)=(E_{\rm TL}/N-\epsilon_{b}/2)/E_{\rm FG}\,$; the new quantity $\zeta(\eta)$ can easily be determined by adding $\epsilon_{b}/2$ in units of $E_{\rm FG}$ to $f(\eta)$. Noting that $\epsilon_{F}=(\pi N\hbar^{2})/(mA)$ in 2D (using $k^{2}_{F}=2\pi N/A$), the chemical potential is related to $\zeta(\eta)$ as follows:

[TABLE]

In the fourth line we have evaluated $\partial\eta/\partial N$ using

[TABLE]

Expressing the chemical potential in units of the non-interacting Fermi gas, we find

[TABLE]

Our result is plotted in Fig. 6, which shows the chemical potential with half of the two-body binding energy subtracted. Our DMC determination is the thick blue line, which is qualitatively similar to the chemical potential one gets starting from the earlier DMC results Bertaina:2011 (which we generated using an analogous fitting strategy). The dotted vertical line corresponds to the interaction strength ($\eta\approx 0.65$) where the chemical potential changes sign. Comparing our results with experimentally extracted values of the chemical potential Boettcher:2016 ; Enss:2015 , we find a nice match in the deep BEC regime. Differences become more significant in the BCS limit, probably due to finite-temperature and quasi-2D effects in the experiment.

4 Summary and Conclusion

In this paper we presented a detailed discussion of scattering in the context of a two-body system confined to two dimensions. We then determined the radial wave function and defined its asymptotic form, which in a purely attractive potential always produces a bound state in contrast to the 3D case. Continuing to explore 2D scattering phenomena we illustrated the dependence of the effective range $r_{e}$ and the scattering length $a$ on the attractive potential’s strength. Varying the potential we plotted the divergence of $r_{e}$ and $a$ upon the approach of new bound states. We also demonstrated the radial wave function and associated asymptotic form dependence on the attractive potential strength, where stronger attractive potentials produce more tightly bound states with small positive scattering lengths. Having quantified the 2D two body interacting system in detail, we moved to the many-body problem. Starting with the BCS determinant we looked at its composite pairing functions and specifically the variational parameters. We presented the explicit form of the $\beta$ pairing function, and briefly described the mixed-estimate method. Then we provided QMC ground-state energy results for various interaction strengths within the BEC-BCS crossover range. Fitting to these energies with a 7th-order polynomial an equation of state was determined. The expression for the chemical potential was derived according to the EOS, for which fitting parameters were explicitly provided. To conclude, the careful study of two-dimensional scattering properties in conjunction with a non-perturbative many-body method containing several variational parameters (like DMC), has led to dependable predictions for the properties of strongly correlated physical systems. Overall, we observe that two-dimensional strongly interacting cold Fermi gases constitute an exciting new development, where theory can be confronted by impressive experimental work.

Acknowledgements.

The authors would like to thank G. E. Astrakharchik, T. Enss, J. Thywissen, and E. Vitali for helpful discussions. This work was supported in part by the Natural Sciences and Engineering Research Council (NSERC) of Canada, the Canada Foundation for Innovation (CFI), the Early Researcher Award (ERA) program of the Ontario Ministry of Research, Innovation and Science, the US Department of Energy, Office of Nuclear Physics, under Contract DE-AC52-06NA25396, and the LANL LDRD program. Computational resources were provided by SHARCNET, NERSC, and Los Alamos Open Supercomputing. The authors would like to acknowledge the ECT* for its warm hospitality during the “Superfluidity and Pairing Phenomena” workshop in March 2017, where part of this work was carried out.