Table Of Content

Novel p-wave superfluids of fermionic polar molecules A.K. Fedorov,1,2 S.I. Matveenko,3,2 V.I. Yudson,1,4 and G.V. Shlyapnikov1,2,5,6 1Russian Quantum Center, Skolkovo, Moscow Region 143025, Russia 2LPTMS, CNRS, Univ. Paris-Sud, Université Paris-Saclay, Orsay 91405, France 3L.D. Landau Institute for Theoretical Physics, Russian Academy of Sciences, Moscow 119334, Russia 4Institute for Spectroscopy, Russian Academy of Sciences, Troitsk, Moscow 142190, Russia 5Van der Waals-Zeeman Institute, University of Amsterdam, Science Park 904, 1098 XH Amsterdam, The Netherlands 6Wuhan Institute of Physics and Mathematics, Chinese Academy of Sciences, Wuhan 430071, China (Dated: July 12, 2016) 6 Recentlysuggestedsubwavelengthlatticesofferremarkableprospectsfortheobservationofnovel 1 superfluids of fermionic polar molecules. It becomes realistic to obtain a topological p-wave super- 0 fluid of microwave-dressed polar molecules in 2D lattices at temperatures of the order of tens of 2 nanokelvins, which is promising for topologically protected quantum information processing. An- otherforeseennovelphaseisaninterlayerp-wavesuperfluidofpolarmoleculesinabilayergeometry. l u J 1 Introduction z 1 E dc Non-conventional superconductors and superfluids at- ] s tract a great deal of interest due to their non-trivial y a transport properties and/or topological behavior [1–11]. g - This behavior has been actively discussed in two di- nt mensions (2D) for the px + ipy superfluid of identical x Eac L a fermions, where Cooper pairs have orbital angular mo- u mentum equal to unity [12–17]. Quantized vortices in q this superfluid carry zero-energy Majorana modes on . t their cores [3, 18, 19]. These modes cause the vortices a m to obey non-Abelian exchange statistics, which is a basis fortopologicallyprotectedquantuminformationprocess- - d ing [20, 21]. However, the p-wave topological superfluid n of ultracold atoms is either collisionally unstable near a (a) (b) o Feshbach resonance, or has a vanishingly low superfluid c Figure 1. Setups for p-wave superfluids of polar molecules: [ transition temperature far from the resonance [22–24]. (a)polarmoleculeinanexternalmicrowavefieldE rotating ac Successful experiments on the creation of ground- in the plane perpendicular to the stationary field E (upper 2 dc state ultracold polar molecules [25–36] opened fascinat- part),andmicrowave-dressedpolarmoleculesloadedina2D v 6 ing prospects for obtaining non-conventional superfluids lattice(lowerpart);(b)bilayersystemofpolarmoleculeswith 2 [37–39]. Inparticular,microwave-dressedpolarmolecules dipolemomentsintheupperandlowerlayers,orientedoppo- 0 confined to 2D may acquire an attractive dipole-dipole site to each other. 3 tail in the interaction potential, which ensures the emer- 0 gence of collisionally stable p-wave superfluid with a . where the lattice constant (interlayer spacing in the bi- 1 reachable transition temperature [16, 17]. Another in- 0 layer system) can be as small as about 50 nm. An interesting system concerns fermionic polar molecules in a 6 crease of energy scales in such lattices makes it realistic bilayer geometry. Here they may form interlayer super- 1 toobtainsizeabletransitiontemperaturesoftheorderof : fluids in which Cooper pairs consist of molecules belong- v tens of nanokelvins. ing to different layers [40–43]. i X It this paper we consider novel p-wave superfluids of r fermionic polar molecules in 2D lattice geometries (see a General relations and qualitative arguments Fig. 1). It is shown that a collisionally stable topological p + ip superfluid of identical microwave-dressed x y polarmoleculesmayemergeina2Dlatticeduetoalong- The superfluid pairing of identical fermions is char- range character of the dipole-dipole interaction. We also acterized by the order parameter ∆(r,r(cid:48)) = V(r−r(cid:48)) show how one can get a p-wave interlayer superfluid of ×(cid:104)ψˆ(r(cid:48))ψˆ(r)(cid:105), where V(r − r(cid:48)) is the interaction po- fermionic polar molecules in a bilayer geometry, which tential, the symbol (cid:104)...(cid:105) denotes the statistical average, can be a quantum simulator of superconductivity in lay- and ψˆ(r) is the field operator of fermions. For spin-1/2 eredcondensedmattersystems[7,8]. Itiscrucialtorely fermions one of the field operators in the expression for on the recently proposed subwavelength lattices [44–54], ∆(r,r )isforspin-↑fermions,andtheotheroneforspin- (cid:48) 2 ↓ fermions. In free space the order parameter depends in the tight binding model, if the extension of the par- on the coordinates r and r only through the difference ticle wavefunction in the lattice site greatly exceeds the (cid:48) (r−r ). In 2D the transition temperature T of a Fermi characteristic radius of the interparticle interaction. An (cid:48) c gas from the normal to superfluid regime is set by the increase of the critical temperature for the s-wave su- Kosterlitz-Thouless transition. However, for a weak at- perfluidity by the lattice potential has been indicated in tractive interaction the order parameter and the super- Refs. [57, 58]. fluidtransitiontemperaturecanbefoundintheBCSap- The situation changes for the p-wave pairing of iden- proach [55]. For both spinless and spin-1/2 fermions the tical fermions attractively interacting via a short-range renormalizedgapequationfortheorderparameterinthe potential. This pairing in an optical lattice at very low (cid:82) momentum space, ∆k= d2r∆(r−r(cid:48))exp[ik(r−r(cid:48))], temperatures has been considered in Ref. [59] (more so- reads (see [16, 17] and references therein): phisticated lattice models, where p-wave pairing is con- structed with the use of s-wave pairing at intermediate (cid:90) d2k (cid:26) 1 (cid:27) ∆k=− (2π)(cid:48)2f(k(cid:48),k)∆k(cid:48) K(k(cid:48))−2(Ek(cid:48)−Ek−i0) stitgahgtesb,inwdeirnegremcoendtellytwsuoggsuescthedferinmiRoenfss.ca[6n0,no6t1]b).eIinn tthhee (1) (cid:90) d2k same lattice site unless one of them occupies a higher − (2π)(cid:48)2δV(k(cid:48),k)∆k(cid:48)K(k(cid:48)), Bloch band. Therefore, the main contribution to the scatteringamplitudecomesfromtheinteractionbetween where f(k(cid:48),k) is the off-shell scattering amplitude two fermions sitting in neighboring sites [59]. In partic- and Ek = (cid:126)2k2/2m with m being the particle ular, the fermions undergo quantum tunneling from the mass. The single particle excitation energy is given by centers of their sites and experience the short-range in- (cid:112) (cid:15)k = (Ek−µ)2+|∆(k)|2 where µ is the chemical po- teraction in the spatial region where their wavefunctions tential, and K(k) = tanh((cid:15)k/2T)/2(cid:15)k. For weak inter- are attenuated. This strongly suppresses the interaction actions chemical potential coincides with the Fermi en- amplitude and leads to a very low critical temperature. ergy EF =(cid:126)2kF2/2m (kF is the Fermi momentum). The We however show below that the picture is drastically quantity δV(k(cid:48),k) is a correction to the bare interpar- differentforanattractivelong-rangeinteractionbetween ticle interaction due to polarization of the medium by the fermions. colliding particles. The leading terms of this quantity introduced by Gor’kov and Melik-Barkhudarov [56], are second order in the bare interaction (see Methods). P-wave pairing of microwave-dressed In order to gain insight in what is happening, we first polar molecules in a 2D lattice omit the correction δV(k,k) in Eq. (1). We then put (cid:48) k = k , and notice that the main contribution to the F Wewillconsideridentical fermionicpolarmoleculesin integral over k in Eq. (1) comes from k close to k . (cid:48) (cid:48) F a 2D lattice of period b. Being dressed with a microwave At temperatures T tending to the critical temperature field, they acquire an attractive dipole-dipole tail in the Tc from below, we put (cid:15)k(cid:48)=|Ek(cid:48) − EF| in K(k(cid:48)). For interaction potential [16, 17, 62, 63]: the pairing channel related to the interaction with orbital angular momentum l, this immediately leads to an V(r)=−d2/r3. (3) estimate: (cid:20) (cid:21) 1 Heredisaneffectivedipolemoment,andweassumethat T ∼E exp − . (2) c F ρ(k )|f (k )| Eq. (3) is valid at intermolecular distances r (cid:38) b. This F l F leads to superfluid p-wave pairing of the molecules. In Thequantityρ(k )=m/2π(cid:126)2intheexponentofEq. (2) F free space the emerging ground state is the topological is the density of states on the Fermi surface, and f (k ) l F p +ip superfluid,andtheleadingpartofthescattering x y is the on-shell scattering amplitude. amplitude can be obtained in the first Born approxima- In the lattice with a period b satisfying the condition tion [16, 17]. We assume the weakly interacting regime kFb (cid:28) 1, the superfluid paring of fermions can be con- at a small filling factor in the lattice, k b(cid:28)1. F sideredasthatofparticleswitheffectivemassm∗ >min TheHamiltonianofthesystemisHˆ =Hˆ +Hˆ ,with 0 int free space. The density of states ρ(k ) is then given by F the same expression, with m replaced by m . Thus, the (cid:88) ∗ Hˆ = ε aˆ aˆ , (4) BCS exponent [ρ(kF)|fl(kF)|]−1 in the lattice is smaller 0 q q †q q than in free space at the same k (density) if there is no F significant reduction in the scattering amplitude. Hence, where aˆq, aˆ†q are the annihilation and creation operators although the Fermi energy decreases by the same factor ofamoleculewithquasimomentumq,andεqisthesingle m/m , the critical temperature T in the lattice can be particle energy. In the low momentum limit we have ∗ c muchlargerthaninfreespace. Thisisthecaseforthes- εq = (cid:126)2q2/2m∗, where m∗ > m is the effective mass in wavepairingofshort-rangeinteractingspin-1/2fermions the lowest Bloch band. The quantity Hˆ describes the int 3 interaction between the molecules and is given by here B follows from the fact that the relative wavefunction is zero for r ≤b (perfectly reflecting wall). 1 (cid:88) d2 Hînt =−2rj(cid:54)=r(cid:48)jψˆ†(rj)ψˆ†(r(cid:48)j)|rj −r(cid:48)j|3ψˆ(r(cid:48)j)ψˆ(rj), (5) cspriaIttciceiaslbcetleceaamrupsteehraatthtufeorreBtCihnSetsehaxemploeant2teDinctedieinsnslEiatqryg.en(r8()tahniasdnskmiFna)lfltreheree. where ψˆ(r ) is the field operator of a particle in the lat- However, in ordinary optical lattices one has the lattice tice site jjlocated at r in the coordinate space. At a constant b (cid:38) 200 nm. In this case, for m∗/m ≈ 2 (still j the tight binding case with b/ξ ≈ 3, where ξ is the small filling factor in the low momentum limit, the main 0 0 contribution to the matrix elements of Hˆ comes from extension of the particle wavefunction in the lattice site) int intermolecular distances |r − r | (cid:29) b (see Methods). and at a fairly small filling factor (let say, kFb = 0.35) j (cid:48)j Therefore, wemayreplacethesummationoverr andr theFermienergyforthelightestalkalinepolarmolecules j (cid:48)j NaLi is about 10 nK (n ≈ 2×107 cm 2). Then, even by the integration over d2r and d2r . As a result the − j (cid:48)j for k r approaching unity the critical temperature is Hamiltonian of the system reduces to F e∗ff only of the order of a nanokelvin (for k b = 0.35 and F (cid:90) (cid:126)2 r /b≈3 Fig. 2 in Methods gives κ∼1). Hˆ =− ψˆ (r)∇2ψˆ(r)d2r e∗ff 2m † The picture is quite different in recently introduced ∗ (6) d2 (cid:90) 1 subwavelengthlattices[44–52],wherethelatticeconstant − ψˆ†(r)ψˆ†(r(cid:48))ψˆ(r(cid:48))ψˆ(r)d2r, can be as small as b (cid:39) 50 nm. This strongly increases 2 |r−r|3 (cid:48) all energy scales, and even for a small filling factor the where the first term in the right hand side is Hˆ (4) Fermi energy may become of the order of hundreds of 0 nanokelvins. Subwavelength lattices can be designed us- rewritten in the coordinate space. We thus see that ing adiabatic dressing of state-dependent lattices [44], the problem becomes equivalent to that of particles with multi-photon optical transitions [45, 46], spin-dependent mass m in free space. ∗ opticallatticeswithtime-dependentmodulations[47],as The scattering amplitude at k =k , which enters the F well as nanoplasmonic systems [48], vortex arrays in su- exponential factor in Eq. (2), is obtained from the solu- perconductingfilms[49],periodicallypatternedgraphene tionofthescatteringprobleminthelattice. Forparticles monolayers [50], magnetic-film atom chips [51], and pho- thathavemassm (seeMethods), theamplitudeiswrit- ∗ tonic crystals [52–54]. These interesting proposals al- ten as follows readystimulatedstudiesrelatedtomany-bodyphysicsin 8 (cid:126)2 π (cid:126)2 such lattices, in particular the analysis of the Hubbard f(k )=− k r + (k r )2ln(Bk r ), (7) F 3m F e∗ff 2m F e∗ff F e∗ff model and engineering of spin-spin Hamiltonians [52]. ∗ ∗ In the considered case of p +ip pairing in the 2D wherek r (cid:28)1,andBisanumericalcoefficientcoming x y F e∗ff lattice, putting b = 50 nm, for the same k b as above fromshort-rangephysics. Sinceforweakinteractionstwo F the Fermi energy for NaLi molecules exceeds 200 nK fermions practically do not get to the same lattice site, (n ≈ 4×108 cm 2). Then, for the same κ ∼ 1 and for calculating B we may introduce a perfectly reflect- − k r approaching unity we have T ∼ 20 nK, which is ingwallatintermoleculardistancesr ∼b(seeMethods). F e∗ff c twice as high as in free space. An additional advantage For the superfluid pairing the most important are parti- ofthelatticesystemistheforeseenquantuminformation cle momenta ∼k . Therefore, the low-momentum limit F processing, since addressing qubits in the lattice is much requires the inequality k b(cid:28)1. F easier than in free space. The solution of the gap equation (1) then leads to Note that there is a (second-order) process, in which thep +ip superfluidwiththecriticaltemperature(see x y the interaction between two identical fermions belonging Methods): to the lowest Bloch band provides a virtual transfer of κ (cid:20) 3π (cid:21) one of them to a higher band. Then, the two fermions T =E exp − , (8) c F(k r )9π2/64 4k r may get to the same lattice site and undergo the inelas- F e∗ff F e∗ff tic process of collisional relaxation. The rate constant where the coefficient κ is related to B and depends on of this second-order process is roughly equal to the rate the ratio r /b (see Methods). There are two important constantinfreespace,multipliedbytheratioofthescat- e∗ff differencesofequation(8)fromasimilarequationinfree teringamplitude(dividedbytheelementarycellarea)to space obtained in Ref. [16]. First, the Fermi energy E the frequency of the potential well in a given lattice site F is smaller by a factor of m/m , and the effective dipole- (the difference in the energies of the Bloch bands). This ∗ dipole distance r is larger than the dipole-dipole dis- ratio originates from the virtual transfer of one of the e∗ff tance in free space by m /m. Second, the coefficient B fermions to a higher band and does not exceed (ξ/b)2. ∗ and, hence, κ in free space is obtained from the solution Even in not a deep lattice, where m /m is 2 or 3, we ∗ oftheSchrödingerequationinthefullmicrowave-induced have (ξ/b)2 <0.1. Typical values of the rate constant of potential of interaction between two molecules, whereas inelasticrelaxationinfreespaceare∼10 8−10 9 cm2/s − − 4 [16], and hence in the lattice it will be lower than 10 9 Thepotentialofinteractionbetweentwomoleculesbe- − or even 10 10 cm2/s. Thus, the rate of this process is longing to different layers has the form: − rather low and for densities approaching 109 cm 2 the − decay time will be on the level of seconds or even tens of VL(r)=−d2(r2−2L2)/(r2+L2)5/2, (9) seconds. whereListheinterlayerspacing,r isthein-layersepara- tionbetweenthemolecules,and−d2isthescalarproduct of the average dipole moments of these molecules. The Interlayer p-wave superfluid √ of fermionic polar molecules in a bilayer system potential VL(r) is repulsive for r < 2L and attractive atlargerr. Thepotentialwellismuchmoreshallowthan in the case of all dipoles oriented in the same direction, Anotherinterestingnovelsuperfluidoffermionicpolar which was considered in Refs. [40–42]. We have checked molecules is expected in a bilayer system, where dipoles that s-wave interlayer dimers, which exist at any r /L, areorientedperpendicularlytothelayersandinopposite are weakly bound even for r /L ≈ 3. Their binding∗en- directions in different layers. ergy at r /L(cid:46)3 is much sm∗aller than the Fermi energy Such a bilayer configuration, but with all dipoles ori- ∗ at least for k L > 0.1. For such r /L, interlayer dimers entedinthesamedirection, hasbeenconsideredinRefs. F with orbital angular momenta |l| ≥∗ 1 do not exist. We [40–43]. As found, it should form an interlayer s-wave thus are dealing with purely fermionic physics. superfluid, where Cooper pairs are formed by dipoles of For the analysis of the superfluid pairing we are inter- different layers due to the s-wave dipolar interaction be- estedinparticlemomenta k ∼k . Aswellasinthecase tween them. F ofalldipolesorientedinthesamedirection[40–43],under For the dipoles of one layer that are opposite to the theconditionk r (cid:28)1(wherer =md2/(cid:126)2)theamplitude dipoles of the other one, the picture of interlayer pairing F ∗ ∗ of interlayer interaction is obtained in the Born approxi- is different. The s-wave pairing is practically impossi- mation. The Fourier transform of the potential (9) is ble, and the system may form p-wave and higher partial wave superfluids. This type of bilayer systems can be V (k,k)=(2π(cid:126)2/m)r |k−k|exp[−|k−k|L], (10) L (cid:48) ∗ (cid:48) (cid:48) created by putting polar molecules with rotational moment J = 0 in one layer, and molecules with J = 1 in and in the first Born approximation the on-shell ampli- the other. Then, applying an electric field (perpendicu- tude of the l-wave scattering at k =k reads (see Meth- F lar to the layers) one gets a field-induced average dipole ods): moment of J = 0 molecules parallel to the field, and the dipole moment of J = 1 molecules oriented in the 2(cid:126)2k r (cid:90) 2π F ∗ f (k )= dφcos(lφ) opposite direction. One should also prevent a flip-flop l F m (11) 0 process in which the dipole-dipole interaction between ×|sin(φ/2)|exp[−2k L|sin(φ/2)|]. given J = 1 and J = 0 molecules reverses their dipoles, F thus inducing a rapid three-body decay in collisions of a The s-wave amplitude is positive, i.e. the s-wave chan- dipole-reversedmoleculewithtwooriginalones. Thiscan nel corresponds to repulsion. Note that for extremely be done by making the electric field inhomogeneous, so lowcollisionenergiescomparablewiththedimerbinding thatitislargerinthelayerwithJ =0moleculesandthe energy, where the Born approximation is not accurate, flip-flop process requires an increase in the Stark energy. the s-wave scattering amplitude can be negative. This, This process will be suppressed if the difference in the however, does not lead to superfluid s-wave pairing. Stark energies of molecules in the layers significantly ex- The channels with |l|≥1 correspond to attraction. A ceeds the Fermi energy, which is a typical kinetic energy straightforwardcalculationshowsthatfork L(cid:46)0.7the F of the molecules (∼ 100 nK for the example considered largestisthep-waveamplitudeand,hence,atsufficiently below). This is realistic for present facilities. lowtemperaturesthesystemwillbeaninterlayerp-wave For the dipole moment close to 1 Debye and the in- superfluid. As for d-wave and higher partial wave super- terlayer spacing of 50 nm, one thus should have the field fluids, they are possible only at extremely low tempera- gradient (perpendicularly to the layers) significantly ex- tures. Thus, we confine ourselves to the p-wave pairing ceeding 0.5 kV/cm2. This could be done by using elec- and employ the BCS approach. trodesconsistingoffourrods,andevenahighergradient A detailed analysis of the gap equation (1), which in- ∼30kV/cm2 shouldbeachievable[64,65]. Bychanging cludesfirstandsecondordercontributionstothescatter- thepositionsoftherodsonecanobtainthefieldgradient ing amplitude and Gor’kov-Melik-Barkhudarov correc- exceeding 0.5 kV/cm2 in the direction perpendicular to tions, is given in Methods. The critical temperature for the layers of the bilayer system. The field itself will not the p-wave superfluidity proves to be (see Methods): beexactlyperpendiculartothelayersandtherewillalso (cid:20) (cid:21) bethefieldgradientparalleltothelayers. This,however, F(k L) T =E β(k L)exp − F , (12) does not essentially influence the physics. c F F k r F ∗ 5 and for not very small k r the validity of the pertur- Superfluidityitselfcanbedetectedinthesamewayas F ∗ bative treatment of the Gor’kov-Melik-Barkhudarov cor- in the case of s-wave superfluids [69, 70]. Rotating the rections requires k L (cid:38) 0.15 (see Methods). The func- p +ip superfluid and inducing the appearance of vor- F x y tionsF(k L)andβ(k L)aregiveninMethods. Fork L tices one can find signatures of Majorana modes on the F F F ranging from 0.15 to 0.3 the function F increases from vortexcoresintheRFabsorptionspectrum[71]. Eventu- 3.4 to 5, and the coefficient β is fairly large, being about ally, onecanthinkofrevealingthestructureoftheorder 80 at k L=0.15 (see Fig. 3 and Methods). parameter by visualizing vortex-related dips in the den- F Creating the bilayer system by using a 1D subwave- sity profile on the approach to the strongly interacting length lattice we may have L ≈ 50 nm. In this case, regime, where these dips should be pronounced at least for k L = 0.15 the Fermi energy of NaLi molecules is in time-of-flight experiments. F close to 100 nK, and the critical temperature for k r F ∗ approaching 0.5 is about 10 nK. Forcompleteness,wealsoconsidertheregimeofstrong Methods interactionswithinasinglelayer. Assumingthatthecou- pling between the layers is still fairly weak, we have su- Scattering problem and superfluid pairing of perfluid (interlayer) pairing between quasiparticles. Re- microwave-dressed polar molecules in a 2D lattice latedproblemshavebeendiscussedforcoupled2DFermi liquids as models for layered superconductors [8]. In this As we concluded in the main text, in the low momen- case,wereplacethebaremassmbytheeffectivemassm tum limit at a small filling factor the system of lattice ∗ and account for renormalization of the fermionic Green polar molecules is equivalent to that of molecules with functions by a factor Z < 1 [66]. Then, the expression effective mass m in free space. We now demonstrate ∗ for the transition temperature takes the form: this explicitly by the calculation of the off-shell scattering amplitude f(k,k). For our problem the main part (cid:20) (cid:21) (cid:48) T ∼E exp −F(kFL) m 1 , (13) of the scattering amplitude can be obtained in the Born c F kFr∗ m∗Z2 approximation [16]. In the lattice the scattering amplitude is, strictly where we can not determine the pre-exponential coeffi- speaking, the function of both incoming quasimomenta cient. Therefore, Eq. (13) only gives an order of mag- q ,q and outgoing quasimomenta q ,q . However, in 1 2 (cid:48)1 (cid:48)2 nitude of Tc. For kFL = 0.3 and L ≈ 50 nm the Fermi the low-momentum limit where qb (cid:28) 1, taking into ac- energy of NaLi molecules is about 400 nK, and for, let count the momentum conservation law the amplitude say, kFr ≈ 2 the dimer physics is still not important. becomes the function of only relative momenta k = ∗ Then, using the effective mass and factor Z from the (q −q )/2 and k =(q −q )/2. For the off-shell scat- 1 2 (cid:48) (cid:48)1 (cid:48)2 Monte Carlo calculations [67] one may think of super- tering amplitude the first Born approximation gives: fluid transition temperatures of the order of several tens (cid:90) of nanokelvins. f(k,k)=S χ (r )χ (r )V(r −r ) (cid:48) ∗q(cid:48) 1 ∗q(cid:48) 2 1 2 1 2 ×χ (r )χ (r )d2r d2r Conclusions q1 1 q2 2 1 2 (14) =−d2b4 (cid:88)exp[i(q1−q(cid:48)1)rj+i(q2−q(cid:48)2)r(cid:48)j], S |r −r |3 We have demonstrated the emergence of the topolog- rj,r(cid:48)j j (cid:48)j ical p +ip superfluid for identical microwave-dressed x y fermionic polar molecules in a 2D lattice. Another where V(r −r ) is given by Eq. (3) of the main text, 1 2 novel p-wave superfluid is found to emerge for fermionic and S is the surface area. The last line of Eq. (14) is molecules in a bilayer system, with dipoles of one layer obtained assuming the tight-binding regime, where the opposite to the dipoles of the other one. In both cases single particle wavefunction is the use of subwavelength lattices with a period b (cid:39) 50 1 (cid:88) nm (creation of the bilayer system with the interlayer χ (r)= √ Φ (r−r )exp[iqr ]. (15) q 0 j j spacingL(cid:39)50nm)allowsonetoobtainsuperfluidtran- N j sition temperature of the order of tens of nanokelvins. This opens interesting prospects for topologically pro- Here, the index j labels the lattice sites located at the tected quantum information processing with p +ip su- points r , and N =S/b2 is the total number of the sites. x y j perfluids in 2D lattices. The interlayer p-wave superfluid The particle wavefunction in a given site j has extension √ in bilayer systems, together with the earlier proposed s- ξ and is expressed as Φ (r−r )=(1/ πξ )exp[−(r− 0 0 j 0 waveinterlayersuperfluid[40–43]andsuperfluidsinmul- r )2/2ξ2]. In the low-momentum limit we may replace j 0 tilayerfermionicsystems[68], canbeastartingpointfor thesummationoverj andj bytheintegrationoverd2r (cid:48) j the creation of more sophisticated layered structures. and d2r taking into account that b2(cid:80) transforms into j(cid:48) j 6 (a)  (b) B 10 100 1 10 1 0.10 0.10 0.01 2 3 4 5 2.0 2.5 3.0 3.5 4.0 4.5 re⇤↵/b re⇤↵/b Figure 2. Coefficients B and κ as functions of r∗ /b. eff (a) (b) F(k L) �(k L) (cid:82) d2r F. Th5.0is immediately yields F 110fermionic polar molecules in a bilayer system j 100 4f.5(k,k)=−d2(cid:90) exp[i(k−k)r]d2r, (16) For t90he interlayer interaction potential VL(r) given by (cid:48) (cid:48) r3 equatio80n (9) in the main text, the scattering amplitude 4.0 for kFr70∗(cid:28)1 can be calculated in the Born approxima- and the p-wave part of the scattering amplitude is ob- tion [40]. The p-wave part of the first order contribution 60 tainedmultiplyingEq. (16)byexp(−iφ)andintegrating 3.5 to the off-shell amplitude is over dφ/2π, where φ is the angle between the vectors k 50 (cid:90) fa[1n(6kd]))k.=(cid:48)T.−hTe(3h8.00oi.(cid:126)1sn52i-/ss3htmehle∗l)askamrme∗pff0el.2,i0trwuesdhueeklrtF(ekaLrs=e∗ffink0=.(cid:48)f2)r5emcea∗sndp2ba/ce(cid:126)e2w(irss0ie.3tet0th,eene.egaf.s-, f1(4k0(cid:48),k)= 2dπ(ϕ(cid:126)02.2(cid:48)k20πd)ϕ2kkeFi(Lϕ(cid:48)k−ϕ0.2k5)VL(r)ei(k(cid:48)−0.30k)r (19) fectivedipole-dipoledistanceinthelattice. Theapplica- =− (kr∗)F1(k(cid:48),k,L), m bilityoftheBornapproximationassumesthatkr (cid:28)1, e∗ff whichisclearlyseenbycalculatingthesecondordercor- where rection to the scattering amplitude. (cid:90) 1 ∞ Uptotheterms∼(kre∗ff)2, theon-shellscatteringam- F1(k(cid:48)L,kL)= kL xdxJ1(k(cid:48)Lx)J1(kLx) plitudefollowingformthesolutionofthescatteringprob- 0 (20) x2−2 lem for particles with mass m∗, is given by [16]: × , (x2+1)5/2 8(cid:126)2 π(cid:126)2 f(k)=− kr∗+ (kr∗)2ln(Bkr∗), (17) andJ1istheBesselfunction. Regardingthesecondorder 3m 2 m contribution,forthesolutionofthegapequationweonly need the on-shell p-wave part, which is given by wherethenumericalcoefficientBcomesformshort-range physics. For calculating B we introduce a perfectly re- 2π(cid:126)2 flectingwallatintermoleculardistancesr∼b,whichtakes f2(k)= (kr∗)2F2(kL), (21) m intoaccountthattwofermionspracticallycannotgetto one and the same lattice site. The coefficient B depends where o2na.the ratio re∗ff/b, and we show this dependence in Fig. F (kL)= π (cid:90) ∞xdxJ2(kLx) x2−2 The treatment of the superfluid pairing is the same 2 (kL)2 0 1 (x2+1)5/2 (22) asinRef. [16], includingtheGorkov-Melik-Barkhudarov (cid:90) ∞ y2−2 × ydyJ (kLy)N (kLy) . correction. We should only replace the mass m with m . 1 1 (y2+1)5/2 ∗ x Theexpressionforthecriticaltemperaturethenbecomes: Infact,thetruep-wavescatteringamplitudefollowsfrom κ (cid:20) 3π (cid:21) the exact relation T =E exp − , (18) c F(k r )9π2/64 4k r (cid:90) F e∗ff F e∗ff f(k ,k)= ∞J (k r)V (r)ψ(k,r)2πrdr, (23) (cid:48) 1 (cid:48) L where κ (cid:39) 0.19B 9π2/64, and it is displayed in Fig. 2b 0 − as a function of r /b. whereψ(k,r)isthetruewavefunctionofthep-waverela- e∗ff tive motion with momentum k, normalized such that for r →∞wehaveψ(k,r)=J (kr)−(imf(k)/4(cid:126)2)H(1)(kr) 1 1 Superfluid pairing of with H(1) being the Hankel function. This amplitude is 1 7 complexanditisrelatedtotherealamplitudef˜=f +f ∆(k ) from Eq. (26) and find 1 2 (cid:48) given by equations (19)–(22) as ∆(2)(k )=−∆(k )ρ(k )f (k ) F F F 1 F f(k(cid:48),k)= 1+if˜m(kf˜(cid:48)(,kk))/4(cid:126)2. (24) ×(cid:20)ln(cid:18)EωF(cid:19)−η(kFL)(cid:21), (28) Inordertocalculatethesuperfluidtransitiontempera- where tureweusetheBCSapproachalongthelinesofRef.[16]. We consider temperatures T tending to Tc from below η(k L)=−2(cid:90) EF−ω k(cid:48)dk(cid:48) (cid:20)f12(kF,k) −1(cid:21) and rely on the renormalized gap equation (1). For the F 0 (kF2 −k(cid:48)2) f12(kF) p-wave pairing the order parameter is ∆k = ∆(k)eiϕk, (cid:90) 1 ydy (cid:40)(cid:20)F (k L,k Ly)(cid:21)2 (cid:41) (29) andwethenmultiplyEq.(1)bye−iϕk andintegrateover =−2 1−y2 1FF(k LF) −1 , dϕk(cid:48) and dϕk. As a result, we obtain the same equation 0 1 F (1) in which ∆k and ∆k(cid:48) are replaced with ∆(k) and and F (k L)≡F (k L,k L) ∆(k ), the off-shell scattering amplitude f(k,k) with 1 F 1 F F (cid:48) (cid:48) Then, we consider the Gor’kov-Melik-Barkhudarov its p-wave part, and δV(k,k) with its p-wave part (cid:82) (cid:48) corrections to the bare interaction of the molecules in δV(k(cid:48),k)= δV(k(cid:48),k)exp[i(ϕk(cid:48) −ϕk)]dϕk(cid:48)dϕk/4π2. the bilayer. These many-body corrections are second or- Calculating the contribution of the pole in the second der in (k r ) and are described by four diagrams (for term in square brackets and using Eq. (24) we obtain F ∗ details, see [16, 41, 56]). For the case of p-wave super- ∆(k)=−P(cid:82) mdEk(cid:48)f˜(k ,k)∆(k )(cid:104)K(k )− 1 (cid:105) fluidity of identical fermionic polar molecules they have 2π(cid:126)2 (cid:48) (cid:48) (cid:48) 2(Ek(cid:48)−Ek) beenconsideredinRef.[16]. Theyhavebeenalsostudied −(cid:82) mdEk(cid:48)δV(k ,k)∆(k )K(k ), (25) for the interlayer s-wave superfluidity of dipoles oriented 2π(cid:126)2 (cid:48) (cid:48) (cid:48) in the same direction in Ref. [41]. where the symbol P stands for the principal value of Weareinterestedinthecaseofsufficientlysmallk L. F the integral. In the first term in the right hand side Following the same treatment as in Refs. [16, 41], in the of Eq. (25) we divide the region of integration into limit of k L→0 we obtain: F two parts: |Ek(cid:48) − EF| < ω and |Ek(cid:48) − EF|>ω, where (cid:126)2 ∆(kF),Tc (cid:28) ω (cid:28) EF. The contribution to the p- δV(kF,kF)=−α (kFr∗)2, (30) wave order parameter from the first region we denote m as ∆(1)(k), and the contribution from the second region whereα(cid:39)10.57. Thedominantcontributiontothisresult as ∆(2)(k). The contribution of the second term in right comesfromthediagramcontainingabubbleintheinter- hand side of equation (25) is denoted as ∆(3)(k). action line (diagram a) in Refs. [16, 41]). This contribu- We first notice that the main contribution to ∆(k) tion strongly decreases with increasing k L. In particu- F comes from k close to k . Retaining only f , which is (cid:48) F 1 lar, for k L (cid:39) 0.15 we have α (cid:39) 2.8, and α (cid:39) 2.2 when F proportional to kr , in the off-shell scattering amplitude ∗ increasing k L to 0.2. Comparing δV with the scat- F and omitting the second term in the right hand side of tering amplitude f (k ) we see that for not very small Eq. (25) (which is proportional to (kr )2) we obtain 1 F ∗ k r the perturbative treatment of the Gor’kov-Melik- F ∗ Barkhudarovcorrectionsisadequatefork L(cid:38)0.15. We ∆(k)=∆(k )f (k ,k)/f (k ). (26) F F 1 F 1 F therefore confine ourselves to these values of k L. F Performing the integration in the second term of Putting k = k and performing the integration in the F Eq.(25)weobtainthecontributionoftheGor’kov-Melik- first region in the first term of Eq (25), where E −ω < Ek(cid:48) <EF +ω, we may put ∆(k(cid:48))=∆(kF) and f˜(Fk(cid:48),k)= Barkhudarov corrections to the order parameter: f˜(kF) = f1(kF)+f2(kF). Then, putting (cid:15)k(cid:48) = |ξk(cid:48)| in α(k L) (cid:18)E (cid:19) K(k(cid:48)) and taking into account that the contribution of ∆(3)(kF)=∆(kF) 2Fπ (kFr∗)2ln TF , (31) c the second term in square brackets is zero, we obtain: The sum of Eqs. (27), (28) and (31) yields (cid:18)2eCω(cid:19) ∆(1)(kF)=−∆(kF)ρ(kF)f˜(kF)ln πTc , (27) ∆(k )=∆(k )(cid:20)(cid:26)(k r )F (k L)ln(cid:16)2eC−η(kFL)EF(cid:17)(cid:27) F F F ∗ 1 F π Tc with C=0.577 being the Euler constant, and ρ(kF) = (cid:110)(cid:16) (cid:17) (cid:16) (cid:17)(cid:111)(cid:21) m/2π(cid:126)2 the density of states. − F2(kFL)−α(k2FπL) (kFr∗)2ln ETFc , (32) Inthesecondregion,whereEk(cid:48)>EF+ωorEk(cid:48) <EF− ω,weputK=1/2|ξk(cid:48)|andkeeponlyf1 inthescattering where we put ω ∼ EF in the terms proportional to amplitude. Fork =kF theintegraloverEk(cid:48) fromEF+ω (kFr∗)2. We should also recall that the bare mass m to ∞ vanishes. In the integral from 0 to E −ω we use should be replaced with the effective mass m = m[1− F ∗ (a)  (b) B 10 100 1 10 1 0.10 0.10 0.01 2 3 4 5 2.0 2.5 3.0 3.5 4.0 4.5 re⇤↵/b re⇤↵/b 8 (a) (b) F(kFL) �(kFL)80 5.0 70 4.5 60 4.0 50 3.5 40 3.0 0.15 0.20 k L 0.25 0.30 0.15 0.20 k L 0.25 0.30 F F Figure 3. The dependence of F and β on k L. F (4/3π)k r ]whichhasbeenfoundinRefs.[41,72]. Since Competing financial interests: The authors declare F ∗ therelativedifferencebetweenm andmissmallask r , no competing financial interests. ∗ F ∗ it is sufficient to replace m with m∗ only in the multiple Corresponding author:Correspondenceshouldbead- r∗∼m in the first term of Eq. (32). This leads to the dressed to A.K. Fedorov ([email protected]). appearance of a new term (cid:18) (cid:19) 4 E −∆(k ) F (k L)(k r )2ln F (33) F 1 F F ∗ 3π T c [1] V.P.MineevandK.V.Samokhin,IntroductiontoUncon- in the right hand side of equation (32). Then, dividing ventional Superconductivity (CRC Press, 1999). bothsidesofEq.(32)by∆(kF)weobtainforthecritical [2] G.E.Volovik,ExoticPropertiesofSuperfluid3He(World temperature: Scientific, Singapore, 1992). (cid:20) (cid:21) [3] N. Read and D. Green, Phys. Rev. B 61, 10267 (2000). F(k L) T =E β(k L)exp − F , (34) [4] K.D. Nelson, Z.Q. Mao, Y. Maeno, and Y. Liu, Science c F F kFr∗ 306, 1151 (2004). [5] T.M. Rice, Science 306, 1142 (2004). where [6] S. Das Sarma, C. Nayak, and S. Tewari, Phys. Rev. B F(kFL)=[F1(kFL)]−1 (35) 73, 220502 (2006). [7] L.N. Bulaevskii, Sov. Phys. JETP 37, 1133 (1973). and [8] M.J. Graf, D. Rainer, and J.A. Sauls, Phys. Rev. B 47, (cid:20) (cid:18) (cid:19) 12089 (1993). 2 β(k L)=exp C+ln −η(k L) [9] P. Fulde and R.A. Ferrell, Phys. Rev 135, A550 (1964). F F π [10] A.I.LarkinandY.N.Ovchinnikov,Sov.Phys.JETP20, (36) F (k L) 4 1 α(k L) (cid:21) 762 (1965). − F22(kFL) − 3πF (k L) + 2Fπ F2(kFL) . [11] For a review see R. Casalbuoni and G. Nardulli, Rev. 1 F 1 F Mod. Phys. 76, 263 (2004). The dependence of F and β on k L is shown in Fig. 3. [12] For a review, see V. Gurarie and L. Radzihovsky, Ann. F Phys. (Amsterdam) 322, 2 (2007). We stop at k L = 0.3 because for larger values of this F [13] V. Gurarie, L. Radzihovsky, and A.V. Andreev, Phys. parameter the function F is so large that the critical Rev. Lett. 94, 230403 (2005). temperature will be negligible. [14] Y. Nishida, Ann. Phys. (N.Y.) 324, 897 (2009). [15] M.S.Foster,V.Gurarie,M.Dzero,andE.A.Yuzbashyan, Phys. Rev. Lett. 113, 076403 (2014). Acknowledgments [16] N.R. Cooper and G.V. Shlyapnikov, Phys. Rev. Lett. 103, 155302 (2009). We acknowledge fruitful discussions with Jun Ye, Ig- [17] J. Levinsen, N.R. Cooper, and G.V. Shlyapnikov, Phys. Rev. A 84, 013603 (2011). nacio Cirac, and Martin Zwierlein. The research lead- [18] A. Stern, Ann. Phys. 323, 204 (2008). ing to these results has received funding from the Eu- [19] G. Mo¨ller, N.R. Cooper, and V. Gurarie, Phys. Rev. B ropean Research Council under European Community’s 83, 014513 (2011). Seventh Framework Programme (FR7/2007-2013 Grant [20] For a review, see C. Nayak, S.H. Simon, A. Stern, M. Agreement no. 341197). Freedman,andS.DasSarma,Rev.Mod.Phys.80,1083 (2008). Author contributions: A.K. Fedorov, S.I. Matveenko, [21] S. Das Sarma, M. Freedman, and C. Nayak, npj Quant. and V.I. Yudson equally contributed to the obtained re- Inform. 1, 15001 (2015). sults. G.V.Shlyapnikovsupervisedtheproject,andA.K. [22] J. Levinsen, N.R. Cooper, and V. Gurarie, Phys. Rev. Fedorov and G.V. Shlyapnikov have written the text. Lett. 99, 210402 (2007). 9 [23] J.Levinsen,N.R.Cooper,andV.Gurarie,Phys.Rev.A [48] M. Gullans, T.G. Tiecke, D.E. Chang, J. Feist, J.D. 78, 063616 (2008). Thompson, J.I. Cirac, P. Zoller, and M. D. Lukin, Phys. [24] M. Jona-Lasinio, L. Pricoupenko, and Y. Castin, Phys. Rev. Lett. 109, 235309 (2012). Rev. A 77, 043611 (2008). [49] O. Romero-Isart, C. Navau, A. Sanchez, P. Zoller, and [25] For a review, see L.D. Carr, D. DeMille, R.V. Krems, J.I. Cirac, Phys. Rev. Lett. 111, 145304 (2013). and J. Ye, New J. Phys. 11, 055049 (2009). [50] C.Huang,F.Ye,Z.Sun,andX.Chen,Opt.Express22, [26] Forareview,seeO.DulieuandC.Gabbanini,Rep.Prog. 30108 (2014). Phys. 72, 086401 (2009). [51] V.Y.F. Leung, A. Tauschinsky, N.J. Druten, and R.J.C. [27] K.-K. Ni, S. Ospelkaus, M.H.G. de Miranda, A. Pe’er, Spreeuw, Quantum Inf. Process. 10, 955 (2011). B. Neyenhuis, J.J. Zirbel, S. Kotochigova, P.S. Julienne, [52] A.Gonza´lez-Tudela,C.-L.Hung,D.E.Chang,J.I.Cirac, D.S. Jin, and J. Ye, Science 322, 231 (2008). and H.J. Kimble, Nat. Photonics 9, 320 (2015). [28] S. Ospelkaus, K.-K. Ni, D. Wang, M.H.G. de Miranda, [53] J.D.Thompson,T.G.Tiecke,N.P.deLeon,J.Feist,A.V. B. Neyenhuis, G. Quéméner, P.S. Julienne, J.L. Bohn, Akimov, M. Gullans, A.S. Zibrov, V. Vuletic, and M.D. D.S. Jin, and J. Ye, Science 327, 853 (2010). Lukin, Science 340, 1202 (2013). [29] K.-K. Ni, S. Ospelkaus, D. Wang, G. Quéméner, B. [54] J.S.Douglas, H.Habibian, C.-L.Hung, A.V.Gorshkov, Neyenhuis,M.H.G.deMiranda,J.L.Bohn,J.Ye,and H.J. Kimble, and D.E. Chang, Nat. Photonics 9, 326 D.S. Jin, Nature (London) 464, 1324 (2010). (2015). [30] M.H.G.deMiranda,A.Chotia,B.Neyenhuis,D.Wang, [55] K. Miyake, Prog. Theor. Phys. 69, 1794 (1983). G. Quéméner, S. Ospelkaus, J.L. Bohn, J. Ye, and D.S. [56] L.P. Gor’kov and T.K. Melik-Barkhudarov, Sov. Phys. Jin, Nature Phys. 7, 502 (2011). JETP 13, 1018 (1961). [31] A.Chotia,B.Neyenhuis,S.A.Moses,B.Yan,J.P.Covey, [57] W.Hofstetter,J.I.Cirac,P.Zoller,E.Demler,andM.D. M. Foss-Feig, A.M. Rey, D.S. Jin, and J. Ye, Phys. Rev. Lukin, Phys. Rev. Lett. 89, 220407 (2002). Lett. 108, 080405 (2012). [58] J.K.Chin,D.E.Miller,Y.Liu,C.Stan,W.Setiawan,C. [32] J. Deiglmayr, A. Grocola, M. Repp, K. Morlbauer, C. Sanner, K. Xu, and W. Ketterle, Nature (London) 443, Gluck, J. Lange, O. Dilieu, R. Wester, and M. Wei- 961 (2006). demuller, Phys. Rev. Lett. 101, 133004 (2008). [59] M.IskinandC.A.R.Sa´deMelo,Phys.Rev.B72,224513 [33] M.-S.Heo,T.T.Wang,C.A.Christensen,T.M.Rvachov, (2005). D.A.Cotta,J.H.Choi,Y.R.Lee,andW.Ketterle,Phys. [60] A.Bu¨hler, N.Lang, C.V.Kraus, G.Mo¨ller, S.D.Huber, Rev. A 86, 021602 (2012). and H.P. Bu¨chler, Nat. Commun. 5, 4504 (2014). [34] J.W. Park, S.A. Will, and M.W. Zwierlein, Phys. Rev. [61] B. Liu, X. Li, B. Wu, and W.V. Liu, Nat. Commun. 5, Lett. 114, 205302 (2015). 5064 (2014). [35] S.A. Moses, J.P. Covey, M.T. Miecnikowski, B. Yan, B. [62] A.Micheli,G.Pupillo,H.P.Bu¨chler,andP.Zoller,Phys. Gadway, J. Ye, and D.S. Jin, Science 350, 659 (2015). Rev. A 76, 043604 (2007). [36] J.P. Covey, S.A. Moses, M. Garttner, A. Safavi-Naini, [63] A.V.Gorshkov,P.Rabl,G.Pupillo,A.Micheli,P.Zoller, M.T. Miecnikowski, Z. Fu, J. Schachenmayer, P.S. Juli- M.D. Lukin, and H.P. Bu¨chler, Phys. Rev. Lett. 101, enne, A.M. Rey, D.S. Jin, and J. Ye, Nat. Commun. 7, 073201 (2008). 11279 (2016). [64] M. Lara, B.L. Lev, and J.L. Bohn, Phys. Rev. A 78, [37] For a review, see M.A. Baranov, Phys. Rep. 464, 71 033433 (2008). (2008). [65] S. Hoekstra, V. Metsa¨lä, P.C. Zieger, L. Scharfenberg, [38] For a review, see T. Lahaye, C. Menotti, L. Santos, M. J.J. Gilijamse, G. Meijer, and S.Y.T. van de Meerakker, Lewenstein, and T. Pfau, Rep. Prog. Phys. 72, 126401 Phys. Rev. A 76, 063408 (2007). (2009). [66] A.A. Abrikosov, L.P. Gorkov, and I.E. Dzyaloschinskii, [39] For a review, see M.A. Baranov, M. Delmonte, G. Methods of Quantum Field Theory in Statistical Physics Pupillo, and P. Zoller, Chem. Rev. 112, 5012 (2012). (Prentice-Hall, 1965). [40] A.Pikovski,M.Klawunn,G.V.Shlyapnikov,andL.San- [67] N. Matveeva and S. Giorgini, Phys. Rev. Lett. 109, tos, Phys. Rev. Lett. 105, 215302 (2010). 200401 (2012). [41] M.A.Baranov,A.Micheli,S.Ronen,andP.Zoller,Phys. [68] A.C.Potter,E.Berg,D.-W.Wang,B.I.Halpern,andE. Rev. A 83, 043602 (2011). Demler, Phys. Rev. Lett. 105, 220406 (2010). [42] M. Klawunn, A. Pikovski, and L. Santos, Phys. Rev. A [69] For a review, see S. Giorgini, L.P. Pitaevskii, and S. 82, 044701 (2010). Stringari, Rev. Mod. Phys. 80, 1215 (2008). [43] N.T. Zinner, B. Wunsch, D. Pekker, and D.-W. Wang, [70] M.W. Zwierlein, J.R. Abo-Shaeer, A. Schirotzek, C.H. Phys. Rev. A 85, 013603 (2012). Schunck, and W. Ketterle, Nature (London) 435, 1047 [44] W.Yi,A.Daley,G.Pupillo,andP.Zoller,NewJ.Phys. (2005). 10, 073015 (2008). [71] E. Grosfeld, N.R. Cooper, A. Stern, and R. Ilan, Phys. [45] Phys. Rev. A 66, 045402 (2002). Rev. B 76, 104516 (2007). [46] B. Dubetsky and P. Berman, Laser Phys. 12, 1161 [72] Z.-K.LuandG.V.Shlyapnikov,Phys.Rev.A85,023614 (2002). (2012). [47] S. Nascimbene, N. Goldman, N.R. Cooper, and J. Dal- ibard, Phys. Rev. Lett. 115, 140401 (2015).