Fermi-Bose mixture in mixed dimensions

TL;DR

This paper proposes a novel method to realize a topological p-wave superfluid in a 2D fermionic gas mediated by a 3D BEC, overcoming previous instability issues and predicting achievable critical temperatures.

Contribution

It introduces a phonon-mediated pairing mechanism in mixed dimensions, calculating the critical temperature with advanced diagrammatic corrections, and outlines an experimental implementation.

Findings

Critical temperature is significantly boosted by higher-order corrections.

The proposed scheme makes p-wave superfluidity experimentally feasible.

The approach overcomes instability issues of previous p-wave superfluid attempts.

Abstract

One of the challenging goals in the studies of many-body physics with ultracold atoms is the creation of a topological $p_{x} + ip_{y}$ superfluid for identical fermions in two dimensions (2D). The expectations of reaching the critical temperature $T_c$ through p-wave Feshbach resonance in spin-polarized fermionic gases have soon faded away because on approaching the resonance, the system becomes unstable due to inelastic-collision processes. Here, we consider an alternative scenario in which a single-component degenerate gas of fermions in 2D is paired via phonon-mediated interactions provided by a 3D BEC background. Within the weak-coupling regime, we calculate the critical temperature $T_c$ for the fermionic pair formation, using Bethe-Salpeter formalism, and show that it is significantly boosted by higher-order diagramatic terms, such as phonon dressing and vertex corrections. We describe in detail an experimental scheme to implement our proposal, and show that the long-sought p-wave superfluid is at reach with state-of-the-art experiments.

Tables

Table 1. Table 1: Parameters of 171,173 Yb- 7 Li mixture. The elastic scattering rate Γ el , B subscript Γ el B \Gamma_{\rm el,B} is given for thermal atoms at a temperature of T = T c p 𝑇 superscript subscript 𝑇 𝑐 𝑝 T=T_{c}^{p} colliding with a BEC at density n B subscript 𝑛 𝐵 n_{B} . Γ 3 − body , B = − N ˙ B / N B subscript Γ 3 body B subscript ˙ 𝑁 𝐵 subscript 𝑁 𝐵 \Gamma_{\rm 3-body,B}=-\dot{N}_{B}/N_{B} is the initial 3-body loss rate of the BEC Gross  et al. ( 2009 ); Pollack  et al. ( 2009a ) .

$n_B$	$6\times {10}^{14}$ atoms/cm3
$a_B$	$8a_0$
$a_{\mathrm{FB}}$	$200a_0$
$a_{\mathrm{eff}}$	$\sqrt{m_B/m_{FB}}{a}_{\mathrm{FB}}=204a_0$
$\alpha$	$v_F/c_s=1.5{\alpha }_{\max }=0.22$
${\gamma }_{BEC}$	$a_B{n}_B^{1/3}=0.004$
${\gamma }_{\mathrm{eff}}$	$a_{\mathrm{eff}}{n}_B^{1/3}=0.1$
$\xi$	$1/\sqrt{8\pi n_B{a}_B}=0.4\mu$m
$X$	$\xi k_F=\xi \sqrt{4\pi n_F}=3.8$
$v_F$	$\hslash k_F/m_F=0.4$cm/s
$c_s$	$\sqrt{n_B{g}_B/m_B}=1.6$cm/s
${\Gamma }_{\mathrm{el},B}$	$21$s-1
${\Gamma }_{3-\mathrm{body},B}$	$0.002$s-1
${\mu }_{BEC}$	$g_B{n}_B=k_B\times 221$nK$=h\times 4.6$kHz
$T_c^{\mathrm{BEC}}$	16.4 $\mu$K
$n_F$	$720$atoms/(10 $\mu$m)2
$E_F$	$k_B\times 130$nK$=h\times 2.7$kHz$=0.6{\mu }_{BEC}$
$T_c^p$	$0.07T_F=5\times {10}^{-4}{T}_c^{\mathrm{BEC}}=9.5$nK
$T_{\mathrm{KT}}$	$0.09T_F=12$nK

Table 2. Table 2: Optical dipole trap configuration. λ L ​ i subscript 𝜆 L 𝑖 \lambda_{{\rm L}i} is the wavelength of dipole-trap beam L ​ i L 𝑖 {\rm L}i , with i = 1 , 2 𝑖 1 2 i=1,2 . w 𝑤 w are the 1 / e 1 𝑒 1/e beam radii. The vertical trap depth for 7 Li, U ⟂ , B subscript 𝑈 perpendicular-to 𝐵 U_{\perp,B} , takes the effect of gravity into account. α L ​ i subscript 𝛼 L 𝑖 \alpha_{{\rm L}i} is the angle between lattice beams L ​ i ​ a L 𝑖 a {\rm L}i{\rm a} and L ​ i ​ b L 𝑖 b {\rm L}i{\rm b} . Δ ​ z Δ 𝑧 \Delta z is the lattice spacing. n 2 ​ D , B subscript 𝑛 2 D B n_{\rm 2D,B} is the density of bosons integrated over the vertical direction. τ B , F = 1 / ∑ i Γ i , B , F subscript 𝜏 B F 1 subscript 𝑖 subscript Γ 𝑖 B F \tau_{\rm B,F}=1/\sum_{i}\Gamma_{i,{\rm B,F}} are limits to the lifetimes of bosons and fermions, where Γ i , B , F subscript Γ 𝑖 B F \Gamma_{i,{\rm B,F}} is the off-resonant scattering rate of photons calculated at peak intensity of dipole trap beam L i 𝑖 i , with i 𝑖 i running over all beams Grimm  et al. ( 2000 ); McKay and DeMarco ( 2011 ); Gordon and Ashkin ( 1980 ) .

${\lambda }_{\mathrm{Lv}}$	1064 nm	$w_{\mathrm{Lv}}$	6 $\mu$m
$z_R$	100 $\mu$m
$U_{⟂,B}$	$k_B\times 0.27\mu$K	${\nu }_{⟂,B}$	1.1 kHz
${\lambda }_{\mathrm{Lh}}$	300 nm or 554 nm
$w_{\mathrm{Lh},⟂}$	2 $\mu$m	$w_{\mathrm{Lh},\parallel }$	200 $\mu$m
${\lambda }_{\mathrm{L1}}$	1064 nm	${\alpha }_{\mathrm{L1}}$	60°
${\lambda }_{\mathrm{L2}}$	470 nm	${\alpha }_{\mathrm{L2}}$	25.5°
$\Delta z$	1064 nm
$U_{⟂,F}$	$h\times 16$kHz	${\nu }_{⟂,F}$	4.1 kHz=1.5 $E_F$
${\tau }_B$	296 s	${\tau }_F$	79 s
$n_{2D,B}$	$3\times {10}^5$atoms/(10 $\mu$m)2

Equations

$$\displaystyle\hat{H}_{B}=\int dz\int d^{2}x\hat{\phi}^{{\dagger}}(t,\mathbf{x},z)\Big{[}-\frac{\nabla^{2}}{2m_{B}}$$

$$\displaystyle\qquad\qquad\quad+\frac{g_{B}}{2}\hat{\phi}^{{\dagger}}(t,\mathbf{x},z)\hat{\phi}(t,\mathbf{x},z)-\mu_{B}\Big{]}\hat{\phi}(t,\mathbf{x},z),$$

$$\displaystyle\hat{H}_{F}=\int d^{2}x\hskip 4.26773pt\hat{{\psi}}^{{\dagger}}(t,\mathbf{x})\Big{[}-\frac{\nabla^{2}}{2m_{F}}-\mu_{F}\Big{]}\hat{\psi}(t,\mathbf{x}),$$

$${\hat{H}}_{FB}=g_{FB}\int dz\int d^{2}x\delta (z){\hat{\psi }}^{†}(t,\mathbf{x}){\hat{\varphi }}^{†}(t,\mathbf{x},z)\hat{\varphi }(t,\mathbf{x},z)\hat{\psi }(t,\mathbf{x}),$$

$$\hat{\varphi }(t,\mathbf{x},z)=\frac{1}{\sqrt{V}}\sum _{\mathbf{q}}{e}^{i\mathbf{q}\cdot \mathbf{r}}{\hat{b}}_{\mathbf{q}}(t),$$

$${\hat{H}}_B(t)=\sum _{\mathbf{q}}\left(\frac{q^2}{2m_B}-{\mu }_B\right){\hat{b}}_{\mathbf{q}}^{†}(t){\hat{b}}_{\mathbf{q}}(t)$$

$$+\frac{g_B}{2V}\sum _{\mathbf{q},{\mathbf{q}}^{\prime },{\mathbf{q}}^{\prime \prime }}{\hat{b}}_{\mathbf{q}+{\mathbf{q}}^{\prime \prime }}^{†}(t){\hat{b}}_{{\mathbf{q}}^{\prime }-{\mathbf{q}}^{\prime \prime }}^{†}(t){\hat{b}}_{\mathbf{q}}(t){\hat{b}}_{{\mathbf{q}}^{\prime }}(t).$$

$$\displaystyle\hat{H}_{B}(t)=\frac{g_{B}N_{0}^{2}}{2V}+\sum_{\mathbf{q}}\Big{(}\frac{q^{2}}{2m_{B}}+n_{B}g_{B}\Big{)}\hat{b}^{\dagger}_{\mathbf{q}}(t)\hat{b}_{\mathbf{q}}(t)$$

$$+\frac{g_B{n}_B}{2}\sum _{\mathbf{q}}\left[{\hat{b}}_{\mathbf{q}}^{†}(t){\hat{b}}_{-\mathbf{q}}^{†}(t)+{\hat{b}}_{\mathbf{q}}(t){\hat{b}}_{-\mathbf{q}}(t)\right].$$

$${\hat{H}}_B(t)=\frac{g_B{n}_B{N}_0}{2}+\sum _{\mathbf{q}(\mathbf{q}\ne 0)}{\omega }_q{\hat{\beta }}_{\mathbf{q}}^{†}(t){\hat{\beta }}_{\mathbf{q}}(t)-\frac{1}{2}\sum _{\mathbf{q}(\mathbf{q}\ne 0)}({\xi }_q-{\omega }_q),$$

$${\xi }_q=\frac{q^2}{2m_B}+g_B{n}_B.$$

$${\hat{H}}_{FB}(t)=g_{FB}{n}_B{N}_F$$

$$+\frac{g_{FB}\sqrt{N_0}}{V}\sum {}_{\mathbf{p},{\mathbf{q}}_{\perp },q_z}^{{}^{\prime }}{V}_q{\hat{a}}_{\mathbf{p}}^{†}(t)\left[{\hat{\beta }}_{\mathbf{q}}(t)+{\hat{\beta }}_{-\mathbf{q}}^{†}(t)\right]{\hat{a}}_{\mathbf{p}-{\mathbf{q}}_{\perp }}(t),$$

$$V_q={\left(\frac{q^2}{q^2+4m_B{g}_B{n}_B}\right)}^{1/4}.$$

$$\Gamma$$

$${\Gamma }^{(2)}$$

$${\lambda }_{\text{eff}}=g_{FB}^2{n}_B{\int }_{-\infty }^{\infty }\frac{d{q}_z}{2\pi }{\left(\frac{\frac{q^2}{2m_B}}{\frac{q^2}{2m_B}+2g_B{n}_B}\right)}^{1/2}\frac{2{\omega }_q}{{\omega }^2-{\omega }_q^2+i\delta }.$$

$${\lambda }_{\text{eff}}$$

$$V_{\text{eff}}(\mathrm{R})=\int d^{2}Q{e}^{i\mathbf{Q}\cdot \mathbf{R}}{V}_{\text{eff}}(Q)$$

$$=-2\pi \frac{g_{FB}^2}{g_B}\frac{1}{{\xi }^2}\frac{1}{\mathrm{R}}{e}^{-\frac{\sqrt{2}}{\xi }\mathrm{R}},$$

$$\displaystyle\hat{H}^{\prime}_{F}=\int\frac{d^{2}p}{(2\pi)^{2}}\bigg{\{}\left(\frac{p^{2}}{2m_{F}}-\mu\right)\hat{a}^{{\dagger}}(\mathbf{p})\hat{a}(\mathbf{p})+\frac{1}{2}\int\frac{d^{2}kd^{2}k^{\prime}}{(2\pi)^{4}}V_{\textrm{eff}}(\mathbf{p},\mathbf{k})$$

$$\displaystyle\times\hat{a}^{{\dagger}}\left({\mathbf{k^{\prime}}}/{2}+\mathbf{k}\right)\hat{a}^{{\dagger}}\left({\mathbf{k^{\prime}}}/{2}-\mathbf{k}\right)\hat{a}\left({\mathbf{k^{\prime}}}/{2}-\mathbf{p}\right)\hat{a}\left({\mathbf{k^{\prime}}}/{2}+\mathbf{p}\right)\bigg{\}},$$

$$V_{\text{eff}}(\mathbf{p},\mathbf{k})=-V_0\frac{1}{\sqrt{|\mathbf{p}-\mathbf{k}{|}^2+2{\xi }^{-2}}},$$

$${\hat{H}}_F^{BCS}=\sum _{\mathbf{p}}{E}_p{\hat{\alpha }}_{\mathbf{p}}^{†}{\hat{\alpha }}_{\mathbf{p}}+$$

$$\displaystyle+\frac{1}{2}{\sum_{\mathbf{p}}}\bigg{\{}\frac{|{\triangle}_{\mathbf{p}}|^{2}}{E_{p}}\Big{[}1-2n_{F}(E_{p})\Big{]}+\left(\epsilon_{p}-E_{p}\right)\bigg{\}},$$

$$\displaystyle\triangle_{\mathbf{p}}=-\int\frac{d^{2}k}{(2\pi)^{2}}V_{\textrm{eff}}(\mathbf{p},\mathbf{k})\frac{{\triangle}_{\mathbf{k}}}{2E_{k}}\Big{[}1-2n_{F}(E_{k})\Big{]}.\quad$$

$$V_{\text{eff}}(\mathbf{p},\mathbf{k})=\sum _{\ell }{V}_{\text{eff}}^{(\ell )}(p,k)\cos [\ell (\theta -\phi )],$$

$$V_{\text{eff}}^{(1)}(k_F)$$

$$\mathcal{F}(X)=\frac{E[-2X^2]-(1+X^2)K[-2X^2]}{X^2},$$

$${△}^{(\ell )}{e}^{i\ell {\theta }_{\hat{\mathbf{p}}}}$$

$$1$$

$$1=-\frac{1}{(2\pi {)}^2}\frac{\pi }{4}{V}_{\text{eff}}^{(\ell )}(k_F)\int kdk\frac{1}{\sqrt{{ϵ}_k^2+|{△}^{(\ell )}{|}^2}}$$

$$1=-\frac{1}{2\pi }{V}_{\text{eff}}^{(\ell )}(k_F)\frac{\pi }{4}\frac{m_F}{2\pi }{\int }_0^{{\Lambda }_{\epsilon }}d\epsilon \frac{1}{\sqrt{{\epsilon }^2+|{△}^{(\ell )}{|}^2}},$$

$$\displaystyle\triangle^{\textrm{Max}}=\triangle^{(1)}=2\Lambda_{\varepsilon}\;\exp\bigg{(}\frac{1}{\rho_{2D}\tilde{V}_{\textrm{eff}}^{(1)}(k_{F})}\bigg{)},$$

$${\tilde{V}}_{\text{eff}}^{(1)}(k_F)=2\sqrt{2}\pi \frac{n_B{a}_{\text{eff}}^2\xi }{m_{FB}}\mathcal{F}(k_F\xi ).$$

$$\xi k_F=\frac{1}{\sqrt{2}}\sqrt{\frac{n_F}{a_B{n}_B}}.$$

Full text

I Introduction

The quest for the experimental realization of a chiral $p_{x}+ip_{y}$ superconductor in two dimensions (2D) is gathering increasing attention because this phase exhibits Majorana modes, which are relevant for constructing fault-tolerant topological quantum computers Kitaev (2001); Tewari et al. (2007). Although a chiral $p$-wave superfluid has been shown to occur in the A-phase of 3He at high pressure vol and experiments have revealed that Strontium ruthenate (Sr2RuO4) is a $p$-wave superconductor Kallin (2012), the manipulation of the Majorana modes in these systems remains difficult. Therefore, the prospect to create a $p$-wave superfluid using ultracold atoms is very appealing because these systems allow for great control of the degrees of freedom.

Several possibilities to generate chiral superfluids have been proposed in the context of ultracold atoms in optical lattices: by using orbital degrees of freedom Soltan-Panahi et al. (2012); Ölschläger et al. (2013), spin-orbit coupling Galitski and Spielman (2013); Zhang et al. (2008a) or dipolar interaction Dutta and Lewenstein (2010); Levinsen et al. (2011). However, these methods either bring new problems to the experimental implementation, such as heating and ultracold chemical-reactions, or require a sophisticated optical-lattice setup and further manipulations to populate the $p$-orbitals.

Here, we adopt a completely different, but feasible route to produce $p$-wave superfluids, which consists of inducing the pairing among the 2D polarized fermionic atoms through a 3D bath of bosonic excitations. The dimensional mismatch between the fermions and the excitations that mediate their interaction leads to a huge increase of the superconducting gap, and consequently of the critical temperature for the observation of the chiral superfluid. The main advantage of our proposal is that it avoids three-body losses and dynamical instabilities (phase separation), which constitute major problems in a strongly-interacting Fermi-Bose mixture.

Mixed-dimension mixtures of two-species fermions with weak interaction were investigated previously Nishida (2009, 2010), with the coupling between polarized fermions in 2D mediated by the particle-hole excitations of a 3D Fermi-sea background. In spite of the high stability of the Fermi-Fermi mixture, the Fermi-Bose mixture, with phonon excitations, provides much higher magnitude for the $p$-wave coupling between fermions. Recently, a 2D-3D mixture of fermions and bosons was considered, and the Berezinskii-Kosterlitz-Thouless (BKT) critical temperature was determined accounting for effects of retardation Wu and Bruun (2016). However, many-body effects were neglected. We argue here that the proximity between the Fermi and sound velocities requires the inclusion of many-body corrections, namely the vertex ladder-diagrams and the RPA dressing of the phonon propagator sch ; Roy et al. (2014).

We calculate these higher-order contributions, which are usually disregarded in the BCS treatment of conventional superconductors, and show that they significantly contribute to increase the magnitude of the anomalous $p$-wave gap in the Fermi-Bose mixture in mixed dimensions. In this calculation, however, we do not consider the renormalization of the pole of the Green’s function, nor take into account retardation effects (the influence of the frequency of the irreducible vertex). The fermions self-energy due to the scattering of the background excitations can be neglected due to the small value of the Fermi-Bose coupling $g_{FB}$, and retardation effects should not provide a relevant contribution to the vertex Efremov et al. (2000) because the singularity for pair formation must come from scattering in the Fermi-surface (Cooper instability sch ; abr ). The simultaneous analysis of both these effects, i.e., retardation and high-order vertex correction, is a tremendous task. Since our calculations are performed in the small momentum limit, if we would consider retardation, it should enhance the positive region of the vertex because correlation between the fermions leads to an even higher prediction to the critical temperature for $p$-wave superfluid formation ($T_{c}^{p}$) Pietronero et al. (1995). Hence, the very high value of $T_{c}^{p}$ that we found due to the vertex correction is actually a lower bound, given the approximations performed.

This paper is structured as follows: Sec. II presents the system Hamiltonian for bosonic and fermionic species, whereas in Sec. III the interaction between the fermions, mediated by the bosonic excitations, is characterized. In sections IV and V, we build the BCS Hamiltonian for the 2D system and solve the associated gap equation, respectively. Higher-order corrections for the gap magnitude are evaluated in Sec. VI, and the experimental feasibility, conclusions and implications of this work are discussed respectively in Sec. VII and Sec. VIII.

II System Hamiltonian

We start by defining the Hamiltonian $\hat{H}=\hat{H}_{B}+\hat{H}_{F}+\hat{H}_{FB}$, where the boson-field operators $\hat{\phi}$ live in 3D, whereas the polarized fermions ${\hat{\psi}}$ live in 2D, (assuming $\hbar=1$)

[TABLE]

with the mass of the bosonic and fermionic species given by $m_{B}$ and $m_{F}$, and their chemical potentials by $\mu_{B}$ and $\mu_{F}$, respectively. The intra- and interspecies contact repulsive interactions are characterized by the coupling constants $g_{B}$ and $g_{FB}$, respectively. We can express the boson-field operators in terms of a discrete set of bosonic modes $\hat{b}_{\mathbf{q}}$, with $V$ the volume of the 3D space,

[TABLE]

which allows us to rewrite the bosonic part of the Hamiltonian in momentum space,

[TABLE]

To characterize the Bose-Einstein condensate, we now use Bogoliubov theory to deal with the macroscopic occupation of the zero-momentum state, that is $\hat{b}_{0}=\hat{b}_{0}^{\dagger}=\sqrt{N_{0}}$. Neglecting higher-order fluctuations, we obtain

[TABLE]

After symmetrizing the above expression, with a sum covering half of the momentum space, and performing a Bogoliubov canonical transformation $\hat{b}_{\mathbf{q}}=u_{q}\hat{\beta}_{\mathbf{q}}-v_{q}\hat{\beta}^{\dagger}_{-\mathbf{q}}$ and $\hat{b}_{-\mathbf{q}}=u_{q}\hat{\beta}_{-\mathbf{q}}-v_{q}\hat{\beta}^{\dagger}_{\mathbf{q}}$, where we select the real parameters $u_{q},v_{q}$ in order to have diagonal-base operators ($\hat{\beta},\hat{\beta}^{\dagger}$) for $H_{B}$, we find

[TABLE]

with the energy spectrum for the free Bogoliubov-modes excitation $\omega_{q}=\sqrt{\xi_{q}^{2}-(g_{B}n_{B})^{2}}$, where

[TABLE]

Applying the same set of transformations for the interspecies-interaction Hamiltonian ($H_{FB}$), and considering $u_{q}=\sqrt{{\xi_{q}}/{\omega_{q}}+1}/{\sqrt{2}}$ and $v_{q}=\sqrt{{\xi_{q}}/{\omega_{q}}-1}/{\sqrt{2}}$, with $\hat{\psi}(t,\mathbf{x})=({1}/{\sqrt{S}})\sum_{\mathbf{p}}e^{i\mathbf{p}\cdot\mathbf{x}}\hat{a}_{\mathbf{p}}(t)$, where $S$ denotes the 2D surface, we get

[TABLE]

with

[TABLE]

In Eq. (LABEL:eq9), the prime symbol in the sum indicates that $\mathbf{q}\neq 0$, and we separate the components of $\mathbf{q}=(\mathbf{q}_{\bot},q_{z})$, to account for momentum conservation in the plane.

III Effective Interaction

As expressed in Eq. (2), there is no direct interaction between the polarized fermions in $H_{F}$, due to the Pauli exclusion principle. We show here, however, how an indirect interaction between fermions arises from $H_{FB}$. For that, we define the effective coupling constant $\lambda_{\textrm{eff}}$ from the four-point function $\Gamma=\Gamma(\mathbf{p},\mathbf{p^{\prime}},\mathbf{k},\mathbf{k^{\prime}};\varepsilon,\varepsilon^{\prime},\nu,\nu^{\prime})$ as follows

[TABLE]

with $G_{0}$ corresponding to the free-fermion propagator and $\omega=\varepsilon-\varepsilon^{\prime}=\nu^{\prime}-\nu$.

Considering the weak-coupling regime, to second order in the interaction (see Fig. 1), we obtain

[TABLE]

where $D_{0}(\mathbf{q},\omega)$ denotes the free-phonon propagator and $\mathbf{q}_{\bot}=\mathbf{p}-\mathbf{p^{\prime}}=\mathbf{k^{\prime}}-\mathbf{k}$. Comparing Eq. (11) and Eq. (12), we find

[TABLE]

For low-energy processes, where the scattered fermions are kept around the 2D Fermi surface, we can assume $\omega\sim 0$, and Eq. (III) can be simplified as

[TABLE]

Hence, an effective potential $\lambda_{\textrm{eff}}=V_{\textrm{eff}}(q_{\bot}=|\mathbf{p^{\prime}}-\mathbf{p}|)$ is generated between the fermions, as a function of the momentum exchange $\mathbf{Q}$ between the scattered particles. In 2D real space, with coordinate $\mathbf{R}$, this yields an attractive Yukawa potential between the fermionic particles in the plane,

[TABLE]

with range given by the healing length $\xi=1/\sqrt{2m_{B}g_{B}n_{B}}$ of the BEC.

IV BCS Hamiltonian

We consider the generalized BCS-type Hamiltonian in momentum space for the fermions in the plane,

[TABLE]

with a momentum-dependent mediated interaction $V_{\textrm{eff}}(\mathbf{p},\mathbf{k})$ and $\mu=\mu_{F}-n_{B}g_{FB}$. According to Eq. (14), we consider the interaction potential

[TABLE]

with $V_{0}=2g_{FB}^{2}n_{B}m_{B}$. After symmetrizing the BCS Hamiltonian properly, we apply the Bogoliubov transformation and find a new basis of operators (see App. A for details) to build the diagonal form

[TABLE]

with the energy dispersion $E_{p}=\sqrt{\epsilon_{p}^{2}+|{\triangle}_{\mathbf{p}}|^{2}}$ and the occupation function $n_{F}(E_{p})=[\exp(\beta E_{p})+1]^{-1}$ of the Bogoliubov modes, where $\beta=(k_{B}T)^{-1}$. As shown in App. A, now we can also write the gap in terms of the mean value over this new basis, to obtain

[TABLE]

V GAP Equation

To solve the integral equation for a momentum-dependent pairing gap in Eq. (19), it is convenient to use the 2D partial-wave expansion of the effective potential Anderson and Morel (1961); Chubukov (1993),

[TABLE]

with $\ell$ integer, $p=|\mathbf{p}|$, $k=|\mathbf{k}|$, and where we associated the angles ${\theta}_{\mathbf{\hat{p}}}=\theta$ and ${\theta}_{\mathbf{\hat{k}}}=\varphi$. Because we are assuming low-energy processes, with the scattered momentum close to the Fermi surface, it is reasonable to consider $p\sim k=k_{F}$ in the coefficients of Eq. (20). For $\ell=1$, considering the even parity of the potential, we have

[TABLE]

where

[TABLE]

with $E[X]$ the complete elliptic integral, $K[X]$ the complete elliptic integral of the first kind, and $X=k_{F}\xi$ (see the inset of Fig. 2).

Since in the weak-coupling limit one expects that the mixing of different angular momentum $\ell$ will be small, we are in a position to solve the gap equation by applying the pure $\ell$-type ansatz $\triangle_{\mathbf{p}}=\triangle^{(\ell)}e^{i\ell\theta_{\hat{\mathbf{p}}}}$ in Eq. (19). That gives

[TABLE]

Analytical solutions for $\triangle^{\textrm{Max}}$ and $T_{c}$ can be obtained in two limiting cases: 1)$\;T\rightarrow 0$, where we have the maximum gap value, and 2) $\;T\rightarrow T_{c}$, where the gap goes to zero. For the first limit, we find $E_{k}=\sqrt{\epsilon_{k}^{2}+|\triangle^{(\ell)}|^{2}}$ and $n_{F}(E_{k})\rightarrow 0$. Then, applying the orthogonality condition given by the angular integral of equation (V), we eliminate the sum in $\ell^{\prime}$ to obtain

[TABLE]

where we can identify the density of states in the Fermi surface $\rho_{2D}=m_{F}/2\pi$ and the cut-off energy scale given by the Fermi energy of the 2D system $\Lambda_{\varepsilon}\sim k_{F}^{2}/2m_{F}$. Since we consider the small-momentum regime, the fermions are scattered to states around the Fermi level. As can be seen from Table 1 in the experimental section, $k_{F}$ is very close to the healing length ($\xi^{-1}$), which characterizes the range of the interaction potential.

One can show that the induced attraction Eq. (17) is strongest in the $p$-wave channel. That means that the dominant pairing instability is in the channel with orbital angular momentum $\ell=1$, and the most stable low-temperature phase, or with highest critical temperature, has $p_{x}+ip_{y}$ symmetry Anderson and Morel (1961); Nishida (2009). We can then solve Eq. (LABEL:eq8) for the maximum gap

[TABLE]

with $\;\tilde{V}_{\textrm{eff}}^{(1)}(k_{F})={V}_{\textrm{eff}}^{(1)}(k_{F})/8$.

The vertex renormalization for two particles in vacuum allows us to express the bare coupling parameter as $g_{FB}\rightarrow-2\pi a_{\textrm{eff}}/\sqrt{m_{B}m_{FB}}$ Nishida and Tan (2008), with the reduced mass $m_{FB}=m_{B}m_{F}/(m_{B}+m_{F})$ and the effective two-body scattering length $a_{\textrm{eff}}$ for a 2D-3D scattering. The latter will be a function of the original 3D scattering length $a_{FB}$ and of the axial confinement. That gives

[TABLE]

Considering $k_{F}=\sqrt{4\pi n_{F}}$ and $\xi=1/\sqrt{8\pi n_{B}a_{B}}$, we get the variable

[TABLE]

Thus, we estimate the gap in Eq. (25) using

[TABLE]

For $a_{B}n_{B}^{1/3}\sim 0.01$ and $a_{\textrm{eff}}k_{F}\sim 0.1$, we consider the maximum value for $\rho_{2D}|\tilde{V}_{\textrm{eff}}^{(1)}(k_{F})|$ with $\mathcal{F}(X)/X\sim-0.15$, restricting $X$ in the interval $[0.5-1.5]$ (see Fig. 2), to determine cal

[TABLE]

VI Higher order correction to the Effective 2D-3D interaction

The previous section shows how to optimize the gap value by manipulating the condensate density, which controls the magnitude and range of the induced potential. In addition, the importance of choosing an appropriate combination of the fermion and boson atomic masses (lighter bosonic species) to maximize the gap became clear. This issue will be further explored in Sec. VII.

By choosing the Fermi wavelength and the healing length such that $\xi k_{F}\sim 1$, the Bogoliubov-sound ($c_{s}$) and the Fermi velocities ($v_{F}$) will also have close values. That requires the inclusion of higher-order diagrammatic terms in our ultracold-atoms model, which are usually disregarded in BCS studies.

In the following, we calculate the four-point function to $4th$ order in the interaction constant $g_{FB}$

[TABLE]

We start with the interaction between the fermions in 2D and the “phonons” of the BEC in 3D as given by Eq. (LABEL:eq9) and Eq. (10). Using the finite temperature formalism with the Matsubara Green’s functions, the effective interaction between the fermions in 2D is given by

[TABLE]

with the free-fermion propagator $\mathcal{G}_{0}$. As seen before, the second-order expansion in the coupling $g_{FB}$ provides

[TABLE]

where $\mathbf{q}\equiv(\mathbf{k}_{1}-\mathbf{k}_{4},q_{z})$ and we applied the static limit to the Bogoliubov-mode propagator $\mathcal{D}_{0}$.

Within a higher-order expansion, we obtain the self-energy bubble diagram (see the details of the calculation in App. B)

[TABLE]

where we identify the static polarization-bubble diagram in 2D

[TABLE]

For $|\mathbf{k}_{1}-\mathbf{k}_{4}|<2k_{F}$, i.e., the external momenta in the Fermi surface, we can easily calculate the RPA series, which yields

[TABLE]

where we defined $\lambda_{0}=-{V_{0}}/{\sqrt{|\mathbf{k}_{1}-\mathbf{k}_{4}|^{2}+2\xi^{-2}}}$ and $P_{0}=-{m_{F}}/{2\pi}=-\rho_{2D}$. For $\lambda_{0}P_{0}<1$, we find

[TABLE]

Replacing Eq. (17) by the effective potential coming from the RPA correction in Eq. (VI), we obtain an increase in the gap magnitude, as predicted by Eq. (25) (see also App. B and Fig. 3).

Since we consider $\lambda_{0}P_{0}$ smaller than one, we do not expect any phase instability driven by a divergence of $\lambda_{\textrm{eff}}^{RPA}$ caused by the vanishing of the denominator of Eq. (VI).

The critical condition given by Eq. (25) can be obtained alternatively through the singularity in the effective interaction, which appears when the total vertex function is calculated in the Fermi surface, considering small total momentum of the colliding particles abr ; pit ; Chubukov (1993). In this case, the $\ell$-th harmonic in the exponent of Eq. (25) will be associated with the irreducible part of the vertex. Here, we determined its $\ell=1$ projection solving the Bethe-Salpeter integral equation for the ladder-series contribution. To build the series, we start with the $4th$-order vertex-correction, which reads

[TABLE]

with $\omega_{\mathbf{q}}=\frac{q}{2m_{B}}\sqrt{q^{2}+2\xi^{-2}}$ and $\mathbf{q}\equiv(\mathbf{k}_{3}-\mathbf{p},q_{z})$. The first term of Eq. (LABEL:vertx1) is related to single-particle behavior, i.e., the scattering of real phonons, whereas the second term corresponds to virtual phonon processes. Only the latter will be relevant in our calculation, which deals with the many-body effects with the 2D momentum integration performed near the Fermi surface.

To evaluate the irreducible-vertex part around the Fermi surface, perturbation theory turns out to be insufficient and we must sum the whole ladder series of diagrams, with terms proportional to the ratio $c_{s}/v_{F}$. The resulting self-consistent vertex equation is presented and solved in the App. B, after performing a partial expansion of the effective interaction $\lambda_{\textrm{eff}}^{V}$ in terms of the angular components $\lambda(|\mathbf{k}_{4}-\mathbf{k}_{1}|)=\sum_{\ell}\lambda^{(\ell)}(k_{F})\cos[\ell(\theta_{4}-\theta_{1})]$ Chubukov (1993); pit , which breaks the integral equation for the total pairing vertex to a set of decoupled algebraic equations for its partial components. Finally, we obtain the vertex correction for the component $\ell=1$

[TABLE]

where we defined $\mathcal{J}[X]=(1+2X^{2})E\left[1-\frac{1}{1+2X^{2}}\right]-(1+X^{2})K\left[1-\frac{1}{1+2X^{2}}\right]$. Remarkably, $\frac{\mathcal{J}[X]}{\mathcal{F}[X]X^{2}\sqrt{1+2X^{2}}}=1!$ Including the correction given by Eq. (38) into the gap equation, according to Eq. (25), we get

[TABLE]

This is the main result of this paper: the inclusion of higher-order diagrams, usually neglected due to their complexity, actually increases the $p$-wave gap by one order of magnitude and brings it to the verge of experimental possibilities.

VII Experimental implementation

We now discuss the experimental feasibility of our proposal. We first examine which quantum gas mixtures are suitable to implement it, then present a scheme for a mixed-dimensional trap, and finally we summarize the experimental proposals to detect a $p$-wave superfluid.

VII.1 Mixture choice

The most important criterion to choose the mixture is that the critical temperature for $p$-wave superfluidity $T_{c}^{p}$ has to be experimentally reachable McKay and DeMarco (2011). As guidance, we note that BECs have been evaporatively cooled to $T=0.02T_{c}^{\rm BEC}=1\,$nK Olf et al. (2015) and Fermi gases with $T/T_{F}\leq 0.05$ have been reached Navon et al. (2010). We maximize $T_{c}^{p}/T_{F}=\gamma\triangle_{V}^{\rm Max}/T_{F}$ sch under constraints imposed by the validity of our theory and experimental constraints ($\gamma$: Euler’s constant $\simeq 0.57$). The static approximation requires that $\alpha=v_{F}/c_{s}\lesssim 1$ Migdal (1958); Roy et al. (2014). In addition, since the effective potential has been obtained within a perturbative treatment, it is necessary that $\gamma_{eff}^{2}<(8\pi\gamma_{BEC})^{1/2}$. Hence, the boundaries of validity of our theoretical studies request $\gamma_{BEC}=a_{B}n_{B}^{1/3}\gtrsim 10^{-3}$ Wu and Bruun (2016) and $\gamma_{\rm eff}=a_{\rm eff}n_{B}^{1/3}\lesssim(8\pi\gamma_{BEC})^{1/4}\approx 0.4$. To be in the superfluid regime we finally require $T^{p}_{c}<T_{\rm KT}$, where $T_{\rm KT}$ is the Kosterlitz-Thouless transition temperature Fisher and Hohenberg (1988); Desbuquois et al. (2012). Since $T_{c}^{p}/T_{F}=8.42\exp(-1/|\rho_{2D}\tilde{V}_{\textrm{eff}}^{(1)}|)$ increases monotonically with $Y=|\rho_{2D}\tilde{V}_{\textrm{eff}}^{(1)}|$ it is sufficient to maximize $Y$, which can be expressed as

[TABLE]

with $X=\alpha(m_{F}/m_{B})/\sqrt{2}$. For large $Y$, a high mass ratio $m_{F}/m_{B}$ should be selected, provided that $\alpha$ is chosen close to $\alpha_{\rm max}=3.56\,m_{B}/m_{F}$, which maximizes $|\mathcal{F}(X)|$. Since $T_{F}=(2\pi\hbar^{2}/k_{B})(m_{F}/m^{2}_{B})n_{B}^{2/3}\alpha^{2}\gamma_{BEC}\propto\alpha^{2}$, we chose in the following a slightly higher value, $\alpha=1.5\,\alpha_{\rm max}$, which barely decreases $|\mathcal{F}(X)|$, but more than doubles $T_{F}$. Furthermore, a low value of $\gamma_{BEC}$ is desired and we chose a value close to its minimum. Finally, a high value of $\gamma_{\textrm{eff}}$ has to be achieved. In order to increase $\gamma_{\rm eff}$, we opt for the rather high value of $n_{B}=6\times 10^{14}\,$atoms/cm*-3* and the relatively low value of $a_{\rm eff}=204\,a_{0}$, where $a_{0}$ is the Bohr radius. The motivation for choosing a large density is that $T_{F}$ increases with $n_{B}$. On the other hand, low values of $a_{\rm eff}$ are more likely available in experiments than large values, and they can be reached without Feshbach or confinement induced resonances. Far from the resonances, the scattering length is given approximately by $a_{\rm eff}\sim\sqrt{m_{B}/m_{FB}}\,a_{\rm FB}$ Massignan and Castin (2006); Nishida and Tan (2008); Lamporesi et al. (2010).

Further limitations arise from experimental constraints. In our scheme, a few thousand fermions will be sympathetically cooled by a much larger bath of evaporatively cooled bosons. To effectively implement evaporative and sympathetic cooling, a sufficient rate of elastic collisions and low rates of heating and loss are required. These conditions limit the range of suitable interaction properties, the gas densities, and the trap designs. An upper limit on $n_{B}$ is imposed by the requirement to keep the BEC in the 3D regime for the finite number of bosons available. A lower limit on $a_{B}$ is imposed by the requirement of a sufficient elastic collision rate between bosons $\Gamma_{\rm el,B}\propto n_{B}a_{B}^{2}$. Together, these requirements lead to an additional, experimental, lower limit on $\gamma_{BEC}$. Attention has also to be given to the rate of 3-body losses involving one fermion and two bosons ($\Gamma_{\rm FBB}\propto n_{B}^{2}a_{\rm FB}^{4}$ D’Incao and Esry (2008); Laurent et al. (2017)), even considering the important role played by the mixed dimensionality in inhibiting the interspecies molecular formation Nishida (2009).

We now discuss possible choices of elements for the mixture. Since $m_{F}/m_{B}$ should be large, we limit our choice of bosons to the lightweight isotopes that have been Bose condensed, 4He∗, 7Li, and 23Na. Among those, 7Li has the great advantage of possessing a broad Feshbach resonance, with which $a_{B}$ can be tuned Chin et al. (2010); Pollack et al. (2009b); Gross et al. (2009); Pollack et al. (2009a). Feshbach resonances in 4He∗ and 23Na are expected or known to be accompanied by strong losses Vassen et al. (2012); Goosen et al. (2010); Borbely et al. (2012); Inouye et al. (1998); Stenger et al. (1999). In the following, we use the triplet-scattering length for 4He∗ and 23Na Tiesinga et al. (1996); Moal et al. (2006). Considering BEC densities for which inelastic collisions limit the BEC lifetime to 10 s Robert et al. (2001); Pereira Dos Santos et al. (2001); Stamper-Kurn et al. (1998), fermion masses up to the mass of the heaviest naturally occurring fermionic isotope (235U) and $a_{\rm eff}=600\,a_{0}$, we find that $T_{c}^{p}/T_{F}<10^{-2}$ for these bosons. Only larger values of $a_{\rm eff}$ might make them suitable for our purposes.

We therefore limit our considerations to 7Li. This choice makes it possible to decrease $a_{B}$ and thereby increase $T_{c}^{p}/T_{F}$. To choose the fermionic element we plot in Fig. 4a) $T_{c}^{p}/T_{F}$ and $T_{c}^{p}$ as a fuction of $m_{F}$. Fermionic isotopes that have been cooled to quantum degeneracy and for which the experimentally relevant regime $T_{c}^{p}/T_{F}>0.05$ can be reached are 171,173Yb, 161Dy, and 167Er DeMarco and Jin (1999); Naylor et al. (2015); DeSalvo et al. (2010); Taie et al. (2010); Lu et al. (2012); Aikawa et al. (2014). A drawback of having to choose such heavy elements could be that they are not well sympathetically cooled by the lightweight Li because during each elastic collision, the energy transfer from the fermion to the boson is suppressed by $4m_{F}m_{B}/(m_{F}+m_{B})^{2}\sim 0.15$ Mudrich et al. (2002). A benefit of Dy and Er compared to Yb is that several interspecies Feshbach resonances will likely be available across the broad 7Li Feshbach resonance, making it possible to tune $a_{B}$ and $a_{\rm FB}$ somewhat independently and to access large values of $a_{\rm FB}$, which would also make tuning of $a_{\rm eff}$ by confinement induced resonances possible.

Nevertheless, since 173,174Yb-6Li mixtures are already available in the lab Hara et al. (2011); Hansen et al. (2013), we concentrate our discussion now on 171,173Yb-7Li. Adapting the existing machines to operate with 7Li instead of 6Li should be straightforward. There are two fermionic Yb isotopes readily available, each providing a chance of possessing suitable interspecies interaction properties with 7Li. Figure 4b) shows the dependence of $T_{c}^{p}/T_{F}$ and $T_{c}^{p}$ on $n_{B}$ and $a_{\rm eff}$. Choosing $a_{B}=8\,a_{0}$ leads to the system parameters given in Table 1. The dotted lines in Fig. 4 are an estimation of the Kosterlitz-Thouless transition temperature, which is given by Fisher and Hohenberg (1988); Prokof’ev et al. (2001)

[TABLE]

where $m$ and $n$ are the mass and density of the superfluid species, while $a$ characterizes the range of the interaction. In particular, for our case of fermionic-pair formation, the interaction between fermions that will form the Cooper pairs is proportional to $a_{FB}^{2}$, with $m=2m_{F}$ and $n\sim n_{F}/2$. Eq. (41) is valid for small interaction parameters $a_{B}$ and $a_{FB}$ - the first makes the range of the potential long enough, such that the superfluid fraction achieves its maximum value Nishida (2010); Wu and Bruun (2016).

The critical temperature $T_{c}^{p}=0.07\,T_{F}=9.5\,$nK is in the regime of temperatures that have already been achieved experimentally, albeit in systems with larger elastic scattering length. However, $T_{c}^{p}/T_{c}^{\rm BEC}=5\times 10^{-4}$ is more than one order of magnitude lower than what has been reached so far. To enhance evaporative cooling, it might be useful to first evaporate at a scattering length above 100 $a_{0}$ and to tune the scattering length to a lower value only when approaching the required low temperature, while compressing the gas at the same time. In doing so, one could even profit from a Li 3-body recombination minimum at $a_{B}=119\,a_{0}$ Pollack et al. (2009a).

VII.2 Trap configuration

Next, we consider suitable trap configurations for the mixture. Whereas the bosons explore a 3D trap, the fermions have to be effectively confined in 2D by a harmonic trap of frequency $\nu_{\perp,F}$, which requires $h\nu_{\perp,F}-E_{F}\gg k_{B}T$. The sample should be as homogeneous as possible to avoid inhomogeneous broadening of $p$-wave superfluidity signals, especially because the number of fermions will be low. Efficient evaporative cooling of the bosons should be possible in order to reach low temperatures. We now take these requirements into account to design an optical dipole trap for the mixture, where we orient the 2D plane of the fermions in the horizontal direction, see Fig. 5a.

The bosonic lithium surrounds the fermions and can be confined by a Gauss-beam dipole trap using a wavelength of 1064 nm. To reach a temperature $T$ by evaporation, the trap depth in the vertical direction $U_{\perp,B}$ should be $\mu_{\rm BEC}+\eta k_{B}T$, where $\mu_{\rm BEC}$ is the chemical potential of the BEC, and $\eta\sim 5$ Ketterle and van Druten (1996). In order to provide a homogeneous vertical trap frequency across the cloud, the horizontal waist should be much larger than the cloud and the vertical Rayleigh length $z_{R}$ much longer than the horizontal sample size. The latter requirement and the additional requirement $h\nu_{\perp,B}\ll\mu_{\rm BEC}$ are only fulfilled if the vertical waist is larger than a minimum size. At the same time, the vertical waist should not be too large in order to limit the size of the 7Li sample in the vertical direction, thereby reducing the required number of 7Li atoms. Gravitational sag of the bosonic cloud is compensated by placing the focus of the Gauss beam slightly above the plane of the fermions. The Gaussian-beam trap creates a nearly constant potential on the fermions, since they explore only a small region in the centre of the trap. A constant potential offset is irrelevant and we can therefore ignore the influence of the Gauss-beam dipole trap on the fermions.

To provide homogeneous confinement for bosons and fermions in the horizontal plane, repulsive dipole trap walls can be erected around the sample using vertically propagating Gauss beams Meyrath et al. (2005a); Gaunt et al. (2013). Four such beams can form a rectangular box with a size of $\sim 10\,\mu$m around the sample, if the waist of the beams is elongated along the sides of the rectangle ($w_{\rm Lh,\parallel}$ of a few 10 $\mu$m) and is narrow orthogonal to that direction ($w_{\rm Lh,\perp}\sim 2\,\mu$m). This rectangular potential box also serves to select the most homogeneous central region of the traps that are used to confine bosons and fermions vertically. The sample density can easily be changed by moving the vertical walls towards each other, which is useful to do while $a_{B}$ is reduced to a low value. If in further studies a cylindrically symmetric system is required, for example to enable the creation of vortices Madison et al. (2000), a Laguerre-Gaussian beam can be used to confine the atoms horizontally Kaplan et al. (2002); Jaouadi et al. (2010); Gaunt et al. (2013).

The confinement of the fermions in quasi-2D is most conveniently done using optical lattices. In comparison to other trap configurations, such as a Hermite-Gaussian beam Meyrath et al. (2005b, a), it is easier to create a more homogeneous confinement in the 2D plane by increasing the diameter of the lattice beams. In order to populate only a single plane of the lattice with fermions, one can use the techniques of Refs. Gemelke et al. (2009); Sherson et al. (2010); Yamamoto et al. (2016); Ville et al. (2016).

The deep dipole potential used to confine the fermions in 2D may only have a negligible effect on the bosons. The parasitic potential on the bosons $U_{\rm lattice,B}$ must be much smaller than $\mu_{\rm BEC}$. This challenge has been met by species-specific dipole traps using a "tune-out" wavelength, for which the AC polarizability of one species is zero Massignan and Castin (2006); LeBlanc and Thywissen (2007); Catani et al. (2009); Lamporesi et al. (2010). Unfortunately, this technique does not work for 7Li because its "tune-out" wavelength is too close to an atomic transition, leading to detrimental off-resonant scattering for the required trap depths LeBlanc and Thywissen (2007). Another option is to use a "tune-in" wavelength, close to an Yb transition and far detuned from any Li transition LeBlanc and Thywissen (2007). In this situation, the potential on Yb $U_{\rm lattice,F}$ can exceed the potential on Li many times. This technique is suitable for our situation, but will limit the lifetime of the fermionic cloud to a few seconds by off-resonant scattering. If this limit is significant depends on the other factors limiting the lifetime of the system, especially the unknown 3-body loss rate $\Gamma_{\rm FBB}$.

If the lifetime limit imposed by a "tune-in" lattice is too severe, a bichromatic dipole trap can be used, consisting of two optical lattices that both confine Yb, but compensate each other for Li. This technique overcomes the possibly excessive off-resonant scattering and replaces it by the technical challenge of creating two lattices with very well controlled intensity profiles. We will explore this scheme in the following. We chose optical lattices with wavelengths of 470 nm and 1064 nm, which are both attractive for Yb. In contrast, for Li only the 1064-nm lattice is attractive, the other is repulsive, see Fig. 5b. In order for the lattice potentials to add up for Yb and cancel sufficiently for Li, the intensity profile of both lattices need to be nearly identical in the region of the atomic clouds. The lattice-well spacing must be the same, and the intensity maxima need to overlap. The lattice spacing can be adjusted by the angle between the two lattice beams of each wavelength. Using an angle of 60∘ between the two beams forming the 1064-nm lattice leads to a lattice spacing of 1064 nm. The same spacing is reached for the 470-nm lattice if the two corresponding beams intersect at an angle of 25.5∘, see Fig. 5a. The position of the intensity maxima along the lattice direction (the vertical direction) depends on the phase difference between the two beams forming a lattice. This phase difference has to be stabilized interferometrically for each lattice to a common reference, combining methods from Refs. Fölling et al. (2007); Wirth et al. (2010). In order for the two lattice potentials to cancel for the bosons, the intensity of the 470-nm lattice beams has to be 1.8 times the intensity of the 1064-nm lattice beams. For Yb the two lattice potentials add up, giving a total potential that is 1.2 times larger than the potential of the 470-nm lattice alone. This total potential needs to confine Yb in quasi-2D and be also deep enough to suppress tunneling of Yb to neighboring lattice planes, see Fig. 6. The cancelation of the lattice potential for the bosons will not be perfect because of intensity and phase fluctuations leading to deviations from the ideal configuration. Phase fluctuations of 90 mrad or intensity imbalances of 9% lead to a residual potential on the order of 10% of $\mu_{\rm BEC}$. This parasitic potential would be tolerable if the timescale of fluctuations is large enough to avoid heating of the sample. In principle, we could have chosen a wavelength for L2 that is further away from the Yb transition, e.g. 532 nm, which would reduce off-resonant scattering and simplify phase locking of the laser sources used for L1 and L2. All the same, we chose 470-nm because at that wavelength we are profiting from less parasitic potential of L2 on Li, reducing the amount of compensation needed from L1. As a result, the overall parasitic potential created for a given intensity or phase mismatch between L1 and L2 is reduced.

Example parameters for the bichromatic dipole trap and important results of using this trap for the Li-Yb mixture are given in Table 2. The 7Li atom number available in current experiments ($3\times 10^{5}$ atoms Pollack et al. (2009b)) is sufficient for a square sample of 10 $\mu$m size. A sample of this size contains about 700 fermions. If this proposal is realisable depends to a large extend on the unknown elastic and inelastic scattering properties of Li-Yb. Similar schemes can be applied to other mixtures, such as Li-Dy or Li-Er, for which some interspecies interaction tuning should be possible.

VII.3 Detection of $p$-wave superfluidity

There are some predictable signatures for the experimental detection of the $p_{x}+ip_{y}$ superfluid phase. Particularly, the density of state (rf absorption spectrum) of a rotating weak pairing $p_{x}+ip_{y}$ phase is expected to exhibit a set of gapless modes Grosfeld et al. (2007), which are a direct consequence of the zero-energy Majorana modes on the vortices. The rf-spectroscopy can be also applied to detect Majorana edge states of the topological superfluid in a 2D square lattice Midtgaard et al. (2016). On the other hand, the time-reversal symmetry broken signature of the chiral $p_{x}+ip_{y}$ fermionic superfluid can be detected with time-of-flight image of the atomic density distribution: an external effective electric field (i.e., dipole interaction between the neutral atoms in the superfluid and the laser field) brings a nonzero antisymmetric transverse mass current in the velocity distribution of the atoms Zhang et al. (2008b).

VIII Conclusion

In the present work, we explored the feasibility of a $p$-wave superfluid by using a Fermi-Bose mixture in a mixed-dimension configuration, where $p$-wave interaction between spin-polarized degenerate fermions in 2D is induced indirectly, through the scattering of the Bogoliubov modes of condensed bosons moving in 3D. We have shown that, even in the weak-coupling regime, the appropriate renormalization of the phonon propagator (BEC modes) with particle-hole fluctuations and the vertex correction significantly increase the gap and the predicted critical temperature for the fermion-pair formation.

It is important to remark that we adopt a minimum value for $\gamma_{BEC}\sim a_{B}n_{B}^{1/3}$, which yields $\upsilon_{F}/c_{s}\leq 1$, thus allowing to disregard retardation effects. According to Wu and Bruun Wu and Bruun (2016), who performed calculations including retardation but no vertex correction to determine $T_{MF}$, in the limit $\upsilon_{F}/c_{s}\leq 1$, it holds that $T_{MF}\sim T_{BCS}$ (see Fig.2 in the cited reference), which confirms the validity of our approximation.

We neglected decay of the BEC phonons, like the Beliaev damping and the lifetime due to the scattered particle-hole pairs of the degenerate fermionic sample. The Beliaev damping is given by the boson-boson scattering potential, resulting in a phonon lifetime proportional to $g_{B}$ Ozeri et al. (2005); Matveenko and Shlyapnikov (2011). In the small-momentum regime, however, the Beliaev decay mechanism is strongly suppressed Ozeri et al. (2005). On the other hand, if we consider the phonon dressed by particle-hole fluctuations of the Fermi sea, it will have a lifetime proportional to $g_{FB}^{2}$. In the static limit considered in the paper, however, the lifetime is infinite (see App. B for details). Hence, we conclude that there is no damping mechanism that could hamper the stability of the BEC in the chosen regime of parameters.

Exploiting the difference in polarizability and mass of the atomic species, and by optimizing the density $n_{B}$ and the scattering length $a_{B}$ of the bosons, our work sets the boundary for the experimental realization of a $p$-wave superfluid within the reachable limit of $T^{p}_{c}=0.05T_{F}$. It identifies a realistic route and provides the details to the accomplishment and manipulation of this long-sought fascinating chiral-superfluid phase in the realm of ultracold atoms in optical lattices.

Acknowledgments

We thank Rodrigo G. Pereira, Servaas Kokkelmans, Frédéric Chevy, and Subhadeep Gupta for discussions and insightful comments. This work was supported by CNPq (Brazil) through the Brazilian government project Science Without Borders. The work of C.M.S. is part of the DITP consortium, a program of the Netherlands Organisation for Scientific Research (NWO) that is funded by the Dutch Ministry of Education, Culture and Science (OCW). F.S. gratefully acknowledges funding from the European Research Council (ERC) under Project No. 615117 QuantStro and by NWO through Vici grant No. 680-47-619.

Appendix A Bogoliubov transformation in the BCS Hamiltonian

Starting with the definition

[TABLE]

we can apply a mean-field approach in Eq. (LABEL:eq4) and replace the pair operator $\hat{A}(\mathbf{k},\mathbf{p})$ by $\langle\hat{A}(\mathbf{k},\mathbf{p})\rangle+\delta\hat{A}(\mathbf{k},\mathbf{p})$ (similar expression for its conjugate), with $\langle\hat{A}(\mathbf{k},\mathbf{p})\rangle=\delta_{\mathbf{k},\mathbf{0}}\langle\hat{a}(-\mathbf{p})\hat{a}(\mathbf{p})\rangle$ and \langle\hat{A}^{\dagger}(\mathbf{k},\mathbf{p})\rangle=\delta_{\mathbf{k},\mathbf{0}}\langle\hat{a}^{\dagger}\big{(}\mathbf{p}\big{)}\hat{a}^{\dagger}(-\mathbf{p})\rangle. Holding terms up to the first order in the fluctuations of this field (neglecting $\mathcal{O}[(\delta\hat{A})^{n}],\;n>1)$, we find

[TABLE]

with $\epsilon_{p}={p^{2}}/{2m_{F}}-\mu$ and the order parameter (or momentum-dependent gap) expressed as

[TABLE]

where we consider the interaction potential

[TABLE]

with $V_{0}=2g_{FB}^{2}n_{B}m_{B}$. Before applying the Bogoliubov transformation, let us first symmetrize this BCS Hamiltonian properly. It is easier to go further with this process in the discrete-momentum space, summing over half of the k-space $\sum_{\mathbf{k}}\rightarrow\sum_{\mathbf{k}}^{\prime}$

[TABLE]

where we used the property ${\triangle}_{-\mathbf{p}}=-{\triangle}_{\mathbf{p}}$, which is simple to prove if we consider that $V_{\textrm{eff}}(-\mathbf{p},\mathbf{k})=V_{\textrm{eff}}(\mathbf{p},-\mathbf{k})$ and $V_{\textrm{eff}}(-\mathbf{p},-\mathbf{k})=V_{\textrm{eff}}(\mathbf{p},\mathbf{k})$, as can be promptly verified from Eq. (45).

Now, we apply the canonical transformation

[TABLE]

with $|u_{\mathbf{p}}|^{2}+|v_{\mathbf{p}}|^{2}=1$. To diagonalize the transformed Hamiltonian, we set the coefficients of the off-diagonal terms to zero, $2\epsilon_{p}u_{\mathbf{p}}\bar{v}_{\mathbf{p}}-{\triangle}^{\ast}_{\mathbf{p}}u_{\mathbf{p}}^{2}+{\triangle}_{\mathbf{p}}\bar{v}_{\mathbf{p}}^{2}=0$. Multiplying this equation by ${\triangle}_{\mathbf{p}}/u_{\mathbf{p}}^{2}$, we get

[TABLE]

and then

[TABLE]

with the energy dispersion $E_{p}=\sqrt{\epsilon_{p}^{2}+|{\triangle}_{\mathbf{p}}|^{2}}$. Using the conjugate of Eq. (49), we can prove that $\frac{|v_{\mathbf{p}}||{\triangle}_{\mathbf{p}}|}{|u_{\mathbf{p}}|}=E_{p}-\epsilon_{p}$. Now, with the previous relation for the parameters $u_{\mathbf{p}}$ and $v_{\mathbf{p}}$, we find

[TABLE]

Finally, we can build the diagonal form

[TABLE]

Considering \big{\langle}\hat{a}_{-\mathbf{p}}\hat{a}_{\mathbf{p}}\rangle=-u_{\mathbf{p}}v_{\mathbf{p}}\big{\langle}\hat{\alpha}^{{\dagger}}_{\mathbf{p}}\hat{\alpha}_{\mathbf{p}}\big{\rangle}+u_{\mathbf{p}}v_{\mathbf{p}}\langle\hat{\alpha}_{-\mathbf{p}}\hat{\alpha}^{{\dagger}}_{-\mathbf{p}}\big{\rangle}, with \big{\langle}\hat{\alpha}^{{\dagger}}_{\mathbf{p}}\hat{\alpha}_{\mathbf{p}}\big{\rangle}=n_{F}(E_{p})=[\exp(\beta E_{p})+1]^{-1}, where $\beta=(k_{B}T)^{-1}$, we obtain the final result

[TABLE]

Appendix B Higher-order correction to the effective 2D-3D interaction

Starting with the interaction between the fermions in 2D and the “phonons" of the BEC in 3D (see the main text)

[TABLE]

where $\mathbf{q}\equiv(\mathbf{p}_{1}-\mathbf{p}_{2},q_{z})$ and

[TABLE]

In the following we calculate the four-point function to $4th$ order in the interaction constant $g_{FB}$

[TABLE]

which corresponds to the Feynman diagrams shown in Fig. 8. We consider the effective interaction between the fermions in 2D, with the free-fermion propagator given by $\mathcal{G}_{0}$

[TABLE]

The second-order expansion in the coupling $g_{FB}$ provides [$\mathbf{q}\equiv(\mathbf{k}_{1}-\mathbf{k}_{4},q_{z})$]

[TABLE]

where we used static approximation to the Bogoliubov-mode propagator $\mathcal{D}_{0}$.

B.1 RPA correction

At higher-order expansion, we obtain for the diagram in Fig. 8(b)

[TABLE]

with $\mathbf{q}_{2}=(\mathbf{k}_{2}-\mathbf{k}_{3},q_{2z})$ and $\mathbf{q}_{3}=(\mathbf{k}_{4}-\mathbf{k}_{1},q_{3z})$, which eventually leads to

[TABLE]

Now we will solve the “polarization bubble” in 2D

[TABLE]

Before we integrate in momentum space, we simplify the above expression by changing the variable in the second term to $\mathbf{p}^{\prime}=\mathbf{p}+\mathbf{k}$. We then obtain

[TABLE]

Since we are interested in the zero-temperature limit, we consider the analytic continuation $i\nu\rightarrow\nu+i\delta$, with $n_{F}(\epsilon_{\mathbf{p}})\rightarrow\Theta(\mu-\varepsilon_{\mathbf{p}})$. Then, we focus on the real part of Eq. (61)

[TABLE]

Starting with the angular integral in Eq. (62) (for $|{k}/{2k_{F}}\pm{\nu m_{F}}/{kk_{F}}|>1$), after changing the variable $p\rightarrow\varepsilon={p^{2}}/{2m_{F}}$ in the resulting integral, we obtain (see Ref.Stern (1967))

[TABLE]

remembering that $\mu={k_{F}^{2}}/{2m_{F}}$. Particularly, in the static limit $\nu=0$, we will have

[TABLE]

and

[TABLE]

Assuming $|\mathbf{k}_{1}-\mathbf{k}_{4}|<2k_{F}$, we can easily calculated the RPA series, which gives

[TABLE]

where we defined $\lambda_{0}=-{2g_{FB}^{2}n_{B}m_{B}}/{\sqrt{|\mathbf{k}_{1}-\mathbf{k}_{4}|^{2}+2\xi^{-2}}}$ and $P_{0}=-{m_{F}}/{2\pi}$. For $\lambda_{0}P_{0}<1$, we find

[TABLE]

Now, we consider the RPA correction to calculate the projected component $\ell=1$ of the potential $V_{\textrm{eff}}^{(1)}$, i.e.,

[TABLE]

and then

[TABLE]

with $Y=V_{0}\rho_{2D}\xi/\sqrt{2}$, and

[TABLE]

where $\Pi[X,Y]$ is the complete elliptic integral of the third kind. One can estimate the RPA gap correction comparing the minima in Fig. 3, which shows the profile of $\mathcal{F}(X)$ and $\mathcal{I}(X,Y)$ in a broad range of $X$ ($Y\sim 0.05/X$, since we consider $n_{B}$ as the only tunable parameter).

B.2 Phonon lifetime

The phonon lifetime ($\tau$) due to particle-hole excitation is

[TABLE]

where

[TABLE]

as shown in Fig.8,

with the polarization bubble

[TABLE]

Since we have

[TABLE]

with the Fermi energy $\mu=k_{F}^{2}/2m_{F}$, then $\tau=\infty$ for $\nu=0$ (static limit considered in the paper).

B.3 Vertex correction

We still have to deal with the $4th$-order vertex-correction in Fig. 8(c)

[TABLE]

with $\mathbf{q}_{2}=(\mathbf{k}_{2}-\mathbf{k}_{3},q_{2z})$ and $\mathbf{q}_{4}=(\mathbf{k}_{3}-\mathbf{p},q_{4z})$.

That leads to

[TABLE]

with $\omega_{\mathbf{q}}=\frac{q}{2m_{B}}\sqrt{q^{2}+2\xi^{-2}}$ and $\mathbf{q}\equiv(\mathbf{k}_{3}-\mathbf{p},q_{z})$.

B.4 Self-consistent vertex equation

Summing the ladder series as shown in Fig. 9, we derive the self-consistent vertex equation

[TABLE]

After considering $\epsilon_{\mathbf{p}}=\varepsilon_{\mathbf{p}}-\mu\sim 0$ and $\lambda=\lambda(|\mathbf{k}_{4}-\mathbf{k}_{1}|)$, again for zero external frequencies $\nu_{i}=0$, we can deal with the remaining sum

[TABLE]

with the additional external momenta constraint $\theta_{2}-\theta_{1}=\pi$ and $\theta_{4}-\theta_{3}=\pi$, and $|\mathbf{k}_{4}-\mathbf{k}_{1}|^{2}=2k_{F}^{2}[1-\cos(\theta_{1}-\theta_{4})]$.

We finally obtain the vertex correction after substituting the angular momentum expansion

[TABLE]

in Eq. (B.4), to obtain the decoupled equation for the projection $\ell=1$

[TABLE]

where $\lambda_{0}^{(1)}=V_{\textrm{eff}}^{(1)}$, as calculated in the main text, and

[TABLE]

After considering $p=k_{F}$ in the integrant of Eq. (81), we have to deal with the angular integrals

[TABLE]

with

[TABLE]

Then, from Eq. (80) we finally get

[TABLE]