Pairing in two-dimensional Fermi gases with a coordinate-space potential

TL;DR

This paper investigates pairing phenomena in two-dimensional Fermi gases with finite-range interactions, deriving a mean-field gap equation, numerically solving it, and using Monte Carlo simulations to evaluate energies and pairing gaps, providing benchmarks for future theories.

Contribution

It introduces a mean-field pairing gap equation for finite-range potentials in 2D Fermi gases and uses Monte Carlo methods to evaluate energies and pairing gaps, offering new benchmarks.

Findings

Short effective ranges yield results identical to zero-range models.

Numerical solutions of the gap equation across interaction strengths.

Monte Carlo energies used to extract pairing gaps via odd-even staggering.

Abstract

In this work we theoretically study pairing in two-dimensional Fermi gases, a system which is experimentally accessible using cold atoms. We start by deriving the mean-field pairing gap equation for a coordinate-space potential with a finite interaction range, and proceed to solve this numerically. We find that for sufficiently short effective ranges the answer is identical to the zero-range one. We then use Diffusion Monte Carlo to evaluate the total energy for many distinct particle numbers; we employ several variational parameters to produce a good ground-state energy and then use these results to extract the pairing gap across a number of interaction strengths in the strongly interacting two-dimensional crossover. Extracting the gap via the odd-even energy staggering, our microscopic results can be used as benchmarks for other theoretical approaches.

Equations

$$e_b^0=-\frac{4{\hslash }^2}{m{a}^2{e}^{2\gamma }}$$

$$E_F=\frac{{\hslash }^2{k}_F^2}{2m}$$

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

$$\eta =\ln (k_{F}a).$$

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

$${\Delta }_{\mathrm{gap}}=\{\begin{array}{l}\sqrt{2E_F|e_b|},\mu >0\\ E_F+|e_b|/2,\mu <0\end{array}$$

$$\mu =E_F+e_b/2$$

$$\Delta (\mathbf{k})=-\sum _{{\mathbf{k}}^{\prime }}⟨\mathbf{k}|V|{\mathbf{k}}^{\prime }⟩\frac{\Delta ({\mathbf{k}}^{\prime })}{2E({\mathbf{k}}^{\prime })}$$

$$E(\mathbf{k})=\sqrt{{\left(\frac{{\hslash }^2\text{absolutevalue}{\mathbf{k}}^2}{2m}-\mu \right)}^2+{\Delta }^2(\mathbf{k})}$$

$$⟨\mathbf{k}|V|{\mathbf{k}}^{\prime }⟩$$

$$=\frac{1}{L^2}\int d\mathbf{r}\int d{\mathbf{r}}^{\prime }{e}^{i\mathbf{k}\cdot \mathbf{r}}⟨\mathbf{r}|V|{\mathbf{r}}^{\prime }⟩e^{-i{\mathbf{k}}^{\prime }\cdot {\mathbf{r}}^{\prime }}$$

$$=\frac{1}{L^2}\int d\mathbf{r}\int d{\mathbf{r}}^{\prime }\sum _{n,m}{i}^n{J}_n(kr)e^{in({\theta }_k-{\theta }_r)}V(r)$$

$$\times \delta (\mathbf{r}-{\mathbf{r}}^{\prime })(-i{)}^m{J}_m(k^{\prime }{r}^{\prime })e^{-im({\theta }_{k^{\prime }}-{\theta }_{r^{\prime }})}$$

$$=\frac{1}{L^2}\sum _{n,m}{i}^n(-i{)}^m{e}^{in{\theta }_k}{e}^{-im{\theta }_{k^{\prime }}}$$

$$\times {\int }_0^{\infty }drr{J}_n(kr)V(r)J_m(k^{\prime }r)$$

$$\times {\int }_0^{2\pi }d{\theta }_r{e}^{i(m-n){\theta }_r}$$

$$=\frac{2\pi }{L^2}\sum _{n=-\infty }^{\infty }{e}^{in({\theta }_k-{\theta }_{k^{\prime }})}{V}_n(k,k^{\prime }).$$

$$e^{i\mathbf{k}\cdot \mathbf{r}}=\sum _{n=-\infty }^{\infty }{i}^n{J}_n(kr)e^{in\theta }$$

$$\Delta (\mathbf{k})\equiv \sum _{n=-\infty }^{\infty }{\Delta }_n(k)e^{in{\theta }_k}$$

$$-\frac{2\pi }{L^2}\sum _{{\mathbf{k}}^{\prime }}\sum _{n,m}{e}^{im{\theta }_{k^{\prime }}}{e}^{in({\theta }_k-{\theta }_{k^{\prime }})}{V}_n(k,k^{\prime })\frac{{\Delta }_m(k^{\prime })}{2E(k^{\prime })}$$

$$=-\frac{2\pi }{(2\pi {)}^2}\int d{\mathbf{k}}^{\prime }\sum _{n,m}{e}^{i(m-n){\theta }_{k^{\prime }}}{e}^{in{\theta }_k}{V}_n(k,k^{\prime })\frac{{\Delta }_m(k^{\prime })}{2E(k^{\prime })}$$

$$=-\frac{1}{2\pi }\sum _{n,m}{e}^{in{\theta }_k}{\int }_0^{\infty }d{k}^{\prime }{k}^{\prime }{V}_n(k,k^{\prime })\frac{{\Delta }_m(k^{\prime })}{2E(k^{\prime })}$$

$$\times {\int }_0^{2\pi }d{\theta }_{k^{\prime }}{e}^{i(m-n){\theta }_{k^{\prime }}}$$

$$=-\sum _n{e}^{in{\theta }_k}{\int }_0^{\infty }d{k}^{\prime }{k}^{\prime }{V}_n(k,k^{\prime })\frac{{\Delta }_n(k^{\prime })}{2E(k^{\prime })}$$

$${\Delta }_l(k)=-{\int }_0^{\infty }d{k}^{\prime }{k}^{\prime }{V}_l(k,k^{\prime })\frac{{\Delta }_l(k^{\prime })}{2E(k^{\prime })}$$

$$\rho ={\int }_0^{\infty }dkk\left(1-\frac{\frac{{\hslash }^2{k}^2}{2m}-\mu }{E(k)}\right)$$

$$\displaystyle\hat{H}=\frac{-\hbar^{2}}{2m}\Bigg{[}\sum_{i=1}^{N_{\uparrow}}\nabla_{i}^{2}+\sum_{j^{\prime}=1}^{N_{\downarrow}}\nabla_{j^{\prime}}^{2}\Bigg{]}+\sum_{i,j^{\prime}}V(r_{ij^{\prime}})\,,$$

$${\mathbf{k}}_{\mathbf{n}}=\frac{2\pi }{L}(n_x,n_y)$$

$$E_v=\frac{\int {\Psi }_T^{\ast }(\text{R})\hat{H}{\Psi }_T(\text{R})d\text{R}}{\int {\Psi }_T^{\ast }(\text{R}){\Psi }_T(\text{R})d\text{R}}\ge E_0$$

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

$$\Psi (\tau \to \infty )=\underset{\tau \to \infty }{\lim }{e}^{-(\mathcal{H}-E_T)\tau }{\Psi }_T\to {\alpha }_0{e}^{-(E_0-E_T)\tau }{\Psi }_0$$

$$\Delta (N)=E(N+1)-\frac{1}{2}\left[E(N+2)+E(N)\right],$$

Full text

I Introduction

The study of ultracold atomic gases has witnessed tremendous progress over the last two decades Bloch:2008 ; Giorgini:2008 ; Levinsen:2014 . Historically, the interactions between different species of particles were fixed by nature; when the binding was strong enough to produce bosonic dimers, one was faced with Bose-Einstein-Condensates (BEC), whereas when the interaction was weaker (but still attractive) one dealt with Bardeen-Cooper-Schrieffer (BCS) Theory. This all changed with the development of Feschbach resonances, which use magnetic fields to tune particle interactions and allow direct experimental measurements of cold Fermi gases in the intermediate BEC-BCS crossover region.

This has enabled a detailed study of the BEC-BCS crossover with direct comparison between theory and experiment, leading to the discovery of the existence of a unitary regime displaying scale invariant properties deep within the crossover region. The BCS theory is not expected to give quantitatively reliable results for ground-state properties in this crossover region. Instead many-body approaches that describe the system from first principals are used. Some of the most successful such approaches applied to three-dimensional Fermi gases are Quantum Monte Carlo (QMC) methods, used at both zero and finite temperature Carlson:2003 ; Astrakharchik:2004 ; Chang:2004 ; Forbes:2011 ; Gandolfi:2011 ; Forbes:2012 ; Carlson:2012 ; Pilati:2014 ; Schonenberg2:2017 ; Dawkins:2017 .

Reduced dimensionality – specifically 2D systems – are a rich area which has been at the forefront of research more recently, as they display distinct and unique properties compared to that of the 3D gas systems. Experimentally, cold Fermi gases are produced in quasi-2D pancake-shape gas clouds Gunter:2005 ; Martiyanov:2010 ; Frohlich:2011 ; Feld:2011 ; Orel:2011 ; Dyke:2011 ; Sommer:2012 ; Makhalov:2014 ; Boettcher:2015 ; Ong:2015 ; Murthy:2015 ; Ries:2015 ; Fenech:2016 ; Boettcher:2016 ; Martiyanov:2016 ; Cheng:2016 ; Luciuk:2017 ; Toniolo:2017 . Both preceding and following the experimental breakthroughs, a number of theoretical approaches have been brought to bear on the subject of lower-dimensional Fermi gases Nishida:2008 ; Liu:2010 ; Valiente:2011 ; Bertaina:2011 ; Werner:2012 ; Bauer:2014 ; Shi:2015 ; Anderson:2015 ; He:2015 ; Rammelmuller:2015 ; He:2016 ; Galea:2016 ; Klawunn:2016 ; Galea:2017 ; Vitali:2017 ; Schonenberg:2017 ; Drut:2018 ; jeszenszki:2019 .

In addition to being experimentally accessible and therefore interesting in their own right, cold Fermi atomic gases can also help shed light on the physics of superfluidity in neutron-star crusts Gezerlis:2008 ; Gandolfi:2015 ; Buraczynski:2016 . In a sense, neutron matter is an attractive Fermi gas where one needs to take the effective range into account, in addition to the scattering length. While the neutron-neutron interaction is fixed by nature, one can envision cold-atomic experiments which probe the regime of relevance to neutron stars Horikoshi:2017 . More generally, two-dimensional cold gases can also help us understand nuclear “pasta” phases, where nucleons end up populating highly deformed quasi-1D or quasi-2D systems.

In earlier works, we have employed Diffusion Monte Carlo Galea:2016 ; Galea:2017 (DMC) to evaluate two-dimensional cold-gas properties such as the equation of state and the pairing gap. These involved, first, the use of a finite-range two-particle interaction and, second, the use of only $N=26$ particles in a periodic area. In the present paper, we investigate the effects of both of these choices in more detail. Specifically, we start from investigating the effect of a finite effective range in the context of another many-body approach, namely mean-field BCS theory; generally speaking, this allows us to use a simpler theory to study effects which would be much harder to tackle in the context of DMC. Then, we turn to our new DMC results for the energy and the pairing gap, which result from using larger-$N$ values; as part of this process, we try to systematize our understanding of the finite-size effects along the 2D crossover.

II Two-body binding energy

For the two-particle problem in two dimensions Miyake:1983 ; Randeria:1989 , the binding energy is:

[TABLE]

where $a$ is the 2D scattering length, $m$ is the particle mass, and $\gamma\approx 0.577215$ also known as the Euler number. The Fermi energy is:

[TABLE]

where $k_{F}$ is the Fermi wave number. In two dimensions the free-gas energy per particle is simply $E_{fg}=E_{F}/2$.

The expression in Eq. (1) involves the scattering length, $a$, but not the effective range, $r_{e}$; this is a zero-effective-range result, i.e., it is valid for $k_{F}r_{e}\rightarrow 0$, where $r_{e}$ is the effective range of the two-particle interaction. On the other hand, many Quantum Monte Carlo calculations are carried out using coordinate-space potentials of finite (but tiny) range. A popular interaction is of the (modified) Pöschl-Teller form:

[TABLE]

where $v_{0}$ and $\nu$ are parameters which we tune to match the desired scattering length $a$ and effective range $r_{e}$, and $m_{r}$ is the reduced mass. As is standard for this problem, we quantify the interaction strength by defining:

[TABLE]

While the $k_{F}$ which appears here and in Eq. (2) is a many-body quantity (which therefore does not arise at the two-body level), we find it convenient to employ it to produce dimensionless quantities throughout.

As further discussed in Ref. Galea:2017 , in two dimensions a partial-wave expansion of the time-independent Schrödinger equation leads to

[TABLE]

for the $l=0$ partial wave. Here $u_{0}(r)=\sqrt{r}R_{0}(r)$ and $R_{0}(r)$ is (proportional to) the wave function. We have numerically solved Eq. (5) to find the finite-effective-range two-body binding energy, $e_{b}$, for the two cases of $\eta=-0.5$ and $\eta=0.5$. The result is shown in Fig. 1: we see that for $\eta=0.5$ the range is largely irrelevant, since we basically get the same binding energy for all effective ranges up to $k_{F}r_{e}\approx 0.1$. The case of $\eta=-0.5$ is different, in that the binding energy changes by a few percent for such effective ranges. For both $\eta$ values, the figure clearly shows that the effective range we employ in our DMC calculations, $k_{F}r_{e}\approx 0.006$, is short enough that we don’t have to worry about finite-range effects. As a result, in what follows we will be employing $e_{b}^{0}$ and $e_{b}$ interchangeably, under the assumption that the effective range is short enough.

III Mean-Field BCS Theory

Bardeen-Cooper-Schrieffer theory for the two-dimensional case was developed before the recent explosion of activity in low-dimensional cold-atom systems Miyake:1983 ; Randeria:1989 . For $s$-wave pairing, which is our focus in this paper, the problem is analytically solvable, leading to simple expressions for the pairing gap and chemical potential. These were:

[TABLE]

and

[TABLE]

respectively. Crucially, the derivation in Ref. Randeria:1989 involved the two-body $T$-matrix, in an attempt to avoid issues arising from the use of a coordinate-space potential which involves a hard-core repulsion. The expressions above are valid in the limit of $k_{F}r_{e}\rightarrow 0$.

In our Quantum Monte Carlo calculations we employ Eq. (3); this is a purely attractive potential, so it does not give rise to problems in momentum space. Motivated by such studies, in this work we first derive the BCS gap equation in 2D for a finite-range potential; to our knowledge, this derivation and the resulting expression have not been published before. We then numerically solve the gap equation, comparing to the $k_{F}r_{e}\rightarrow 0$ solution, thereby testing if the effective range to be used in later sections (on DMC) is sufficiently small.

The BCS gap equation has the form:

[TABLE]

Where $\Delta(\mathbf{k})$ is the pairing gap function, $V$ is the potential energy, and $E(k)$ is the quasiparticle energy, given by the relationship:

[TABLE]

We are here interested only in qualitative features resulting from a small but finite effective range; thus, we are employing a free single-particle spectrum, as shown in Eq. (9); a more complete study would also involve normal-state interactions. In these equations, the gap function and potential matrix element are given in momentum space. We would like to see how that relates to the coordinate-space potential, which could be of the form of Eq. (3). We have:

[TABLE]

In the first equality we introduced two resolutions of the identity. In the second equality we plugged in the appropriate plane waves, normalized to fit inside a square of area $L^{2}$. In the third equality we employed the locality of the interaction and used (twice) the two-dimensional expansion of a plane wave:

[TABLE]

where $J_{n}$ is the Bessel function of order $n$, and $\theta$ is the angle between the $\mathbf{k}$ and $\mathbf{r}$ vectors, or $\theta=\theta_{k}-\theta_{r}$. In the fourth equality we merely rearranged terms and in the fifth equality we carried out the integral over $\theta_{r}$ and used the result to eliminate one of the sums; we also took the opportunity to introduce $V_{n}(k,k^{\prime})$ for the Bessel-transformed potential. The result is analogous to a Fourier transform of $V(r)$, in that when the potential in real space is narrow and deep, the potential in $k$-space is wide and shallow, and vice versa. This is important because in what follows we will be interested precisely in extremely deep, narrow real-space potentials, making $V(k,k^{\prime})$ very wide.

We now carry out an analogous expansion for the gap function:

[TABLE]

Inserting this, as well as the equation for potential matrix element, Eq. (10), the right-hand side of the gap equation in Eq. (8) turns into:

[TABLE]

In the first equality, we cancelled the $L^{2}$ when converting from sum to integral. In the second equality we wrote out the integral over $\mathbf{k}^{\prime}$ and in the third eqaulity we used the Kronecker delta to eliminate one of the sums.

If we now recall that the left-hand side of the BCS gap equation of Eq. (8) will also take the form Eq. (12), we can multiply both sides by $e^{-il\theta_{k}}$ and integrate, to get:

[TABLE]

In what follows, we will be interested in taking $l=0$, corresponding to the $s$-wave problem. Thus, we have arrived at the two-dimensional BCS gap equation in the thermodynamic limit. In practice, when carrying out the integral over $k^{\prime}$ one must be careful to go up to sufficiently high momenta where $V_{0}$ has died off.

In addition to the BCS gap equation, one must also self-consistently solve the constant-average-density equation:

[TABLE]

In practice, we solve Eq. (14) together with Eq. (15) iteratively: this is fixed-point iteration for an entire function. Given that we are dealing with a purely attractive potential, the process converges; as part of this, we tune the chemical potential $\mu$ to get the desired density $\rho=k_{F}^{2}/2\pi$.

Results of our mean-field calculations are given in Fig. 2 where the solid line of $k_{F}r_{e}=0$ is the analytical solution from Eq. (6). The other two curves are new results, which match very closely with the analytical solution for $k_{F}r_{e}=0.006$, and even for the larger $k_{F}r_{e}=0.1$ for positive interaction strengths. These results inticate that the choices of $k_{F}r_{e}$’s used in our DMC calculations (in Refs. Galea:2016 ; Galea:2017 and below) are small enough to get meaningful results. This is fully consistent with what we discovered at the two-body level in Fig. 1.

IV Quantum Monte Carlo

The problem we are faced with is how to tackle the following two-component Hamiltonian in two dimensions:

[TABLE]

where $m$ is the particle mass and the particle number is $N_{\uparrow}+N_{\downarrow}=N$. The $V(r_{ij^{\prime}})$ is taken to be of the modified Pöschl-Teller form, Eq. (3). In what follows, we work in a periodic simulation area, as is standard in continuum Quantum Monte Carlo approaches. Specifically, we work with $N$ particles in a square of length $L$; allowed momentum states are governed by the equation:

[TABLE]

where $n_{x}$ and $n_{y}$ are integers.

The first step in our study is always a Variational Monte Carlo (VMC) calculation: this provides an upper bound on the ground-state energy. We use R will denote a walker which represents the coordinates of all of the particles in the simulation box, and is $2N$ dimensional. Using a trial wave function, $\Psi_{T}(\textbf{R})$, we determine the variational energy, $E_{v}$, which is the upper bound on the ground-state energy $E_{0}$:

[TABLE]

This is nothing other than the Rayleigh-Ritz principle applied to a many-particle system. For a “good” wave function, this will provide a decent estimate of the ground-state energy of an interacting system. The many-dimensional integrals involved are carried out stochastically, using a large number of walkers.

The question of what constitutes a good wave function naturally arises at this point. Following Refs. Galea:2016 ; Galea:2017 , we employ a Jastrow-BCS trial wave function:

[TABLE]

where $f_{J}(r)$ is a nodeless Jastrow term, $\cal A$ is an antisymmetrizing operator, and we took $N_{\uparrow}=N_{\downarrow}$. We express the pairing function $\phi({\bf r})$ as a sum of 10 plane wave terms, each multiplied with an unknown parameter, and we capture higher-momentum contributions by a spherically symmetric function. Then, we employ VMC to determine the best possible set of these 10 parameter values; crucially, this is a process which we repeat for each new particle number $N$. As a technical aside we note that in two-dimensional physics, where the $k_{F}a$ spans multiple orders of magnitude, the Jastrow term is implemented using a table with substantially more points than in the three-dimensional problem.

The next step in our calculations is to use Diffusion Monte Carlo (DMC): this also approximates the ground-state energy using a trial wave function, but is superior in that it propagates in imaginary time to reduce the contribution of the parts of the trial wave function that are not the ground state (i.e., excited states). This is done by solving the Schrödinger equation in imaginary time. Schematically:

[TABLE]

where $E_{T}$ is the so-called “trial energy” which helps us ensure the walker number doesn’t get out of hand. In practice, we propagate up to sufficiently large imaginary time, such that there is no longer a decay going on.

Apart from technical practicalities like the trial energy and the Jastrow term, the only other consideration is the fermion sign problem. This arises due to the presence of both positive and negative nodal pockets in the complicated many-particle wave function. As is standard in zero-temperature continuum QMC algorithms, we employ the fixed-node approximation, which does not allow walkers to change sign in the wave function. Assuming the trial wave function has the appropriate physical content (e.g., pairing properties for a superfluid system), the effect of the fixed-node approximation ends up being tiny. Fortunately, both our VMC and our DMC results obey a variational property, meaning that one can keep trying to improve the answers by using a better wave function.

Even so, it’s worth highlighting that one of our goals in this work is to calculate the pairing gap via the odd-even staggering in the energy:

[TABLE]

Here $E(N)$ is the total energy of a closed shell $N=10,18,26,42,50,58$, $E(N+1)$ is the total energy of the system of a closed shell with one extra particle that does not have a partner particle to pair with, and $E(N+2)$ is one more particle added to the system which again causes all particles to be paired. We have explored a couple higher closed shells, but the statistical errors were so large that it became impossible to trust the energy differences; note that QMC methods employ determinants, so the computational cost typically scales as $N^{3}$. That means that a calculation for $N=58$ (the largest particle number employed here) is roughly 11 times more demanding than one for $N=26$ (which was the particle number used in earlier works). It’s worth emphasizing that, while a variational property applies to the total energy, it does not apply to energy differences such as in Eq. (20). That being said, similar QMC predictions for the problem of the three-dimensional pairing gap Gezerlis:2008 turned out to be verified experimentally Horikoshi:2017 .

With the systematic errors under control, the only things left to investigate are the finite-size effects (i.e., the dependence on $N$) as well as the placement of the $(N+1)$-th particle. We discuss the former in detail in the following section; as for the latter, we note that in Eq. (20) one must determine the pairing gap as the minimum. For fermions in a non-interacting gas, or in a normal gas, the placement of the next particle is trivial: it goes onto the next available momentum state. For the present case, where pairing correlations play a major role, the choice of momentum states is not quite so intuitive. This is illustrated in Fig. 3, in which the interacting system’s optimal $k$-state is not simply the next unoccupied shell. This is analogous to what we saw in Eq. (9). The VMC values can be quite close and within error so we chose to include a DMC run in the determination of the $k$-states. (Incidentally, we have encountered one case where the 27-particle energies of Ref. Galea:2016 gave different answers in VMC and DMC; we have corrected this here, always employing the lowest-possible values.)

V Results

We start our discussion of the finite-size effects by looking at the energy per particle in the non-interacting gas, see Fig. 4. The red horizontal line is the thermodynamic limit value (infinite volume, infinite particle number, constant density). We see that for small particle numbers there are large deviations from the thermodynamic-limit answer, up to 5% or more, but these go away as the particle number is increased. Of course, the question arises whether one should expect the interacting system finite-size effects to match those of the non-interacting gas.

The results of our DMC calculations within the strong coupling crossover for various interaction strengths in descending order from $\eta=0.5$ to $\eta=3.0$ are shown in Fig. 5. Looking at the general shape of the DMC energies, the bottom three panels follow roughly the same trend where the $N=10$ is large then a local minimum occurs at $N=18$ and a local maximum at $N=26$ which proceeds to have a slight dip at $N=42$ then small increase to $N=58$; this trend closely follows the shape of the free gas plot of Fig. 4. These results imply that the pairing correlations on the BCS side of the crossover region ($\eta$ greater than or equal to 1.5) have a small effect on the energy dependence on particle number. This suggests that even for other observables one could get guidance from the non-interacting gas behavior.

Turning to the stronger interaction strengths (the top two panels in Fig. 5), we find fewer similarities with the non-interacting gas trend. Going from $N=10$ to $N=18$ there is a drop, but the commonalities largely end there. The panel that is the most different from Fig. 4 is that corresponding to $\eta=0.5$ , where $N=18$ does not display a local minimum. Overall, we can say that the trend (as a function of $N$) is “flatter” for the top two panels, i.e., the results for $N=50$ and $N=58$ are very close to each other; this is a result of stronger pairing correlations. It’s crucial to note that for small $\eta$ the finite-size effects are less significant: for $\eta=0.5$ the fluctuation from the largest to the smallest value is roughly 5%, to be compared with the corresponding change of approximately 15% in the non-interacting gas. In short, for even smaller $\eta$ (say, [math] or $-1$) the finite-size effects become irrelevant.

Since we will be computing the pairing gap from the odd-even staggering of Eq. (20), we now spend some time discussing the energies of systems with $N$, $N+1$, and $N+2$ particles. Similarly to what we did in Fig. 4, we first go over the finite-size effects in the non-interacting gas, to possibly some guidance about how to interpret our DMC results. In Fig. 6 we show the result of applying the odd-even staggering formula to the non-interacting gas. Since there is no pairing in the absence of interactions, one would expect the pairing gap to be zero. However, this figure clearly shows that at lower paticle numbers the $\Delta_{fg}$ is largest and decreases as $N$ increases. This tells us that from the particle numbers we are employing $N=50$ and $N=58$ should have the least amount of error. Of course, this fact should be combined with our earlier finding, namely that the finite-size effects of the interacting gas are insignificant for small $\eta$ (but track those of the non-interacting gas for large $\eta$).

We now turn to the energies of strongly interacting systems with $N$, $N+1$, and $N+2$ particles; we show results for $\eta=1$ and $\eta=3$ in Fig. 7 and Fig. 8, respectively (statistical errors are smaller than the symbols shown). In both cases, we observe the same general pattern: the three adjacent-in-$N$ values are quite far apart from each other for small $N$, but then get closer as $N$ is increased, leading to the characteristic odd-even staggering pattern. One should not confuse this trend with the conclusion that the gap decreases as the $N$ increases, since the even partners are not placed symmetrically for small-$N$ values. This is not surprising given what we saw in Fig. 4. Looking back at the odd-even staggering prescription of Eq. (20), it is important to realize that the energies plotted are per particle, so these need to be multiplied by the appropriate particle numbers before the difference is taken. Another similarity between the two plots is that the results around $N=50$ and $N=58$ behave in roughly the same way. It’s also worth discussing the differences between Fig. 7 and Fig. 8. In Fig. 7 the result for $N+1$ is always higher than its neighbors, a fact which is not always true in Fig. 8.

To re-cap, it is certainly plausible to use a small $N$ and then employ a correction on the closed-shell energies as per, say, Fig. 4, i.e., from the non-interacting gas; this would allow one to make the interacting BCS-side results closer to the thermodynamic limit. We’ve already seen in Fig. 5 that the closed-shell interacting results bear this interpretation out. Something analogous could be done when trying to extract a thermodynamic limit value for the pairing gap, in which case Fig. 6 would be relevant. Of course, that would still leave open the question of how to handle finite-size effects away from the deep BCS regime. Instead, we have chosen to approximate the thermodynamic-limit value, throughout the crossover, as the average of our best two sets of points, namely $N=50$ and $N=58$ and their neighbors. Beyond simply following from Fig. 6, this is an attempt to estimate finite-size effects more generally. Our overall DMC thermodynamic-limit for the pairing gap is shown in Fig. 9. We can see that the energy of the gap sharply increases as the BEC regime is approached and correspondingly a much reduced pairing gap in the BCS regime.

Another choice of dependent variable, which is commonly made in the literature, is to show $\Delta_{gap}+e_{b}/2$ vs the interaction strength; this is the same plot as in Fig. 2 and effectively produces a “zoomed-in” version of what’s going on. In Fig. 10 we compare our new results to the mean-field BCS value from Eq. (6). Overall, we find a consistent suppression with respect to the mean-field BCS result. Crucially, this is different from our pairing gap values in Ref. Galea:2016 (or those of Ref. Bertaina:2011 ), which tended to become larger than the mean-field value on the BCS side. In Ref. Galea:2016 , our $N=26$ DMC pairing gaps were explicitly described as “finite-size uncorrected” (the same could be said about the pioneering calculations of Ref. Bertaina:2011 ): in the present work we carried out considerably more DMC calculations to approximate the thermodynamic limit better. Our consistent suppression is similar to what one finds using the theory of Gorkov and Melik-Barkhudarov Gorkov:1961 ; Petrov:2003 : in two dimensions this gives suppression by a factor of $e$ with respect to the mean-field BCS result. We also compare with the values produced in a recent work which employs imaginary-time Green’s functions and analytic continuation Vitali:2017 . For the intermediate region, where these results were produced, the two sets of points are qualitatively similar, though the agreement is not perfect.

VI Conclusion

In summary, we have discussed pairing in two dimensions employing a number of formalisms. We started from mean-field BCS theory, which was solved (for the first time) for a finite interaction range; this confirmed that the effective ranges employed in the rest of this work were sufficiently short so as not to impact the final extraction of pairing gaps. We then proceeded to discuss finite-size effects in the strongly interacting crossover, and their impact on the extraction of the pairing gap via odd-even staggering. More specifically, we used DMC to determine the ground-state energy at various interaction strengths within the BEC-BCS crossover. We then determined the energy pairing gap for closed-shell values. As part of this process, we carried out an investigation of the effect of finite simulation sizes at various interaction strengths. We saw that the finite-size effects closer to the BCS regime closely followed the trend predicted by the non-interacting gas model, and deep into the BEC regime as the attractive interaction became increasingly strong-finite size effects decreased in relative magnitude and followed less closely the trend of the non-interacting gas. The main result is that as we go toward the BCS regime we find a suppression with respect to the mean-field BCS pairing gap prediction, qualitatively similar to that found using the theory of Gorkov and Melik-Barkhudarov.

Acknowledgements.

The authors would like to thank F. Diakonos, S. Giorgini, G. Palkanoglou, E. Vitali, and S. Zhang for several discussions. This work was supported in part by the Natural Sciences and Engineering Research Council (NSERC) of Canada, the Canada Foundation for Innovation (CFI), and the Early Researcher Award (ERA) program of the Ontario Ministry of Research, Innovation and Science. Computational resources were provided by SHARCNET and NERSC.