P-wave superfluidity of atomic lattice fermions A.K. Fedorov,1,2,3 V.I. Yudson,4,1 and G.V. Shlyapnikov1,2,3,5,6,7 1Russian Quantum Center, Skolkovo, Moscow 143025, Russia 2LPTMS, CNRS, Univ. Paris-Sud, Universit´e Paris-Saclay, Orsay 91405, France 3Russian Quantum Center, National University of Science and Technology MISIS, Moscow 119049, Russia 4LaboratoryforCondensedMatterPhysics, NationalResearchUniversityHigherSchoolofEconomics, Moscow101000, Russia 5SPEC, CEA, CNRS, Universit´e Paris-Saclay, CEA Saclay, Gif sur Yvette 91191, France 6Van der Waals-Zeeman Institute, Institute of Physics, University of Amsterdam, Science Park 904, 1098 XH Amsterdam, The Netherlands 7 7Wuhan Institute of Physics and Mathematics, Chinese Academy of Sciences, 430071 Wuhan, China 1 (Dated: March 27, 2017) 0 2 Wediscusstheemergenceofp-wavesuperfluidityofidenticalatomicfermionsinatwo-dimensional optical lattice. The optical lattice potential manifests itself in an interplay between an increase in r a the density of states on the Fermi surface and the modification of the fermion-fermion interaction M (scattering)amplitude. Thedensityofstatesisenhancedduetoanincreaseoftheeffectivemassof atoms. Indeeplatticesthescatteringamplitudeisstronglyreducedcomparedtofreespaceduetoa 4 small overlap of wavefunctions of fermions sitting in the neighboring lattice sites, which suppresses 2 the p-wave superfluidity. However, for moderate lattice depths the enhancement of the density of states can compensate the decrease of the scattering amplitude. Moreover, the lattice setup ] significantly reduces inelastic collisional losses, which allows one to get closer to a p-wave Feshbach s resonance. Thisopenspossibilitiestoobtainthetopologicalp +ip superfluidphase,especiallyin a x y g therecentlyproposedsubwavelengthlattices. Wedemonstratethisforthetwo-dimensionalversion - of the Kronig-Penney model allowing a transparent physical analysis. t n PACS numbers: 67.85.De, 03.65.Vf, 03.67.Lx, 03.75.Ss a u q I. INTRODUCTION spin 1/2 spinless . � t a) b) a m P-wavepairingoffermionsisabasisofsuperfluidityin 3He [1], and it provides superconductivity in unconven- - d tional superconductors [2]. Presently, the p-wave super- n fluid pairing attracts a great deal of interest in ultracold o atomic gases [3–9]. One of the reasons is the search for b c topologicalp +ip superfluidofidenticalfermionsinthe c) [ x y two-dimensional(2D)geometry.Topologicalpropertiesof 2 this phase emerge from zero-energy Majorana modes on U0 v the vortex cores [10], and Non-Abelian statistics of the 0 vortices forms a basis for the implementation of topolog- 2 icallyprotectedquantuminformationprocessing[11–15]. Figure 1. Superfluid pairing of lattice fermions in various 5 setups. In(a)twocomponent(spin-1/2)latticefermionswith 8 Despite a significant progress in theory [3–9],the px+ a short-range interaction. The two spin components are la- 0 ip superfluid has not been observed. The crucial obsta- y belled by filled and unfilled circles. In (b) single component . cle to achieve this phase for spinless short-range inter- 1 (spinless) short-range interacting lattice fermions. In (c) 1D 0 acting fermions comes from a small value of the p-wave projectionofatomicfermionsloadedinthe2DKronig-Penney 7 interaction. Therefore, in order to obtain a sizable tran- lattice. 1 sition temperature one has to approach a p-wave Fesh- : bach resonance. The p-wave resonances have been stud- v i ied in experiments with fermionic potassium [16–18] and intermolecularinteraction[28,29]leadstosimilarresults X lithium [19–24] atoms. Close to the resonance the rate regarding the critical temperature as in free space [30], r of inelastic collisional losses becomes very large [25–27]. at least in subwavelength lattices. For short-range in- a Thus, the superfluid of short-range interacting atomic teracting atomic fermions the situation is different. The fermions is characterized either by vanishingly low criti- effectofthelatticepotentialontheformationofasuper- caltemperatureorbyinstabilityduetocollisionallosses. fluidphaseofatomicfermionshasbeenactivelydiscussed The creation of p +ip atomic or molecular topolog- [31–38]. In particular, for the s-wave pairing of spin-1/2 x y ical superfluids in 2D optical lattices can be a promis- fermions an increase in the depth of the optical poten- ing path for future prospects, since addressing qubits tial results in a stronger atom localization and hence in in the lattice should be much easier than in the gas increasing the on-site interaction. At the same time, the phase. For microwave-dressed polar molecules the long- tunneling becomes weaker. The combined effect of these range character of the acquired attractive dipole-dipole twofactorsisastrongincreaseinthecriticaltemperature 2 [31–33]. This has been observed in the MIT experiment II. GENERAL RELATIONS [34]. For the lattice filling somewhat smaller than unity, the physical picture can be rephrased as follows. An in- Let us first present a general framework for the inves- crease in the lattice depth increases an effective mass of tigation of superfluid pairing of weakly interacting lat- atoms and, hence, makes the density of states (DOS) tice fermions. We will do this for 2D identical (spinless) larger. The effective fermion-fermion scattering ampli- fermions, having in mind that the approach for spin-1/2 tude is also increasing. The critical temperature in the fermions is very similar. The grand-canonical Hamilto- BtioCnSalatpoprtohaechpriosdTucct∝ofexthpe[−(m1/oλdcu],luwshoefr)etλhce isscaptrtoeprionrg- pnaiarnticolef tphaertsiysstgeimvenisbHyˆ(h=erHeˆin0a+fteHrˆiwnte, pauntd(cid:126)th=e 1sinagnlde amplitude and the DOS on the Fermi surface. Thus, an set the normalization volume (surface) equal to unity): increase in the lattice potential increases T . c latOtincesth(eticgohnt-tbrainryd,infgormidoednetl)icathlefefremrmioinosn-ifnerfmaiirolny sdceaetp- Hˆ0 =(cid:90) d2rψˆ†(r)(cid:20)−2∇m2 +U(r)−µ(cid:21)ψˆ(r), (1) tering amplitude is strongly reduced. In the lowest band approach two fermions do not occupy the same lattice withµbeingthechemicalpotential,mtheparticlemass, site, and the amplitude is proportional to a very small U(r) the 2D periodic lattice potential, and ψˆ(r) the overlap of the wavefunctions of fermions sitting in the fermionic field operator. neighboring sites. This suppresses the p-wave super- The term Hˆint describes the interaction between par- fluid pairing for fairly small filling factors in deep lat- ticles: tices, which is consistent with numerical calculations of (cid:90) 1 Ref. [33]. Nevertheless, there remains a question about Hˆ = d2rd2r ψˆ (r)ψˆ (r)V(r r)ψˆ(r)ψˆ(r), (2) int (cid:48) † † (cid:48) (cid:48) (cid:48) 2 − an interplay between an increase of the DOS and the modificationofthefermion-fermionscatteringamplitude whereV(r r)isthepotentialofinterparticleinteraction (cid:48) for moderate lattice depths. However, in sinusoidal opti- − of radius r . 0 callatticessingleparticlestatesaredescribedbycompli- In the absence of interactions, fermions in the peri- catedMathieufunctions,whichcomplicatesthequestion. odic potential U(r) fill single particle energy levels ε (k) ν In this paper we study identical fermionic atoms in a determined by the Schr¨odinger equation: 2DversionoftheKronig-Penneymodelallowingatrans- (cid:20) 2 (cid:21) parentphysicalanalysisformoderatelatticedepths. The ∇ +U(r) χ (r)=ε (k)χ (r). (3) νk ν νk 2DversionoftheKronig-Penneymodelisasuperposition −2m of two Kronig-Penney potentials (in the x and y direc- Here ν = 0,1,2,... numerates energy bands, the wave tions,respectively). Withtheeigenfunctionsbeingpiece- vector k = k ,k takes values within the Brillouin wiseplanewaves, theKronig-Penneypotentialisusedin { x y} zone: π/b<k <π/b;i = x,y , and b is the lattice pe- cold atom theory (see, e.g. [39–41]) to mimic sinusoidal {− i } riod. Theeigenfunctionsχ (r)obeytheperiodicitycon- potentials of common optical lattices. In particular, this νk dition model allows us to investigate two important questions. Thefirstoneisaboutaninterplaybetweenanincreaseof χ (r+R )=χ (r)exp[ikR ], (4) νk n νk n theDOSandthemodificationofthefermion-fermionin- teractioninlatticesofmoderatedepths. Wedemonstrate where n = (nx,ny) is the index of the lattice site, with thatthereductionofthescatteringamplitudestilldomi- integer nx,ny. In the described Bloch basis the field op- nates over the enhance of the DOS. The second question erator reads: is about the stability of the system with respect to colli- (cid:88) ψˆ(r)= aˆ χ (r), (5) sional losses. We show that the lattice setup reduces in- ν,k νk νk elastic collisional losses compared to free space, and one withaˆ beingtheannihilationoperatoroffermionswith can approach the Feshbach resonance without a strong νk quasimomentum k in the energy band ν. collisional instability. This opens a possibility to observe We assume a dilute regime where the 2D density n is thelatticep +ip 2Dsuperfluidandmaybeotherinter- x y such that nb2(cid:46)1, and all fermions are in the lowest Bril- esting many-body phases. louin zone (hereinafter we omit the corresponding index The paper is organized as follows. In Sec. II we de- ν =0). In the low momentum limit (small filling factor) scribe a general approach for studying superfluidity of that we consider, their Fermi energy E is small com- F 2Dlatticefermions(Fig.1). Sec.IIIcontainsthedemon- pared to the energy bandwidth E . The lattice poten- B stration of how the ordinary tight-binding optical lattice tial amplitude U is assumed to be sufficiently large, so 0 promotes the s-wave superfluidity of spin-1/2 fermionic that both E and E are smaller than the gap between F B atomsandsuppressesthep-wavesuperfluidityofspinless the first and second lattice bands. The single particle fermions. In Sec. IV we develop a theory of p-wave su- dispersion relation then takes the form: perfluidity of spinless fermions in the 2D Kronig-Penney lattice. In Sec. V we discuss inelastic decay processes in k2 the lattice and in Sec. VI we conclude. Ek = , (6) 2m ∗ 3 where m >m is the effective mass. The renormalized gap equation for the function ∆(k) ∗ In2DthetransitionofaFermigasfromthenormalto then takes the form similar to that in free space (see superfluid state is set by the Kosterlitz-Thouless mech- Ref. [29] and references therein): anism. However, in the weakly interacting regime the Kosterlitz-Thouless transition temperature is very close ∆(k)=(cid:90) d2k(cid:48) f(k,k)∆(k) to Tc calculated in the Bardeen-Cooper-Schrieffer (BCS) (2π)2 (cid:48) (cid:48) (16) approach[42]. WethenreducetheHamiltoniangivenby (cid:26) (cid:27) 1 Eqs. (1) and (2) to the standard BCS form: (k(cid:48)) . × K − 2(Ek(cid:48) Ek) (cid:88) (cid:110) − HˆBCS = k (Ek−µ)aˆ†kaˆk In the weakly interacting regime the chemical potential 1(cid:104) (cid:105)(cid:27) (7) coincides with the Fermi energy EF = kF2/2m∗, where +2 aˆ†kaˆ†k∆(k)+h.c. , kF =√4πn is the Fermi momentum. Note that we omit − a correction to the bare interparticle interaction due to where the momentum-space order parameter ∆(k) is polarization of the medium by colliding particles [43]. given by We will see below that the scattering amplitude and the corresponding critical temperature of the superfluid (cid:88) ∆(k)= V(k,k(cid:48)) aˆ k(cid:48)aˆk(cid:48) , ∆(k)= ∆( k), (8) transition of lattice fermions depend drastically on the k(cid:48) (cid:104) − (cid:105) − − presence or absence of spin and on the pairing angular withV(k,k(cid:48))beingthematrixelementoftheinteraction momentum. Beforeanalyzingvariousregimes,wediscuss potential between the corresponding states. the situation in general. The Hamiltonian (7) is then decomposed in a set of The efficiency of superfluid pairing first of all de- independent quadratic Hamiltonians and the anomalous pends on the symmetry of the order parameter. For averages are determined by the standard BCS expres- the pairing with orbital angular momentum l we have sions: ∆(k) ∆ (k)exp[ilφ ], where φ is the angle of the l k k → vector k with respect to the quantization axis. Integrat- aˆ aˆ = ∆(k) (k), (9) (cid:104) −k k(cid:105) − K ingEq.(16)overφk andφk(cid:48) weobtainthesameequation where (k)=tanh[ (k)/2T]/2 (k), and in which ∆(k) and ∆(k(cid:48)) are replaced with ∆l(k) and K E E ∆ (k ), and f(k,k) is replaced with its l-wave part (cid:112) l (cid:48) (cid:48) (k)= (E µ)2+ ∆ (k)2 (10) k ν E − | | (cid:90) istheenergyofexcitationwithquasimomentumk. From fl(k(cid:48),k)= dφ(2kπd)φ2k(cid:48)f(k(cid:48),k)exp[ilφk−ilφk(cid:48)]. (17) Eqs.(8)and(9)wehaveanequationfor∆(k)(gapequa- tion): Alternatively, we can write (cid:88) ∆(k)= V(k,k) (k )∆(k). (11) (cid:90) (cid:48) (cid:48) (cid:48) − k(cid:48) K fl(k(cid:48),k)= d2r1d2r2Φ(lk0(cid:48))∗(r1,r2) (18) Eq. (11) can be expressed [29] in terms of the effective V(r r )Φ (r ,r ). off-shell scattering amplitude f(k(cid:48),k) of a fermion pair × | 1− 2| lk 1 2 with momenta k and k defined as − where the l-wave parts of the wavefunctions, Φ(0) and (cid:90) lk(cid:48) f(k(cid:48),k)= d2r1d2r2Φ(k0(cid:48))∗(r1,r2) (12) Φlk, are given by (cid:90) ×V(r1−r2)Φk(r1,r2). Φ(lk0(cid:48))(r1,r2)= d2φπk(cid:48)Φ(k0(cid:48))(r1,r2)exp[ilφk(cid:48)], (19) Here (cid:90) dφ Φ (r ,r )= kΦ (r ,r )exp[ilφ ]. (20) lk 1 2 k 1 2 k Φ(0)(r ,r )=χ (r )χ (r ), (13) 2π k 1 2 k 1 k 2 − is the wavefunction of a pair of non-interacting fermions As well as in free space (see Ref. [29]), we turn from with quasimomenta k and k. The quantity Φk(r1,r2) fl(k(cid:48),k) to the (real) function − is the true (i.e., accounting for the interaction) wave- function, which develops from the incident wavefunction f˜l(k(cid:48),k)=fl(k(cid:48),k)[1 itanδ(k)], (21) − Φ(0)(r ,r ) of a free pair. The wavefunction Φ (r ,r ) k 1 2 k 1 2 whereδ(k)isthescatteringphaseshift. Thisleadstothe satisfies the Schr¨odinger equation gap equation: [Hˆ 2E ]Φ (r ,r )=0, (14) 12− k k 1 2 (cid:90) d2k ∆ (k)= P (cid:48) f˜(k ,k)∆ (k ) with the two-particle Hamiltonian: l − (2π)2 l (cid:48) l (cid:48) (22) (cid:26) (cid:27) 2+ 2 1 Hˆ = ∇1 ∇2 +U(r )+U(r )+V(r r ).(15) (k(cid:48)) , 12 − 2m 1 2 1− 2 × K − Ek(cid:48) Ek − 4 where the symbol P denotes the principal value of the Using a general formula for the effective mass from integral. Ref. [44], for a deep potential of the form (27) one ob- In order to estimate the critical temperature T , we tains: c firstputk =k andnoticethatthemaincontributionto F m ξ2 (cid:20) 2 b2(cid:21) theintegraloverk(cid:48) inEq.(22)comesfromk(cid:48) closetokF. ∗ π 0 exp . (30) AttemperaturesT tendingtothecriticaltemperatureTc m0 (cid:39) b2 π2ξ02 frombelow,weputE(k(cid:48))=|Ek(cid:48)−EF|inK(k(cid:48)). Thenfor Wewillconsiderfermionicatomsinteractingwitheach thepairingchannelrelatedtotheinteractionwithorbital other via a short-range potential V(r) of radius r and angular momentum l, we have the following estimate: 0 assume the following hierarchy of length scales: (cid:20) (cid:21) 1 Tc ∼EF exp −λc , λc =ρ(kF)|fl(kF)|. (23) r0 (cid:28)ξ0 <b<1/kF. (31) The quantity ρ(k )=m /2π is the effective density of We first discuss the s-wave pairing of spin-1/2 fermions F ∗ states on the Fermi surface, and f (k ) is the on-shell with attractive intercomponent interaction (l=0). l F l-wave scattering amplitude of lattice fermions. The Turning to Eq. (18) for l=0, we notice that the main derivationforspin-1/2fermionswithattractiveintercom- contribution to the s-wave scattering amplitude in the ponent interaction leads to the same gap equations (16), lattice comes from the interaction between spin-up and (22) and estimate (23) in which spin-down fermions sitting in one and the same lattice ∆(k)=(cid:88)V(k,k(cid:48)) aˆ k(cid:48)aˆ k(cid:48) (24) site. The wavefunctions Φ(00k)(cid:48) and Φ0k can be written as k(cid:48) (cid:104) ↓− ↑ (cid:105) Φ(0)(r ,r )=χ (r )χ (r ), (32) 0k(cid:48) 1 2 0 1 0 2 andf(k,k),f (k ,k)aretheamplitudesoftheintercom- Φ (r ,r )=χ (r )χ (r )ζ (r r ), (33) (cid:48) l (cid:48) 0k 1 2 0 1 0 2 0 | 1− 2| ponent interaction. where the function ζ (r r ) is a solution of the Eq. (23) shows that compared to free space we have 0 | 1 − 2| Schr¨odinger equation for the s-wave relative motion of an additional pre-exponential factor m/m <1. Assum- ∗ twoparticlesinfreespaceatzeroenergy,anditistending ing that the lattice amplitude f (k ) and the free-space l F tounityforinteratomicseparationsgreatlyexceedingr . amplitude f0(k ) are related to each other as 0 l F We put l = 0 in Eq. (18) and integrate over r = r r 1 2 f (k )= f0(k ), (25) and r = (r + r )/2. Then, owing to the inequa−lity l F Rl l F + 1 2 r ξ , this equation is reduced to weseethattheexponentialfactorλ inEq.(23)becomes 0 (cid:28) 0 c (cid:90) (cid:90) λc =Rlmm∗λ0c, (26) f0(k(cid:48),k)= d2rV(r)ζ(r) d2r+|χ0(r+)|4. (34) where 1/λ0 is the BCS exponent in free space. Below Recalling that in the low momentum limit the free space c we compare T in various lattice setups with the critical scattering amplitude is given by c temperature in free space. (cid:90) f0 = V(r)ζ(r)d2r (35) 0 III. SHORT-RANGE INTERACTING and using Eq. (28) for the function χ (r), we obtain for FERMIONIC ATOMS IN A DEEP 2D LATTICE 0 the ratio of the lattice to free space amplitude: Westartwiththeanalysisofsuperfluidpairingindeep 1 b2 2D lattices. As an example, we consider a quadratic lat- Rl=0 = 2πξ2, (36) tice with the lattice potential of the form: 0 (cid:20) (cid:18) (cid:19) (cid:18) (cid:19)(cid:21) where we made a summation over the lattice sites and 2π 2π U(r)=U0 cos x +cos y . (27) put = 1/b2 as the normalization volume is set to be b b N unity. Thus, according to Eqs. (26) and (30) the BCS Forsufficientlydeeplattices,thesingleparticlewavefunc- exponent λ 1 becomes smaller than in free space by the −c tion has the Wannier form: following factor: 1 (cid:88) χk(r)= √ φ0(r−Rj)exp[ikRj], (28) m∗ 1exp(cid:20) 2 b2(cid:21). (37) N j Rl=0 m (cid:39) 2 π2ξ2 0 0 where isthenumberoflatticesites. Thegroundstate N For example, taking b/ξ = 4 the BCS exponent λ 1 wavefunction in the lattice cell has an extention ξ and 0 −c 0 decreases by a factor of 0.08, whereas the effective mass is given by becomes higher by a factor of 5 compared to the bare 1 (cid:20) r2 (cid:21) mass m (see Fig. 2). Then, for 6Li atoms at density 108 φ (r)= exp . (29) 0 √πξ0 −2ξ02 cm−2 (b(cid:39)250 nm, kFb(cid:39)0.5) we have the Fermi energy 5 ��0cc100 p�wave Rmm⇤0 5tshiteescjontdhiattioanrsenkse(cid:48)�bar(cid:28)wesat1vneaenigdhbro∼urrs0o(cid:28)fthξe02s/pib�te(cid:28)wi.aξvA0e,ssfuomr tinhge 10 4p-wave part of this wavefunction equation (19) at l = 1 gives: 1 3 Φ(0)(r,r ,φ )= k(cid:48)rb2 (cid:88)exp(cid:8) r2 /ξ2 b2/4ξ2(cid:9) 0.100 2 1k(cid:48) + r N8πξ04 i,j −+j 0− 0 (40) 0.010 s�wave 1 ×[exp(iφr)+exp(−iφr+2iφj)], where φ and φ are the angles of the vectors r and b r j j 0.001 0with respect to the quantization axis. The p-wave part 4.6 4.8 5.0 5.2 5.4 5.6 5.8 6.0 2.0 2.5 3.0 3.5 4.0 ofthetruerelative-motionwavefunctionΦ (r ,r )under b/⇠0 b/⇠0 k 1 2 Figur2e.02. TheratiooftheBCSexponentinthetightbinding 1.0the same conditions is given by adsitansuthhseoe1d.i5dscaaumlrlveaetdtseihcneoswittoys λtah0cn/edλBschCaoSsrtae-rxfapunongnceteincotonuinopflifntrmhegee⇤s/lstamprteatncicgeet,hλp.0ce/rTiλohcde, 0.8Φ1k(r,r+,φr)= N4bπ2ξ04ζ1(r)(cid:88)i(cid:48),j(cid:48)exp(cid:8)−r+2j(cid:48)/ξ02−b2/4ξ02(cid:9) (in units of the harmonic oscillator length ξ �)0f/o�r�the s-wave 0 c c 0.6 [exp(iφ )+exp( iφ +2iφ )]. pairin1g.0of spin-1/2 fermions, and the solid curve is λ0c/λc for × r − r j (41) the p-wave pairing of identical fermions. (2)(⌘) 0.4 F The function ζ (r) is a solution of the Schr¨odinger equa- 1 0.5 tion for the p-wave relative motion of two particles at 40 nK. Assuming that the free space BCS exponent is 0.2 ∼ energytendingtozeroinfreespace. Su(3ffi)(⌘c)ientlyfarfrom about 30 and the related critical temperature is practi- F resonance, where the on-shell scattering amplitude satis- cally0z.0e0ro,int2helatt4icewe6obtain8Tc ∼130nK.W12ethu1s4see 0.0fi0es the2 inequ4ality m6 f1(k)8 11,0the f1u2nction14ζ1(r) be- that the lattice setup may sGtrongly promote the s-wave |G | (cid:28) comes kr/2 at distances r r . 0 superfluidity of spin-1/2 fermions. (cid:29) Looking at the product of the free and true relative- The situation with p-wave superfluidity of identical motion wavefunctions we notice that the main contri- fermions is drastically different. In the single band ap- bution to the scattering amplitude (18) comes from the proximation (tight binding model) two such fermions terms in which Ri +Rj = Ri(cid:48) +Rj(cid:48), i.e. r+j = r+j(cid:48). can not occupy one and the same lattice site. This This is realized for i=i, j =j or i =j, j =i. Then, (cid:48) (cid:48) (cid:48) (cid:48) is clearly seen using the functions χ (r ) and χ (r ) k 1 k 2 recalling that for k r 1 and kr 1 the free space from Eq. (28) at the same R , so that the wavefu−nction (cid:48) 0 (cid:28) 0 (cid:28) j on-shell scattering amplitude is Φ(0)(r ,r ) becomes independent of k. Therefore, the (cid:90) k(cid:48) 1 2 (cid:48) p-wave part of this wavefunction Φ(0) and the p-wave f10(k(cid:48),k)= V(r)(k(cid:48)r/2)ζ1(r)d2r, (42) 1k(cid:48) scattering amplitude f (k ,k) following from Eqs. (19) 1 (cid:48) wefirstintegrateeachtermofthesumoveri,j,ij inthe and (18) at l=1 are equal to zero. (cid:48) (cid:48) The main contribution to the interaction amplitude product Φ(10k)(cid:48)∗Φ1k over d2r and d2r+ in Eq. (18). After then comes from the overlap of the wavefunctions of that we make a summation over the neighbouring sites j fermions sitting in the neighbouring sites. We then use andoverthesitesiandtakeintoaccountthat =1/b2. N Eqs. (28) and (29) and write: Eventually, this gives for the ratio of the lattice to free space p-wave amplitude: (cid:88)Φ(ke0(cid:48)x)p(r(cid:8)1−,r(2r)1−=Rχik)(cid:48)2(/r21ξ)02χ−−(kr(cid:48)(2r−2)R=j)(2/12/ξN02−πξi02k)(cid:48)×bj(cid:9), (38) Rl=1 = 21π (cid:18)ξb0(cid:19)6exp(cid:20)−2bξ202(cid:21). (43) i,j Thus,withthehelpofEq.(30)theinverseBCSexponent in the lattice becomes: with b =R R and R ,R being the coordinates of j j i i j thesitesiand−j. Fortheshort-rangeinteractionbetween m λ0 (cid:18) b (cid:19)4 (cid:20) cb2(cid:21) λ = ∗λ0 = c exp , (44) particles the main contribution to the scattering ampli- c Rl=1 m c 2 ξ − ξ2 0 0 tude comes from distances r ,r that are very close to 1 2 where c 0.3. each other, and for given i,j both coordinates should be (cid:39) WenowclearlyseethattheinverseBCSexponentλ in closeto(R +R )/2. Therefore, Eq.(38)isconveniently c j i the lattice is exponentially smaller compared to its value rewritten as infreespace. Inparticular,alreadyforb/ξ =5theratio (cid:110) 0 Φ(k0(cid:48))(r1,r2)=(1/Nπξ02)(cid:80)i,jexp −ik(cid:48)bj −r+2j/ξ02− λsu0cp/λercfliutidisitayboofuitde6n,twichailchferpmraiocntisca(lsleyesFuipgp.r2e)s.seHsopw-wevaevre, (cid:111) r2/4ξ02−b2/4ξ02−rbj/2ξ02 , (39) this ratio rapidly reduces with decreasing the ratio b/ξ0 and becomes 1 for b/ξ =4. It is therefore interesting 0 ∼ where r = r r , r = (r (R + R )/2, r = to analyze more carefully the case of moderate lattice 1 2 +j + i j + − − (r +r )/2, and the summation is performed over the depths. 1 2 �0 10�4cc p�wave Rmm⇤0 5 s�wave p�wave 1000 10 4 100 3 1 10 2 1 0.10 0.10 s wave 1 � 0.010.01 60 4.0 44..22 4.44.4 4.6 4.6 4.8 4.8 5.0 5.0 2.0 2.5 3.0 3.5 4.0 b/⇠ b/⇠ 0 0 IV. SUPERFLUID P-WAVE PAIRING IN THE 2.0 1.0 2D KRONIG-PENNEY LATTICE m⇤/m 0.8 1.5 Wewilldosousinga2DversionoftheKronig-Penney �0/�� c c 0.6 model,namelyasuperpositionoftwo1DKronig-Penney 1.0 lattices (in the x and y directions, respectively), with a (2)(⌘) 0.4 F δ-functional form of potential barriers: 0.5 0.2 + (3)(⌘) (cid:88)∞ F U(x,y)=U b [δ(x jb)+δ(y jb)]. (45) 0 − − 0.0 0.0 j=−∞ 0 2 4 6 G 8 10 12 14 0 2 4 6 G 8 10 12 14 Figure 3. The ratio of the BCS exponent in the 2D δ- Withtheeigenfunctionsbeingpiecewiseplanewaves,the functionalKronig-PenneylatticetotheBCSexponentinfree 1D Kronig-Penney potential is used in ultracold atom space, λ0/λδ, at the same density and short-range coupling theory (see, e.g. [39–41]) to mimic sinusoidal potentials. c c strength. The solid blue curve shows λ0/λδ as a function of The model (45) catches the key physics and allows for c c the lattice depth G, and the dashed red curve the effective transparent calculations. The latter circumstance is a massm∗/mversusG. Thedottedpartsofthesecurvesshow great advantage compared to sinusoidal lattices where our expectation at G(cid:46)1, where the single-band approxima- singleparticlestatesaredescribedbycomplicatedMath- tion used in our calculations does not work. ieufunctions. Theconsideredmodelallowsustoinvesti- gatetwoimportantquestions. Thefirstquestionisabout an interplay between an increase of the DOS and the closeto4E ,andtheratioE /E issignificantlysmaller B F B modification of the fermion-fermion interaction for mod- than unity if k b < 0.5. This justifies the single-band F erate lattice depths. The second one is the stability of approximation and the use of the quadratic dispersion the system with respect to collisional losses. relation (6). Single-particle energies in the periodic potential (45) Single-particle wavefunctions χ (r) are of the form k are represented as χ (r)=χ (x)χ (y), where k kx ky Ek =E(kx)+E(ky), (46) √2sin(η/2) j=(cid:88)+∞ χ (x)= A (x)exp[ik jb] where E(kx,y) > 0 is the dispersion relation for the 1D kx (cid:112)1+sinη/η j= j x (50) Kronig–Penneymodel. Itfollowsfromtheequation(see, −∞ (cid:26)eiqbeiq(x jb) e iqbe iq(x jb)(cid:27) e.g., Ref. [44]): − − − − × eiqb−eikxb − e−iqb−eikxb sin(qb) cos(qb)+G =cos(kb), (47) qb is the exact eigenfunction of the 1D Kronig-Penney model, with A (x) = 1 for (j 1)b < x < jb and zero j (cid:112) − where q= 2mE(k) > 0, and G = mU0b2. As well as in otherwise. The function χky(y) has a similar form. For the previous section, we consider a dilute regime where k(cid:48)b 1 and kb 1 the p-wave parts of the wavefunc- (cid:28) (cid:28) the filling factor is ν = nb2(cid:46)1 and the fermions fill only tions, Φ(0) and Φ , following from Eqs. (13), (19), and 1k(cid:48) 1k a small energy interval near the bottom of the lowest (20) at l=1 turn out to be Brillouinzone. Thentheenergycountedfromthebottom of the zone is given by Eq. (6) and for the effective mass Φ(0) =ik r ηcot(η/2) (cid:88)∞ A (x )A (y ) (51) Eq. (47) yields: 1k(cid:48) (cid:48) [1+sinη/η]2 jx + jy + m∗ tan(η/2)(cid:20)1+ sinη(cid:21), (48) (cid:26)cosφ cos2(cid:18)q y j b+bj(cid:19)x,+jyi=s−in∞φ cos2(cid:18)q x j b+b(cid:19)(cid:27), m ≈ η η × r 0 +− y 2 r 0 +− x 2 with η being the smallest root of the equation: ηcot(η/2) (cid:88)∞ Φ =2iζ (r) A (x )A (y ) (52) 1k 1 [1+sinη/η]2 jx + jy + ηtan(η/2)=G. (49) jx,jy=−∞ (cid:26) (cid:18) (cid:19) (cid:18) (cid:19)(cid:27) b b cosφ cos2q y j b+ +isinφ cos2q x j b+ , Actually,η =q0bwhereq0followsfromEq.(47)atk =0. × r 0 +− y 2 r 0 +− x 2 For m m we have m /m = G/π2, which means ∗ ∗ (cid:29) thatthequantityGshouldbeverylarge. Thenthewidth where the function ζ (r) is defined after equation (41). 1 of the lowest Brillouin zone is E = 2/m b2 and it is Fortheratioofthelatticetofreespacescatteringampli- B ∗ much larger than the Fermi energy E = k2/2m for tude we then obtain: F ∗ k b<0.5. Thegapbetweenthelowestandsecondzones isFE = 3π2/2mb2 and it greatly exceeds E and E . η2cot2(η/2) (cid:20)3 2sinη sin2η(cid:21) G B F = + + , (53) Note that even for m∗ 1.3m (G 5) we have EG Rl=1 [1+sinη/η]4 2 η 4η (cid:39) (cid:39) �0 10�4cc p�wave Rmm⇤0 5 s�wave p�wave 1000 10 4 100 3 1 10 2 1 0.10 0.10 s wave 1 � 0.010.01 0 7 4.0 44..22 4.44.4 4.6 4.6 4.8 4.8 5.0 5.0 2.0 2.5 3.0 3.5 4.0 b/⇠ b/⇠ 0 0 2.0 1.0 withHˆ beingtheHamiltonianofelasticinteraction,and 0 m⇤/m 0.8 (cid:90) 1.5 Hˆ (0)= d(cid:126)r d(cid:126)r V ((cid:126)r (cid:126)r ) (cid:48) 1 2 r 1 2 − �0/�� (57) c c 0.6 (cid:104) (cid:105) 1.0 Φˆ†((cid:126)r2)Φˆ†((cid:126)r1)ψˆ((cid:126)r1)ψˆ((cid:126)r2)+h.c. . (2)(⌘) × 0.4 F Here (cid:126)r and (cid:126)r are the 3D coordinates of the atoms, 0.5 1 2 0.2 (3)(⌘) ψˆ((cid:126)r) is the field operator of the initial-state atoms, and F V ((cid:126)r (cid:126)r ) is the interaction potential causing the in- r 1 2 0.00 2 4 6 8 10 12 14 0.00 2 4 6 8 10 12 14 elastic−relaxation. The field operator of atoms in the G G final (ground) internal state is Figure 4. Coefficients F(2) and F(3) as functions of the lattice depth G. The dotted parts of the curves show our (cid:88) expectation at G (cid:46) 1, where the single-band approximation Φ((cid:126)r)= aˆ(cid:126)qexp(i(cid:126)q(cid:126)r), (58) (cid:126)q used in our calculations does not work. andinitiallythesestatesarenotoccupied. Wethushave: atincde uissienxgpEreqs.se(d48t)htrhoeuginhvetrhseeiBnCveSrseexpBoCnSenetxipnotnheentlaitn- W2 =(cid:82)−∞∞dt(cid:82) d(cid:126)r1d(cid:126)r2dr(cid:126)(cid:48)1dr(cid:126)(cid:48)2Vr((cid:126)r1−(cid:126)r2)Vr(r(cid:126)(cid:48)1−r(cid:126)(cid:48)2) free space as ×exp{i(cid:126)q1((cid:126)r1−r(cid:126)(cid:48)1)+i(cid:126)q2((cid:126)r2−r(cid:126)(cid:48)2)−i[2E0−(q12+q22)/2m]t} ψˆ ((cid:126)r ,0)ψˆ ((cid:126)r ,0)ψˆ(r(cid:126) ,t)ψˆ(r(cid:126) ,t) . (59) λδ = m∗λ0 = ×(cid:104) † 1 † 2 (cid:48)2 (cid:48)1 (cid:105) c Rl=1 m c The momenta q and q are large but the center of (cid:20) (cid:21) (54) 1 2 ηcot(η/2) 3 2sinη sin2η mass momentum (cid:126)q +(cid:126)q is almost zero. The energy = + + λ0. | 1 2| [1+sinη/η]3 2 η 4η c conservation law then reads: IntheextremelimitofG 1wehaveη (π 2π/G), p2 2E = , (60) so that π4/G2 and λ(cid:29)δ/λ0 π2/G (cid:39) 1.−We thus 0 m Rl=1 (cid:39) c c (cid:39) (cid:28) arrive at the same conclusion as in the previous section where (cid:126)p=((cid:126)q (cid:126)q )/2 is the relative momentum. From for sinusoidal lattices: in a very deep lattice the p-wave 1− 2 thesummationover(cid:126)q ,(cid:126)q weturntotheintegrationover pairingofidenticalfermionsissuppressed. However,even 1 2 (cid:126)p and ((cid:126)q +(cid:126)q ). The coordinate-dependent part of the for G 20 the BCS exponent in the lattice exceeds the 1 2 (cid:39) exponent in Eq. (59) takes the form: exponentinfreespaceonlybyafactorof1.7atthesame density and short-range coupling strength (see Fig. 3). exp i(cid:126)q ((cid:126)r r(cid:126) )+i(cid:126)q ((cid:126)r r(cid:126) ) It is thus crucial to understand what happens with the { 1 1− (cid:48)1 2 2− (cid:48)2 } (61) rates of inelastic decay processes in the lattice setup. =exp i((cid:126)q +(cid:126)q )(R(cid:126) R(cid:126) )+i(cid:126)p((cid:126)r r(cid:126)) , { 1 2 − (cid:48) − (cid:48) } where R(cid:126) = ((cid:126)r +(cid:126)r )/2, R(cid:126) = (r(cid:126) +r(cid:126) )/2,(cid:126)r =(cid:126)r (cid:126)r , V. INELASTIC DECAY PROCESSES 1 2 (cid:48) (cid:48)1 (cid:48)2 1− 2 and r(cid:126) =r(cid:126) r(cid:126) . The integration over ((cid:126)q +(cid:126)q ) yields: (cid:48) (cid:48)1− (cid:48)2 1 2 We first consider the two-body relaxation, assuming (cid:90) d((cid:126)q +(cid:126)q ) 1 2 exp i((cid:126)q +(cid:126)q )(R(cid:126) R(cid:126) ) that both colliding atoms are in an excited (internal en- (2π)3 { 1 2 − (cid:48) } (62) ergy E ) hyperfine state and they relax to the ground 0 =δ(R(cid:126) R(cid:126) ), state. The released hyperfine-state energy 2E0 goes to − (cid:48) the kinetic energy of the atoms. It greatly exceeds the and the correlation function becomes: Fermienergyandthelatticepotentialdepth, sothatthe relativemotionoftheatomsinthefinalstateisdescribed ψˆ ((cid:126)r ,0)ψˆ ((cid:126)r ,0)ψˆ(r(cid:126) ,t)ψˆ(r(cid:126) ,t) = (63) by a 3D plane wave with a high momentum and they es- (cid:104) † 1 † 2 (cid:48)2 (cid:48)1 (cid:105) cape from the system. Then the number of relaxation (cid:104)ψˆ†(R(cid:126) +r(cid:126)(cid:48)/2)ψˆ†(R(cid:126) −r(cid:126)(cid:48)/2)ψˆ(R(cid:126) −(cid:126)r/2,t)ψˆ(R(cid:126) +(cid:126)r/2,t)(cid:105). events per unit time can be written in the form (see, e.g. Characteristic times t on which the correlation function Ref. [45]): changes are of the order of the inverse Fermi energy or W2 =(cid:90) ∞ dt(cid:88)ρi Hˆ(cid:48)(0)Hˆ(cid:48)(t) , (55) tehveantdlaormgeinr.atTehtehyeainretemgruaclhovloenrgdetrinthEanq.t(h5e9)t.imTehse∼refEo0r−e1, (cid:104)| |(cid:105) −∞ i we may put t = 0 in the correlation function, which re- duces Eq. (59) to whereρ istheequilibriumdensitymatrix, andHˆ isthe i (cid:48) Hamiltonian responsible for the relaxation process: W =(cid:82) W˜ ((cid:126)r,r(cid:126))dR(cid:126)d(cid:126)rdr(cid:126) 2 2 (cid:48) (cid:48) Hˆ (t)=exp iHˆ t Hˆ (0)exp iHˆ t , (56) ψˆ (R(cid:126) +(cid:126)r/2)ψˆ (R(cid:126) (cid:126)r/2)ψˆ(R(cid:126) (cid:126)r/2)ψˆ(R(cid:126) +(cid:126)r/2) ,(64) (cid:48) 0 (cid:48) 0 † (cid:48) † (cid:48) { } {− } ×(cid:104) − − (cid:105) 8 with have been already extracted from the wavefunctions. In the leading (linear) order in small r we have: (cid:90) (cid:18) p2(cid:19) d(cid:126)p W˜2= Vr((cid:126)r)Vr(r(cid:126)(cid:48))exp{i(cid:126)p((cid:126)r−r(cid:126)(cid:48))}δ 2E0−m (2π)2.(65) (cid:16) r(cid:17) (cid:26) 1 (cid:27) χ R χ (R) 1 r ln[χ (R)] , (72) k k R k ± 2 ≈ ± 2 ·∇ In the quasi-2D geometry the field operator can be writtenasψˆ((cid:126)r )=ψ (z )ψ(r ),wherer isthe2D where the Slater determinant takes the form: 1,2 0 1,2 1,2 1,2 vector in the x,y plane, and (r,R;k ,k )=χ (R)χ (R) D 1 2 k1 k2 (73) 1 (cid:34) z2 (cid:35) ×r·∇R{ln[χk1(R)]−ln[χk2(R)]}. ψ (z )= exp 1,2 (66) 0 1,2 (πl2)1/4 −2l2 Astheleadingcontributiontothescatteringofslowiden- 0 0 tical fermions comes from the p-wave scattering channel, the expression in the curly brackets in Eq. (73) is linear is the wavefunction in the tightly confined z-direction (in the leading order) in the difference (k k ). For [46]. As the inelastic relaxation occurs at interparticle 1 − 2 instance, the “x-component” of this expression has the distances much smaller than the confinement length l , 0 form: the product of four field operators in Eq. (64) becomes ln[χ (X)] ln[χ (X)]= ψˆ†(R+r(cid:48)/2)ψˆ†(R r(cid:48)/2) i(k k1x k )b− sink2[xq(X nb)] − (67) 1x− 2x − , (74) ψˆ(R r/2)ψˆ(R+r/2)ψ4(Z), 2sin(η/2) cos[q(X nb)+b/2)] × − 0 − where X varies from (n 1)b to nb (see Eq. (50)). In where R,r and r(cid:48) are 2D vectors in the x,y plane, and the considered low densit−y limit (kb 1) we may put Z =(z1+z2)/2. Integrating over Z in Eq. (64) we then k = k = 0 in the product χ (r )χ(cid:28)(r ) in Eq. (73). have: 1 2 k1 1 k2 1 As a result we transform Eq. (68) to (cid:82) W = w˜ (r,r)dRdrdr (cid:90) 2 2 (cid:48) (cid:48) (η) ψˆ†(R+r(cid:48)/2)ψˆ†(R r(cid:48)/2)ψˆ(R r/2)ψˆ(R+r/2) ,(68) W2=F22 rr(cid:48)w2(r,r(cid:48))drdr(cid:48) ×(cid:104) − − (cid:105) (cid:90) dk dk (75) where × Nk1Nk2(k12+k22) (21π)42. (cid:90) dzdz w˜2(r,r(cid:48))= W˜2((cid:126)r,r(cid:126)(cid:48))√2πl(cid:48)0, (69) T2Dhelaqtutiacnet:ity F2 is determined by the integral over the (cid:90) (cid:90) witUhsizng=ezx1pa−nzsi2onan(d5)zo(cid:48)n=ecza1(cid:48)n−exz2p(cid:48).resstheaveragedprod- F2 =8sin12η/2(cid:88)∞ ∞ ∞ dxdyAn(x)Am(y) uct of four 2D field operators in terms of the standard n,m −∞ −∞ Slater determinants (r,R;k ,k ): χ0(x)4 χ0(y)4(cid:2)P2(x)+P2(y)(cid:3), (76) 1 2 ×| | | | D where the functions A (x) are defined below Eq. (50), ψˆ (R+r/2)ψˆ (R r/2)ψˆ(R r/2)ψˆ(R+r/2) j † (cid:48) † (cid:48) andthefunctionP resultsfromthedifferentiationofthe (cid:104) − − (cid:105) 1 (cid:88) curly brackets in Eq. (73) with the use of Eq. (74): = N N (r,R;k ,k ) (r,R;k ,k ),(70) 2! k1 k2D∗ (cid:48) 1 2 D 1 2 k1,k2 P(u) d sin[qu] . (77) ≡ du cos[q(u+b/2)] where N is the Fermi distribution function, and k Performing the integration in Eq. (76) we find (r,R;k ,k ) 1 2 D=Det(cid:18)χk1(R+r/2) χk1(R−r/2)(cid:19). (71) F2(η)=Rl=1(η), (78) χk2(R+r/2) χk2(R−r/2) with the lattice factor Rl=1 given by Eq. (53). In the absence of the 2D lattice (i.e. in free 2D space) The distance r between relaxing particles is small com- wealsoarriveatEq.(68). Then,usingχ (r)=exp(ikr), k pared to the lattice period and particle wavelengths. the Slater determinant becomes: Therefore,allthewavefunctionsenteringEq.(71)should (r,R;k ,k ) i(k k )rexp[i(k +k )R].(79) be taken within the same lattice cell (n,m) of the con- 1 2 1 2 1 2 D (cid:39) − sidered 2D lattice, so that Eq. (68) will contain only one PerformingintegrationswegetEq.(75)with replaced 2 double lattice summation over n and m. The Slater de- F by unity. Thus, we obtain that in the lattice the two- terminant (71) within a given cell (n,m) contains a fac- body inelastic relaxation is reduced by a factor of 2 tor exp[i(k +k )bn+i(k +k )bm] [see Eq. (50)], F 1x 2x 1y 2y compared to free space: which does not contribute to the product . Below ∗ we will imply that the corresponding exponeDntDial factors Wlat = (η)Wfree. (80) 2 F2 2 9 The function (η) following from Eqs. (53) and (78) is (56) and Hˆ (0) (81)) takes the form: 2 (cid:48) F displayed in Fig. 4 versus the lattice depth G, which is (cid:90) (cid:90) related to η by Eq. (49). W = ∞ dt dR(cid:126)dR(cid:126) d(cid:126)rd(cid:126)rd(cid:126)ud(cid:126)uV((cid:126)r,(cid:126)u)V((cid:126)r,(cid:126)u) 3 (cid:48) (cid:48) (cid:48) (cid:48) (cid:48) We complete this section with the discussion of three- (cid:88) −(cid:110)∞ (cid:104) (cid:105)(cid:111) body recombination, assuming that the binding energy exp i (cid:126)p(R(cid:126) R(cid:126) +(cid:126)u (cid:126)u)+(cid:126)q(R(cid:126) R(cid:126) ) (cid:48) (cid:48) (cid:48) of the molecule formed in this process greatly exceeds − − − − (cid:126)p,(cid:126)q,s the Fermi energy and the lattice depth. In this case the (cid:20) (cid:18) p2 q2 (cid:19) (cid:21) kinetic energies of the molecule and atom in the output χ ((cid:126)r)χ ((cid:126)r)exp i + E t channel of the recombination are very high and they es- × s (cid:48) ∗s 2m 4m − s cape from the system. The results for the ratio of the ψˆ (R(cid:126) +(cid:126)r/2,0)ψˆ (R(cid:126) (cid:126)r/2,0)ψˆ (R(cid:126) +(cid:126)u,0) † (cid:48) (cid:48) † (cid:48) (cid:48) † (cid:48) (cid:48) three-body recombination rate in the lattice to the rate ×(cid:104) − ψˆ (R(cid:126) +(cid:126)u,t)ψˆ(R(cid:126) (cid:126)r/2,t)ψˆ(R(cid:126) +(cid:126)r/2,t) , (85) in free space are obtained in a way similar to that for × † − (cid:105) the two-body relaxation. The number of recombination where V(r(cid:126),u(cid:126) ) V(R(cid:126) +(cid:126)r/2,R(cid:126) (cid:126)r/2,R(cid:126) +(cid:126)u) and (cid:126)u= events per unit time, W3, is given by Eq. (55) in which (cid:126)r R(cid:126). O(cid:48) m(cid:48)itt≡ing a small diff−erence between q and Hˆ(cid:48)(t) follows from Eq. (56), and the Hamiltonian Hˆ(cid:48)(0) p3 i−n the time-dependent exponent transforms it to is given by exp[i(3p2/4m E )t] and after putting t=0 in the cor- s − relationfunctiontheintegrationovertyieldsδ(3p2/4m (cid:90) − Hˆ(cid:48)(0)= d(cid:126)r1d(cid:126)r2(cid:126)r3V((cid:126)r1,(cid:126)r2,(cid:126)r3)× Es). Thesummationover(cid:126)qgivesδ(R(cid:126) −R(cid:126)(cid:48)). Asaresult, (81) Eq. (85) reduces to (cid:104) (cid:105) Bˆ ((cid:126)r ,(cid:126)r )ψˆ ((cid:126)r )ψˆ((cid:126)r )ψˆ((cid:126)r )ψˆ((cid:126)r )+h.c. † 1 2 † 3 3 2 1 (cid:90) W = W˜ ((cid:126)r,(cid:126)r,(cid:126)u,(cid:126)u)dR(cid:126)d(cid:126)rd(cid:126)rd(cid:126)ud(cid:126)u 3 3 (cid:48) (cid:48) (cid:48) (cid:48) × with V((cid:126)r1,(cid:126)r2,(cid:126)r3) being the sum of three pair interac- ψˆ (R(cid:126) +(cid:126)r/2)ψˆ (R(cid:126) (cid:126)r/2)ψˆ (R(cid:126) +(cid:126)u) (86) tion potentials, and Bˆ ((cid:126)r ,(cid:126)r ) the field operator of the (cid:104) † (cid:48) † − (cid:48) † (cid:48) × † 1 2 molecules. The latter can be written as ψˆ(R(cid:126) +(cid:126)u)ψˆ(R(cid:126) (cid:126)r/2)ψˆ(R(cid:126) +(cid:126)r/2) , − (cid:105) Bˆ†((cid:126)r1,(cid:126)r2)=(cid:88)(cid:126)q,sexp−i(cid:126)qR(cid:126) χ∗s((cid:126)r)ˆb(cid:126)†qs, (82) with (cid:90) d(cid:126)p W˜ ((cid:126)r,(cid:126)r,(cid:126)u,(cid:126)u)=V((cid:126)r,(cid:126)u)V((cid:126)r,(cid:126)u) 3 (cid:48) (cid:48) (cid:48) (cid:48) (2π)2 where ˆb(cid:126)†qs is the creation operator of the molecule with (cid:88) (cid:18)3p2 (cid:19) momentum (cid:126)q in the internal state s, χs((cid:126)r) is the wave- ×exp[i(cid:126)p((cid:126)u−(cid:126)u(cid:48))] δ 4m−Es χ∗s((cid:126)r)χs((cid:126)r(cid:48)).(87) function of this state, and the notations for coordinates s are the same as in the above discussion of two-body re- Integratingoutthemotionofparticlesinthetightlycon- laxation. finedz-directioninawaysimilartothatforthetwo-body Initially molecules are not present in the system and, relaxation, we transform Eq. (86) to hence,fortheaverageofthemolecularfieldoperatorswe (cid:90) have: W = dRdrdrduduw˜ (r,r,u,u) ψˆ (R+r/2)ψˆ (R r/2) 3 (cid:48) (cid:48) 3 (cid:48) (cid:48) † (cid:48) † (cid:48) (cid:104) − (cid:104)Bˆ((cid:88)(cid:126)r(cid:48)1,(cid:126)r(cid:48)2,(cid:26)0)Bˆ†((cid:126)r1,(cid:126)r2,t)(cid:105)=χ(cid:18)s((cid:126)rq2(cid:48))χ∗s((cid:126)r)(cid:19) (cid:27) (83) ×ψˆ†(R+u(cid:48))ψˆ(R+u)ψˆ(R−r/2)ψˆ(R+r/2)(cid:105), (88) exp i(cid:126)q(R(cid:126) R(cid:126) )+i E t , where R,r,u and R,r,u are 2D vectors in the x,y (cid:48) s (cid:48) (cid:48) (cid:48) × − − 4m − (cid:126)q,s plane and (cid:90) dzdz du du wstiattheEs;sR(cid:126)be=in(g(cid:126)rth+e(cid:126)rbi)n/d2i;n(cid:126)rg=en(cid:126)rergy(cid:126)rof(atnhde msimolielcaurllyefionrtR(cid:126)he w˜3(r,r(cid:48),u,u(cid:48))= W˜3((cid:126)r,(cid:126)r(cid:48),(cid:126)u,(cid:126)u(cid:48)) √(cid:48)3πlz02 (cid:48)z.(89) 1 2 1 2 (cid:48) − and(cid:126)r(cid:48)). The momentum p of the atom in the outgoing Similarly to Eq. (70), the averaged product of six recombination channel is very high, and the states with fermionic field operators is represented as such momenta are not initially occupied. Therefore, we get S ψˆ (r )ψˆ (r )ψˆ (r )ψˆ(r )ψˆ(r )ψˆ(r ) 3 ≡(cid:104) † (cid:48)1 † (cid:48)2 † (cid:48)3 3 2 1 (cid:105) 1 (cid:88) = N N N (90) ψˆ((cid:126)r ,0)ψˆ ((cid:126)r ,t) 3! k1 k2 k3 (cid:104) (cid:48)3 =†(cid:88)3 ex(cid:105)p(cid:26)−i(cid:126)p((cid:126)r3−(cid:126)r(cid:48)3)+i2pm2 t(cid:27). (84) ×D∗(kr1(cid:48)1,k,r2,(cid:48)2k,1r(cid:48)3;k1,k2,k3)D(r1,r2,r3;k1,k2,k3), (cid:126)p where (r ,r ,r ;k ,k ,k ) is the (Slater) determinant 1 2 3 1 2 3 D ofthe3 3matrix χ (r ) . Usingtheexpansionofthe Thus, the initial expression for W (Eq. (55) with Hˆ (t) wavefun×ctions in (s{maklil) rjel}ative coordinates r=r r 3 (cid:48) 1 2 − 10 andu=r (r +r )/2wefindthat isbilinearinthe Using (97) in Eqs. (90) and (88) we arrive at the re- 3 1 2 − D D components of these quantities: combination rate given by Eq. (94) without the factor in the right hand side. Thus, the relation between D(R+r/2,R−r/2,R+u;k1,k2,k3) Fth3e recombination decay rate in free space and the one 1 (cid:88) = χ (R)χ (R)χ (R) (r u r u ) in the 2D lattice reads: 2 k1 k2 k3 α β − β α α,β Wlat =Wfree (η). (98) ×{∇α[lnχk1(R)−lnχk2(R)]∇β[lnχk2(R)−lnχk3(R)] 3 3 F3 −∇β[lnχk1(R)−lnχk2(R)]∇α[lnχk2(R)−lnχk3(R)]},(91) The function F3 is shown in Fig. 4 versus the lattice depth G related to η by Eq. (49). The results obtained where α,β = x,y . Using Eqs. (74) and (77), in the { } in this section indicate that both two-body and three- leading order in small relative wavevectors equation (91) body inelastic collisions are significantly suppressed in takes the form: the lattice setup even at moderate depths. [χ (R)]3b2 (R+r/2,R r/2,R+u;k ,k ,k ) 0 For usual sinusoidal optical lattices used in experi- D − 1 2 3 (cid:39) 4sin2(η/2) ments with ultracold atoms, one can proceed along the (cid:88) samelinesasinthecaseofthe2DKronig-Penneymodel. (r u r u )(k k ) (k k ) P(R )P(R ).(92) α β β α 1 2 α 3 2 β α β × − − − For the two-body relaxation, Eqs. (55)-(71) remain the α,β same. Then, for fairly deep lattices (b/ξ (cid:38) 4) where 0 Substituting the result of Eq. (92) into equation (90) we the function χ (r) can be still used in the form (28), k find for the correlation function: we obtain the ratio of the lattice to free space relax- S3 = 2|7χs0i(nR4)(|η6/b24) agtivioenn brayteEqW. 2(s4l/3W).2frTeehe(cid:39)caRlcul=la1t,iownisthforthtehefatchtroere-bRold=y1 recombination are more involved. The estimate using 1(cid:90) d2k d2k d2k N N N [k2k2+k2k2+k2k2] 1 2 3 an analogy with the Kronig-Penney model at large G, ×3 k1 k2 k3 1 2 1 3 1 3 (2π)6 leads to the ratio of the lattice to free space recombina- ×(cid:88)(rαuβ −rβuα)(rα(cid:48) u(cid:48)β −rβ(cid:48)u(cid:48)α)P2(Rα)P2(Rβ).(93) t(mion/rmate W5)3slt/hWe3tfrweeo∼-bRo2ld=y1.relIanxaptairotniciuslasru,pfporresbs/eξd0b=y 4a ∗ α,β (cid:39) factor of 5 and the three-body recombination by about a Having in mind that only the terms with β = α con- factor of 25. (cid:54) tributetothesummationover2DCartesianindices,from Eqs. (88), (90), and (93) we obtain for the decay rate: (η)(cid:90) VI. CONCLUSIONS AND OUTLOOK W3 = F3 drdr(cid:48)dudu(cid:48)w˜3(r,r(cid:48),u,u(cid:48))[(cid:126)r (cid:126)u]z[(cid:126)r(cid:48) (cid:126)u(cid:48)]z 12 × × (cid:90) d2k d2k d2k Theresultsofthepresentpaperindicatethatthereare × Nk1Nk2Nk3[k12k22+k22k32+k12k32] 1(2π)26 3, (94) possibilities to create the superfluid topological px+ipy phase of atomic lattice fermions. In deep lattices the p- whereweexpressedthecombinationr1u2 r2u1 interms wave superfluid pairing is suppressed and even for mod- − of 3D vectors(cid:126)r and(cid:126)u. The quantity 3(η) in Eq. (94) is erate lattice depths the BCS exponent is larger than in F given by free space at the same density and short-range coupling b2 (cid:90) strength. However,thelatticesetupsignificantlyreduces (η)= dRχ (R)6P2(X)P2(Y).(95) the inelastic collisional losses, so that one can get closer F3 16sin4(η/2) | 0 | to the p-wave Feshbach resonance and increase the in- HeretheintegrationoverRisonlyinthe2Dlatticecell, teraction strength without inducing a rapid decay of the while the summation over all lattice cells resulted in the system. multiplicationoftheresultbythecellnumber1/b2. Per- For ultracold 6Li the p-wave resonance is observed for forming the integration we obtain: atoms in the lowest hyperfine state (1/2,1/2) [19–24], andtheonlydecaychannelisthree-bodyrecombination. (cid:20) (cid:21) 4 3(η)=η4cot4(η/2) 1+ sinη − . (96) In the 2D Kronig-Penney lattice with the depth G (cid:39) 12 F η andtheperiodb 200nm(m /m 2and m /m ∗ l=1 ∗ (cid:39) (cid:39) R ≈ 0.7), at k b 0.5 the Fermi energy is close to 100 nK F LetusnowcomparetheresultofEq.(94)withthatin and the 2D d(cid:39)ensity is about 0.5 108 cm 2. Slightly − free space. Taking the wavefunction χ = exp(ikr) the × k awayfromtheFeshbachresonance(atthescatteringvol- expression for becomes: ume V 8 10 15 cm3) we are still in the weakly D sc − (cid:39) × (R+r/2,R r/2,R+u;k ,k ,k ) interacting regime, and the 3D recombination rate con- 1 2 3 D − stant is α3D 10 24 cm6/s [19]. Then, using Eq. (23) =exp[i(k1+k2+k3)R] (97) andthequreacsi∼2Dsc−atteringamplitudeexpressedthrough (cid:88) (cid:112) i (rαuβ rβuα)(k1 k2)α(k3 k2)β. Vsc and the tight confinement length l0 = 1/mω0 [47], × − − − for the confinement frequency ω 100 kHz we obtain α,β 0 (cid:39)