Introduction

Non-conventional superconductors and superfluids attract a great deal of interest due to their non-trivial transport properties and/or topological behavior1,2,3,4,5,6,7,8,9,10,11. This behavior has been actively discussed in two dimensions (2D) for the px + ipy superfluid of identical fermions, where Cooper pairs have orbital angular momentum equal to unity12,13,14,15,16,17. Quantized vortices in this superfluid carry zero-energy Majorana modes on their cores3,18,19. These modes cause the vortices to obey non-Abelian exchange statistics, which is a basis for topologically protected quantum information processing20,21. However, the p-wave topological superfluid of ultracold atoms is either collisionally unstable near a Feshbach resonance, or has a vanishingly low superfluid transition temperature far from the resonance22,23,24.

Successful experiments on the creation of ground-state ultracold polar molecules25,26,27,28,29,30,31,32,33,34,35,36 opened fascinating prospects for obtaining non-conventional superfluids37,38,39. In particular, microwave-dressed polar molecules confined to 2D may acquire an attractive dipole-dipole tail in the interaction potential, which ensures the emergence of collisionally stable p-wave superfluid with a reachable transition temperature16,17. Another interesting system concerns fermionic polar molecules in a bilayer geometry. Here they may form interlayer superfluids in which Cooper pairs consist of molecules belonging to different layers40,41,42,43.

It this paper we consider novel p-wave superfluids of fermionic polar molecules in 2D lattice geometries (see Fig. 1). It is shown that a collisionally stable topological px + ipy superfluid of identical microwave-dressed polar molecules may emerge in a 2D lattice due to a long-range character of the dipole-dipole interaction. We also show how one can get a p-wave interlayer superfluid of fermionic polar molecules in a bilayer geometry, which can be a quantum simulator of superconductivity in layered condensed matter systems7,8. It is crucial to rely on the recently proposed subwavelength lattices44,45,46,47,48,49,50,51,52,53,54, where the lattice constant (interlayer spacing in the bilayer system) can be as small as about 50 nm. An increase of energy scales in such lattices makes it realistic to obtain sizeable transition temperatures of the order of tens of nanokelvins.

Figure 1
figure 1

Setups for p-wave superfluids of polar molecules: (a) polar molecule in an external microwave field Eac rotating in the plane perpendicular to the stationary field Edc (upper part), and microwave-dressed polar molecules loaded in a 2D lattice (lower part); (b) bilayer system of polar molecules with dipole moments in the upper and lower layers, oriented opposite to each other.

Full size image

General Relations and Qualitative Arguments

The superfluid pairing of identical fermions is characterized by the order parameter Δ(r, r′) = V(r − r′) , where V(r − r′) is the interaction potential, the symbol 〈...〉 denotes the statistical average, and 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 Tc 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 approach55. For both spinless and spin-1/2 fermions the renormalized gap equation for the order parameter in the momentum space, , reads (see16,17 and references therein):

where f(k′, k) is the off-shell scattering amplitude and Ek =k2/2m with m being the particle mass. The single particle excitation energy is given by where μ is the chemical potential, and . For weak interactions chemical potential coincides with the Fermi energy (kF 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-Barkhudarov56, 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 = kF, and notice that the main contribution to the integral over k′ in Eq. (1) comes from k′ close to kF. At temperatures T tending to the critical temperature Tc from below, we put in . For the pairing channel related to the interaction with orbital angular momentum l, this immediately leads to an estimate:

The quantity ρ(kF) = m/2π in the exponent of Eq. (2) is the density of states on the Fermi surface, and fl(kF) is the on-shell scattering amplitude.

In the lattice with a period b satisfying the condition kFb ≪ 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 ρ(kF) is then given by the same expression, with m replaced by m*. Thus, the BCS exponent [ρ(kF)| fl(kF)|]−1 in the lattice is smaller than in free space at the same kF (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 Tc 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 sites59. 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 potential16,17,62,63:

Here d is an effective dipole moment, and we assume that Eq. (3) is valid at intermolecular distances . This leads to superfluid p-wave pairing of the molecules. In free space the emerging ground state is the topological px + ipy superfluid, and the leading part of the scattering amplitude can be obtained in the first Born approximation16,17. We assume the weakly interacting regime at a small filling factor in the lattice, kFb ≪ 1.

The Hamiltonian of the system is , with

where , 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 = q2/2m*, where m* > m is the effective mass in the lowest Bloch band. The quantity describes the interaction between the molecules and is given by

where is the field operator of a particle in the lattice site j located at rj in the coordinate space. At a small filling factor in the low momentum limit, the main contribution to the matrix elements of comes from intermolecular distances (see Methods). Therefore, we may replace the summation over rj and by the integration over d2rj and . As a result the Hamiltonian of the system reduces to

where the first term in the right hand side is (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 = kF, 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

where , 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 ~kF. Therefore, the low-momentum limit requires the inequality kFb ≪ 1.

The solution of the gap equation (1) then leads to the px + ipy superfluid with the critical temperature (see Methods):

where the coefficient κ is related to B and depends on the ratio (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 EF is smaller by a factor of m/m*, and the effective dipole-dipole distance 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 kF) 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 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, kFb = 0.35) the Fermi energy for the lightest alkaline polar molecules NaLi is about 10 nK (n ≈ 2 × 107 cm−2). Then, even for approaching unity the critical temperature is only of the order of a nanokelvin (for kFb = 0.35 and Fig. 2 in Methods gives κ ~ 1).

Figure 2
figure 2

Coefficients B and κ as functions of .

Full size image

The picture is quite different in recently introduced subwavelength lattices44,45,46,47,48,49,50,51,52, where the lattice constant can be as small as 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 lattices44, multi-photon optical transitions45,46, spin-dependent optical lattices with time-dependent modulations47, as well as nanoplasmonic systems48, vortex arrays in superconducting films49, periodically patterned graphene monolayers50, magnetic-film atom chips51, and photonic crystals52,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 Hamiltonians52.

In the considered case of px + ipy pairing in the 2D lattice, putting b = 50 nm, for the same kFb as above the Fermi energy for NaLi molecules exceeds 200 nK (n ≈ 4 × 108 cm−2). Then, for the same κ ~ 1 and approaching unity we have Tc ~ 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 cm2/s16, and hence in the lattice it will be lower than 10−9 or even 10−10 cm2/s. Thus, the rate of this process is rather low and for densities approaching 109 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, 41, 42, 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/cm2. This could be done by using electrodes consisting of four rods, and even a higher gradient ~30 kV/cm2 should be achievable64,65. By changing the positions of the rods one can obtain the field gradient exceeding 0.5 kV/cm2 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 −d2 is the scalar product of the average dipole moments of these molecules. The potential VL(r) is repulsive for 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, 41, 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 is much smaller than the Fermi energy at least for kFL > 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 ~ kF. As well as in the case of all dipoles oriented in the same direction40,41,42,43, under the condition kFr* ≪ 1 (where r* = md2/) the amplitude of interlayer interaction is obtained in the Born approximation. The Fourier transform of the potential (9) is

and in the first Born approximation the on-shell amplitude of the l-wave scattering at k = kF 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 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):

and for not very small kFr* the validity of the perturbative treatment of the Gor’kov-Melik-Barkhudarov corrections requires (see Methods). The functions F(kFL) and β(kFL) are given in Methods. For kFL ranging from 0.15 to 0.3 the function F increases from 3.4 to 5, and the coefficient β is fairly large, being about 80 at kFL = 0.15 (see Fig. 3 and Methods).

Figure 3
figure 3

The dependence of F and β on kFL.

Full size image

Creating the bilayer system by using a 1D subwavelength lattice we may have L ≈ 50 nm. In this case, for kFL = 0.15 the Fermi energy of NaLi molecules is close to 100 nK, and the critical temperature for kFr* 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 superconductors8. 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 < 166. Then, the expression for the transition temperature takes the form:

where we can not determine the pre-exponential coefficient. Therefore, Eq. (13) only gives an order of magnitude of Tc. For kFL = 0.3 and L ≈ 50 nm the Fermi energy of NaLi molecules is about 400 nK, and for, let say, kFr* ≈ 2 the dimer physics is still not important. Then, using the effective mass and factor Z from the Monte Carlo calculations67 one may think of superfluid transition temperatures of the order of several tens of nanokelvins.

Conclusions

We have demonstrated the emergence of the topological px + ipy 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 nm (creation of the bilayer system with the interlayer spacing 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 px + ipy superfluids in 2D lattices. The interlayer p-wave superfluid in bilayer systems, together with the earlier proposed s-wave interlayer superfluid40,41,42,43 and superfluids in multilayer fermionic systems68, 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 superfluids69,70. Rotating the px + ipy superfluid and inducing the appearance of vortices one can find signatures of Majorana modes on the vortex cores in the RF absorption spectrum71. 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.

Methods

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 approximation16.

In the lattice the scattering amplitude is, strictly speaking, the function of both incoming quasimomenta q1, q2 and outgoing quasimomenta . 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 = (q1 − q2)/2 and . For the off-shell scattering amplitude the first Born approximation gives:

where V(r1 − r2) 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 rj, and N = S/b2 is the total number of the sites. The particle wavefunction in a given site j has extension ξ0 and is expressed as . In the low-momentum limit we may replace the summation over j and j′ by the integration over d2rj and taking into account that b2 ∑j transforms into ∫d2rj. This immediately yields

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 (see16). The on-shell amplitude (k = k′) can be written as , where is the effective dipole-dipole distance in the lattice. The applicability of the Born approximation assumes that , which is clearly seen by calculating the second order correction to the scattering amplitude.

Up to the terms , the on-shell scattering amplitude following from the solution of the scattering problem for particles with mass m*, is given by16:

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 , 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:

where , and it is displayed in Fig. 2b as a function of .

Superfluid pairing of fermionic polar molecules in a bilayer system

For the interlayer interaction potential VL(r) given by equation (9) in the main text, the scattering amplitude for kFr* ≪ 1 can be calculated in the Born approximation40. The p-wave part of the first order contribution to the off-shell amplitude is

where

and J1 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

In fact, the true p-wave scattering amplitude follows from the exact relation

where ψ(k, r) is the true wavefunction of the p-wave relative motion with momentum k, normalized such that for r → ∞ we have with being the Hankel function. This amplitude is complex and it is related to the real amplitude given by equations (19)–(22), as

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 Tc from below and rely on the renormalized gap equation (1). For the p-wave pairing the order parameter is , and we then multiply Eq. (1) by 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: |Ek′ − EF| < ω and |Ek′ − EF| > ω, where Δ(kF), Tc ≪ ω ≪ EF. 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 kF. Retaining only f1, 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 = kF and performing the integration in the first region in the first term of Eq (25), where EF − ω < Ek′ < EF + ω, we may put Δ(k′) = Δ(kF) and . Then, putting in 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 ρ(kF) = m/2π the density of states.

In the second region, where Ek′ > EF + ω or Ek′ < EF − ω, we put and keep only f1 in the scattering amplitude. For k = kF the integral over Ek′ from EF + ω to ∞ vanishes. In the integral from 0 to EF − ω we use Δ(k′) from Eq. (26) and find

where

and .

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 (kFr*) and are described by four diagrams (for details, see16,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 kFL. Following the same treatment as in refs 16 and 41, in the limit of kFL → 0 we obtain:

where . 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 kFL. In particular, for we have , and when increasing kFL to 0.2. Comparing δV with the scattering amplitude f1(kF) we see that for not very small kFr* the perturbative treatment of the Gor’kov-Melik-Barkhudarov corrections is adequate for . We therefore confine ourselves to these values of kFL.

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:

the sum of Eqs (27), (28) and (31) yields

where we put ω ~ EF in the terms proportional to (kFr*)2. We should also recall that the bare mass m should be replaced with the effective mass m* = m[1 − (4/3π)kFr*] which has been found in refs 41 and 72. Since the relative difference between m* and m is small as kFr*, 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 Δ(kF) we obtain for the critical temperature:

where

and

The dependence of F and β on kFL is shown in Fig. 3. We stop at kFL = 0.3 because for larger values of this parameter the function F is so large that the critical temperature will be negligible.

Additional Information

How to cite this article: Fedorov, A. K. et al. Novel p-wave superfluids of fermionic polar molecules. Sci. Rep. 6, 27448; doi: 10.1038/srep27448 (2016).