Theory of strongly paired fermions with arbitrary short-range interactions

TL;DR

This paper develops an effective field theory for strongly interacting fermions with short-range attractions, incorporating Gaussian pair fluctuations to predict equations of state and non-universal energy dependence, validated against quantum Monte Carlo results.

Contribution

It extends existing theories by including non-trivial phase shift effects and Gaussian fluctuations, providing a more accurate description of strongly paired fermion systems.

Findings

Good agreement with quantum Monte Carlo simulations on effective-range dependence.

The theory predicts non-universal ground-state energy variations.

Proposes experimental probing via dark-state optical control of Feshbach resonances.

Abstract

We develop an effective field theory to describe the superfluid pairing in strongly interacting fermions with arbitrary short-range attractions, by extending Kaplan's idea of coupling fermions to a fictitious boson state in Nucl. Phys. B \textbf{494}, 471 (1997). This boson field is assigned with unconventional kinetic term to recover the exact scattering phase shift obtained either from scattering data or model calculations. The theory works even if the explicit form of the interaction potential has not been constructed from scattering data. The contact boson-fermion coupling allows us to go beyond mean-field to include Gaussian pair fluctuations, yielding reliable predictions on equations of state. As an application, we use our theory to address the non-univerisal ground-state energy of strongly paired fermions, due to the non-trivial momentum dependence of the phase shift characterized, for example, by effective range. We find a good agreement between our predictions and recent quantum Monte Carlo simulations on the effective-range dependence in both three and two spatial dimensions. We propose that in cold-atom experiments, the non-universal dependence in thermodynamics can be probed using dark-state optical control of Feshbach resonances.

Equations

$${\Delta }_{\mathbf{k}}=-\sum _{{\mathbf{k}}^{\prime }}{V}_{\mathbf{k}{\mathbf{k}}^{\prime }}\frac{{\Delta }_{{\mathbf{k}}^{\prime }}}{2E_{{\mathbf{k}}^{\prime }}},$$

$$\mathcal{A}(E)=\frac{4\pi }{m}\frac{1}{k\cot \delta (k)-ik}.$$

$$k\cot \delta (k)=-\frac{4\pi }{m}\sum _{n=0}^{\infty }{c}_n{E}^n.$$

$${\mathcal{L}}_{\text{eff}}=\sum _{\sigma =\uparrow ,\downarrow }{\psi }_{\sigma }^{†}{\hat{K}}_{\mathrm{F}}{\psi }_{\sigma }^{}+{\varphi }^{†}\mathcal{F}({\hat{K}}_{\mathrm{B}})\varphi -\left({\varphi }^{†}{\psi }_{\downarrow }{\psi }_{\uparrow }+\text{h.c.}\right),$$

$$\mathcal{A}(E)$$

$$\mathcal{B}(E)=\sum _{\mathbf{p}}\frac{1}{E+iϵ-2{\epsilon }_{\mathbf{p}}},$$

$$\mathcal{B}(E)=-\frac{m}{4\pi }ik+\mathcal{D}(\Lambda ),$$

$$\mathcal{D}(\Lambda )=\frac{m\Lambda }{2{\pi }^2}=\sum _{\mathbf{p}}\frac{1}{2{\epsilon }_{\mathbf{p}}}.$$

$$k\cot \delta (k)=-\frac{4\pi }{m}{\mathcal{F}}_{\mathrm{R}}(E).$$

$${\mathcal{F}}_{\mathrm{R}}({\hat{K}}_{\mathrm{B}})=\sum _{n=0}^{\infty }{c}_n({\hat{K}}_{\mathrm{B}}{)}^n.$$

$$\mathcal{A}(E)=\frac{4\pi }{m}\frac{1}{\pi \cot \delta (k)-i\pi }.$$

$$\pi \cot \delta (k)=\ln \left(\frac{E}{{\epsilon }_{2D}}\right)-\frac{4\pi }{m}\sum _{n=1}^{\infty }{c}_n{E}^n,$$

$$\mathcal{B}(E)=-\frac{m}{4\pi }\left[\ln \left(\frac{{\Lambda }^2}{mE}\right)+i\pi \right],$$

$$\pi \cot \delta (k)=\ln \left(\frac{E}{{\epsilon }_0}\right)-\frac{4\pi }{m}{\mathcal{F}}_{\mathrm{R}}(E),$$

$$\mathcal{D}(\Lambda )=\frac{m}{4\pi }\ln \left(\frac{{\Lambda }^2}{m{\epsilon }_0}\right)=\sum _{\mathbf{p}}\frac{1}{2{\epsilon }_{\mathbf{p}}+{\epsilon }_0}.$$

$${\mathcal{F}}_{\mathrm{R}}({\hat{K}}_{\mathrm{B}})=\sum _{n=1}^{\infty }{c}_n({\hat{K}}_{\mathrm{B}}{)}^n.$$

$$\mathcal{Z}=\int [d\psi ][d{\psi }^{†}][d\varphi ][d{\varphi }^{†}]\exp \left[\int dx{\mathcal{L}}_{\text{eff}}\right],$$

$$\displaystyle{\cal Z}=\int[d\phi][d\phi^{\dagger}]\exp{\Big{\{}-{\cal S}_{{\rm eff}}[\phi,\phi^{\dagger}]\Big{\}}},$$

$${\mathcal{S}}_{\mathrm{eff}}=-\int dx{\varphi }^{†}\mathcal{F}({\hat{K}}_{\mathrm{B}}){\varphi }^{}-\mathrm{Tr}\ln M_{\mathrm{F}}[\varphi ,{\varphi }^{†}],$$

$$\displaystyle M_{\rm F}[\phi,\phi^{\dagger}]=\left(\begin{array}[]{cc}\hat{K}_{{\rm F}}&\phi\\ \phi^{\dagger}&-\hat{K}_{{\rm F}}^{*}\end{array}\right).$$

$$\mathcal{Z}=\int [d\varphi ][d{\varphi }^{†}]\det M_{\mathrm{F}}\exp \left[\int dx{\varphi }^{†}\mathcal{F}({\hat{K}}_{\mathrm{B}}){\varphi }^{}\right],$$

$${\Omega }_{\mathrm{MF}}=\sum _{\mathbf{k}}\left({\xi }_{\mathbf{k}}-E_{\mathbf{k}}+\frac{|\Delta {|}^2}{2{\epsilon }_{\mathbf{k}}}\right)-|\Delta {|}^2{\mathcal{F}}_{\mathrm{R}}(2\mu ).$$

$$\sum _{\mathbf{k}}\left(\frac{1}{2E_{\mathbf{k}}}-\frac{1}{2{\epsilon }_{\mathbf{k}}}\right)=-{\mathcal{F}}_{\mathrm{R}}(2\mu ).$$

$${\mathcal{A}}^{-1}(2\mu =-{\epsilon }_{\mathrm{B}})=0,$$

$$\Delta \simeq \frac{8}{e^2}{\epsilon }_{\mathrm{F}}\exp \left(\frac{\pi }{2k_{\mathrm{F}}a}-\frac{\pi }{4}{k}_{\mathrm{F}}{r}_{\mathrm{e}}\right),$$

$${\mathcal{S}}_{\mathrm{eff}}=\beta V{\Omega }_{\mathrm{MF}}+{\mathcal{S}}_{\mathrm{FL}}\left[\phi ,{\phi }^{†}\right]$$

$$\mathcal{Z}=e^{-\beta V{\Omega }_{\mathrm{MF}}}\int \left[d\phi \right]\left[d{\phi }^{†}\right]e^{-{\mathcal{S}}_{\mathrm{FL}}}.$$

$${\mathcal{S}}_{\mathrm{GF}}=\frac{1}{2}\sum _Q{\Phi }^{†}(Q)\mathbf{M}(Q)\Phi (Q),$$

$${\mathbf{M}}_{11}(Q)=\sum _{\mathbf{k}}\left(\frac{u_{+}^2{u}_{-}^2}{Y_{-}}-\frac{{\upsilon }_{+}^2{\upsilon }_{-}^2}{Y_{+}}+\frac{1}{2{\epsilon }_{\mathbf{k}}}\right)-{\mathcal{F}}_{\mathrm{R}}\left(Z\right),$$

$${\mathbf{M}}_{12}(Q)=\sum _{\mathbf{k}}{u}_{+}{\upsilon }_{+}{u}_{-}{\upsilon }_{-}\left(\frac{1}{Y_{+}}-\frac{1}{Y_{-}}\right),$$

$${\Omega }_{\mathrm{GF}}=-\sum _{\mathbf{q}}{\int }_{-\infty }^{\infty }\frac{d\omega }{\pi }\frac{1}{e^{\beta \omega }-1}{\delta }_{\mathrm{M}}(\omega ,\mathbf{q}),$$

Full text

I Introduction

In quantum many-body Fermi systems, attractive inter-particle interaction leads to Cooper pairing and superfluidity Review01 ; Review02 ; Review03 ; Exciton . Increasing the attraction strength will induce a crossover from a Bardeen-Cooper-Schrieffer (BCS) superfluid state with largely overlapping Cooper pairs to a Bose-Einstein condensate (BEC) of tightly bound molecules Eagles ; Leggett ; NSR ; BCSBEC1 ; BCSBEC2 ; BCSBEC3 ; BCSBEC4 . Ultracold atomic Fermi gases near magnetic-field-tuned Feshbach resonances provide clean systems to demonstrate the BCS-BEC crossover Regal2004 ; Zwierlein2004 ; Kinast2004 and explore novel many-body phenomena such as universal thermodynamics Ho2004 ; Hu2007 ; Thomas2009 ; Nascimbene2010 ; Navon2010 ; Horikoshi2010 ; Ku2012 .

The universal properties of an ultracold Fermi gas stem from its simple form of the scattering phase shift in the dilute limit, characterized by a large $s$-wave scattering length and negligible effective range TH01 ; TH02 ; TH03 ; TH04 ; TH05 ; TH06 ; TH07 ; TH08 ; TH09 ; TH10 ; TH11 ; TH12 ; EOSmc1 ; EOSmc2 ; EOSmc3 ; EOSmc4 . In realistic systems, however, non-universal effects could be important due to the non-trivial momentum dependence in the phase shift, which leads to, for example, nonzero effective range. For instance, neutron superfluid with a large $s$-wave scattering length $a_{{\rm nn}}\simeq-18.5$ fm and a sizable effective range $r_{{\rm nn}}\simeq 2.7$ fm may exist in the inner crust of neutron stars Schwenk2005 ; Gezerlis2010 ; Wyk2018 . Certainly, it is of great interest to use cold atoms to simulate and understand any non-universal properties associated with realistic short-range interaction potential, particularly in thermodynamics.

In this work, we aim to establish a genuine and elegant theory to describe strongly paired fermions with arbitrary short-range interaction $V(r)$. The conventional way to handle this problem is technically difficult. For example, within mean-field theory one needs to solve the gap equation,

[TABLE]

where $\Delta_{{\bf k}}$ is the gap function, $V_{{\bf k}{\bf k}^{\prime}}$ is the Fourier transform of the interaction potential, and $E_{{\bf k}}$ is the BCS-type single-particle dispersion. Even at the mean-field level, this integral equation is not easy to solve MRE-01 ; MRE-02 ; MRE-03 . Apart from numerically expensive quantum Monte Carlo (QMC) simulations, going beyond mean-field seems impossible, since $V_{{\bf k}{\bf k}^{\prime}}$ is generally not separable. Here, our strategy is to develop an effective field theory following the pioneering work by Kaplan Kaplan01 , who introduced a fictitious boson state and coupled it to fermions via a contact interaction. In this way, the effective-range expansion of the scattering phase shift is recovered in the two-body limit Kaplan01 ; Kaplan02 . Interestingly, the boson state introduced by Kaplan is no longer fictitious with the recent realization of magnetic-field-tuned Feshbach resonances Chin2010 : it is a real dimer state in the closed channel.

The key advantage of our effective field theory is that it is constructed to precisely reproduce the full two-body scattering phase shift of the potential $V(r)$ of interest, which is assumed to be known, either directly from scattering data (i.e., for nucleon superfluids Gezerlis2010 ) or from model calculations (i.e., for quasi-two-dimensional gases Petrov2001 ). Thus, all the information of the interaction potential is retained, beyond the effective-range expansion adopted earlier Tajima01 ; Tajima02 ; Hu-2D-01 ; Hu-2D-02 . The contact boson-fermion interaction then allows us to include crucial quantum fluctuations beyond mean-field and hence to provide a reliable description of strongly paired fermions.

As a simple application, we predict the effective-range dependence of the ground-state energy near $s$-wave resonances in both three (3D) and two dimensions (2D). Our results are in good agreement with existing QMC calculations. To demonstrate the potential of using cold atoms to understand the non-universal properties of strongly paired fermions due to the effective-range effect and beyond, we propose dark-state optical control of Feshbach resonances, which leads to a non-trivial momentum dependence of the phase shift and consequently non-universal thermodynamics. Our results pave the way to using cold atoms to simulate realistic many-body Fermi systems, which exist ubiquitously in all fields of physics.

The paper is organized as follows. In Sec. II, we construct a general effective Lagrangian for an arbitrary short-range interaction. In Sec. III, we formulate the many-body theory of strongly paired fermionic superfluids based on the effective Lagrangian. The effective-range dependence of strongly paired fermions is studied in Sec. IV, and an experimental scheme to probe the non-universal thermodynamics is proposed in Sec. V. We summarize in Sec. VI.

II Effective Lagrangian

For the sake of simplicity, we assume that the short-range two-body interaction $V(r)$ is around an $s$-wave resonance, and thus we neglect the contributions from higher partial waves. The $s$-wave scattering amplitude ${\cal A}(E)$ can be expressed in terms of the $s$-wave scattering phase shift $\delta(k)$, where $E=k^{2}/m$ is the scattering energy and $m$ is the mass of fermions ($\hbar=1$ hereafter). In 3D, we have

[TABLE]

For a short-range interaction, $k\cot\delta(k)$ is an analytical function of $E$, leading to the expansion

[TABLE]

Truncating to the first two terms, we obtain the so-called effective range expansion, with $c_{0}=m/(4\pi a)$ and $c_{1}=-m^{2}r_{{\rm e}}/(8\pi)$, where $a$ and $r_{{\rm e}}$ are the scattering length and effective range, respectively.

A low-energy effective Lagrangian including only fermion fields $\psi_{\sigma}^{\phantom{{\dagger}}}(x)$ ($\sigma=\uparrow,\downarrow$) can be constructed order by order according to the expansion (3) Kaplan03 , where $x=(t,{\bf r})$ with $t$ being the time and ${\bf r}$ the spatial coordinates. Here we aim to construct an effective Lagrangian that recovers exactly the scattering phase shift $\delta(k)$. Following Kaplan Kaplan01 ; Kaplan02 , we introduce a boson field $\phi(x)$ that couples to the fermions. Because of the Galilean invariance, the general effective Lagrangian takes the form

[TABLE]

where $\hat{K}_{{\rm F}}=i\partial_{t}+\mu+\nabla^{2}/(2m)$ and $\hat{K}_{{\rm B}}=i\partial_{t}+2\mu+\nabla^{2}/(4m)$ are the Galilean invariant kinetic operators for fermion and boson, respectively, with $\mu$ being the chemical potential of fermions. The crucial point of our construction is that the boson field has an unconventional kinetic term, represented by the function ${\cal F}(\hat{K}_{{\rm B}})$, which can be designed to recover precisely the phase shift $\delta(k)$.

To see this, let us recall that the scattering amplitude ${\cal A}(E)$ is given by the ladder summation with an interaction vertex ${\cal F}^{-1}(E)$, which gives

[TABLE]

Here, the two-particle bubble diagram reads

[TABLE]

with $\varepsilon_{{\bf p}}={\bf p}^{2}/(2m)$ and $\sum_{\bf p}=\int d^{3}{\bf p}/(2\pi)^{3}$ for 3D. The ultraviolet (UV) divergence in ${\cal B}$ can be regularized via a hard cutoff $\Lambda$ for $|{\bf p}|$, leading to

[TABLE]

where the divergent part is given by

[TABLE]

Thus, the scattering amplitude takes the same form of Eq. (2), with the phase shift given by

[TABLE]

Here ${\cal F}_{{\rm R}}(E)={\cal F}(E)+{\cal D}(\Lambda)$ is the renormalized version of the ${\cal F}$-function. It is readily seen that the expansion Eq. (3) allows us to determine the ${\cal F}$-function as a polynomial in $\hat{K}_{{\rm B}}$, i.e.,

[TABLE]

Since all terms in Eq. (3) are included, our construction of the effective Lagrangian Eq. (4) is valid beyond the radius of convergence of the effective-range expansion.

The effective Lagrangian (4) also applies to 2D systems, where $\nabla^{2}=\partial_{x}^{2}+\partial_{y}^{2}$. In 2D, the the $s$-wave scattering amplitude is given by scatter2D

[TABLE]

For a short-range interaction, we have

[TABLE]

where $\varepsilon_{\rm 2D}=1/(ma_{\rm 2D}^{2})$, with $a_{\rm 2D}$ being the 2D scattering length. The 2D effective range can be defined as $R_{\rm 2D}=-4\pi c_{1}/m^{2}$, which has units of length2. Even though the leading term is nonanalytical in $E$, it is purely from the two-particle bubble diagram ${\cal B}(E)$. Using the same cutoff regularization we obtain

[TABLE]

Direct ladder summation shows that

[TABLE]

where the the renormalized ${\cal F}$-function reads ${\cal F}_{\rm R}(E)={\cal F}(E)+{\cal D}(\Lambda)$, with the counter term

[TABLE]

The energy scale $\varepsilon_{0}$ can be chosen arbitrarily, and we set $\varepsilon_{0}=\varepsilon_{\rm 2D}$ for convenience. Thus in 2D, the ${\cal F}$-function is given by

[TABLE]

III Many-Body Theory

We consider a Fermi gas with a short-range two-body interaction $V(r)$. Now we can solve the many-body problem based on the effective Lagrangian (4), instead of directly using the interaction potential $V(r)$. In the imaginary time path integral formalism, the partition function ${\cal Z}$ is given by

[TABLE]

where $x=(\tau,{\bf r})$ and $\int dx=\int_{0}^{\beta}d\tau\int d{\bf r}$ after the replacement $t\rightarrow-i\tau$, with $\tau$ being the imaginary time and $\beta=1/(k_{B}T)$ the inverse temperature. The fermions can be directly integrated out and we obtain

[TABLE]

where the bosonic effective action reads

[TABLE]

with the fermion matrix

[TABLE]

The partition function can be alternatively expressed as

[TABLE]

Thus the many-body problem can be simulated using the Lattice Monte Carlo method since the fermion determinant $\det M_{\rm F}$ is positive, which has been applied to the zero-range interaction case Lattice01 .

In this work, we aim to develop an analytical theory. One advantage of the effective Lagrangian (4) is that the saddle point or classical part of the boson field $\phi$, directly serves as the superfluid order parameter. The mean-field theory amounts to searching for the static and uniform saddle-point solution $\phi(x)=\Delta$ that minimizes the effective action ${\cal S}_{{\rm eff}}$. In 3D and at zero temperature, the mean-field contribution to the grand potential $\Omega=-(T/V)\ln\mathcal{Z}$ can be evaluated as

[TABLE]

Here $\xi_{{\bf k}}=\varepsilon_{{\bf k}}-\mu$ and $E_{{\bf k}}=(\xi_{{\bf k}}^{2}+|\Delta|^{2})^{1/2}$. Without loss of generality, we can set $\Delta$ to be real and positive. At $T=0$, the gap equation determining $\Delta(\mu)$ is given by

[TABLE]

In the strong attraction limit, the system forms a BEC of tightly bound dimers and we have $\Delta\ll|\mu|$. The gap equation thus reduces to a two-body equation determining the negative-energy pole of the scattering amplitude,

[TABLE]

where $\varepsilon_{{\rm B}}$ is precisely the binding energy of the dimer state determined by solving the Schroedinger equation with the interaction potential $V(r)$. In the weak attraction limit, where $\mu\simeq\varepsilon_{{\rm F}}$, with $\varepsilon_{{\rm F}}=k_{{\rm F}}^{2}/(2m)$ being the Fermi energy, the gap equation provides a reasonable effective-range dependence of the pairing gap. By approximating the ${\cal F}$-function as ${\cal F}_{{\rm R}}(E)\simeq c_{0}+c_{1}E$, the pairing gap reads

[TABLE]

indicating that a positive (negative) effective range suppresses (enhances) the pairing gap.

The mean-field theory is only qualitatively correct for strongly paired fermions. To have a more quantitative description, we consider quantum fluctuations around the saddle point by writing $\phi(x)=\Delta+\varphi(x)$. The effective action can be expressed as

[TABLE]

and the partition function becomes

[TABLE]

An exact analytical treatment of the fluctuation contribution ${\cal S}_{{\rm FL}}$ is impossible. Here we consider only the Gaussian fluctuations, i.e., the contributions that are quadratic in $\varphi(x)$ and $\varphi^{\dagger}(x)$, corresponding to the contributions from collective modes. In the momentum space, this Gaussian fluctuation contribution, ${\cal S}_{{\rm GF}}$, can be expressed as

[TABLE]

where $\Phi(Q)=[\varphi(Q),\varphi^{\dagger}(-Q)]^{{\rm T}}$. The inverse Green’s function of collective bosonic modes, ${\bf M}(Q)$, is a $2\times 2$ matrix, with elements satisfying the relations ${\bf M}_{11}(Q)={\bf M}_{22}(-Q)$ and ${\bf M}_{12}(Q)={\bf M}_{21}(Q)$. Here and in the following, we use the notation $Q=(iq_{l},{\bf q})$, with $q_{l}=2\pi lT$ ($l\in\mathbb{Z}$) being the bosonic Matsubara frequency. At $T=0$, the elements of ${\bf M}$ can be explicitly evaluated as

[TABLE]

where $Z=iq_{l}+2\mu-{\bf q}^{2}/(4m)$, $Y_{\pm}=iq_{l}\pm(E_{+}+E_{-})$, and the BCS distribution functions are defined as $u_{{\bf k}}^{2}=1-v_{{\bf k}}^{2}=(1+\xi_{{\bf k}}/E_{{\bf k}})/2$. The plus and minus signs denote the momenta ${\bf k}+{\bf q}/2$ and ${\bf k}-{\bf q}/2$, respectively.

Within the Gaussian pair fluctuation (GPF) approximation TH03 ; TH04 ; TH05 ; TH06 , i.e., ${\cal S}_{{\rm FL}}\simeq{\cal S}_{{\rm GF}}$, the path integral over the fluctuations can be carried out, and the grand potential is given by $\Omega\simeq\Omega_{{\rm MF}}+\Omega_{{\rm GF}}$, where the contribution from Gaussian fluctuations reads

[TABLE]

and the phase shift $\delta_{{\rm M}}(\omega,{\bf q})=-{\rm Im}{\cal W}(\omega+i\epsilon,{\bf q})$, with the ${\cal W}$-function given by

[TABLE]

The grand potential $\Omega(\mu)$ in the GPF theory can be determined, by solving $\Delta(\mu)$ from the gap equation (25). The density equation of state is then calculated using $n(\mu)=-\partial\Omega(\mu)/\partial\mu$.

It is worth noting that as the full scattering phase-shift is reproduced by our theory in the two-body limit, we recover correctly the virial expansion of the equation of state at high temperature, i.e., the Beth-Uhlenbeck formalism for the second-order virial coefficient can be derived. At high temperature, the system is a normal gas with large but negative $\mu$. The equation of state can be expanded in powers of the fugacity $z=e^{\beta\mu}\ll 1$ Liu2013 . In 3D we have

[TABLE]

where $\lambda_{T}=\sqrt{2\pi\beta/m}$ is the thermal wavelength. Our theory recovers the correct virial equation of state up to the order $O(z^{2})$. To see this, we write $b_{2}=b_{2}^{(1)}+b_{2}^{(2)}$, where the one-body contribution $b_{2}^{(1)}=-2^{-5/2}$. To find the two-body contribution $b_{2}^{(2)}$, it is sufficient to use Eq. (32), with the phase shift $\delta_{{\rm M}}(\omega,{\bf q})$ replaced by the two-body one $\delta_{{\rm 2B}}(\omega,{\bf q})$. Using a new variable $E=\omega+2\mu-{\bf q}^{2}/(4m)$, we obtain

[TABLE]

where $\delta_{{\rm 2B}}(E)=-{\rm Im}\ln[{\cal A}^{-1}(E)]$. Using Eq. (2), we recover the elegant Beth-Uhlenbeck formalism,

[TABLE]

where $\varepsilon_{{\rm B}}$ is the exact binding energy and $\delta(k)$ is the exact phase shift.

The many-body theory for 2D is quite similar to the 3D case. Due to the energy scale $\varepsilon_{0}=\varepsilon_{{\rm 2D}}$, the mean-field thermodynamics and gap equation are modified to,

[TABLE]

and

[TABLE]

respectively. Moreover, for the Gaussian fluctuations, the matrix element $\mathbf{M}_{11}(Q)$ takes the form,

[TABLE]

All the notations, i.e., $u_{\pm}^{2}$, $v_{\pm}^{2}$ and $Y_{\pm}$, are the same as in the 3D case.

The consistency and validity of the GPF theory has been studied for both 3D TH03 ; TH04 and 2D TH05 ; TH06 cases with zero-range interactions. The truncation of the pair fluctuations at the Gaussian level provides a quantitatively good description of the BCS-BEC crossover at $T=0$, since the most important fluctuation contribution, the Goldstone mode fluctuation, is taken into account properly. The missing Fermi liquid correction in the BCS limit of the mean-field theory is naturally recovered by the GPF contribution TH04 . In the BEC limit, the GPF contribution is significant to give a quantitatively good boson-boson interaction TH03 ; TH04 . In particular, in 2D, the boson-boson interaction is missing in the mean-field theory, leading to a qualitatively incorrect equation of state in the BCS-BEC crossover. The correct boson-boson interaction is naturally recovered by the GPF contribution, leading to a correct equation of state TH05 . In this work, the GPF theory has been generalized to finite-range interactions. We will see that the GPF theory also provides a better description of the effective-range dependence than the mean-field theory.

IV Effective-Range Dependence

We now consider the zero-temperature equation of state of a Fermi gas with fixed density $n=k_{{\rm F}}^{3}/(3\pi^{2})$ in 3D. We focus on the effective-range dependence of the ground-state energy at resonance, where the ${\cal F}$-function is approximated as

[TABLE]

While the theory in Sec. III applies for both positive and negative effective ranges, here we focus on a negative effective range since it is relevant to cold atom systems and the computational cost is small (see the Appendix). For a negative effective range, $r_{{\rm e}}<0$, this truncation is equivalent to the two-channel model description of the Feshbach resonance BCSBEC4 . At large negative effective range, the model can be treated perturbatively according to a small parameter $(k_{{\rm F}}r_{{\rm e}})^{-1}$ BCSBEC4 . The mean-field theory provides an accurate description for $k_{{\rm F}}r_{{\rm e}}\rightarrow-\infty$. As shown in Fig. 1(a), we find that the mean-field and the GPF results converge at large $k_{{\rm F}}|r_{{\rm e}}|$, as anticipated. At small and moderate effective range, the GPF result agrees well with the QMC data Conduit3D . For small $k_{{\rm F}}r_{{\rm e}}$, the ground-state energy at resonance can be expressed as

[TABLE]

where the Bertsch parameter reads $\xi=0.591$ in the mean-field theory and $\xi=0.401$ in the GPF theory. The GPF result agrees well with the latest experimental Thomas2009 ; Nascimbene2010 ; Navon2010 ; Horikoshi2010 ; Ku2012 and QMC EOSmc1 ; EOSmc2 ; EOSmc3 ; EOSmc4 ; Conduit3D results, which lie in the range $0.36-0.42$. The coefficient $\zeta$ can also be determined. It reads $\zeta=0.273$ in the mean-field theory and $\zeta=0.105$ in the GPF theory. We note that our GPF result $\zeta=0.105$ is in good agreement with the result $\zeta=0.087(1)$ Conduit3D or $\zeta=0.127(4)$ Forbes2012 from the diffusion QMC and $\zeta=0.11(3)$ EOSmc4 from the auxiliary-field QMC.

In the BEC limit, the effective interaction between two composite dimers may be deduced. As $\mu\rightarrow-\infty$ and $\Delta\ll|\mu|$, we expand the gap equation (10) in powers of $\Delta/|\mu|$ and obtain,

[TABLE]

The solution can be expressed as $2\mu=-\varepsilon_{{\rm B}}+\mu_{{\rm B}}$, with $\mu_{{\rm B}}\ll\varepsilon_{{\rm B}}$ being the dimer chemical potential. We find

[TABLE]

where ${\cal C}=1+8\pi\sqrt{m\varepsilon_{{\rm B}}}{\cal F}_{{\rm R}}^{\prime}(-\varepsilon_{{\rm B}})/m^{2}$, with ${\cal F}_{{\rm R}}^{\prime}(x)=\partial{\cal F}_{{\rm R}}(x)/\partial x$. Meanwhile, the number equation becomes

[TABLE]

where $\alpha\sim O(1)$ comes from the Gaussian-fluctuation contribution and depends on the details of the ${\cal F}$-function. Thus, we recover the Bogoliubov equation of state for weakly interacting bosons, $\mu_{{\rm B}}=4\pi a_{{\rm dd}}n_{{\rm B}}/m_{{\rm B}}$, where $m_{{\rm B}}=2m$ is the mass of the dimers and $n_{{\rm B}}=n/2$ is the density of the dimers. The dimer-dimer scattering length $a_{{\rm dd}}$ is then given by

[TABLE]

Note that the above discussions are valid only for the case ${\cal C}>0$, i.e., for a repulsive dimer-dimer interaction.

For zero-range interaction, we have $a_{{\rm dd}}(0)=2a$ from the mean-field theory and $a_{{\rm dd}}(0)\simeq 0.57a$ from the GPF theory. The GPF result is close to the exact result $a_{{\rm dd}}(0)=0.6a$ from solving the four-body problem Petrov2004 . Considering only the effective range effect, i.e., ${\cal F}_{{\rm R}}(E)\simeq c_{0}+c_{1}E$, and neglecting the weak dependence of $\alpha$ on $r_{\textrm{e}}$, we obtain

[TABLE]

For a large negative effective range, $|r_{{\rm e}}|\gg a$, we have $a_{{\rm dd}}\simeq a^{2}/(2|r_{{\rm e}}|)\ll a$. For a positive effective range, $a_{{\rm dd}}$ is enhanced. The divergence at $r_{{\rm e}}=a/2$ is artificial due to our simple truncation to effective range and is likely cured by the inclusion of the shape term $O(E^{2})$. We also observe that the quantity ${\cal C}$ vanishes at $r_{{\rm e}}=a/2$, indicating that the dimer-dimer interaction turns to be attractive for a larger effective range. The qualitative change around $r_{{\rm e}}\sim a$ indicates that the present analysis fails and the ground state in this regime remains to be explored. Recent few-body calculation shows that two dimers may form a cluster state for $r_{{\rm e}}>0.46a$ Yin2019 .

We also calculated the ground-state energy of a 2D Fermi gas with fixed density $n=k_{{\rm F}}^{2}/(2\pi)$ as a function of the effective-range parameter $k_{{\rm F}}^{2}R_{{\rm 2D}}$. In Fig. 1(b), we show the energy in the strongly interacting regime where the 2D scattering length $a_{{\rm 2D}}$ is determined by requiring $\mu=0$ within mean-field theory. If the effective-range parameter is not large, our GPF result shows an excellent agreement with the QMC result Conduit2D . However, both in 2D and 3D, our GPF predictions with the simple truncation ${\cal F}_{{\rm R}}(E)\simeq c_{0}+c_{1}E$ notably deviate from the QMC results at high density or a large effective range, indicating that the higher-order contributions beyond the effective range expansion may become important. These corrections could depend sensitively on the model potentials used in QMC simulations.

V Probing Non-universal Thermodynamics

In cold atom experiments, the $s$-wave scattering length $a$ is tuned by magnetic field Timmermans1999 ; Chin2010 and the effective range is

[TABLE]

where $a_{{\rm bg}}$ is the open-channel background scattering length, $\gamma$ is the difference of the magnetic moment between the open and the closed channels, and $B_{\Delta}$ is the resonance width. For broad resonances in experiments, such as 6Li at $832$G and 40K at $202$G, the effective range $k_{{\rm F}}r_{0}$ is negligible.

Here we consider dark-state optical control of the Feshbach resonance Thomas2012 ; Thomas2016 ; Thomas2018 , where two laser beams couple the molecular state $|1\rangle$ responsible for the resonance and another molecular state $|2\rangle$ in the closed channel to an excited molecular state $|{\rm e}\rangle$. Near the resonance, the scattering phase shift is modified to Thomas2012 ; Zhang2017 ; He2018

[TABLE]

where $\Sigma(E)$ is the optically induced Stark shift. In the dark-state regime, we have $\Sigma(0)=0$ and $\Sigma^{\prime}(0)=-\Omega_{1}^{2}/\Omega_{2}^{2}$, where $\Omega_{1}$ and $\Omega_{2}$ are the Rabi frequencies for the transitions $|1\rangle\leftrightarrow|{\rm e}\rangle$ and $|2\rangle\leftrightarrow|{\rm e}\rangle$. Thus the resonance does not shift but the effective range can be tuned by changing the ratio $\Omega_{1}/\Omega_{2}$ Thomas2012 :

[TABLE]

Figure 2 shows the ground-state energy as a function of $k_{{\rm F}}r_{{\rm e}}$ for 6Li and 40K atom gases, at resonances $832$ G and $202$ G, respectively. The density $n$ is chosen as $2\times 10^{14}$ cm*-3*, a typical value realized in cold atom experiments. The effective range parameter without optical control is $k_{{\rm F}}r_{0}=-4.5\times 10^{-4}$ for 6Li and $k_{{\rm F}}r_{0}=-2.7\times 10^{-2}$ for 40K. We find that the result from a 40K atom gas agrees well with the effective range expansion. This is because the resonance for 40K has a relatively large intrinsic effective range $r_{0}$. To reach $k_{{\rm F}}r_{{\rm e}}\sim O(1)$, we need $\Omega_{1}/\Omega_{2}\sim 10$ for 40K and a much larger value $\Omega_{1}/\Omega_{2}\sim 100$ for 6Li. As a result, for 6Li the higher-order terms beyond the effective range expansion of Eq. (48) become important.

VI Summary

We have proposed a convenient way to describe strongly paired fermions without knowing the details of the short-range interaction potential. Within our effective field theory, it is easy to go beyond mean field to include quantum fluctuations for arbitrary potential. The predicted effective-range dependence of the ground-state energy is in good agreement with recent quantum Monte Carlo simulations. At a large effective range, the effect beyond the effective-range expansion may become significant. Our theory can be readily generalized to include the contributions from higher partial waves Yu2016 ; Yu2017 and extended to study bosonic systems You2017 ; Pan2019 . It is also interesting to study the case of a positive effective range, which is relevant to the equation of state at low density neutron matter of astrophysical interest Schwenk2005 ; Gezerlis2010 ; Wyk2018 .

Acknowledgements.

We thank Xiangyu Yin and Ren Bi for useful discussions. This research was supported by the National Natural Science Foundation of China, Grant Nos. 11775123 and 11890712 (L.H.), National Key Research and Development Program of China, Grant No. 2018YFA0306503 (L.H.), and Australian Research Council’s (ARC) Discovery Program, Grant No. FT140100003 (X.-J.L), Grant No. DP180102018 (X.-J.L), and Grant No. DP170104008 (H.H.).

Appendix A Evaluating the fluctuation contribution $\Omega_{\rm GF}$ at $T=0$

At $T=0$, there is a much easier approach to evaluate $\Omega_{\rm GF}$ for some cases proposed by Diener, Sensarma, and Randeria TH04 . We define two functions ${\bf M}^{\rm C}_{11}$ and ${\bf M}^{\rm C}_{22}$, which in 3D are given by

[TABLE]

and ${\bf M}^{\rm C}_{22}(Q)={\bf M}^{\rm C}_{11}(-Q)$. In 2D, we need to replace $1/(2\varepsilon_{\bf k})$ with $1/(2\varepsilon_{\bf k}+\varepsilon_{\rm 2D})$. At $T=0$, the boson Matsubara sum is converted to a continuous integral over an imaginary frequency $i\omega$, i.e.,

[TABLE]

Therefore, if the Matsubara sum $\sum_{l}\ln{\bf M}^{\rm C}_{11}(Q)$ vanishes even without the convergent factor $e^{iq_{l}0^{+}}$ at $T=0$, we arrive at an alternative convergent expression for $\Omega_{\rm GF}$ TH04 :

[TABLE]

Here we have used the fact that the integrand is even in $\omega$.

If we consider only the effective-range effect and approximate the ${\cal F}$-function as ${\cal F}_{{\rm R}}(E)\simeq c_{0}+c_{1}E$, the counter function in 3D can be expressed as

[TABLE]

For negative effective range ($r_{\rm e}<0$), the function ${\bf M}_{11}^{\rm C}(z,{\bf q})$ has no zeros in the left half plane (${\rm Re}z<0$) for an arbitrary value of ${\bf q}$. Therefore, the Matsubara sum $\sum_{l}\ln{\bf M}^{\rm C}_{11}(Q)$ vanishes at $T=0$ and the above trick applies. After some manipulations, we obtain

[TABLE]

The functions $A,B,C,D,$ and $F$ are defined as

[TABLE]

This expression for $\Omega_{\rm GF}$ leads to a rapidly convergent result. However, for positive effective range ($r_{\rm e}>0$), the function ${\bf M}_{11}^{\rm C}(z,{\bf q})$ has zeros in the left half plane for sufficiently large values of $|{\bf q}|$. In this case, the above trick fails and we need to use the phase shift expression (32). The convergence of the integral over the real frequency is rather slow and the computational cost becomes large. Similar discussions also apply to the 2D case.