Novel p-wave superfluids of fermionic polar molecules

Recently suggested subwavelength lattices offer remarkable prospects for the observation of novel superfluids of fermionic polar molecules. It becomes realistic to obtain a topological p-wave superfluid of microwave-dressed polar molecules in 2D lattices at temperatures of the order of tens of nanokelvins, which is promising for topologically protected quantum information processing. Another foreseen novel phase is an interlayer p-wave superfluid of polar molecules in a bilayer geometry.


General Relations and Qualitative Arguments
The superfluid pairing of identical fermions is characterized by the order parameter Δ (r, r′) = V(r − r′) ψ ψ ′ ×ˆr r ( ) ( ) , where V(r − r′) is the interaction potential, the symbol 〈 ...〉 denotes the statistical average, and ψ r ( ) is the field operator of fermions. For spin-1/2 fermions one of the field operators in the expression for Δ (r, r′) is for spin-↑ fermions, and the other one for spin-↓ fermions. In free space the order parameter depends on the coordinates r and r′ only through the difference (r − r′ ). In 2D the transition temperature T c of a Fermi gas from the normal to superfluid regime is set by the Kosterlitz-Thouless transition. However, for a weak attractive interaction the order parameter and the superfluid transition temperature can be found in the BCS approach 55 . For both spinless and spin-1/2 fermions the renormalized gap equation for the order parameter in the momentum space, Scientific RepoRts | 6:27448 | DOI: 10.1038/srep27448 where f(k′ , k) is the off-shell scattering amplitude and E k =  2 k 2 /2m with m being the particle mass. The single particle excitation energy is given by where μ is the chemical potential, and K ε ε = k T ( ) tanh( /2 )/2 k k . For weak interactions chemical potential coincides with the Fermi energy = E k m /2 F F 2 2  (k F is the Fermi momentum). The quantity δV(k′, k) is a correction to the bare interparticle interaction due to polarization of the medium by 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).
In order to gain insight in what is happening, we first omit the correction δV(k′, k) in Eq. (1). We then put k = k F , and notice that the main contribution to the integral over k′ in Eq. (1) comes from k′ close to k F . At temperatures T tending to the critical temperature T c from below, we put  = − For the pairing channel related to the interaction with orbital angular momentum l, this immediately leads to an estimate: The quantity ρ(k F ) = m/2π 2 in the exponent of Eq. (2) is the density of states on the Fermi surface, and f l (k F ) is the on-shell scattering amplitude.
In the lattice with a period b satisfying the condition k F b ≪ 1, the superfluid paring of fermions can be considered as that of particles with effective mass m * > m in free space. The density of states ρ(k F ) is then given by the same expression, with m replaced by m * . Thus, the BCS exponent [ρ(k F )| f l (k F )|] −1 in the lattice is smaller than in free space at the same k F (density) if there is no significant reduction in the scattering amplitude. Hence, although the Fermi energy decreases by the same factor m/m * , the critical temperature T c in the lattice can be much larger than in free space. This is the case for the s-wave pairing of short-range interacting spin-1/2 fermions in the tight binding model, if the extension of the particle wavefunction in the lattice site greatly exceeds the characteristic radius of the interparticle interaction. An increase of the critical temperature for the s-wave superfluidity by the lattice potential has been indicated in refs 57 and 58.
The situation changes for the p-wave pairing of identical fermions attractively interacting via a short-range potential. This pairing in an optical lattice at very low temperatures has been considered in ref. 59 (more sophisticated lattice models, where p-wave pairing is constructed with the use of s-wave pairing at intermediate stages, were recently suggested in refs 60 and 61). In the tight binding model two such fermions can not be in the same lattice site unless one of them occupies a higher Bloch band. Therefore, the main contribution to the scattering amplitude comes from the interaction between two fermions sitting in neighboring sites 59 . In particular, the fermions undergo quantum tunneling from the centers of their sites and experience the short-range interaction in the spatial region where their wavefunctions are attenuated. This strongly suppresses the interaction amplitude and leads to a very low critical temperature. We however show below that the picture is drastically different for an attractive long-range interaction between the fermions.

P-wave Pairing of Microwave-dressed Polar Molecules in a 2D Lattice
We will consider identical fermionic polar molecules in a 2D lattice of period b. Being dressed with a microwave field, they acquire an attractive dipole-dipole tail in the interaction potential 16,17,62,63 Here d is an effective dipole moment, and we assume that Eq. (3) is valid at intermolecular distances  r b. This leads to superfluid p-wave pairing of the molecules. In free space the emerging ground state is the topological p x + ip y superfluid, and the leading part of the scattering amplitude can be obtained in the first Born approximation 16,17 . We assume the weakly interacting regime at a small filling factor in the lattice, k F b ≪ 1.
The Hamiltonian of the system is = +ˆĤ H q q q q 0 where â q , ˆ † a q are the annihilation and creation operators of a molecule with quasimomentum q, and ε q is the single particle energy. In the low momentum limit we have ε q = 2  q 2 /2m * , where m * > m is the effective mass in the lowest Bloch band. The quantity Ĥ int describes the interaction between the molecules and is given by where ψ r ( ) j is the field operator of a particle in the lattice site j located at r j in the coordinate space. At a small filling factor in the low momentum limit, the main contribution to the matrix elements of Ĥ int comes from inter- . Therefore, we may replace the summation over r j and ′ r j by the integration over d 2 r j and ′ d r j 2 . As a result the Hamiltonian of the system reduces to where the first term in the right hand side is Ĥ 0 (4) rewritten in the coordinate space. We thus see that the problem becomes equivalent to that of particles with mass m * in free space.
The scattering amplitude at k = k F , which enters the exponential factor in Eq. (2), is obtained from the solution of the scattering problem in the lattice. For particles that have mass m * (see Methods), the amplitude is written as follows , and B is a numerical coefficient coming from short-range physics. Since for weak interactions two fermions practically do not get to the same lattice site, for calculating B we may introduce a perfectly reflecting wall at intermolecular distances r ~ b (see Methods). For the superfluid pairing the most important are particle momenta ~k F . Therefore, the low-momentum limit requires the inequality k F b ≪ 1.
The solution of the gap equation (1) then leads to the p x + ip y superfluid with the critical temperature (see Methods): where the coefficient κ is related to B and depends on the ratio ⁎ r b / eff (see Methods). There are two important differences of equation (8) from a similar equation in free space obtained in ref. 16. First, the Fermi energy E F is smaller by a factor of m/m * , and the effective dipole-dipole distance ⁎ r eff is larger than the dipole-dipole distance in free space by m * /m. Second, the coefficient B and, hence, κ in free space is obtained from the solution of the Schrödinger equation in the full microwave-induced potential of interaction between two molecules, whereas here B follows from the fact that the relative wavefunction is zero for r ≤ b (perfectly reflecting wall).
It is clear that for the same 2D density n (and k F ) the critical temperature in the lattice is larger than in free space because the BCS exponent in Eq. (8) is smaller. However, in ordinary optical lattices one has the lattice constant  b 200 nm. In this case, for m * /m ≈ 2 (still the tight binding case with b/ξ 0 ≈ 3, where ξ 0 is the extension of the particle wavefunction in the lattice site) and at a fairly small filling factor (let say, k F b = 0.35) the Fermi energy for the lightest alkaline polar molecules NaLi is about 10 nK (n ≈ 2 × 10 7 cm −2 ). Then, even for The picture is quite different in recently introduced subwavelength lattices 44-52 , where the lattice constant can be as small as  b 50 nm. This strongly increases all energy scales, and even for a small filling factor the Fermi energy may become of the order of hundreds of nanokelvins. Subwavelength lattices can be designed using adiabatic dressing of state-dependent lattices 44 , multi-photon optical transitions 45,46 , spin-dependent optical lattices with time-dependent modulations 47 , as well as nanoplasmonic systems 48 , vortex arrays in superconducting films 49 , periodically patterned graphene monolayers 50 , magnetic-film atom chips 51 , and photonic crystals [52][53][54] . These interesting proposals already stimulated studies related to many-body physics in such lattices, in particular the analysis of the Hubbard model and engineering of spin-spin Hamiltonians 52 .
In the considered case of p x + ip y pairing in the 2D lattice, putting b = 50 nm, for the same k F b as above the Fermi energy for NaLi molecules exceeds 200 nK (n ≈ 4 × 10 8 cm −2 ). Then, for the same κ ~ 1 and ⁎ k r F eff approaching unity we have T c ~ 20 nK, which is twice as high as in free space. An additional advantage of the lattice system is the foreseen quantum information processing, since addressing qubits in the lattice is much easier than in free space.
Note that there is a (second-order) process, in which the interaction between two identical fermions belonging to the lowest Bloch band provides a virtual transfer of one of them to a higher band. Then, the two fermions may get to the same lattice site and undergo the inelastic process of collisional relaxation. The rate constant of this second-order process is roughly equal to the rate constant in free space, multiplied by the ratio of the scattering amplitude (divided by the elementary cell area) to the frequency of the potential well in a given lattice site (the difference in the energies of the Bloch bands). This ratio originates from the virtual transfer of one of the fermions to a higher band and does not exceed (ξ/b) 2 . Even in not a deep lattice, where m * /m is 2 or 3, we have (ξ/b) 2 < 0.1. Typical values of the rate constant of inelastic relaxation in free space are ~10 −8 -10 −9 cm 2 /s 16 , and hence in the lattice it will be lower than 10 −9 or even 10 −10 cm 2 /s. Thus, the rate of this process is rather low and for densities approaching 10 9 cm −2 the decay time will be on the level of seconds or even tens of seconds.

Interlayer p-wave Superfluid of Fermionic Polar Molecules in a Bilayer System
Another interesting novel superfluid of fermionic polar molecules is expected in a bilayer system, where dipoles are oriented perpendicularly to the layers and in opposite directions in different layers.
Such a bilayer configuration, but with all dipoles oriented in the same direction, has been considered in refs 40-43. As found, it should form an interlayer s-wave superfluid, where Cooper pairs are formed by dipoles of different layers due to the s-wave dipolar interaction between them.
For the dipoles of one layer that are opposite to the dipoles of the other one, the picture of interlayer pairing is different. The s-wave pairing is practically impossible, and the system may form p-wave and higher partial wave superfluids. This type of bilayer systems can be created by putting polar molecules with rotational moment J = 0 in one layer, and molecules with J = 1 in the other. Then, applying an electric field (perpendicular to the layers) one gets a field-induced average dipole moment of J = 0 molecules parallel to the field, and the dipole moment of J = 1 molecules oriented in the opposite direction. One should also prevent a flip-flop process in which the dipole-dipole interaction between given J = 1 and J = 0 molecules reverses their dipoles, thus inducing a rapid three-body decay in collisions of a dipole-reversed molecule with two original ones. This can be done by making the electric field inhomogeneous, so that it is larger in the layer with J = 0 molecules and the flip-flop process requires an increase in the Stark energy. This process will be suppressed if the difference in the Stark energies of molecules in the layers significantly exceeds the Fermi energy, which is a typical kinetic energy of the molecules (~100 nK for the example considered below). This is realistic for present facilities.
For the dipole moment close to 1 Debye and the interlayer spacing of 50 nm, one thus should have the field gradient (perpendicularly to the layers) significantly exceeding 0.5 kV/cm 2 . This could be done by using electrodes consisting of four rods, and even a higher gradient ~30 kV/cm 2 should be achievable 64,65 . By changing the positions of the rods one can obtain the field gradient exceeding 0.5 kV/cm 2 in the direction perpendicular to the layers of the bilayer system. The field itself will not be exactly perpendicular to the layers and there will also be the field gradient parallel to the layers. This, however, does not essentially influence the physics.
The potential of interaction between two molecules belonging to different layers has the form: where L is the interlayer spacing, r is the in-layer separation between the molecules, and − d 2 is the scalar product of the average dipole moments of these molecules. The potential V L (r) is repulsive for < r L 2 and attractive at larger r. The potential well is much more shallow than in the case of all dipoles oriented in the same direction, which was considered in refs 40-42. We have checked that s-wave interlayer dimers, which exist at any r * /L, are weakly bound even for r * /L ≈ 3. Their binding energy at  * r L / 3 is much smaller than the Fermi energy at least for k F L > 0.1. For such r * /L, interlayer dimers with orbital angular momenta |l| ≥ 1 do not exist. We thus are dealing with purely fermionic physics.
For the analysis of the superfluid pairing we are interested in particle momenta k ~ k F . As well as in the case of all dipoles oriented in the same direction [40][41][42][43] , under the condition k F r * ≪ 1 (where r * = md 2 / 2  ) the amplitude of interlayer interaction is obtained in the Born approximation. The Fourier transform of the potential (9) is L 2  and in the first Born approximation the on-shell amplitude of the l-wave scattering at k = k F reads (see Methods): The s-wave amplitude is positive, i.e. the s-wave channel corresponds to repulsion. Note that for extremely low collision energies comparable with the dimer binding energy, where the Born approximation is not accurate, the s-wave scattering amplitude can be negative. This, however, does not lead to superfluid s-wave pairing.
The channels with |l| ≥ 1 correspond to attraction. A straightforward calculation shows that for  . k L 0 7 F the largest is the p-wave amplitude and, hence, at sufficiently low temperatures the system will be an interlayer p-wave superfluid. As for d-wave and higher partial wave superfluids, they are possible only at extremely low temperatures. Thus, we confine ourselves to the p-wave pairing and employ the BCS approach.
A detailed analysis of the gap equation (1), which includes first and second order contributions to the scattering amplitude and Gor'kov-Melik-Barkhudarov corrections, is given in Methods. The critical temperature for the p-wave superfluidity proves to be (see Methods):  Fig. 3 and Methods). Creating the bilayer system by using a 1D subwavelength lattice we may have L ≈ 50 nm. In this case, for k F L = 0.15 the Fermi energy of NaLi molecules is close to 100 nK, and the critical temperature for k F r * approaching 0.5 is about 10 nK.
For completeness, we also consider the regime of strong interactions within a single layer. Assuming that the coupling between the layers is still fairly weak, we have superfluid (interlayer) pairing between quasiparticles. Related problems have been discussed for coupled 2D Fermi liquids as models for layered superconductors 8 . In this case, we replace the bare mass m by the effective mass m * and account for renormalization of the fermionic Green functions by a factor Z < 1 66 . Then, the expression for the transition temperature takes the form:

Conclusions
We have demonstrated the emergence of the topological p x + ip y superfluid for identical microwave-dressed fermionic polar molecules in a 2D lattice. Another novel p-wave superfluid is found to emerge for fermionic molecules in a bilayer system, with dipoles of one layer opposite to the dipoles of the other one. In both cases the use of subwavelength lattices with a period  b 50 nm (creation of the bilayer system with the interlayer spacing  L 50 nm) allows one to obtain superfluid transition temperature of the order of tens of nanokelvins. This opens interesting prospects for topologically protected quantum information processing with p x + ip y superfluids in 2D lattices. The interlayer p-wave superfluid in bilayer systems, together with the earlier proposed s-wave interlayer superfluid [40][41][42][43] and superfluids in multilayer fermionic systems 68 , can be a starting point for the creation of more sophisticated layered structures.
Superfluidity itself can be detected in the same way as in the case of s-wave superfluids 69,70 . Rotating the p x + ip y superfluid and inducing the appearance of vortices one can find signatures of Majorana modes on the vortex cores in the RF absorption spectrum 71 . Eventually, one can think of revealing the structure of the order parameter by visualizing vortex-related dips in the density profile on the approach to the strongly interacting regime, where these dips should be pronounced at least in time-of-flight experiments.

Scattering problem and superfluid pairing of microwave-dressed polar molecules in a 2D lattice.
As we concluded in the main text, in the low momentum limit at a small filling factor the system of lattice polar molecules is equivalent to that of molecules with effective mass m * in free space. We now demonstrate this explicitly by the calculation of the off-shell scattering amplitude f(k′, k). For our problem the main part of the scattering amplitude can be obtained in the Born approximation 16 .
In the lattice the scattering amplitude is, strictly speaking, the function of both incoming quasimomenta q 1 , q 2 and outgoing quasimomenta ′ ′ q q , 1 2 . However, in the low-momentum limit where qb ≪ 1, taking into account the momentum conservation law the amplitude becomes the function of only relative momenta k = (q 1 − q 2 )/2 and ′ = ′ − ′ k q q ( ) /2 1 2 . For the off-shell scattering amplitude the first Born approximation gives: where V(r 1 − r 2 ) is given by Eq. (3) of the main text, and S is the surface area. The last line of Eq. (14) is obtained assuming the tight-binding regime, where the single particle wavefunction is Here, the index j labels the lattice sites located at the points r j , and N = S/b 2 is the total number of the sites. The particle wavefunction in a given site j has extension ξ 0 and is expressed as 2 . In the low-momentum limit we may replace the summation over j and j′ by the integration over d 2 r j and ′ d r j 2 taking into account that b 2 ∑ j transforms into ∫ d 2 r j . This immediately yields 3 and the p-wave part of the scattering amplitude is obtained multiplying Eq. (16) by exp(− iφ) and integrating over dφ/2π, where φ is the angle between the vectors k and k′. This is the same result as in free space (see 16 ). The on-shell amplitude (k = k′ ) can be written as  is the effective dipole-dipole distance in the lattice. The applicability of the Born approximation assumes that  ⁎ kr 1 eff , which is clearly seen by calculating the second order correction to the scattering amplitude.
Up to the terms ⁎ kr ( ) eff 2 , the on-shell scattering amplitude following from the solution of the scattering problem for particles with mass m * , is given by 16 : where the numerical coefficient B comes from short-range physics. For calculating B we introduce a perfectly reflecting wall at intermolecular distances r ~ b, which takes into account that two fermions practically can not get to one and the same lattice site. The coefficient B depends on the ratio ⁎ r b / eff , and we show this dependence in Fig. 2a. The treatment of the superfluid pairing is the same as in ref. 16, including the Gorkov-Melik-Barkhudarov correction. We should only replace the mass m with m * . The expression for the critical temperature then becomes: , and it is displayed in Fig. 2b as a function  Superfluid pairing of fermionic polar molecules in a bilayer system. For the interlayer interaction potential V L (r) given by equation (9) in the main text, the scattering amplitude for k F r * ≪ 1 can be calculated in the Born approximation 40 . The p-wave part of the first order contribution to the off-shell amplitude is and J 1 is the Bessel function. Regarding the second order contribution, for the solution of the gap equation we only need the on-shell p-wave part, which is given by where ψ(k, r) is the true wavefunction of the p-wave relative motion with momentum k, normalized such that for In order to calculate the superfluid transition temperature we use the BCS approach along the lines of ref. 16. We consider temperature T tending to T c from below and rely on the renormalized gap equation (1). For the p-wave pairing the order parameter is ∆ = ∆ ϕ k e ( ) k i k , and we then multiply Eq. (1) by ϕ − e i k and integrate over dϕ k′ and dϕ k . As a result, we obtain the same equation (1) in which Δ k and Δ k′ are replaced with Δ (k) and Δ (k′ ), the off-shell scattering amplitude f(k′, k) with its p-wave part, and δV(k′, k) with its p-wave part Calculating the contribution of the pole in the second term in square brackets and using Eq. (24) we obtain where the symbol P stands for the principal value of the integral. In the first term in the right hand side of Eq. (25) we divide the region of integration into two parts: The contribution to the p-wave order parameter from the first region we denote as Δ (1) (k), and the contribution from the second region as Δ (2) (k). The contribution of the second term in right hand side of equation (25) is denoted as Δ (3) (k). We first notice that the main contribution to Δ (k) comes from k′ close to k F . Retaining only f 1 , which is proportional to kr * , in the off-shell scattering amplitude and omitting the second term in the right hand side of Eq. (25) (which is proportional to (kr * ) 2 ) we obtain Putting k = k F and performing the integration in the first region in the first term of Eq (25), where E F − ω < E k′ < E F + ω, we may put Δ (k′ ) = Δ (k F ) and and taking into account that the contribution of the second term in square brackets is zero, we obtain: with C = 0.577 being the Euler constant, and ρ(k F ) = m/2π 2  the density of states. In the second region, where E k′ > E F + ω or E k′ < E F − ω, we put  ξ = | ′ 1/2 k and keep only f 1 in the scattering amplitude. For k = k F the integral over E k′ from E F + ω to ∞ vanishes. In the integral from 0 to E F − ω we use Δ (k′ ) from Eq. (26) and find . Then, we consider the Gor'kov-Melik-Barkhudarov corrections to the bare interaction of the molecules in the bilayer. These many-body corrections are second order in (k F r * ) and are described by four diagrams (for details, see 16,41,56 ). For the case of p-wave superfluidity of identical fermionic polar molecules they have been considered in ref. 16. They have been also studied for the interlayer s-wave superfluidity of dipoles oriented in the same direction in ref. 41. We are interested in the case of sufficiently small k F L. Following the same treatment as in refs 16 and 41, in the limit of k F L → 0 we obtain: where α .  10 57. The dominant contribution to this result comes from the diagram containing a bubble in the interaction line (diagram a) in refs 16 and 41). This contribution strongly decreases with increasing k F L. In particular, for .  k L 0 15 F we have α .  2 8, and α .  2 2 when increasing k F L to 0.2. Comparing δV with the scattering amplitude f 1 (k F ) we see that for not very small k F r * the perturbative treatment of the Gor'kov-Melik-Barkhudarov corrections is adequate for  . k L 0 15 F . We therefore confine ourselves to these values of k F L. Performing the integration in the second term of Eq. (25) we obtain the contribution of the Gor'kov-Melik-Barkhudarov corrections to the order parameter: where we put ω ~ E F in the terms proportional to (k F r * ) 2 . We should also recall that the bare mass m should be replaced with the effective mass m * = m[1 − (4/3π)k F r * ] which has been found in refs 41 and 72. Since the relative difference between m * and m is small as k F r * , it is sufficient to replace m with m * only in the multiple r * ~ m in the first term of Eq. (32). This leads to the appearance of a new term in the right hand side of equation (32). Then, dividing both sides of Eq. (32) by Δ (k F ) we obtain for the critical temperature: The dependence of F and β on k F L is shown in Fig. 3. We stop at k F L = 0.3 because for larger values of this parameter the function F is so large that the critical temperature will be negligible.